text:
## we implement a method to load the original MAT file of then
## binarized MNIST

using MAT
using MLDataUtils
using Distributions
import Distributions: logpdf
using MLDatasets: MNIST

img(x) = Gray.(permutedims(reshape(x, 28, 28)))

function load_binarized_mnist(batch_size, flatten)
    # we start by loading the training set
    file = open("Original_MNIST_binarized/binarized_mnist_train.amat", "r")
    train_imgs = readlines(file)
    #C =  findall(x->x=="1", vec(train_imgs[2,:]))
    #@show C
    train_imgs = [parse.(Float64, split(train_imgs[i]," ")) for i in 1:size(train_imgs, 1)]
    example = train_imgs[1]
    save("check_while_load_fin.png", img(example))
    #C =  findall(x->x==1.0, vec(train_imgs[2,:][1]))
    #@show C
    #train_imgs = permutedims(reshape(hcat(train_imgs...), (length(train_imgs[1]), length(train_imgs))))
    train_imgs = hcat(train_imgs...)
    @show(size(train_imgs))
    example = train_imgs[:,1]
    save("check_while_load_fin2.png", img(example))
    #C =  findall(x->x==1.0, vec(train_imgs[2,:]))
    #@show C
    @show(size(train_imgs))

    # valid test
    file = open("Original_MNIST_binarized/binarized_mnist_valid.amat", "r")
    valid_imgs = readlines(file)
    valid_imgs = [parse.(Float64, split(valid_imgs[i]," ")) for i in 1:size(valid_imgs, 1)]
    valid_imgs = hcat(valid_imgs...)
    #valid_imgs = permutedims(reshape(hcat(valid_imgs...), (length(valid_imgs[1]), length(valid_imgs))))
    #valid_imgs = permutedims(reshape(hcat(valid_imgs...), (length(valid_imgs), length(valid_imgs[1]))))
    #@show(size(valid_imgs))

    # test set
    file = open("Original_MNIST_binarized/binarized_mnist_test.amat", "r")
    test_imgs = readlines(file)
    test_imgs = [parse.(Float64, split(test_imgs[i]," ")) for i in 1:size(test_imgs, 1)]
    test_imgs = hcat(test_imgs...)
    #test_imgs = permutedims(reshape(hcat(test_imgs...), (length(test_imgs[1]), length(test_imgs))))
    #test_imgs = permutedims(reshape(hcat(test_imgs...), (length(test_imgs), length(test_imgs[1]))))

    #@show(size(test_imgs))

    # flatten
    if flatten
        # we want to reshape them into (L, BS)
        # training set
        x_train = train_imgs
        # validation set
        x_valid = valid_imgs
        # test set
        x_test = test_imgs

    else
        # we want to return the images (H,W,C,n) where n = |set|
        # training set
        x_train = reshape(train_imgs, 28, 28, 1, size(train_imgs, 2))
        # validation set
        x_valid = reshape(valid_imgs, 28, 28, 1, size(valid_imgs, 2))
        # test set
        x_test = reshape(test_imgs, 28, 28, 1, size(test_imgs, 2))

    end

    return x_train, x_valid, x_test
end

#const batch_size = 64
#x_train, x_valid, x_test = load_binarized_mnist(batch_size, false)

## function to load MNIST (credit: Jesse Bettencourt)
function loadMNIST(batch_size)
    # we use MLDataUtils LabelEnc for the one-hot-conversion
    onehot(labels_raw) =  convertlabel(LabelEnc.OneOfK, labels_raw, LabelEnc.NativeLabels(collect(0:9)))
    # load MNIST from the library
    imgs, labels_raw = MNIST.traindata();
    # process the images into (H,W,C,BS) batches
    @show(size(imgs))
    x_train = reshape(imgs, size(imgs, 1), size(imgs, 2), 1, size(imgs, 3))
    @show(size(x_train))
    x_train = batchview(x_train, batch_size)
    @show(size(x_train))
    # one-hot and batch the labels
    y_train = onehot(labels_raw)
    y_train = batchview(y_train, batch_size)

    ## we should do the same with the test set
    imgs, labels_raw = MNIST.testdata();
    # process the images into (H,W,C,BS) batches
    x_test = reshape(imgs, size(imgs, 1), size(imgs, 2), 1, size(imgs, 3))
    x_test = batchview(x_test, batch_size)
    # one-hot and batch the labels
    y_test = onehot(labels_raw)
    y_test = batchview(y_test, batch_size)

    return x_train, y_train, x_test, y_test
end

#const batch_size = 64
#x_train, x_valid, x_test = load_binarized_mnist(batch_size, true)
#@show size(x_train)
#x_train, y_train, x_test, y_test = loadMNIST(batch_size)


text:
using BinaryBuilder

src_version = v"6.1.0"  # also change below in script

# Collection of sources required to build PROJ
sources = [
    "http://download.osgeo.org/proj/proj-$src_version.tar.gz" =>
    "676165c54319d2f03da4349cbd7344eb430b225fe867a90191d848dc64788008",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd proj-6.1.0

# sqlite needed to build proj.db, so this should not be the
# cross-compiled one since it needs to be executed on the host
apk add sqlite

if [[ ${target} == *mingw* ]]; then
    SQLITE3_LIBRARY=$prefix/bin/libsqlite3-0.dll
elif [[ ${target} == *darwin* ]]; then
    SQLITE3_LIBRARY=$prefix/lib/libsqlite3.dylib
else
    SQLITE3_LIBRARY=$prefix/lib/libsqlite3.so
fi

mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain \
      -DSQLITE3_INCLUDE_DIR=$prefix/include \
      -DSQLITE3_LIBRARY=$SQLITE3_LIBRARY \
      -DHAVE_PTHREAD_MUTEX_RECURSIVE_DEFN=1 \
      -DBUILD_LIBPROJ_SHARED=ON \
      ..
cmake --build .
make install

# add proj-datumgrid files directly to the result
wget https://download.osgeo.org/proj/proj-datumgrid-1.8.tar.gz
tar xzf proj-datumgrid-1.8.tar.gz -C $prefix/share/proj/
"""

platforms = supported_platforms()
platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libproj", :libproj),

    ExecutableProduct(prefix, "cct", :cct_path),
    ExecutableProduct(prefix, "cs2cs", :cs2cs_path),
    ExecutableProduct(prefix, "geod", :geod_path),
    ExecutableProduct(prefix, "gie", :gie_path),
    ExecutableProduct(prefix, "proj", :proj_path),
    ExecutableProduct(prefix, "projinfo", :projinfo_path),

    # complete contents of share/proj, must be kept up to date
    FileProduct(prefix, joinpath("share", "proj", "CH"), :ch_path),
    FileProduct(prefix, joinpath("share", "proj", "GL27"), :gl27_path),
    FileProduct(prefix, joinpath("share", "proj", "ITRF2000"), :itrf2000_path),
    FileProduct(prefix, joinpath("share", "proj", "ITRF2008"), :itrf2008_path),
    FileProduct(prefix, joinpath("share", "proj", "ITRF2014"), :itrf2014_path),
    FileProduct(prefix, joinpath("share", "proj", "nad.lst"), :nad_lst_path),
    FileProduct(prefix, joinpath("share", "proj", "nad27"), :nad27_path),
    FileProduct(prefix, joinpath("share", "proj", "nad83"), :nad83_path),
    FileProduct(prefix, joinpath("share", "proj", "null"), :null_path),
    FileProduct(prefix, joinpath("share", "proj", "other.extra"), :other_extra_path),
    FileProduct(prefix, joinpath("share", "proj", "proj.db"), :proj_db_path),
    FileProduct(prefix, joinpath("share", "proj", "world"), :world_path),

    # part of files from proj-datumgrid which are added to the default ones
    # all are added but only the few below are checked if they are added
    # note that none of proj-datumgrid-europe, proj-datumgrid-north-america,
    # proj-datumgrid-oceania, proj-datumgrid-world is added by default,
    # though users are free to add them to the rest themselves
    FileProduct(prefix, joinpath("share", "proj", "alaska"), :alaska_path),
    FileProduct(prefix, joinpath("share", "proj", "conus"), :conus_path),
    FileProduct(prefix, joinpath("share", "proj", "egm96_15.gtx"), :egm96_15_path),
    FileProduct(prefix, joinpath("share", "proj", "ntv1_can.dat"), :ntv1_can_path),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaDatabases/SQLiteBuilder/releases/download/v0.10.0/build_SQLite.v3.28.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "PROJ", src_version, sources, script, platforms, products, dependencies)

ext: jl


text:
function normal_model(t :: Trace, data :: Vector{Float64})
    loc = sample(t, :loc, Normal())
    scale = sample(t, :scale, LogNormal())
    for i in 1:length(data)
        observe(t, (:obs, i), Normal(loc, scale), data[i])
    end
end

@testset "likelihood weighting 1" begin
    true_loc = 2.0
    true_scale = 1.0
    @info "True loc = $true_loc"
    @info "True scale = $true_scale"
    n_dpts = 10
    data = true_loc .+ true_scale .* randn(n_dpts)
    prior_samples = forward_sampling(normal_model; params = (data,), num_iterations = 100)
    @info "Prior E[loc] = $(mean(prior_samples[:loc]))"
    @info "Prior E[scale] = $(mean(prior_samples[:scale]))"

    posterior = likelihood_weighting(normal_model, data; nsamples = 1000)
    log_p_x = log_evidence(posterior)
    @info "Log evidence = $log_p_x"

    post_loc = sample(posterior, :loc, 2000)
    post_scale = sample(posterior, :scale, 2000)
    @info "Posterior E[loc] = $(mean(post_loc))"
    @info "Posterior E[scale] = $(mean(post_scale))"
end

function is_proposal(t :: Trace, data :: Vector{Float64})
    loc_guess = mean(data)
    std_guess = std(data)
    propose(t, :loc, Normal(loc_guess, std_guess * 0.25))
    propose(t, :scale, truncated(Normal(std_guess, std_guess * 0.25), 0.1, Inf))
end

@testset "importance sampling 1" begin
    true_loc = 2.0
    true_scale = 1.0
    @info "True loc = $true_loc"
    @info "True scale = $true_scale"
    n_dpts = 10
    data = true_loc .+ true_scale .* randn(n_dpts)
    prior_samples = forward_sampling(normal_model; params = (data,), num_iterations = 100)
    @info "Prior E[loc] = $(mean(prior_samples[:loc]))"
    @info "Prior E[scale] = $(mean(prior_samples[:scale]))"

    posterior = importance_sampling(normal_model, is_proposal; params = (data,), nsamples = 1000)
    log_p_x = log_evidence(posterior)
    @info "Log evidence = $log_p_x"

    post_loc = sample(posterior, :loc, 2000)
    post_scale = sample(posterior, :scale, 2000)
    @info "Posterior E[loc] = $(mean(post_loc))"
    @info "Posterior E[scale] = $(mean(post_scale))"
end


text:
module ADMM

using LinearAlgebra: norm

"""
    admm(x0, z0, u0, prox_f, prox_g, ϵ=1e-4)

Alternating Direction Method of Multipliers for minimizing the sum of two functions f and g
(also known as Douglas-Rachford algorithm).

# Arguments
- `x0::Vector`: starting estimate for x primal variable.
- `z0::Vector`: starting estimate for z consensus variable.
- `u0::Vector`: starting estimate for u dual variable.
- `prox_f`: proximal operator of function f.
- `prox_g`: proximal operator of function g.

"""
function admm(
    x0::Vector{T}, u0::Vector{T}, prox_f::Function, prox_g::Function, ϵ::T=1e-4
) where {T <: AbstractFloat}
    xk, uk = copy(x0), copy(u0)
    xkk, zkk, ukk = xk .+ ϵ, xk .+ ϵ, uk .+ ϵ

    while norm(xkk - xk) / norm(xkk) + norm(ukk - uk) / norm(ukk) >= ϵ
        xk, uk = copy(xkk), copy(ukk)

        zkk = prox_g(xk + uk)
        xkk = prox_f(zkk - uk)
        ukk = uk + xkk - zkk
    end
    return xkk
end

end

max_stars_repo_name: jzazo/InverseValueProblem

max_issues_repo_name: jzazo/InverseValueProblem

max_forks_repo_name: jzazo/InverseValueProblem


text:
@testset "approx_posterior_gp" begin
    rng = MersenneTwister(123456)
    N_cond = 3
    N_a = 5
    N_b = 6

    # Specify prior.
    f = GP(sin, Matern32Kernel())

    # Sample from prior.
    x = collect(range(-1.0, 1.0; length=N_cond))
    fx = f(x, 1e-15)
    y = rand(rng, fx)

    # Construct posterior.
    f_post = posterior(fx, y)

    # Construct optimal approximate posterior.
    f_approx_post = posterior(VFE(f(x, 1e-12)), fx, y)

    # Verify that approximate posterior ≈ posterior at various inputs.
    x_test = randn(rng, 100)
    @test mean(f_post, x_test) ≈ mean(f_approx_post, x_test)
    @test cov(f_post, x_test) ≈ cov(f_approx_post, x_test)

    # Verify that AbstractGP interface is implemented fully and consistently.
    a = collect(range(-1.0, 1.0; length=N_a))
    b = randn(rng, N_b)
    TestUtils.test_internal_abstractgps_interface(rng, f_approx_post, a, b)

    @testset "update_posterior (new observation)" begin
        rng = MersenneTwister(1)
        X = rand(rng, 10)
        y = rand(rng, 10)
        Z = rand(rng, 4)

        f = GP(SqExponentialKernel())

        # online learning
        p_fx1 = posterior(VFE(f(Z)), f(X[1:7], 0.1), y[1:7])
        u_p_fx1 = update_posterior(p_fx1, f(X[8:10], 0.1), y[8:10])

        # batch learning
        p_fx2 = posterior(VFE(f(Z)), f(X, 0.1), y)

        @test inducing_points(u_p_fx1) ≈ inducing_points(p_fx2) atol = 1e-5
        @test u_p_fx1.data.m_ε ≈ p_fx2.data.m_ε atol = 1e-5
        @test u_p_fx1.data.Λ_ε.U ≈ p_fx2.data.Λ_ε.U atol = 1e-5
        @test u_p_fx1.data.U ≈ p_fx2.data.U atol = 1e-5
        @test u_p_fx1.data.α ≈ p_fx2.data.α atol = 1e-5
        @test u_p_fx1.data.b_y ≈ p_fx2.data.b_y atol = 1e-5
        @test u_p_fx1.data.B_εf ≈ p_fx2.data.B_εf atol = 1e-5
        @test u_p_fx1.data.x ≈ p_fx2.data.x atol = 1e-5
        @test u_p_fx1.data.Σy ≈ p_fx2.data.Σy atol = 1e-5
    end

    @testset "update_posterior (new pseudo-points)" begin
        rng = MersenneTwister(1)
        X = rand(rng, 10)
        y = rand(rng, 10)
        Z1 = rand(rng, 4)
        Z2 = rand(rng, 3)
        Z = vcat(Z1, Z2)

        f = GP(SqExponentialKernel())

        # Sequentially adding pseudo points
        p_fx1 = posterior(VFE(f(Z1)), f(X, 0.1), y)
        u_p_fx1 = update_posterior(p_fx1, f(Z2))

        # Adding all pseudo-points at once
        p_fx2 = posterior(VFE(f(Z)), f(X, 0.1), y)

        @test inducing_points(u_p_fx1) ≈ inducing_points(p_fx2) atol = 1e-5
        @test u_p_fx1.data.m_ε ≈ p_fx2.data.m_ε atol = 1e-5
        @test u_p_fx1.data.Λ_ε.U ≈ p_fx2.data.Λ_ε.U atol = 1e-5
        @test u_p_fx1.data.U ≈ p_fx2.data.U atol = 1e-5
        @test u_p_fx1.data.α ≈ p_fx2.data.α atol = 1e-2 rtol = 1e-2
        @test u_p_fx1.data.b_y ≈ p_fx2.data.b_y atol = 1e-5
        @test u_p_fx1.data.B_εf ≈ p_fx2.data.B_εf atol = 1e-5
        @test u_p_fx1.data.x ≈ p_fx2.data.x atol = 1e-5
        @test u_p_fx1.data.Σy ≈ p_fx2.data.Σy atol = 1e-5
    end

    @testset "elbo / dtc" begin
        x = collect(Float64, range(-1.0, 1.0; length=N_cond))
        f = GP(SqExponentialKernel())
        fx = f(x, 0.1)
        y = rand(rng, fx)

        # Ensure that the elbo is close to the logpdf when appropriate.
        @test elbo(VFE(f(x)), fx, y) isa Real
        @test elbo(VFE(f(x)), fx, y) ≈ logpdf(fx, y)
        @test elbo(VFE(f(x .+ randn(rng, N_cond))), fx, y) < logpdf(fx, y)

        # Ensure that the dtc is close to the logpdf when appropriate.
        @test dtc(VFE(f(x)), fx, y) isa Real
        @test dtc(VFE(f(x)), fx, y) ≈ logpdf(fx, y)
    end

    @testset "Type Stability - $T" for T in [Float64, Float32]
        x = collect(T, range(-1.0, 1.0; length=N_cond))
        f = GP(T(0), SqExponentialKernel())
        fx = f(x, T(0.1))
        y = rand(rng, fx)
        jitter = T(1e-12)

        @test elbo(VFE(f(x, jitter)), fx, y) isa T
        @test dtc(VFE(f(x, jitter)), fx, y) isa T

        post = posterior(VFE(f(x, jitter)), fx, y)
        p_fx = post(x, jitter)

        @test rand(rng, p_fx) isa Vector{T}
        @test logpdf(p_fx, y) isa T
    end
end


text:
using Plots, LaTeXStrings, Statistics

function Deposition_2DG(; len, tot_time, dep_rate, color_step)
    surface = InitialSurf(len)
    for n in 0:tot_time
        randsurf = rand(1:len, dep_rate)
        for i in randsurf
            surface = check_and_depose(surface, i, len, (n%color_step) + 1)
        end
    end
    return surface
end

function check_and_depose(surf, index, len, addval) #checking down and sides and first layer
    i1 , i2 = sides(index, len)
    if surf[end][index] != 0
        newsurf = [0 for j = 1:len]
        newsurf[index] = addval
        push!(surf, newsurf)
    end
    for num in 0:length(surf) -1
        if surf[end - num][index] == 0 && (surf[end - num][i1] != 0 || surf[end - num][i2] != 0)
            surf[end - num][index] = addval
            break
        elseif length(surf) > num + 1 && surf[end - 1 - num][index] != 0
            surf[end - num][index] = addval
            break
        end
    end
    return surf
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function InitialSurf(len)
    pos1 = [0 for i=1:len]
    pos1[floor(Int, len/2)] = 1
    return [pos1]
end

Parameters = Dict(:len => 300,
                    :tot_time => 1000,
                        :dep_rate => 50,
                            :color_step => 100)

surf = Deposition_2DG(;Parameters...)

heatmap(transpose(hcat(surf...)), c = :solar, legend = false)
savefig("../../Figs/Q3/Graphic.pdf")


text:
@testset "deltatest" begin
	p=DCEGM.Param(par=Dict(:delta => 0.9))
	@test p.delta == 0.9
end

@testset "basics" begin
	p = DCEGM.Param()
	@testset "run minimal egm" begin
		m = DCEGM.minimal_EGM(p)
		@test true
	end

	@testset "ccp" begin
		x = rand(5,10)
		c = DCEGM.ccp(x,p)
		@test length(c) == 10
		@test minimum(c) >= 0.0
		@test maximum(c) <= 1.0
	end
	@testset "logsum" begin
		x = rand(5,10)
		c = DCEGM.logsum(x,p)
		@test length(c) == 10
		@test maximum(abs,DCEGM.logsum(zeros(2,2),p) .- p.lambda * [log(2) log(2)]) .< 1e-6
		p.lambda = 1
		@test maximum(abs,DCEGM.logsum(zeros(2,2),p) .- p.lambda * [log(2) log(2)]) .< 1e-6
	end
end


text:
# Index of representation element of type ElemT
"""
    Index{T, ElemT}

Index of an element of type `ElemT` in a `Rep{T}`.
"""
struct Index{T, ElemT}
    value::Int
end

# Type of the value associated to this index
valuetype(::Union{Index{T, ElemT}, Type{Index{T, ElemT}}}) where {T, ElemT} = ElemT

const HyperPlaneIndex{T} = Index{T, <:HyperPlane{T}}
const HalfSpaceIndex{T} = Index{T, <:HalfSpace{T}}
const HIndex{T} = Union{HyperPlaneIndex{T}, HalfSpaceIndex{T}}

#const SymPointIndex{T} = Index{T, <:SymPoint{T}}
const PointIndex{T} = Index{T, <:AbstractVector{T}}
#const PIndex{T} = Union{SymPointIndex{T}, PointIndex{T}}
const PIndex{T} = PointIndex{T}
const LineIndex{T} = Index{T, <:Line{T}}
const RayIndex{T} = Index{T, <:Ray{T}}
const RIndex{T} = Union{LineIndex{T}, RayIndex{T}}
const VIndex{T} = Union{PIndex{T}, RIndex{T}}
islin(::Union{Index{T, ElemT}, Type{Index{T, ElemT}}}) where {T, ElemT} = islin(ElemT)
ispoint(::Union{Index{T, ElemT}, Type{Index{T, ElemT}}}) where {T, ElemT} = ispoint(ElemT)

"""
    Indices{T, ElemT, RepT<:Rep{T}}

Iterator over the indices of the elements of type `ElemT` of the field `rep`.
"""
struct Indices{T, ElemT, RepT<:Rep{T}}
    rep::RepT
    function Indices{T, ElemT}(rep) where {T, ElemT}
        new{T, ElemT, typeof(rep)}(rep)
    end
end

Base.eltype(::Indices{T, ElemT}) where {T, ElemT} = Index{T, ElemT}
valuetype(idxs::Indices) = valuetype(eltype(idxs))

const HyperPlaneIndices{T, RepT} = Indices{T, <:HyperPlane{T}, RepT}
const HalfSpaceIndices{T, RepT} = Indices{T, <:HalfSpace{T}, RepT}
const HIndices{T, RepT} = Union{HyperPlaneIndices{T, RepT}, HalfSpaceIndices{T, RepT}}

#const SymPointIndices{T, RepT} = Indices{T, <:SymPoint{T}, RepT}
const PointIndices{T, RepT} = Indices{T, <:AbstractVector{T}, RepT}
#const PIndices{T, RepT} = Union{SymPointIndices{T, RepT}, PointIndices{T, RepT}}
const PIndices{T, RepT} = PointIndices{T, RepT}
const LineIndices{T, RepT} = Indices{T, <:Line{T}, RepT}
const RayIndices{T, RepT} = Indices{T, <:Ray{T}, RepT}
const RIndices{T, RepT} = Union{LineIndices{T, RepT}, RayIndices{T, RepT}}
const VIndices{T, RepT} = Union{PIndices{T, RepT}, RIndices{T, RepT}}

undouble_it(idx::Nothing) = nothing
double_it(idx::Nothing) = nothing
undouble_it(idx::NTuple{2, Index}) = idx[1]
double_it(idx::Index) = idx, idx
function Base.iterate(idxs::Indices)
    return double_it(startindex(idxs))
end
function Base.iterate(idxs::Indices{T, ElemT},
                      idx::Index{T, ElemT}) where {T, ElemT}
    return double_it(nextindex(idxs.rep, idx))
end

# For polyhedron, redirect to hrep or vrep depending on whether it is an
# HRepElement or VRepElement
repfor(p::Polyhedron, ::Type{<:HRepElement}) = hrep(p)
repfor(p::Polyhedron, ::Type{<:VRepElement}) = vrep(p)
function startindex(idxs::Indices{T, ElemT,
                                  <:Polyhedron{T}}) where {T, ElemT}
    return startindex(Indices{T, ElemT}(repfor(idxs.rep, ElemT)))
end
function nextindex(p::Polyhedron{T}, idx::Index{T, ElemT}) where {T, ElemT}
    return nextindex(repfor(p, ElemT), idx)
end
function Base.length(idxs::Indices{T, ElemT, <:Polyhedron{T}}) where {T, ElemT}
    return length(Indices{T, ElemT}(repfor(idxs.rep, ElemT)))
end
function Base.isempty(idxs::Indices{T, ElemT, <:Polyhedron{T}}) where {T, ElemT}
    return isempty(Indices{T, ElemT}(repfor(idxs.rep, ElemT)))
end
function Base.get(p::Polyhedron{T}, idx::Index{T, ElemT}) where {T, ElemT}
    return get(repfor(p, ElemT), idx)
end

"""
The representation `rep` does not contain any `elem`.
"""
macro norepelem(rep, elem)
    idxs = Symbol(string(elem) * "Indices")
    idx = Symbol(string(elem) * "Index")
    quote
        Base.length(idxs::$idxs{T, <:$rep{T}}) where {T} = 0
        Base.isempty(idxs::$idxs{T, <:$rep{T}}) where {T} = true
        Base.iterate(idxs::$idxs{T, <:$rep{T}}) where {T} = nothing
    end
end

abstract type HAffineSpace{T} <: HRepresentation{T} end
@norepelem HAffineSpace HalfSpace

_promote_reptype(P1::Type{<:HAffineSpace}, ::Type{<:HAffineSpace}) = P1
_promote_reptype(P1::Type{<:HAffineSpace}, ::Type{<:HRep}) = hreptype(P1)

abstract type VPolytope{T} <: VRepresentation{T} end
@norepelem VPolytope Line
@norepelem VPolytope Ray

_promote_reptype(P1::Type{<:VPolytope}, ::Type{<:VPolytope}) = P1
_promote_reptype(P1::Type{<:VPolytope}, ::Type{<:VRep}) = vreptype(P1)

abstract type VSymPolytope{T} <: VPolytope{T} end
@norepelem VSymPolytope Point

abstract type VCone{T} <: VRepresentation{T} end
#@norepelem VCone SymPoint
# See issue #28
Base.length(idxs::PointIndices{T, <:VCone{T}}) where {T} = hasallrays(idxs.rep) ? 1 : 0
Base.isempty(idxs::PointIndices{T, <:VCone{T}}) where {T} = !hasallrays(idxs.rep)
Base.get(L::VCone{T}, ::PointIndex{T}) where {T} = origin(vvectortype(typeof(L)), fulldim(L))
startindex(idxs::PointIndices{T, <:VCone{T}}) where {T} = hasallrays(idxs.rep) ? eltype(idxs)(1) : nothing
nextindex(::VCone{T}, idx::PointIndex{T}) where {T} = nothing

_promote_reptype(P1::Type{<:VCone}, ::Type{<:VCone}) = P1
_promote_reptype(P1::Type{<:VCone}, ::Type{<:VRep}) = vreptype(P1)

abstract type VLinearSpace{T} <: VCone{T} end
@norepelem VLinearSpace Ray

_promote_reptype(P1::Type{<:VLinearSpace}, ::Type{<:VLinearSpace}) = P1
_promote_reptype(P1::Type{<:VLinearSpace}, ::Type{<:VCone}) = conetype(P1)
_promote_reptype(P1::Type{<:VLinearSpace}, ::Type{<:VRep}) = vreptype(P1)

"""
The representation `rep` contain the elements `elem` inside a vector in the field `field`.
"""
macro vecrepelem(rep, elem, field)
    idxs = Symbol(string(elem) * "Indices")
    idx = Symbol(string(elem) * "Index")
    esc(quote
        Base.length(idxs::$idxs{T, <:$rep{T}}) where {T} = length(idxs.rep.$field)
        Base.isempty(idxs::$idxs{T, <:$rep{T}}) where {T} = isempty(idxs.rep.$field)
        function Polyhedra.startindex(idxs::$idxs{T, <:$rep{T}}) where {T}
            if isempty(idxs.rep.$field)
                return nothing
            else
                return eltype(idxs)(1)
            end
        end
        Base.get(rep::$rep{T}, idx::$idx{T}) where {T} = rep.$field[idx.value]
        function Polyhedra.nextindex(rep::$rep{T}, idx::$idx{T}) where {T}
            if idx.value >= length(rep.$field)
                return nothing
            else
                return typeof(idx)(idx.value + 1)
            end
        end
    end)
end

"""
The representation `rep` contain the elements `elem` inside a representation in the field `field`.
"""
macro subrepelem(rep, elem, field)
    idxst = Symbol(string(elem) * "Indices")
    idxs = :(Polyhedra.$idxst)
    idxt = Symbol(string(elem) * "Index")
    idx = :(Polyhedra.$idxt)
    subidxs = :(Polyhedra.Indices{T, Polyhedra.valuetype(idxs)}(idxs.rep.$field))
    esc(quote
        Base.length(idxs::$idxs{T, <:$rep{T}}) where {T} = length($subidxs)
        Base.isempty(idxs::$idxs{T, <:$rep{T}}) where {T} = isempty($subidxs)
        Polyhedra.startindex(idxs::$idxs{T, <:$rep{T}}) where {T} = Polyhedra.startindex($subidxs)
        Base.get(rep::$rep{T}, idx::$idx{T}) where {T} = get(rep.$field, idx)
        Polyhedra.nextindex(rep::$rep{T}, idx::$idx{T}) where {T} = Polyhedra.nextindex(rep.$field, idx)
    end)
end


text:
@doc doc"""
compute_ps(Z::Array{Int, 2})

Estimate of Posterior Similarity Matrix

Assuming Z to be [Trails x Observations]

""" ->
function compute_psm(Z::Array{Int, 2})

	(N, M) = size(Z)

	valuerange = collect(1:M)
	inRange = true

	for n in 1:N
		for m in 1:M
				if !(Z[n, m] in valuerange)
					inRange = false
				end
		end
	end

	@assert inRange "Values are not in range: All values must be larger equal 1 and smaller equal "

	psm = zeros(N, N)

	for m in 1:M
		for i in 1:N
			for j in 1:N
				if Z[i, m] == Z[j, m]
					psm[i, j] += 1
				end
			end
		end
	end

	return psm ./ M

end


text:
# download input files from: https://gist.github.com/isaacsas/7648fe065f1884307a2906a0a48ea12b
# then add ReactionNetworkImporter package (unregistered)
#]
# add https://github.com/isaacsas/ReactionNetworkImporter.jl.git

using DiffEqBase, DiffEqBiological, Plots, OrdinaryDiffEq, DataFrames, CSVFiles, LinearAlgebra
using ReactionNetworkImporter

# parameters
doplot = false
networkname = "testbcrbng"
tf = 4*10^4

# BNG simulation data
datadir  = joinpath(@__DIR__,"../data/repressilator")
fname    = joinpath(datadir, "Repressilator.net")
gdatfile = joinpath(datadir, "Repressilator.gdat")
print("getting gdat file...")
gdatdf = DataFrame(load(File(format"CSV", gdatfile), header_exists=true, spacedelim=true))
println("done")

# load the BNG reaction network in DiffEqBio
prnbng = loadrxnetwork(BNGNetwork(),string(networkname,"bng"), fname)
rnbng = prnbng.rn; u₀ = prnbng.u₀; p = prnbng.p;
u0 = u₀
addodes!(rnbng; build_jac=false, build_symfuncs=false)
boprob = ODEProblem(rnbng, u0, (0.,tf), p)

# BNG simulation data testing
pTetRid = prnbng.groupstoids[:pTetR][1]
bngsol = gdatdf[:pTetR]

# note solvers run _much_ faster the second time
bsol = solve(boprob, Tsit5(), abstol=1e-12, reltol=1e-12, saveat=tf/(4*10^2));

if doplot
    plotlyjs()
    plot(gdatdf[:time], gdatdf[:pTetR], label=:BNGPTetR)
    plot!(bsol.t, bsol[pTetRid,:], label=:DEBPTetR)
end

# test the error
@assert all(bsol.t .== gdatdf[:time])
#norm(gdatdf[:pTetR] - bsol[pTetRid,:], Inf)

@assert maximum(abs.(gdatdf[:pTetR][2:end] - bsol[pTetRid,2:end]) ./ abs.(gdatdf[:pTetR][2:end])) < 1e-6

max_issues_repo_name: isaacsas/ParseRxNetwork


text:
### Day 12 - Problem 1 & 2
## Leonardo's Monorail
## Author: Thanasis Georgiou

# A CPU with registers, a program counter and program ROM
type CPU
  gpreg::Dict{String, Int}
  pc::Int
  program::Array{String}
end

function cpy(cpu::CPU, from, to)
  # Check if origin is register or number
  from = get(tryparse(Int, from), String(from))

  # Value to copy
  value = isa(from, String)? cpu.gpreg[from] : from

  # Set target register
  cpu.gpreg[to] = value
end

function inc(cpu::CPU, reg)
  # Increment the register
  cpu.gpreg[reg] += 1
end

function dec(cpu::CPU, reg)
  # Decrement the register
  cpu.gpreg[reg] -= 1
end

function jnz(cpu::CPU, cond, jump)
  # Check if condition is a number or a register
  cond = get(tryparse(Int, cond), String(cond))

  # Value to check
  value = isa(cond, String)? cpu.gpreg[cond] : cond

  # Jump if the value is not zero
  if value != 0
    cpu.pc += parse(Int, jump)
  end
end

# Go forward one clock cycle
function clock(cpu::CPU)::Bool
  # Get current instruction from program memory
  instr = split(cpu.program[cpu.pc], ' ')

  # Store program counter
  pc = cpu.pc

  # Find out which instruction it is and execute it
  if instr[1] == "cpy"
    cpy(cpu, instr[2], instr[3])
  elseif instr[1] == "inc"
    inc(cpu, instr[2])
  elseif instr[1] == "dec"
    dec(cpu, instr[2])
  elseif instr[1] == "jnz"
    jnz(cpu, instr[2], instr[3])
  else
    error("Invalid opcode $(instr[1])! Halting...")
  end

  # If the program counter didn't change by some instruction, we should move it to the next line in ROM
  if cpu.pc == pc
    cpu.pc += 1
  end

  # Check if there are more instructions to execute
  return cpu.pc <= length(cpu.program)
end

function createCPU(rom::Array{String})::CPU
  # Initialize registers
  registers = Dict{String, Int}()
  registers["a"] = 0
  registers["b"] = 0
  registers["c"] = 0
  registers["d"] = 0

  # Create CPU
  return CPU(registers, 1, rom)
end

# Read input from file
f = open("cpu.asm")
rom = map(line -> chomp(line), readlines(f))
close(f)

# Create CPU
cpu = createCPU(rom)

# Problem 2: Initialize register c
cpu.gpreg["c"] = 1

# Execute program
i = 0
while clock(cpu)
  i += 1
end

# Print cpu state
println("Execution finished in $i clocks.")
render(cpu)


text:
using BlockMatching
using Distances
using ImageCore, ImageTransformations, TestImages, TestImages, ImageDistances
using OffsetArrays
using Test

CUDA_INSTALLED = false
try
    global CUDA_INSTALLED
    # This errors with `IOError` when nvidia driver is not available,
    # in which case we don't even need to try `using CUDA`
    run(pipeline(`nvidia-smi`, stdout=devnull, stderr=devnull))
    push!(LOAD_PATH, "@v#.#") # force using global CUDA installation

    @eval using CUDA
    CUDA.allowscalar(false)
    CUDA_INSTALLED = true
catch e
    e isa Base.IOError || @warn e LOAD_PATH
end
CUDA_FUNCTIONAL = CUDA_INSTALLED && CUDA.functional()
if CUDA_FUNCTIONAL
    @info "CUDA test: enabled"
else
    @warn "CUDA test: disabled"
end

@testset "BlockMatching.jl" begin
    include("utils.jl")
    include("strategies/full_search.jl")

    if CUDA_FUNCTIONAL
        include("strategies/full_search_cuda.jl")
    end
end


text:
using Pkg
Pkg.activate(".")

using DelimitedFiles

print("Using $blas_num_threads BLAS thread")
blas_num_threads > 1 ? println("s") : println()

println("Using $omp_num_threads OpenBLAS threads")
println("Maximum DMRG bond dimension is set to $maxdim.")

seperator = "#"^70

Base.run(`make`)

N = length(omp_num_threads)
for j in 1:N
  omp_num_thread = omp_num_threads[j]
  println(seperator)
  println("Run C++ benchmark in path $(@__DIR__) with $omp_num_thread OpenMP threads.")
  println()

  open("run.sh", "w") do io
    write(io, """#!/bin/bash
                 MKL_NUM_THREADS=$blas_num_threads OMP_NUM_THREADS=$omp_num_thread ./run $maxdim""")
  end
  chmod("run.sh", 0o777)
  time = @elapsed Base.run(`./run.sh`)
  rm("run.sh")
  @show time
  println()

  # Printing results
  # TODO: add version number to data file name
  filename = joinpath(@__DIR__, "data", "data_maxdim_$(maxdim)_omp_num_threads_$(omp_num_thread).txt")
  println("Writing results to $filename")
  mkpath(dirname(filename))
  writedlm(filename, time)
end


text:
# Tests for Cost Additive DEA Models
@testset "CostAdditiveDEAModel" begin

    # Test using Book data
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];
    Y = [1; 1; 1; 1; 1; 1; 1; 1];
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # With :Ones
    deacostaddones = deacostadd(X, Y, W, :Ones)

    @test typeof(deacostaddones) == CostAdditiveDEAModel

    @test nobs(deacostaddones) == 8
    @test ninputs(deacostaddones) == 2
    @test noutputs(deacostaddones) == 1
    @test ismonetary(deacostaddones) == false

    @test efficiency(deacostaddones) == efficiency(deacostaddones, :Economic)
    @test efficiency(deacostaddones, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 6.0; 3.0; 3.0; 5.6]
    @test efficiency(deacostaddones, :Technical)  ≈ [0.0; 0.0; 0.0; 3.0; 6.0; 2.0; 3.0; 5.2]
    @test efficiency(deacostaddones, :Allocative) ≈ [0.0; 1.0; 1.0; 0.0; 0.0; 1.0; 0.0; 0.4]
    @test normfactor(deacostaddones) ≈ [1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]

    # With :Ones and CRS (Same result in this example)
    deacostaddonescrs = deacostadd(X, Y, W, :Ones, rts = :CRS)

    @test efficiency(deacostaddonescrs, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 6.0; 3.0; 3.0; 5.6]
    @test efficiency(deacostaddonescrs, :Technical)  ≈ [0.0; 0.0; 0.0; 3.0; 6.0; 2.0; 3.0; 5.2]
    @test efficiency(deacostaddonescrs, :Allocative) ≈ [0.0; 1.0; 1.0; 0.0; 0.0; 1.0; 0.0; 0.4]

    # Default is :Ones
    @test efficiency(deacostadd(X, Y, W)) == efficiency(deacostaddones)

    # With :MIP
    deacostaddmip = deacostadd(X, Y, W, :MIP)

    @test efficiency(deacostaddmip, :Economic)   ≈ [0.0; 1.0; 1.0; 1.0; 1.2; 3.0; 1.5; 3.5] atol = 1e-5
    @test efficiency(deacostaddmip, :Technical)  ≈ [0.0; 0.0; 0.0; 5/6; 1.2; 1/3; 7/10; 0.875] atol = 1e-5
    @test efficiency(deacostaddmip, :Allocative) ≈ [0.0; 1.0; 1.0; 1/6; 0.0; 8/3; 4/5; 21/8] atol = 1e-5
    @test normfactor(deacostaddmip) ≈ [2.0; 1.0; 1.0; 3.0; 5.0; 1.0; 2.0; 1.6]

    # With :MIP and monetary = true
    deacostaddmipmon = deacostadd(X, Y, W, :MIP, monetary = true)

    @test ismonetary(deacostaddmipmon) == true

    @test efficiency(deacostaddmipmon) ≈ efficiency(deacostaddmip) .* normfactor(deacostaddmip)
    @test efficiency(deacostaddmipmon, :Technical) ≈ efficiency(deacostaddmip, :Technical) .* normfactor(deacostaddmip)
    @test efficiency(deacostaddmipmon, :Allocative) ≈ efficiency(deacostaddmip, :Allocative) .* normfactor(deacostaddmip)

    # With custom weights
    deacostaddcustom = deacostadd(X, Y, W, rhoX = 1 ./ X)

    @test efficiency(deacostaddcustom) ≈ efficiency(deacostaddmip)

    # Print
    show(IOBuffer(), deacostaddones)
    show(IOBuffer(), deacostadd(X, Y, W, :Ones, dispos = :Weak))

    # Test errors
    @test_throws DimensionMismatch deacostadd([1; 2 ; 3], [4 ; 5], [1; 1; 1] , :Ones) #  Different number of observations
    @test_throws DimensionMismatch deacostadd([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], :Ones) # Different number of observation in prices
    @test_throws DimensionMismatch deacostadd([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], :Ones) # Different number of input prices and inputs
    @test_throws DimensionMismatch deacostadd([1; 2; 3] , [4; 5; 6], [1; 1; 1], rhoX = [1 1; 2 2; 3 3]) # Different size of inputs weights
    @test_throws ArgumentError deacostadd([1; 2; 3], [4; 5; 6], [1; 2; 3],:Ones, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deacostadd([1; 2; 3], [1; 2; 3], [1; 1; 1], :Error) # Invalid inuts weights
    @test_throws ArgumentError deacostadd([1; 2; 3], [1; 2; 3], [1; 1; 1], :Ones, dispos = :Error) # Invalid inputs disposal
    @test_throws ArgumentError deacostadd([1; 2; 3], [1; 2; 3], [1; 1; 1], :Ones, rhoX= [1; 1; 1]) # Weights not allowed if model != :Custom

end


text:
"""
    OptionsMod

Sub-module of ModelBaseEcon, although it can be used independently.
Implements the [`Options`](@ref) data structure.

## Contents
  * [`Options`](@ref)
  * [`getoption`](@ref) - read the value of an option
  * [`getoption!`](@ref) - if not present, also create an option
  * [`setoption!`](@ref) - create or update the value of an option

"""
module OptionsMod

export Options, getoption, getoption!, setoption!

"""
    Options

A collection of key-value pairs representing the options controlling the
behaviour or the definition of a Model object. The key is the option name and is
always a Symbol, or converted to Symbol, while the value can be anything.

The options can be accessed using dot notation. Functions [`getoption`](@ref)
and [`setoption!`](@ref) are also provided. They can be used for programmatic
processing of options as well as when the option name is not a valid Julia
identifier.

See also: [`Options`](@ref), [`getoption`](@ref), [`getoption!`](@ref),
[`setoption!`](@ref)

# Examples
```jldoctest
julia> o = Options(maxiter=20, tol=1e-7)
Options:
    maxiter=20
    tol=1.0e-7

julia> o.maxiter = 25
25

julia> o
Options:
    maxiter=25
    tol=1.0e-7

```
"""
struct Options
    contents::Dict{Symbol,Any}
    # Options() = new(Dict())
    Options(c::Dict{Symbol,<:Any}) = new(c)
end

############
# Constructors

"""
    Options(key=value, ...)
    Options(:key=>value, ...)

Construct an Options instance with key-value pairs given as keyword arguments or
as a list of pairs. If the latter is used, each key must be a `Symbol`.
"""
Options(; kwargs...) = Options(Dict{Symbol,Any}(kwargs))
# Options(pair::Pair{Symbol, T}) where T = Options(Dict(pair))
Options(pairs::Pair{Symbol,<:Any}...) = Options(Dict(pairs...))
Options(pairs::Pair{<:AbstractString,<:Any}...) = Options(Dict(Symbol(k) => v for (k, v) in pairs))

"""
    Options(::Options)

Construct an Options instance as an exact copy of an existing instance.
"""
Options(opts::Options) = Options(deepcopy(Dict(opts.contents)))

############
# compare

Base.:(==)(opts::Options, opts2::Options) = opts.contents == opts2.contents
Base.:(==)(opts::Dict, opts2::Options) = opts == opts2.contents
Base.:(==)(opts::Options, opts2::Dict) = opts.contents == opts2

############
# merge

"""
    merge(o1::Options, o2::Options, ...)

Merge the given Options instances into a new Options instance.
If the same option key exists in more than one instance, keep the value from
the last one.
"""
Base.merge(o1::Options, o2::Options...) = Options(merge(o1.contents, (o.contents for o in o2)...))

"""
    merge!(o1::Options, o2::Options...)

Update the first argument, adding all options from the remaining arguments. If the same
option exists in multiple places, use the last one.
"""
Base.merge!(o1::Options, o2::Options...) = (merge!(o1.contents, (o.contents for o in o2)...); o1)

############
# Access by dot notation

Base.propertynames(opts::Options) = tuple(keys(opts.contents)...)

Base.setproperty!(opts::Options, name::Symbol, val) = opts.contents[name] = val

Base.getproperty(opts::Options, name::Symbol) =
    name ∈ fieldnames(Options) ? getfield(opts, name) :
    name ∈ keys(opts.contents) ? opts.contents[name]  :
                                 error("option $name not set.");

############
# Pretty printing

function Base.show(io::IO, opts::Options)
    print(io, "Options(")
    str = String[
        sprint(print, on, "=", ov, context=io, sizehint=0)
        for (on, ov) in pairs(opts)
    ]
    print(io, join(str, ", "))
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", opts::Options)
    recur_io = IOContext(io, :SHOWN_SET => opts.contents,
                             :typeinfo => eltype(opts.contents),
                             :compact => get(io, :compact, true))
    print(io, length(opts.contents), " Options:")
    if !isempty(opts)
        for (key, value) in opts.contents
            print(io, "\n    ", key, " = ")
            show(recur_io, value)
        end
    end
    # print(io, "\n")
end
export getoption!, getoption, setoption!

############
# Iteration

Base.iterate(opts::Options) = iterate(opts.contents)
Base.iterate(opts::Options, state) = iterate(opts.contents, state)

############
# getoption, and setoption

"""
    getoption(o::Options; name=default [, name=default, ...])
    getoption(o::Options, name, default)

Retrieve the value of an option or a set of options.  The provided defaults
are used when the option doesn't exit.

The return value is the value of the option requested or, if the option doesn't
exist, the default. In the first version of the function, if there are more than
one options requested, the return value is a tuple.

In the second version, the name could be a symbol or a string, which can be helpful
if the name of the option is not a valid identifier.
"""
function getoption end
function getoption(opts::Options; kwargs...)
    if length(kwargs) == 1
        return get(opts.contents, first(kwargs)...)
    else
        return tuple((get(opts.contents, kv...) for kv in kwargs)...)
    end
end
getoption(opts::Options, name::Symbol, default) = get(opts.contents, name, default)
getoption(opts::Options, name::S where S <: AbstractString, default) = get(opts.contents, Symbol(name), default)

"""
    getoption!(o::Options; name=default [, name=default, ...])
    getoption!(o::Options, name, default)

Retrieve the value of an option or a set of options. If the name does not match
an existing option, the Options instance is updated by inserting the given name
and default value.

The return value is the value of the option requested (or the default). In the
first version of the function, if there are more than one options requested, the
return value is a tuple.

In the second version, the name could be a symbol or a string, which can be
helpful if the name of the option is not a valid identifier.
"""
function getoption! end
function getoption!(opts::Options; kwargs...)
    if length(kwargs) == 1
        return get!(opts.contents, first(kwargs)...)
    else
        return tuple((get!(opts.contents, kv...) for kv in kwargs)...)
    end
end
getoption!(opts::Options, name::Symbol, default) = get!(opts.contents, name, default)
getoption!(opts::Options, name::S where S <: AbstractString, default) = get!(opts.contents, Symbol(name), default)

"""
    setoption!(o::Options; name=default [, name=default, ...])
    setoption!(o::Options, name, default)

Retrieve the value of an option or a set of options. If the name does not match
an existing option, the Options instance is updated by inserting the given name
and default value.

The return value is the value of the option requested (or the default). In the
first version of the function, if there are more than one options requested, the
return value is a tuple.

In the second version, the name could be a symbol or a string, which can be
helpful if the name of the option is not a valid identifier.
"""
function setoption! end
setoption!(opts::Options; kwargs...) = (push!(opts.contents, kwargs...); opts)
setoption!(opts::Options, name::S where S <: AbstractString, value) = (push!(opts.contents, Symbol(name) => value); opts)
setoption!(opts::Options, name::Symbol, value) = (push!(opts.contents, name => value); opts)

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

Base.in(name, o::Options) = Symbol(name) ∈ keys(o.contents)

Base.keys(o::Options) = keys(o.contents)
Base.values(o::Options) = values(o.contents)
end # module

using .OptionsMod
export Options, getoption, getoption!, setoption!


text:
################################################################################
################################### Utilities ##################################
################################################################################

export
	BFSIterator,
	SortedBFSIterator,
	update_iterator!,
	depth_first_search,
	node_iterator,
	filtered_topological_sort,
    get_nodes_of_type,
    transplant!,
    forward_pass!,
	backward_pass!,
	backup_descendants,
	backup_ancestors,
	get_biggest_tree,
	collect_subtree,
	collect_descendants,
	collect_ancestors,

    validate_edge,

    matches_template,
	required_predecessors,
	required_successors,
	num_required_predecessors,
	num_required_successors,
    eligible_successors,
    eligible_predecessors,
	num_eligible_predecessors,
	num_eligible_successors

"""
	matches_template(template,node)

Checks if a candidate `node` satisfies the criteria encoded by `template`.
"""
matches_template(template::Type{T},node::Type{S}) where {T,S} = S<:T
matches_template(template::Type{T},node::S) where {T,S} = S<:T
matches_template(template,node) = matches_template(typeof(template),node)
matches_template(template::Tuple,node) = any(map(t->matches_template(t,node), template))

matches_template(template::Type{T},n::CustomNode) where {T} = matches_template(template,node_val(n))
# matches_template(template::Type{T},n::CustomNode{N,ID}) where {T,N,ID} = matches_template(template,N)

abstract type AbstractBFSIterator end
Base.IteratorSize(::AbstractBFSIterator) = Base.SizeUnknown()
Base.IteratorEltype(::AbstractBFSIterator) = Base.HasEltype()
Base.eltype(::AbstractBFSIterator) = Int
function Base.iterate(iter::AbstractBFSIterator,v=nothing)
	if !(v === nothing)
		update_iterator!(iter,v)
	end
	if isempty(iter)
		return nothing
	end
	vp = pop!(iter)
	return vp, vp
end

function _indicator_vec(N,idxs)
	vec = zeros(Bool,N)
	for i in idxs
		vec[i] = true
	end
	vec
end
@with_kw struct BFSIterator{G} <: AbstractBFSIterator
	graph::G 				= DiGraph()
	frontier::Set{Int} 		= get_all_root_nodes(graph)
	next_frontier::Set{Int}	= Set{Int}()
	explored::Vector{Bool}  = _indicator_vec(nv(graph),frontier)
	replace::Bool			= false # if true, allow nodes to reused
end
BFSIterator(graph) = BFSIterator(graph=graph)
BFSIterator(graph,frontier) = BFSIterator(graph=graph,frontier=frontier)
Base.pop!(iter::BFSIterator) = pop!(iter.frontier)
Base.isempty(iter::BFSIterator) = isempty(iter.frontier)
function update_iterator!(iter::BFSIterator,v)
	iter.explored[v] = true
	for vp in outneighbors(iter.graph,v)
		if iter.replace || !iter.explored[vp]
			push!(iter.next_frontier,vp)
			iter.explored[vp] = true
		end
	end
	if isempty(iter.frontier)
		union!(iter.frontier, iter.next_frontier)
		empty!(iter.next_frontier)
	end
end

@with_kw struct SortedBFSIterator{G} <: AbstractBFSIterator
	graph::G 				= DiGraph()
	frontier::Vector{Int}   = sort(collect(get_all_root_nodes(graph)))
	next_frontier::Vector{Int} = Vector{Int}()
	explored::Vector{Bool}  = _indicator_vec(nv(graph),frontier)
	replace::Bool			= false # if true, allow nodes to reused
end
SortedBFSIterator(graph) = SortedBFSIterator(graph=graph)
SortedBFSIterator(graph,frontier) = SortedBFSIterator(graph=graph,frontier=frontier)
function Base.pop!(iter::SortedBFSIterator)
	v = popfirst!(iter.frontier)
end
Base.isempty(iter::SortedBFSIterator) = isempty(iter.frontier)
function update_iterator!(iter::SortedBFSIterator,v)
	iter.explored[v] = true
	for vp in outneighbors(iter.graph,v)
		if iter.replace || !iter.explored[vp]
			push!(iter.frontier,vp)
			iter.explored[vp] = true
		end
	end
end

"""
	depth_first_search(graph,v,goal_function,expand_function,
		neighbor_function=outneighbors)

Returns the first vertex satisfying goal_function(graph,v). Only expands v if
expand_function(graph,v) == true.
"""
function depth_first_search(graph,v,
		goal_function,
		expand_function,
		neighbor_function=outneighbors,
		explored=falses(nv(graph));
		skip_first=false,
		)
	if goal_function(v)
		if !(skip_first && sum(explored) == 0)
			return v
		end
	end
	explored[v] = true
	if expand_function(v)
		for vp in neighbor_function(graph,v)
			if !explored[vp]
				u = depth_first_search(graph,vp,
					goal_function,
					expand_function,
					neighbor_function,
					explored)
				if has_vertex(graph,u)
					return u
				end
			end
		end
	end
	return -1
end

"""
	dfs_check_cycle(graph,v,neighbor_func=inneighbors)

Check if `graph` contains a cycle with `v` in it
"""
function dfs_check_cycle(graph,v,neighbor_func=inneighbors)
	vtx = depth_first_search(graph,v,
		vtx->vtx==v,
		vtx->true,
		neighbor_func,
		falses(nv(graph)),
		skip_first=true,
		)
	has_vertex(graph,vtx)
end

"""
	has_path(graph,v,v2)

Return true if `graph` has a path from `v` to `v2`
"""
function has_path(graph,v,v2)
	vtx = depth_first_search(graph,v,vtx->vtx==v2,vtx->true,outneighbors)
	return vtx > 0
end

"""
	node_iterator(graph,it)

Wraps an iterator over ids or vertices to return the corresponding node at each
iteration.
"""
node_iterator(graph,it) = transform_iter(v->get_node(graph,v),it)

"""
	filtered_topological_sort(graph,template)

Iterator over nodes that match template.
"""
function filtered_topological_sort(graph,template)
	Base.Iterators.filter(n->matches_template(template,n),
		# transform_iter(v->get_node(graph,v), topological_sort_by_dfs(graph))
		node_iterator(graph, topological_sort_by_dfs(graph))
		)
end

"""
    transplant!(graph,old_graph,id)

Share node with `id` in `old_graph` to `graph` with the same id.
"""
function transplant!(graph,old_graph,id)
    add_node!(graph,get_node(old_graph,id),id)
end

# get_nodes_of_type(g::AbstractCustomNGraph,T) = Dict(id=>get_node(g, id) for id in get_vtx_ids(g) if isa(id,T))
# get_nodes_of_type(g::AbstractCustomNGraph,T) = Dict(id=>get_node(g, id) for id in get_vtx_ids(g) if matches_template(T,id))
# get_nodes_of_type(g::AbstractCustomNGraph,T) = Dict(get_vtx_id(g,v)=>node for (v,node) in enumerate(get_nodes(g)) if matches_template(T,id))
function get_nodes_of_type(g::AbstractCustomNGraph,T)
	d = Dict()
	for (v,node) in enumerate(get_nodes(g))
		id = get_vtx_id(g,v)
		if matches_template(T,id) || matches_template(T,node)
			d[id] = node
		end
	end
	return d
end

function forward_pass!(g::AbstractCustomNGraph,init_function,update_function)
    init_function(g)
    for v in topological_sort_by_dfs(g)
        update_function(g,v)
    end
    return g
end
function backward_pass!(g::AbstractCustomNGraph,init_function,update_function)
    init_function(g)
    for v in reverse(topological_sort_by_dfs(g))
        update_function(g,v)
    end
    return g
end

"""
    backup_descendants(g::AbstractCustomNGraph{G,N,ID},template)

Return a dictionary mapping each node's id to the id of it's closest descendant
matching `template`.
"""
backup_descendants(g,f) = _backup_descendants(g,f,get_graph(g))
backup_ancestors(g,f) = _backup_descendants(g,f,reverse(get_graph(g)))
function _backup_descendants(g::AbstractCustomNGraph{G,N,ID},f,
		graph=get_graph(g),
		) where {G,N,ID}
    descendant_map = Dict{ID,Union{ID,Nothing}}()
    for v in reverse(topological_sort_by_dfs(graph))
		node = get_node(g,v)
        if f(node)
            descendant_map[node_id(node)] = node_id(node)
		elseif is_terminal_node(graph,v)
            descendant_map[node_id(node)] = nothing
		else
			descendant_map[node_id(node)] = nothing
			for vp in outneighbors(graph,v)
				id = get!(descendant_map,get_vtx_id(g,vp),nothing)
				if !(id === nothing)
					descendant_map[node_id(node)] = id
				end
            end
        end
    end
    descendant_map
end
backup_multiple_descendants(g,f) = _backup_multiple_descendants(g,f,get_graph(g))
backup_multiple_ancestors(g,f) = _backup_multiple_descendants(g,f,reverse(get_graph(g)))
function _backup_multiple_descendants(g::AbstractCustomNGraph{G,N,ID},f,
		graph=get_graph(g),
		) where {G,N,ID}
    descendant_map = Dict{ID,Set{ID}}()
    for v in reverse(topological_sort_by_dfs(graph))
		node = get_node(g,v)
		dset = descendant_map[node_id(node)] = Set{ID}()
        if f(node)
			push!(dset,node_id(node))
		elseif is_terminal_node(graph,v)
		else
			for vp in outneighbors(graph,v)
				union!(dset,descendant_map[get_vtx_id(g,vp)])
            end
        end
    end
    descendant_map
end

forward_cover(g,vtxs,f) = _forward_cover(g,vtxs,f,get_graph(g))
backward_cover(g,vtxs,f) = _forward_cover(g,vtxs,f,reverse(get_graph(g)))
function _forward_cover(g,vtxs,f,graph=get_graph(g))
	iter = SortedBFSIterator(graph,Vector(unique(vtxs)))
	cover = Set{_node_type(g)}()
	while !isempty(iter)
		v = pop!(iter)
		n = get_node(g,v)
		if f(n) && !(v in vtxs)
			push!(cover,n)
		else
			update_iterator!(iter,v)
		end
	end
	cover
end
_forward_cover(g,id::AbstractID,args...) = _forward_cover(g,get_vtx(g,id),args...)

"""
	contract_by_node_type(g,f)

Constructs a new graph of the same type as `g`, whose nodes are formed by those
for which `f(n) == true`. Edges are contracted so that, e.g.
	`n1 => n2 => n3 => n4` collapses to `n1 => n4`
	if  `f(n1) && !(f(n2)) && !(f(n3)) && f(n4) `.
"""
function contract_by_predicate(g,f)
    tg = typeof(g)()
    descendant_map = backup_multiple_descendants(g,f)
    ancestor_map = backup_multiple_ancestors(g,f)
    for n in get_nodes(g)
        if f(n)
            add_node!(tg,n,node_id(n))
        end
    end
    for n in get_nodes(g)
        if !f(n)
			for a in ancestor_map[node_id(n)]
				for d in descendant_map[node_id(n)]
					if !has_edge(tg,a,d)
						add_edge!(tg,a,d)
					end
				end
			end
        end
    end
    tg
end

"""
	get_biggest_tree(graph,dir=:in)

Return the root/terminal vertex corresponding to the root of the largest tree in
the graph.
"""
function get_biggest_tree(graph,dir=:in)
	if dir==:in
    	leaves = collect(get_all_terminal_nodes(graph))
	else
    	leaves = collect(get_all_root_nodes(graph))
	end
	v = argmax(map(v->ne(bfs_tree(graph,v;dir=dir)),leaves))
	leaves[v]
end

"""
	collect_subtree(graph,v,dir=:out)

Return a set of all nodes in the subtree of `graph` starting from `v` in
direction `dir`.
"""
function collect_subtree(graph,v,dir=:out,keep=true)
	descendants = Set{Int}()
	for e in edges(bfs_tree(graph,v;dir=dir))
		push!(descendants,e.dst)
	end
	if keep
		push!(descendants,v)
	end
	descendants
end
collect_subtree(graph,vec::AbstractVector{Int},args...) = collect_subtree(graph,Set{Int}(vec),args...)
function collect_subtree(graph,starts::Set{Int},dir=:out,keep=true)
	frontier = Set{Int}(starts)
	explored = Set{Int}()
	f = dir == :out ? outneighbors : inneighbors
	while !isempty(frontier)
		v = pop!(frontier)
		push!(explored,v)
		for vp in f(graph,v)
			if !(vp in explored)
				push!(frontier,vp)
			end
		end
	end
	if !(keep == true)
		setdiff!(explored,starts)
	end
	return explored
end

"""
	collect_descendants(graph,v) = collect_subtree(graph,v,:out)
"""
collect_descendants(graph,v,keep=false) = collect_subtree(graph,v,:out,keep)
"""
	collect_ancestors(graph,v) = collect_subtree(graph,v,:in)
"""
collect_ancestors(graph,v,keep=false) = collect_subtree(graph,v,:in,keep)

"""
	required_predecessors(node)

Identifies the types (and how many) of required predecessors to `node`
Return type: `Dict{DataType,Int}`
"""
function required_predecessors end

"""
	required_successors(node)

Identifies the types (and how many) of required successors to `node`
Return type: `Dict{DataType,Int}`
"""
function required_successors end

"""
	eligible_predecessors(node)

Identifies the types (and how many) of eligible predecessors to `node`
Return type: `Dict{DataType,Int}`
"""
function eligible_predecessors end

"""
	eligible_successors(node)

Identifies the types (and how many) of eligible successors to `node`
Return type: `Dict{DataType,Int}`
"""
function eligible_successors end

"""
	num_required_predecessors(node)

Returns the total number of required predecessors to `node`.
"""
function num_required_predecessors(node)
	n = 1
	for (key,val) in required_predecessors(node)
		n += val
	end
	n
end

"""
	num_required_successors(node)

Returns the total number of required successors to `node`.
"""
function num_required_successors(node)
	n = 1
	for (key,val) in required_successors(node)
		n += val
	end
	n
end

"""
	num_eligible_predecessors(node)

Returns the total number of eligible predecessors to `node`.
"""
function num_eligible_predecessors(node)
	n = 1
	for (key,val) in eligible_predecessors(node)
		n += val
	end
	n
end

"""
	num_eligible_successors(node)

Returns the total number of eligible successors to `node`.
"""
function num_eligible_successors(node)
	n = 1
	for (key,val) in eligible_successors(node)
		n += val
	end
	n
end

illegal_edge_msg(n1,n2) = "Edges from `$(typeof(n1))` → `$(typeof(n2))` are illegal"
illegal_indegree_msg(n,val,lo,hi) = "`$(typeof(n))` nodes should $lo ≦ indegree(n) ≦ $hi, but indegree(n) = val"
illegal_outdegree_msg(n,val,lo,hi) = "`$(typeof(n))` nodes should $lo ≦ outdegree(n) ≦ $hi, but indegree(n) = val"

"""
	validate_edge(n1,n2)

For an edge (n1) --> (n2), checks whether the edge is legal and the nodes
"agree".
"""
# function validate_edge(n1,n2)
#     return false
# end
function validate_edge(a,b)
    valid = false
    for (key,val) in eligible_successors(a)
        if matches_template(key,b) && val >= 1
            valid = true
        end
    end
    for (key,val) in eligible_predecessors(b)
        if matches_template(key,a) && val >= 1
            valid = valid && true
            return valid
        end
    end
    return false
end

validate_edge(n1::CustomNode,n2) = validate_edge(node_val(n1),n2)

validate_edge(n1::CustomNode,n2::CustomNode) = validate_edge(n1,node_val(n2))
for op in [:required_successors,:required_predecessors,:eligible_successors,:eligible_predecessors]
	@eval $op(n::CustomNode) = $op(node_val(n))
end

# extending the validation interface to allow dispatch on graph type
for op in [
	:required_successors,
	:required_predecessors,
	:eligible_successors,
	:eligible_predecessors,
	:num_required_successors,
	:num_required_predecessors,
	:num_eligible_successors,
	:num_eligible_predecessors,
	]
	@eval $op(graph,node) = $op(node)
end
for op in [
	:validate_edge,
	]
	@eval $op(graph,n1,n2) = $op(n1,n2)
end

# indegree_bounds(n,g,v) = (0,0)
# outdegree_bounds(n,g,v) = (0,0)
# function validate_indegree(n,g,v)
#     lo,hi = indegree_bounds(n,g,v)
#     @assert (lo <= indegree(g,v) <= hi) illegal_indegree_msg(n,indegree(g,v),lo,hi)
# end
# function validate_outdegree(n,g,v)
#     lo,hi = outdegree_bounds(n,g,v)
#     @assert (lo <= outdegree(g,v) <= hi) illegal_outdegree_msg(n,outdegree(g,v),lo,hi)
# end

function validate_neighborhood(g,v)
	n = get_node(g,v)
	try
		for (d,list,required,eligible) in [
				(:out,outneighbors(g,v),required_successors(g,n),eligible_successors(g,n)),
				(:in,inneighbors(g,v),required_predecessors(g,n),eligible_predecessors(g,n)),
			]
			for vp in list
				np = get_node(g,vp)
				has_match = false
				for k in keys(required)
					if matches_template(k,np)
						required[k] -= 1
						@assert required[k] >= 0 "Node $v has too many $(string(d))neighbors of type $k"
						has_match = true
						break
					end
				end
				for k in keys(eligible)
					if matches_template(k,np)
						eligible[k] -= 1
						@assert eligible[k] >= 0 "Node $v has too many $(string(d))neighbors of type $k"
						has_match = true
						break
					end
				end
				@assert has_match "Node $vp should not be an $(string(d))neighbor of node $v"
			end
		end
	catch e
		if isa(e,AssertionError)
			bt = catch_backtrace()
            showerror(stdout,e,bt)
		else
			rethrow(e)
		end
		return false
	end
	return true
end

function validate_graph(g::AbstractCustomGraph)
    try
        for e in edges(g)
            node1 = get_node(g,e.src)
            node2 = get_node(g,e.dst)
            @assert(validate_edge(g,node1,node2), string(" INVALID EDGE: ", string(node1), " --> ",string(node2)))
        end
        for v in vertices(g)
			if !validate_neighborhood(g,v)
				return false
			end
        end
    catch e
        if typeof(e) <: AssertionError
            bt = catch_backtrace()
            showerror(stdout,e,bt)
        else
            rethrow(e)
        end
        return false
    end
    return true
end

function print_tree_level(io,tree,v,start,f=summary,spacing=" ")
    println(io,start,f(get_node(tree,v)))
    for vp in outneighbors(tree,v)
        print_tree_level(io,tree,vp,string(start,spacing),f,spacing)
    end
end
function Base.print(io::IO,tree::AbstractCustomTree,f=summary,spacing=" ")
    @assert !is_cyclic(tree)
    println(io,typeof(tree))
    for v in get_all_root_nodes(tree)
        print_tree_level(io,tree,v,"",f,spacing)
    end
end

function log_graph_edges(graph,v0=-1;
        show_upstream = true,
        show_downstream = true,
        show_all = true,
    )
    if has_vertex(graph,v0)
        if show_upstream
            @info "upstream from $(summary(get_node(graph,v0))):"
			bfs_graph = reverse(bfs_tree(graph,v0;dir=:in))
			for v in topological_sort_by_dfs(bfs_graph)
				for vp in outneighbors(bfs_graph,v)
					e = Edge(v,vp)
    	        	@info "$(summary(get_node(graph,e.src))) => $(summary(get_node(graph,e.dst))),"
				end
			end
        end
        if show_downstream
            @info "downstream from $(summary(get_node(graph,v0))):"
			bfs_graph = bfs_tree(graph,v0;dir=:out)
			for v in topological_sort_by_dfs(bfs_graph)
				for vp in outneighbors(bfs_graph,v)
					e = Edge(v,vp)
    	        	@info "$(summary(get_node(graph,e.src))) => $(summary(get_node(graph,e.dst))),"
				end
			end
        end
    end
    if show_all
		for v in topological_sort_by_dfs(graph)
			for vp in outneighbors(graph,v)
				e = Edge(v,vp)
            	@info "$(summary(get_node(graph,e.src))) => $(summary(get_node(graph,e.dst))),"
			end
		end
        # for e in edges(graph)
        #     @info "$(summary(get_node(graph,e.src))) => $(summary(get_node(graph,e.dst))),"
        # end
    end
end


text:
using LinearAlgebra, SparseArrays, Printf
function GenerateGS(Node, Elem, Lvl::Int, ColTol::Float64)

	# connectivity matrix
	Nn, Ne, Nne = size(Node,1), size(Elem,1), size(Elem,2)

	A1 = spzeros(Bool, Nn, Nn)
	for i=1:Ne; A1[Elem[i,:],Elem[i,:]] .= true; end
	A1 = A1 .* .~I(Nn); An = A1

	# Level 1 connectivity
	IdxJ, IdxI = findnz(An)
	Bars = [IdxI IdxJ]

	# Bar length and normalized direction
	D = Node[IdxI,:] - Node[IdxJ,:];
	D = D./sqrt.(sum(D.^2, dims=2))

	# Levels 2 and above
	for i = 2:Lvl
		Aold = An; An = map(x->x>0, An*A1); Gn = An - Aold
		IdxJ, IdxI = findnz(Gn .* .~I(Nn))
		if isempty(IdxJ); @printf("No new bars at Level %d\n", Lvl); Lvl = i-1; break; end

		# RmFlag = [] # considering no domain restrictions for now
		# deleteat!(IdxI, RmFlag)
		# deleteat!(IdxJ, RmFlag)

		newD = Node[IdxI,:] - Node[IdxJ,:]
		newD = newD./sqrt.(sum(newD.^2, dims=2))

		# Collinearity check
		p = 1; RmFlag = zeros(Int, length(IdxI)); Nb = size(Bars,1)
		for j = 1:Nn
			# Find IdxI[p:q] - NEW bars starting @ node 'j'
			for P=p:length(IdxI); p=P; if IdxI[p]>=j; break; end; end
			q = p; for Q=p:length(IdxI); q=Q; if IdxI[q]>j; break; end; end
			if IdxI[q]>j; q = q - 1; end

			if IdxI[p]==j
				# Find Bars[m:n] - OLD bars starting @ node 'j'
				m = 1; for M=1:Nb; m=M; if Bars[m,1]>=j; break; end; end
				n = m; for N=m:Nb; n=N; if Bars[n,1]>j; break; end; end
				if Bars[n,1]>j; n = n - 1; end

				if Bars[n,1]==j
					# Dot products of old vs. new bars, if collinear: mark
					C = vec(maximum(D[m:n,:] * newD[p:q,:]', dims=1))
					RmFlag[(p-1).+findall(C.>ColTol)] .= 1
				end
			end
		end

		# Remove collinear bars and make symmetric again.
		# Bars that have one angle marked as collinear but the other not, will be spared
		Ind = findall(RmFlag.==0)
		H = sparse(IdxI[Ind], IdxJ[Ind], true, Nn, Nn)
		IdxJ, IdxI = findnz(H+H')

		BarsIJ = [Bars; IdxI IdxJ]
		Bars = BarsIJ[sortperm(BarsIJ[:,1]),:]

		D = Node[Bars[:,1],:] - Node[Bars[:,2],:];
		D = D./sqrt.(sum(D.^2, dims=2))
	end

	# Only return bars {i,j} with i<j (no duplicated bars)
	A = sparse(Bars[:,1], Bars[:,2], true, Nn, Nn)
	IdxJ, IdxI = findnz(tril(A)); Bars = [IdxI IdxJ]
	return Bars
end

max_stars_repo_name: ricardoaf/jutop

max_issues_repo_name: ricardoaf/jutop

max_forks_repo_name: ricardoaf/jutop


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

@doc raw"""Represents an empty directory for a pod. Empty directory volumes support ownership management and SELinux relabeling.

    IoK8sApiCoreV1EmptyDirVolumeSource(;
        medium=nothing,
        sizeLimit=nothing,
    )

    - medium::String : What type of storage medium should back this directory. The default is \"\" which means to use the node's default medium. Must be an empty string (default) or Memory. More info: https://kubernetes.io/docs/concepts/storage/volumes#emptydir
    - sizeLimit::IoK8sApimachineryPkgApiResourceQuantity : Total amount of local storage required for this EmptyDir volume. The size limit is also applicable for memory medium. The maximum usage on memory medium EmptyDir would be the minimum value between the SizeLimit specified here and the sum of memory limits of all containers in a pod. The default is nil which means that the limit is undefined. More info: http://kubernetes.io/docs/user-guide/volumes#emptydir
"""
mutable struct IoK8sApiCoreV1EmptyDirVolumeSource <: SwaggerModel
    medium::Any # spec type: Union{ Nothing, String } # spec name: medium
    sizeLimit::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: sizeLimit

    function IoK8sApiCoreV1EmptyDirVolumeSource(;medium=nothing, sizeLimit=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EmptyDirVolumeSource, Symbol("medium"), medium)
        setfield!(o, Symbol("medium"), medium)
        validate_property(IoK8sApiCoreV1EmptyDirVolumeSource, Symbol("sizeLimit"), sizeLimit)
        setfield!(o, Symbol("sizeLimit"), sizeLimit)
        o
    end
end # type IoK8sApiCoreV1EmptyDirVolumeSource

const _property_map_IoK8sApiCoreV1EmptyDirVolumeSource = Dict{Symbol,Symbol}(Symbol("medium")=>Symbol("medium"), Symbol("sizeLimit")=>Symbol("sizeLimit"))
const _property_types_IoK8sApiCoreV1EmptyDirVolumeSource = Dict{Symbol,String}(Symbol("medium")=>"String", Symbol("sizeLimit")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiCoreV1EmptyDirVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1EmptyDirVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1EmptyDirVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EmptyDirVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EmptyDirVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EmptyDirVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1EmptyDirVolumeSource)
    true
end

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


text:
module Callbacks

using KernelAbstractions
using MPI
using Printf
using Requires
using Statistics

using ..GenericCallbacks

_sync_device(::Type{Array}) = nothing
@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    using .CuArrays, .CuArrays.CUDAdrv
    _sync_device(::Type{CuArray}) = synchronize()
end

"""
    wall_clock_time_per_time_step(interval, array_type, comm)

Returns a callback function that gives wall-clock time per time-step statistics
across MPI ranks in the communicator `comm`.  The times are averaged over the
`interval` of time steps requested for output.  The `array_type` is used
to synchronize the compute device.
"""
function wall_clock_time_per_time_step(interval, array_type, comm)
    _sync_device(array_type)
    before = time_ns()

    callback = GenericCallbacks.EveryXSimulationSteps(interval) do
        _sync_device(array_type)
        after = time_ns()

        time_per_timesteps = after - before

        times = MPI.Gather(time_per_timesteps, 0, comm)
        if MPI.Comm_rank(comm) == 0
            ns_per_s = 1e9
            times = times ./ ns_per_s ./ interval

            @info @sprintf(
                """Wall-clock time per time-step (statistics across MPI ranks)
                   maximum (s) = %25.16e
                   minimum (s) = %25.16e
                   median  (s) = %25.16e
                   std     (s) = %25.16e
                """,
                maximum(times),
                minimum(times),
                median(times),
                std(times),
            )
        end

        _sync_device(array_type)
        before = time_ns()

        nothing
    end

    return callback
end

end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
@symbol_func function h_CS(cur_reactor::AbstractReactor)
  cur_h = cur_reactor.kappa_95 * a(cur_reactor)

  cur_h += b(cur_reactor)

  cur_h += c(cur_reactor)

  cur_h *= 2

  cur_h
end


text:
##########################################################################################
# Consts for edge defaults
##########################################################################################

BASEWIDTH = 1

# Currently the colours are a mix of the OSM standard style and the German OSM style.
# They are modified to work with a white background because we don't have any buildings and
# areas on the map. This specifically affects all roads which are white or light grey in
# OSM standard style, e.g. "residential" and "tertiary". #444 was chosen in these cases.
WAYTYPECOLORS = Dict{String, Colorant}(
    # TODO: roads and link roads should have thin outlines in darker colour
    #       currently impossible with linesegments, Makie feature request is open
    # roads
    "motorway" => colorant"#EB7D54",
    "trunk" => colorant"#E27272",
    "primary" => colorant"#E892A2",
    "secondary" => colorant"#FCD6A4",
    "tertiary" => colorant"#F7FABF",
    "unclassified" => colorant"#444",
    "residential" => colorant"#444",
    # link roads
    "motorway_link" => colorant"#EB7D54",
    "trunk_link" => colorant"#E27272",
    "primary_link" => colorant"#E892A2",
    "secondary_link" => colorant"#FCD6A4",
    "tertiary_link" => colorant"#F7FABF",
    # special road types
    "living_street" => colorant"#444",
    "service" => colorant"#444",
    "pedestrian" => colorant"#939894",
    "track" => colorant"#AC8327",
    "bus_guideway" => colorant"#6666FF",
    "escape" => colorant"#444",
    "raceway" => colorant"#FFC0CA",
    "road" => colorant"#444",
    "busway" => colorant"#444",
    # path
    "footway" => colorant"#F49C8B",
    "bridleway" => colorant"#008000",
    "steps" => colorant"#FC7F6F",
    "corridor" => colorant"#444",
    "path" => colorant"#F49C8B",
    # cycleway
    "cycleway" => colorant"#1619FC",
)

WAYTYPEWIDTHS = Dict{String, Float64}(
    # roads
    "motorway" => 6.0,
    "trunk" => 6.0,
    "primary" => 6.0,
    "secondary" => 6.0,
    "tertiary" => 6.0,
    "unclassified" => 3.0,
    "residential" => 3.0,
    # link roads
    "motorway_link" => 3.0,
    "trunk_link" => 3.0,
    "primary_link" => 3.0,
    "secondary_link" => 3.0,
    "tertiary_link" => 3.0,
    # special road types
    "living_street" => 3.0,
    "service" => 3.0,
    "pedestrian" => 3.0,
    "track" => 1.0,
    "bus_guideway" => 2.0,
    "escape" => 3.0,
    "raceway" => 3.0,
    "road" => 3.0,
    "busway" => 3.0,
    # path
    "footway" => 1.0,
    "bridleway" => 1.0,
    "steps" => 2.0,
    "corridor" => 1.0,
    "path" => 1.0,
    # cycleway
    "cycleway" => 1.0,
)

WAYTYPESTYLES = Dict{String, Symbol}(
    # roads
    "motorway" => :solid,
    "trunk" => :solid,
    "primary" => :solid,
    "secondary" => :solid,
    "tertiary" => :solid,
    "unclassified" => :solid,
    "residential" => :solid,
    # link roads
    "motorway_link" => :solid,
    "trunk_link" => :solid,
    "primary_link" => :solid,
    "secondary_link" => :solid,
    "tertiary_link" => :solid,
    # special road types
    "living_street" => :solid,
    "service" => :solid,
    "pedestrian" => :solid,
    "track" => :dashdot,
    "bus_guideway" => :dashed,
    "escape" => :solid,
    "raceway" => :solid,
    "road" => :solid,
    "busway" => :solid,
    # path
    "footway" => :dot,
    "bridleway" => :dash,
    "steps" => :dot,
    "corridor" => :solid,
    "path" => :dot,
    # cycleway
    "cycleway" => :dot,
)

BUILDINGSCOLORS = colorant"#bbb"


text:
""" Functions for computing and manipulating structure factors """

# TODO:
#  1. Many optimizations + clean-ups still possible in this file.
#       In particular, still a ton of allocations?
#  2. Figure out how to best reduce periodic artifacts along the time FFT

"""
    StructureFactor

Type responsible for computing and updating the static/dynamic structure factor
averaged across multiple spin configurations.

Note that the initial `sys` provided does _not_ enter the structure factor,
it is purely used to determine the size of various results.

The full dynamic structure factor is
``𝒮^{αβ}_{jk}(𝐪, ω) = ⟨M^α_j(𝐪, ω) M^β_k(𝐪, ω)^∗⟩``,
which is an array of shape `[3, 3, B, B, Q1, Q2, Q3, T]`
where `B = nbasis(sys.lattice)`, `Qi = max(1, bz_size_i * L_i)` and
`T = dyn_meas`. By default, `bz_size=ones(d)`.

Indexing the `.sfactor` attribute at `(α, β, j, k, q1, q2, q3, w)`
gives ``𝒮^{αβ}_{jk}(𝐪, ω)`` at `𝐪 = q1 * 𝐛_1 + q2 * 𝐛_2 + q3 * 𝐛_3`, and
`ω = maxω * w / T`, where `𝐛_1, 𝐛_2, 𝐛_3` are the reciprocal lattice vectors
of the system supercell.

Allowed values for the `qi` indices lie in `-div(Qi, 2):div(Qi, 2, RoundUp)`, and allowed
 values for the `w` index lie in `0:T-1`.

The maximum frequency sampled is `ωmax = 2π / (dynΔt * meas_rate)`, and the frequency resolution
is set by `dyn_meas` (the number of spin snapshots measured during dynamics). By default,
`dyn_meas=1`, and the static structure factor is computed. However, beyond
increasing the frequency resolution, increasing `dyn_meas` will also make frequencies become
more accurate.

Setting `reduce_basis` performs the phase-weighted sums over the basis/sublattice
indices, resulting in a size `[3, 3, Q1, Q2, Q3, T]` array.

Setting `dipole_factor` applies the dipole form factor, further reducing the
array to size `[Q1, Q2, Q3, T]`.
"""
struct StructureFactor{A1, A2}
    sfactor       :: A1
    _mag_ft       :: Array{ComplexF64, 6}                    # Buffer for FT of a mag trajectory
    _bz_buf       :: A2                                      # Buffer for phase summation / BZ repeating
    lattice       :: Lattice
    reduce_basis  :: Bool                                    # Flag setting basis summation
    dipole_factor :: Bool                                    # Flag setting dipole form factor
    bz_size       :: NTuple{3, Int}                          # Num of Brillouin zones along each axis
    dynΔt         :: Float64                                 # Timestep size in dynamics integrator
    meas_rate     :: Int                                     # Num timesteps between snapshot saving
    dyn_meas      :: Int                                     # Total number of snapshots to FT
    integrator    :: SphericalMidpoint
    plan          :: FFTW.cFFTWPlan{ComplexF64, -1, true, 6, UnitRange{Int64}}
end

Base.show(io::IO, sf::StructureFactor) = print(io, join(size(sf.sfactor), "x"),  " StructureFactor")
Base.summary(io::IO, sf::StructureFactor) = string("StructureFactor: ", summary(sf.sfactor))

function StructureFactor(sys::SpinSystem; bz_size=(1,1,1), reduce_basis=true,
                         dipole_factor=false, dynΔt::Float64=0.01,
                         dyn_meas::Int=1, meas_rate::Int=10,)
    nb = nbasis(sys.lattice)
    spat_size = size(sys)[2:end]
    q_size = map(s -> s == 0 ? 1 : s, bz_size .* spat_size)
    result_size = (3, q_size..., dyn_meas)
    min_q_idx = -1 .* div.(q_size .- 1, 2)

    spin_ft = zeros(ComplexF64, 3, nb, spat_size..., dyn_meas)
    if reduce_basis
        bz_buf = zeros(ComplexF64, 3, q_size..., dyn_meas)
        bz_buf = OffsetArray(bz_buf, OffsetArrays.Origin(1, min_q_idx..., 0))
    else
        bz_buf = zeros(ComplexF64, 3, nb, q_size..., dyn_meas)
        bz_buf = OffsetArray(bz_buf, OffsetArrays.Origin(1, 1, min_q_idx..., 0))
    end

    if reduce_basis
        if dipole_factor
            sfactor = zeros(Float64, q_size..., dyn_meas)
            sfactor = OffsetArray(sfactor, OffsetArrays.Origin(min_q_idx..., 0))
        else
            sfactor = zeros(ComplexF64, 3, 3, q_size..., dyn_meas)
            sfactor = OffsetArray(sfactor, OffsetArrays.Origin(1, 1, min_q_idx..., 0))
        end
    else
        if dipole_factor
            sfactor = zeros(Float64, nb, nb, q_size..., dyn_meas)
            sfactor = OffsetArray(sfactor, OffsetArrays.Origin(1, 1, min_q_idx..., 0))
        else
            sfactor = zeros(ComplexF64, 3, 3, nb, nb, q_size..., dyn_meas)
            sfactor = OffsetArray(sfactor, OffsetArrays.Origin(1, 1, 1, 1, min_q_idx..., 0))
        end
    end

    integrator = SphericalMidpoint(sys)
    plan = plan_spintraj_fft!(spin_ft)

    StructureFactor{typeof(sfactor), typeof(bz_buf)}(
        sfactor, spin_ft, bz_buf, sys.lattice, reduce_basis, dipole_factor,
        bz_size, dynΔt, meas_rate, dyn_meas, integrator, plan
    )
end

"""
    StructureFactor(snaps::Vector{SpinSystem}; kwargs...)

Construct a `StructureFactor` from a list of spin configurations.
All `SpinSystem`s should have the same underlying lattice -- for
all intents, they should be the "same" system only with different
spin configs. `kwargs` are passed onto the default `StructureFactor`
constructor.
"""
function StructureFactor(snaps::Vector{SpinSystem}; kwargs...)
    if length(snaps) == 0
        error("No snapshots provided, cannot construct StructureFactor")
    end
    sf = StructureFactor(snaps[1]; kwargs...)
    for snap in snaps
        update!(sf, snap)
    end
    sf
end

"""
    StructureFactor(snaps::Vector{Array{Vec3, 4}}, crystal; kwargs...)

Construct a `StructureFactor` from a list of spin configurations,
and a `Crystal` specifying the underlying lattice. `kwargs` are passed
onto the default `StructureFactor` constructor.
"""
function StructureFactor(snaps::Vector{Array{Vec3, 4}}, crystal; kwargs...)
    if length(snaps) == 0
        error("No snapshots provided, cannot construct StructureFactor")
    end
    sys = SpinSystem(crystal, Vector{Interaction}(), size(snaps[1])[2:end])
    sf = StructureFactor(sys; kwargs...)
    for snap in snaps
        sys.sites .= snap
        update!(sf, sys)
    end
    sf
end

"""
Updates `M` in-place to hold the magnetization vectors obtained by scaling `s`
 by the appropriate spin magnitudes and g-tensors in `sites_info`.
This function assumes `M` has a first index of length 3, which correspond
 to the magnetization components. (Rather than storing an Array{Vec3}).
"""
function _compute_mag!(M, sys::SpinSystem)
    for b in 1:nbasis(sys)
        gS = sys.sites_info[b].g * sys.sites_info[b].S
        for idx in eachcellindex(sys)
            M[:, b, idx] .= gS * sys[b, idx]
        end
    end
end

"""
    update!(sf::StructureFactor, sys::SpinSystem{3})

Accumulates a contribution to the dynamic structure factor from the spin
configuration currently in `sys`.
"""
function update!(sf::StructureFactor, sys::SpinSystem)
    @unpack sfactor, _mag_ft, _bz_buf = sf
    @unpack reduce_basis, dipole_factor, bz_size = sf

    # Evolve the spin state forward in time to form a trajectory
    # Save off the magnetic moments 𝐦_i(t) = g_i S_i 𝐬_i(t) into _mag_ft
    dynsys = deepcopy(sys)
    sf.integrator.sys = dynsys
    T_dim = ndims(_mag_ft)
    _compute_mag!(selectdim(_mag_ft, T_dim, 1), dynsys)
    for nsnap in 2:sf.dyn_meas
        for _ in 1:sf.meas_rate
            evolve!(sf.integrator, sf.dynΔt)
        end
        _compute_mag!(selectdim(_mag_ft, T_dim, nsnap), dynsys)
    end

    # Fourier transform the trajectory in space + time
    fft_spin_traj!(_mag_ft, plan=sf.plan)

    # Optionally sum over basis sites then accumulate the conjugate outer product into sfactor
    # Accumulate the conjugate outer product into sfactor, with optionally:
    #   1) Doing a phase-weighting sum to reduce the basis atom dimensions
    #   2) Applying the neutron dipole factor to reduce the spin component dimensions
    if reduce_basis
        phase_weight_basis!(_bz_buf, _mag_ft, sys.lattice)
        if dipole_factor
            accum_dipole_factor!(sfactor, _bz_buf, sys.lattice)
        else
            outerprod_conj!(sfactor, _bz_buf, 1)
        end
    else
        expand_bz!(_bz_buf, _mag_ft)
        if dipole_factor
            accum_dipole_factor_wbasis!(sfactor, _bz_buf, sys.lattice)
        else
            outerprod_conj!(sfactor, _bz_buf, (1, 2))
        end
    end
end

"""
    zero!(sf::StructureFactor)

Zeros out the accumulated structure factor.
"""
function zero!(sf::StructureFactor)
    sf.sfactor .= 0
end

"""
    apply_dipole_factor(sf::StructureFactor) :: StructureFactor

Apply the neutron dipole factor to a dynamic structure factor.
"""
function apply_dipole_factor(sf::StructureFactor)
    if sf.dipole_factor == true
        return sf
    end

    dip_sfactor = apply_dipole_factor(sf.sfactor, sf.lattice)
    StructureFactor(
        dip_sfactor, copy(sf._mag_ft), copy(sf._bz_buf), sf.lattice,
        sf.reduce_basis, true, sf.bz_size, sf.dynΔt, sf.meas_rate,
        sf.dyn_meas, sf.integrator, sf.plan
    )
end

function apply_dipole_factor(struct_factor::OffsetArray{ComplexF64}, lattice::Lattice)
    recip = gen_reciprocal(lattice)

    T = size(struct_factor)[end]
    result = zeros(Float64, axes(struct_factor)[3:end])
    for q_idx in CartesianIndices(axes(struct_factor)[3:end-1])
        q = recip.lat_vecs * SVector{3, Float64}(Tuple(q_idx) ./ lattice.size)
        q = q / (norm(q) + 1e-12)
        dip_factor = reshape(I(3) - q * q', 3, 3, 1)
        for t in 0:T-1
            result[q_idx, t] = real(dot(dip_factor, struct_factor[:, :, q_idx, t]))
        end
    end
    return result
end

"""
    dynamic_structure_factor(sys, sampler; nsamples=10, dynΔt=0.01, meas_rate=10,
                             dyn_meas=100, bz_size, thermalize=10, reduce_basis=true,
                             verbose=false)

Measures the full dynamic structure factor tensor of a spin system, for the requested range
of 𝐪-space and range of frequencies ω. Returns ``𝒮^{αβ}(𝐪, ω) = ⟨S^α(𝐪, ω) S^β(𝐪, ω)^∗⟩``,
which is an array of shape `[3, 3, Q1, ..., Qd, T]`
where `Qi = max(1, bz_size_i * L_i)` and `T = dyn_meas`. By default, `bz_size=ones(d)`.

Setting `reduce_basis=false` makes it so that the basis/sublattice indices are not
phase-weighted and summed over, making the shape of the result `[3, 3, B, B, Q1, ..., Qd, T]`
where `B = nbasis(sys)` is the number of basis sites in the unit cell.

`nsamples` sets the number of thermodynamic samples to measure and average
 across from `sampler`. `dynΔt` sets the integrator timestep during dynamics,
 and `meas_rate` sets how often snapshots are recorded during dynamics. `dyn_meas`
 sets the total number snapshots taken. The sampler is thermalized by sampling
 `thermalize` times before any measurements are made.

The maximum frequency sampled is `ωmax = 2π / (dynΔt * meas_rate)`, and the frequency resolution
 is set by `dyn_meas` (the number of spin snapshots measured during dynamics). However, beyond
 increasing the resolution, `dyn_meas` will also make all frequencies become more accurate.

Indexing the result at `(α, β, q1, ..., qd, w)` gives ``S^{αβ}(𝐪, ω)`` at
    `𝐪 = q1 * a⃰ + q2 * b⃰ + q3 * c⃰`, and `ω = maxω * w / T`, where `a⃰, b⃰, c⃰`
    are the reciprocal lattice vectors of the system supercell.

Allowed values for the `qi` indices lie in `-div(Qi, 2):div(Qi, 2, RoundUp)`, and allowed
 values for the `w` index lie in `0:T-1`.
"""
function dynamic_structure_factor(
    sys::SpinSystem, sampler::S; nsamples::Int=10,
    thermalize::Int=10, bz_size=(1,1,1), reduce_basis::Bool=true,
    dipole_factor::Bool=false, dynΔt::Float64=0.01, dyn_meas::Int=100,
    meas_rate::Int=10, verbose::Bool=false
) where {S <: AbstractSampler}

    sf  = StructureFactor(sys; dynΔt=dynΔt, meas_rate=meas_rate, dyn_meas=dyn_meas,
                               bz_size=bz_size, reduce_basis=reduce_basis,
                               dipole_factor=dipole_factor)

    if verbose
        println("Beginning thermalization...")
    end

    # Equilibrate the system by sampling from it `nsamples` times (discarding results)
    thermalize!(sampler, thermalize)

    if verbose
        println("Done thermalizing. Beginning measurements...")
    end

    progress = Progress(nsamples; dt=1.0, desc="Sample: ", enabled=verbose)
    for n in 1:nsamples
        sample!(sampler)
        update!(sf, sys)
        next!(progress)
    end

    return sf
end

"""
    static_structure_factor(sys, sampler; nsamples, dynΔt, meas_rate, dyn_meas
                                          bz_size, thermalize, verbose)

Measures the static structure factor tensor of a spin system, for the requested range
of 𝐪-space. Returns ``𝒮^{αβ}(𝐪) = ⟨S^α(𝐪) S^β(𝐪)^∗⟩``,
which is an array of shape `[3, 3, Q1, ..., Qd]` where `Qi = max(1, bz_size_i * L_i)`.
By default, `bz_size=ones(d)`.

`nsamples` sets the number of thermodynamic samples to measure and average
 across from `sampler`. `dynΔt` sets the integrator timestep during dynamics,
 and `meas_rate` sets how many timesteps are performed between recording snapshots.
 `dyn_meas` sets the total number snapshots taken. The sampler is thermalized by sampling
 `thermalize` times before any measurements are made.

Indexing the result at `(α, β, q1, ..., qd)` gives ``𝒮^{αβ}(𝐪)`` at
    `𝐪 = q1 * a⃰ + q2 * b⃰ + q3 * c⃰`, where `a⃰, b⃰, c⃰`
    are the reciprocal lattice vectors of the system supercell.

Allowed values for the `qi` indices lie in `-div(Qi, 2):div(Qi, 2, RoundUp)`.
"""
function static_structure_factor(sys::SpinSystem, sampler::S; kwargs...) where {S <: AbstractSampler}
    dynamic_structure_factor(sys, sampler; dyn_meas=1, kwargs...)
end

#= Non-exposed internal functions below =#

"""
    plan_spintraj_fft(spin_traj::Array{Vec3})

Prepares an out-of-place FFT plan for a spin trajectory array of
size [B, D1, ..., Dd, T]
"""
function plan_spintraj_fft(spin_traj::Array{Vec3})
    spin_traj = _reinterpret_from_spin_array(spin_traj)
    return plan_fft(spin_traj, 3:ndims(spin_traj))
end

"""
    plan_spintraj_fft!(spin_traj::Array{ComplexF64})

Prepares an in-place FFT plan for a spin trajectory array of
size [3, B, D1, ..., Dd, T].
"""
function plan_spintraj_fft!(spin_traj::Array{ComplexF64})
    return plan_fft!(spin_traj, 3:ndims(spin_traj))
end

"""
    fft_spin_traj!(res, spin_traj; plan=nothing)

In-place version of `fft_spin_traj`. `res` should be an `Array{ComplexF64}` of size
`[3, B, D1, ..., Dd, T]` to hold the result, matching the size `[B, D1, ..., Dd, T]`
of `spin_traj`.
"""
function fft_spin_traj!(res::Array{ComplexF64}, spin_traj::Array{Vec3};
                        plan::Union{Nothing, FFTW.cFFTWPlan}=nothing)
    @assert size(res) == tuplejoin(3, size(spin_traj)) "fft_spins size not compatible with spin_traj size"

    # Reinterpret array to add the spin dimension explicitly
    # Now of shape [3, B, D1, ..., Dd, T]
    spin_traj = _reinterpret_from_spin_array(spin_traj)

    # FFT along the spatial indices, and the time index
    if isnothing(plan)
        res .= spin_traj
        fft!(res, 3:ndims(spin_traj))
    else
        mul!(res, plan, spin_traj)
    end

    return res
end

function fft_spin_traj!(spin_traj::Array{ComplexF64};
                        plan::Union{Nothing, FFTW.cFFTWPlan}=nothing)
    if isnothing(plan)
        fft!(spin_traj, 3:ndims(spin_traj))
    else
        spin_traj = plan * spin_traj
    end
end

"""
    fft_spin_traj(spin_traj; bz_size, plan=nothing)

Takes in a `spin_traj` array of spins (Vec3) of shape `[B, D1, ..., Dd, T]`,
 with `D1 ... Dd` being the spatial dimensions, B the sublattice index,
 and `T` the time axis.
Computes and returns an array of the shape `[3, B, D1, ..., Dd, T]`,
 holding spatial and temporal fourier transforms ``S^α_b(𝐪, ω)``. The spatial
 fourier transforms are done periodically, but the temporal axis is
 internally zero-padded to avoid periodic contributions. *(Avoiding
 periodic artifacts not implemented yet)*
"""
function fft_spin_traj(spin_traj::Array{Vec3};
                       plan::Union{Nothing, FFTW.cFFTWPlan}=nothing)
    fft_spins = zeros(ComplexF64, 3, size(spin_traj)...)
    fft_spin_traj!(fft_spins, spin_traj; plan=plan)
end

"""
    phase_weight_basis(spin_traj_ft, bz_size, lattice)

Combines the sublattices of `spin_traj_ft` with the appropriate phase factors, producing
 the quantity ``S^α(𝐪, ω)`` within the number of Brillouin zones requested by `bz_size`.
Specifically, computes:

``S^α(𝐪, ω) = ∑_b e^{-i𝐫_b ⋅ 𝐪} S^α_b(𝐪, ω)``

where ``b`` is the basis index and ``𝐫_b`` is the associated basis vector.
``S^α_b(𝐪, ω)`` is periodically repeated past the first Brillouin zone,
but the resulting ``S^α(𝐪, ω)`` will not necessarily be periodic.
"""
function phase_weight_basis(spin_traj_ft::Array{ComplexF64},
                            lattice::Lattice, bz_size=nothing)
    if isnothing(bz_size)
        bz_size = ones(ndims(lattice) - 1)
    end

    bz_size = convert(SVector{3, Int}, bz_size)                  # Number of Brilloin zones along each axis
    spat_size = lattice.size                                     # Spatial lengths of the system
    T = size(spin_traj_ft, ndims(spin_traj_ft))                  # Number of timesteps in traj / frequencies in result
    q_size = map(s -> s == 0 ? 1 : s, bz_size .* spat_size)      # Total number of q-points along each q-axis of result
    result_size = (3, q_size..., T)
    min_q_idx = -1 .* div.(q_size .- 1, 2)

    result = zeros(ComplexF64, result_size)
    result = OffsetArray(result, OffsetArrays.Origin(1, min_q_idx..., 0))
    phase_weight_basis!(result, spin_traj_ft, lattice)
end

"""
    phase_weight_basis!(res, spin_traj_ft, lattice)

Like `phase_weight_basis`, but in-place. Infers `bz_size` from `size(res)`.
"""
function phase_weight_basis!(res::OffsetArray{ComplexF64},
                             spin_traj_ft::Array{ComplexF64},
                             lattice::Lattice)
    # Check that spatial size of spin_traj_ft same as spatial size of lattice
    spat_size = size(lattice)[2:end]
    valid_size = size(spin_traj_ft)[3:end-1] == spat_size
    @assert valid_size "`size(spin_traj_ft)` not compatible with `lattice`"
    # Check that q_size is elementwise either an integer multiple of spat_size, or is 1.
    q_size = size(res)[2:end-1]
    valid_q_size = all(map((qs, ss) -> qs % ss == 0 || qs == 1, q_size, spat_size))
    @assert valid_q_size "`size(res)` not compatible with `size(spin_traj_ft)`"

    recip = gen_reciprocal(lattice)

    T = size(spin_traj_ft)[end]

    fill!(res, 0.0)
    for q_idx in CartesianIndices(axes(res)[2:end-1])
        q = recip.lat_vecs * SVector{3, Float64}(Tuple(q_idx) ./ lattice.size)
        wrap_q_idx = modc(q_idx, spat_size) + one(CartesianIndex{3})
        for (b_idx, b) in enumerate(lattice.basis_vecs)
            phase = exp(-im * (b ⋅ q))
            # Note: Lots of allocations here. Fix?
            # Warning: Cannot replace T with 1:end due to Julia issues with end and CartesianIndex
            @. res[:, q_idx, 0:T-1] += spin_traj_ft[:, b_idx, wrap_q_idx, 1:T] * phase
        end
    end

    return res
end

# === Helper functions for outerprod_conj === #

# TODO: Bounds checking
""" Given `size`, compute a new size tuple where there is an extra `1` before each dim in `dims`.
"""
function _outersizeα(size, dims)
    if length(dims) == 0
        return size
    end

    newsize = tuplejoin(size[1:dims[1]-1], 1)
    for i in 2:length(dims)
        newsize = tuplejoin(newsize, size[dims[i-1]:dims[i]-1], 1)
    end
    tuplejoin(newsize, size[dims[end]:end])
end

""" Given `size`, compute a new size tuple where there is an extra `1` after each dim in `dims`.
"""
_outersizeβ(size, dims) = length(dims) == 0 ? size : _outersizeα(size, dims .+ 1)

# ========================================== #

"""
    outerprod_conj(S, [dims=1])

Computes the outer product along the selected dimensions, with a complex
conjugation on the second copy in the product.

I.e. given a complex array of size `[D1, ..., Di, ..., Dd]`, for each
dimension `i` in `dims` this will create a new axis of the same size to make an
array of size `[D1, ..., Di, Di, ..., Dd]` where the new axes are formed by
an outer product of the vectors of the original axes with a complex conjugation
on one copy.
"""
function outerprod_conj(S, dims=1)
    sizeα = _outersizeα(axes(S), dims)
    sizeβ = _outersizeβ(axes(S), dims)
    Sα = reshape(S, sizeα)
    Sβ = reshape(S, sizeβ)
    @. Sα * conj(Sβ)
end

"""
    outerprod_conj!(res, S, [dims=1])

Like `outerprod_conj`, but accumulates the result in-place into `res`.
"""
function outerprod_conj!(res, S, dims=1)
    sizeα = _outersizeα(axes(S), dims)
    sizeβ = _outersizeβ(axes(S), dims)
    Sα = reshape(S, sizeα)
    Sβ = reshape(S, sizeβ)
    @. res += Sα * conj(Sβ)
end

"""
    expand_bz!(res::OffsetArray, S::Array)

Copy S periodically into res, with the periodic boundaries set by the
spatial axes of S. Assumes that S is of shape [3, B, L1, L2, L3, T], and
that res is of shape [3, B, Q1, Q2, Q3, T], with all Qi >= Li.
"""
function expand_bz!(res::OffsetArray{ComplexF64}, S::Array{ComplexF64})
    spat_size = size(S)[end-3:end-1]
    T = size(S, ndims(S))

    for t in 1:T
        for q_idx in CartesianIndices(axes(res)[end-3:end-1])
            wrap_q_idx = modc(q_idx, spat_size) + CartesianIndex(1, 1, 1)
            res[:, :, q_idx, t-1] = S[:, :, wrap_q_idx, t]
        end
    end
end

#= These two "accumulate with dipole factor" functions are so close that it seems
    like they should be joined, but I cannot think of a clever way to do so.
=#

"""
    accum_dipole_factor!(res, S, lattice)

Given complex `S` of size [3, Q1, ..., QD, T] and `res` of size [Q1, ..., QD, T],
accumulates the structure factor from `S` with the dipole factor applied into `res`.
"""
function accum_dipole_factor!(res, S, lattice::Lattice)
    recip = gen_reciprocal(lattice)
    for q_idx in CartesianIndices(axes(res)[end-3:end-1])
        q = recip.lat_vecs * SVector{3, Float64}(Tuple(q_idx) ./ lattice.size)
        q = q / (norm(q) + 1e-12)
        dip_factor = I(3) - q * q'

        for α in 1:3
            for β in 1:3
                dip_elem = dip_factor[α, β]
                @. res[q_idx, :] += dip_elem * real(S[α, q_idx, :] * conj(S[β, q_idx, :]))
            end
        end
    end
end

"""
    accum_dipole_factor_wbasis!(res, S, lattice)

Given complex `S` of size [3, B, Q1, ..., QD, T] and real `res` of size [B, B, Q1, ..., QD, T],
accumulates the structure factor from `S` with the dipole factor applied into `res`.
"""
function accum_dipole_factor_wbasis!(res, S, lattice::Lattice)
    recip = gen_reciprocal(lattice)
    nb = nbasis(lattice)
    Sα = reshape(S, _outersizeα(axes(S), 2))  # Size [3, 1, B, ...]
    Sβ = reshape(S, _outersizeβ(axes(S), 2))  # Size [3, B, 1, ...]

    for q_idx in CartesianIndices(axes(res)[end-3:end-1])
        q = recip.lat_vecs * SVector{3, Float64}(Tuple(q_idx) ./ lattice.size)
        q = q / (norm(q) + 1e-12)
        dip_factor = I(3) - q * q'

        for α in 1:3
            for β in 1:3
                dip_elem = dip_factor[α, β]
                @. res[:, :, q_idx, :] += dip_elem * real(Sα[α, :, :, q_idx, :] * Sβ[β, :, :, q_idx, :])
            end
        end
    end
end


text:
"""
    logreg(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame, Fscore::Bool=true)

Performs a logistic regression on selected variables.

# Arguments

- `topredict::String`: The name of the column of the dataframe of interest that is to be predicted.
- `varnames::Vector{Symbol}`: The explanatory variables to include in the regression.
- `train::DataFrame`: The dataset used to train the model.
- `test::DataFrame`: The dataset used to test the model.
- `Fscore::Bool=true`: If true, returns the F1 score of the model. If false, returns the vector of predictions. Default = true.

"""
function logreg(topredict::String, varnames::Vector{Symbol}, train::DataFrame, test::DataFrame, Fscore::Bool=true)
    threshold = collect(0.1:.01:.6)

    model = glm(Term(Symbol(topredict)) ~ sum(Term.(varnames)), train, Bernoulli(), LogitLink())
    Ŷ = convert(Vector{Float64}, GLM.predict(model, test))
    r = MLBase.roc(test[:, Symbol(topredict)], Ŷ, threshold)

    Fscore ? (return maximum(f1score.(r))) : return Ŷ
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO