text:
##########
## fill ##
##########

function Base.fill(
    value::TrackedReal,
    dims::Vararg{Union{Integer, AbstractUnitRange}},
)
    return track(fill, value, dims...)
end
@grad function fill(v::Real, dims...)
    return fill(value(v), dims...), function(Δ)
        size(Δ) ≢  dims && error("Dimension mismatch")
        return (sum(Δ), map(_->nothing, dims)...)
    end
end

###############
## any & all ##
###############

Base.any(f::Function, x::TrackedArray; dims=:) = any(f, value(x), dims = dims)
Base.all(f::Function, x::TrackedArray; dims=:) = all(f, value(x), dims = dims)

#########
## cat ##
#########

function combinations(xs, n)
    n < 1 && return [[]]
    cs = combinations(xs, n-1)
    [[x, c...] for x in xs, c in cs]
end

for f in [:hcat, :vcat]
    for i = 0:2, c = combinations([:AbstractVector, :TrackedVector, :AbstractMatrix, :TrackedMatrix, :Number, :TrackedReal], i)
        cnames = map(_ -> gensym(), c)
        @eval begin
            Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::TrackedVector) = track($f, $(cnames...), x)
            Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::TrackedMatrix) = track($f, $(cnames...), x)
            Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::TrackedReal) = track($f, $(cnames...), x)
        end
        for T in [
            :AbstractVector,
            :AbstractMatrix,
            :Number,
            :AbstractVecOrMat,
            :(Union{AbstractVector, Number}),
        ]
            @eval begin
                Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::TrackedVector, xs::$T...) = track($f, $(cnames...), x, xs...)
                Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::TrackedMatrix, xs::$T...) = track($f, $(cnames...), x, xs...)
                Base.$f($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::TrackedReal, xs::$T...) = track($f, $(cnames...), x, xs...)
            end
        end
    end
end

@grad function vcat(xs::Union{Number, AbstractVecOrMat}...)
    xs_value = value.(xs)
    out_value = vcat(xs_value...,)
    function back(Δ)
        start = 0
        Δs = map(xs) do xsi
            if xsi isa Number
                d = Δ[start+1]
            else
                d = Δ[start+1:start+size(xsi,1), :]
            end
            start += size(xsi, 1)
            d
        end
        return (Δs...,)
    end
    return out_value, back
end

@grad function hcat(xs::Union{Number, AbstractVecOrMat}...)
    xs_value = value.(xs)
    out_value = hcat(xs_value...,)
    function back(Δ)
        start = 0
        Δs = map(xs) do xsi
            d = if ndims(xsi) == 0
                Δ[start+1]
            elseif ndims(xsi) == 1
                Δ[:, start+1]
            else
                i = map(_ -> :, size(xsi)) |> Base.tail |> Base.tail
                Δ[:, start+1:start+size(xsi,2), i...]
            end
            start += size(xsi, 2)
            d
        end
        return (Δs...,)
    end
    return out_value, back
end

for i = 0:2, c = combinations([:AbstractArray, :TrackedArray, :Number, :TrackedReal], i)
    cnames = map(_ -> gensym(), c)
    @eval Base.cat($([:($x::$c) for (x, c) in zip(cnames, c)]...), x::Union{TrackedArray,TrackedReal}, xs::Union{AbstractArray,Number}...; dims) = track(cat, $(cnames...), x, xs...; dims=dims)
end
@grad function cat(Xs::Union{Number, AbstractArray}...; dims)
    Xs_value = value.(Xs)
    return cat(Xs_value...; dims = dims), Δ -> begin
        start = ntuple(i -> 0, Val(ndims(Δ)))
        Δs = map(Xs) do xs
            if xs isa Number
                d = Δ[start+1]
                start = start .+ 1
            else
                dim_xs = 1:ndims(xs)
                till_xs = ntuple((i -> i in dims ? (i in dim_xs ? size(xs,i) : 1) : 0), Val(ndims(Δ)))
                xs_in_Δ = ntuple(i -> till_xs[i] > 0 ? (start[i]+1:start[i]+till_xs[i]) : Colon(), Val(ndims(Δ)))
                d = reshape(Δ[xs_in_Δ...],size(xs))
                start = start .+ till_xs
            end
            d
        end
        return (Δs...,)
    end
end


text:
module SSLTestDataFrames
using SimpleSDMLayers
using DataFrames
using Test

temperature = SimpleSDMPredictor(WorldClim, BioClim, 1)

df = DataFrame(longitude = [30.0, 31.0], latitude = [0.0, 1.0], values = [42.0, 15.0])

@test eltype(temperature[df]) <: Number

temperature_clip = clip(temperature, df)

@test typeof(temperature_clip) == typeof(temperature)

df1 = DataFrame(temperature_clip)
df2 = DataFrame([temperature_clip, temperature_clip])

@test typeof(df1) == DataFrame
@test eltype(df1.values) == Union{Missing, eltype(temperature_clip)}
@test typeof(df2) == DataFrame
@test eltype(df2.x1) == Union{Missing, eltype(temperature_clip)}

@test typeof(SimpleSDMPredictor(df, :values, temperature_clip)) <: SimpleSDMLayer
@test typeof(SimpleSDMPredictor(df, :values, temperature_clip)) <: SimpleSDMPredictor

l1 = SimpleSDMPredictor(df1, :values, temperature_clip)
l2 = SimpleSDMPredictor(df2, :x2, temperature_clip)
for l in (l1, l2)
    @test isequal(l.grid, temperature_clip.grid)
    @test isequal(longitudes(l), longitudes(temperature_clip))
    @test isequal(latitudes(l), latitudes(temperature_clip))
end

mbool = mask(temperature_clip, df, Bool)
@test eltype(mbool) == Bool

mfloat = mask(temperature_clip, df, Float64)
@test eltype(mfloat) == Float64

@test sum(mfloat) >= sum(mbool)
@test sum(mfloat) == nrow(df)

end

ext: jl


text:
using GLPK

@testset "IO" begin
    @testset "Input parameters" begin
        optimizer = GLPK.Optimizer
        params = create_param( stages                  = 3,
                            model_constructor_grid  = DCPPowerModel,
                            post_method             = PowerModels.post_opf,
                            optimizer                  = optimizer)
        # no hydro data
        @test params["stages"] == 3
        @test params["post_method"] == PowerModels.post_opf
        @test params["optimizer"] == optimizer
        @test params["model_constructor_grid"] == DCPPowerModel

    end
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ f63bf213-75fd-4722-bc7a-8cded43bc53a
begin
	using Parameters
end

# ╔═╡ 1c7e58cb-e4e9-49af-b782-b107e61196d7
# how many iterations?
number_of_iterations = 1000

# ╔═╡ a0bb564e-5c10-4d0e-8685-b4d5af6bbed2
# how many agents?
number_of_agents = 10

# ╔═╡ 624d2477-c19a-49c8-888d-23a580b63842
# how many assets?
number_of_assets = 10

# ╔═╡ 95345a88-747e-4d33-96d2-0d17b81d4518
mutable struct Test
	a::Int64
	update::Function
	Test() = new()
end

# ╔═╡ d11be206-e79b-47de-a05b-d654693eb7f3
function _update!(obj::Test, value::Int64)
	obj.a = value
end

# ╔═╡ 71255cb5-c21d-4ed2-ae10-a430af4b4511
begin

	# build test -
	test_obj = Test()
	test_obj.a = 10
	test_obj.update = _update!

	test_obj
end

# ╔═╡ c445c71a-689f-4449-ba57-82abdf755218
test_obj.update(test_obj,200)

# ╔═╡ ce08a9c6-8638-482c-8f34-9f7ad91ae364
test_obj

# ╔═╡ bcebefc3-d63e-4340-a944-63fe273356b9
function ingredients(path::String; module_name::Symbol = :model)

	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object

	# JV change ... want to
	# name = Symbol(basename(path))
	name = module_name

	# original -
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 31e84c23-0f9a-42f1-ab34-00def17db241
begin

	# load my code into module M -
	M = ingredients("../Include.jl"; module_name = :ABM)

	# alias my types that are declared -
	VLAgentModel = M.VLAgentModel
	VLGameWorld = M.VLGameWorld

	# return -
	nothing
end

# ╔═╡ 8b63e2b6-25e1-40cf-beb6-dc74d930806e
# setup the game world -
game_world = M.build(VLGameWorld; number_of_iterations=number_of_iterations, number_of_assets=number_of_assets)

# ╔═╡ b57bc814-c397-450e-97c9-e3a0a84c03f7
# setup the agents -
begin

	my_agent_array = Array{VLAgentModel,1}()
	for agent_index = 1:number_of_agents

		# let's use the defaults for the various agent settings -
		agent_model = M.build(VLAgentModel; number_of_iterations=number_of_iterations, number_of_assets=number_of_assets)

		# grab -
		push!(my_agent_array, agent_model)
	end

	my_agent_array
end

# ╔═╡ db7b8487-2467-47ad-80b0-413574d02a38
# finally, update the game world, attach the agents that live in the game world -
game_world.agent_array = my_agent_array

# ╔═╡ 6c53c6b9-f870-4fce-8302-093c9ae74e0f
html"""
<style>
main {
    max-width: 1200px;
    margin: auto;
	width: 90%;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"

[compat]
Parameters = "~0.12.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"
"""

# ╔═╡ Cell order:
# ╠═f63bf213-75fd-4722-bc7a-8cded43bc53a
# ╠═31e84c23-0f9a-42f1-ab34-00def17db241
# ╠═1c7e58cb-e4e9-49af-b782-b107e61196d7
# ╠═a0bb564e-5c10-4d0e-8685-b4d5af6bbed2
# ╠═624d2477-c19a-49c8-888d-23a580b63842
# ╠═8b63e2b6-25e1-40cf-beb6-dc74d930806e
# ╠═b57bc814-c397-450e-97c9-e3a0a84c03f7
# ╠═db7b8487-2467-47ad-80b0-413574d02a38
# ╠═95345a88-747e-4d33-96d2-0d17b81d4518
# ╠═d11be206-e79b-47de-a05b-d654693eb7f3
# ╠═71255cb5-c21d-4ed2-ae10-a430af4b4511
# ╠═c445c71a-689f-4449-ba57-82abdf755218
# ╠═ce08a9c6-8638-482c-8f34-9f7ad91ae364
# ╟─bcebefc3-d63e-4340-a944-63fe273356b9
# ╟─6c53c6b9-f870-4fce-8302-093c9ae74e0f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using SafeTestsets

@safetestset "My function testing" begin include("testsets.jl") end

#nn_arr 3  [1, :]392    4    7   63   62    0    2    0   61    0  332
# [end, arr] 0  527    0  508    0  417  396  397  395  526  336  507


text:
# SPDX-FileCopyrightText: 2021 Johannes Hampp <johannes.hampp@zeu.jlug.de>
#
# SPDX-License-Identifier: MIT

# Activate julia environment
using Pkg
Pkg.activate("./workflow/envs")
using GlobalEnergyGIS
import YAML

region_name = snakemake.wildcards["region"]
region_definitions = YAML.load_file(snakemake.input["region_definitions"])

region = region_definitions["regions"][region_name]

# Ensure members are alphabetically sorted by code
region = sort(collect(region))

# GEGIS expects a matrix with name + GADM(...) entry for each
# member of a region. Convert to appropriate shape.
members = [[k GADM(v)] for (k,v) in region]
members = vcat(members...)

# Create and save the GEGIS region
saveregions(region_name, members)

max_stars_repo_name: euronion/synde

max_issues_repo_name: euronion/synde

max_forks_repo_name: euronion/synde


text:
#########################################################################
# Wrapper around Python exceptions

type PyError <: Exception
    msg::String # message string from Julia context, or "" if none

    # info returned by PyErr_Fetch/PyErr_Normalize
    T::PyObject
    val::PyObject
    traceback::PyObject

    # generate a PyError object.  Should normally only be called when
    # PyErr_Occurred returns non-NULL, and clears the Python error
    # indicator.  Assumes Python is initialized!
    function PyError(msg::String)
        exc = Array(PyPtr, 3)
        pexc = convert(Uint, pointer(exc))
        # equivalent of passing C pointers &exc[1], &exc[2], &exc[3]:
        ccall((@pysym :PyErr_Fetch), Void, (Uint,Uint,Uint),
              pexc, pexc + sizeof(PyPtr), pexc + 2*sizeof(PyPtr))
        ccall((@pysym :PyErr_NormalizeException), Void, (Uint,Uint,Uint),
              pexc, pexc + sizeof(PyPtr), pexc + 2*sizeof(PyPtr))
        new(msg, exc[1], exc[2], exc[3])
    end
end

function show(io::IO, e::PyError)
    print(io, "PyError",
          isempty(e.msg) ? e.msg : string(" (",e.msg,")"),
          " ")

    if e.T.o == C_NULL
        println(io, "None")
    else
        println(io, pystring(e.T), "\n", pystring(e.val))
    end

    if e.traceback.o != C_NULL
        o = pycall(format_traceback::PyObject, PyObject, e.traceback)
        if o.o != C_NULL
            for s in PyVector{String}(o)
                print(io, s)
            end
        end
    end
end

#########################################################################
# Conversion of Python exceptions into Julia exceptions

# call to discard Python exceptions
pyerr_clear() = ccall((@pysym :PyErr_Clear), Void, ())

function pyerr_check(msg::String, val::Any)
    # note: don't call pyinitialize here since we will
    # only use this in contexts where initialization was already done
    if ccall((@pysym :PyErr_Occurred), PyPtr, ()) != C_NULL
        throw(PyError(msg))
    end
    val # the val argument is there just to pass through to the return value
end

pyerr_check(msg::String) = pyerr_check(msg, nothing)
pyerr_check() = pyerr_check("")

# Macros for common pyerr_check("Foo", ccall((@pysym :Foo), ...)) pattern.
# (The "i" variant assumes Python is initialized.)
macro pychecki(ex)
    :(pyerr_check($(string(ex.args[1].args[2].args[1])), $ex))
end
macro pycheck(ex)
    quote
        @pyinitialize
        @pychecki $ex
    end
end

# Macros to check that ccall((@pysym :Foo), ...) returns value != bad
# (The "i" variants assume Python is initialized.)
macro pycheckvi(ex, bad)
    quote
        val = $ex
        if val == $bad
            # throw a PyError if available, otherwise throw ErrorException
            pyerr_check($(string(ex.args[1].args[2].args[1])))
            error($(string(ex.args[1].args[2].args[1])), " failed")
        end
        val
    end
end
macro pycheckni(ex)
    :(@pycheckvi $ex C_NULL)
end
macro pycheckzi(ex)
    :(@pycheckvi $ex -1)
end
macro pycheckv(ex, bad)
    quote
        @pyinitialize
        @pycheckvi $ex $bad
    end
end
macro pycheckn(ex)
    quote
        @pyinitialize
        @pycheckni $ex
    end
end
macro pycheckz(ex)
    quote
        @pyinitialize
        @pycheckzi $ex
    end
end

#########################################################################
# Mapping of Julia Exception types to Python exceptions

pyexc = Dict{DataType, PyPtr}()
type PyIOError <: Exception end

function pyexc_initialize()
    global pyexc
    exc = [Exception => :PyExc_RuntimeError,
           ErrorException => :PyExc_RuntimeError,
           SystemError => :PyExc_SystemError,
           TypeError => :PyExc_TypeError,
           ParseError => :PyExc_SyntaxError,
           ArgumentError => :PyExc_ValueError,
           KeyError => :PyExc_KeyError,
           LoadError => :PyExc_ImportError,
           MethodError => :PyExc_RuntimeError,
           EOFError => :PyExc_EOFError,
           BoundsError => :PyExc_IndexError,
           DivideError => :PyExc_ZeroDivisionError,
           DomainError => :PyExc_RuntimeError,
           OverflowError => :PyExc_OverflowError,
           InexactError => :PyExc_ArithmeticError,
           MemoryError => :PyExc_MemoryError,
           StackOverflowError => :PyExc_MemoryError,
           UndefRefError => :PyExc_RuntimeError,
           InterruptException => :PyExc_KeyboardInterrupt,
           PyIOError => :PyExc_IOError]
    for (k,v) in exc
        p = convert(Ptr{PyPtr}, pysym_e(v))
        if p != C_NULL
            (pyexc::Dict)[k] = unsafe_load(p)
        end
    end
end

function pyexc_finalize()
    global pyexc
    pyexc::Dict = Dict{DataType, PyPtr}()
end

function pyraise(e)
    global pyexc
    eT = typeof(e)
    pyeT = haskey(pyexc::Dict, eT) ? (pyexc::Dict)[eT] : (pyexc::Dict)[Exception]
    ccall((@pysym :PyErr_SetString), Void, (PyPtr, Ptr{Uint8}),
          pyeT, bytestring(string("Julia exception: ", e)))
end

function pyraise(e::PyError)
    ccall((@pysym :PyErr_Restore), Void, (PyPtr, PyPtr, PyPtr),
          e.T, e.val, e.traceback)
    e.T.o = e.val.o = e.traceback.o = C_NULL # refs were stolen
end


text:
module Oolong

const OL = Oolong
export OL

include("config.jl")
include("logging.jl")
include("base.jl")
include("core/core.jl")
include("start.jl")

end


text:
# Common Utils

function resid2DLinear(μ, mus, Lambdas; diffop::Function=-)  # '-' exploits EuclideanManifold commutativity a-b = b-a
  dμ = broadcast(diffop, μ, mus)  # mus .- μ  ## μ .\ mus
  # @show round.(dμ, digits=4)
  ret = sum( Lambdas.*dμ )
  return ret
end

function solveresid2DLinear!(res, x, mus, Lambdas; diffop::Function=-)::Nothing
  res[1] = resid2DLinear(x, mus, Lambdas, diffop=diffop)
  nothing
end

# import ApproxManifoldProducts: resid2DLinear, solveresid2DLinear
function solveresid2DLinear(res, x, mus, Lambdas; diffop::Function=-)::Float64
  solveresid2DLinear!(res, x, mus, Lambdas, diffop=diffop)
  return res[1]
end

"""
    $SIGNATURES

A clunky repeat calculation of one product kernel.
"""
function updateProductSample( dest::BallTreeDensity,
                              proposals::Vector{BallTreeDensity},
                              manifolds::Tuple,
                              smplIdx::Int,
                              labels::Vector{Int})
  #

  Ndens = length(proposals)
  Ndim = Ndim(dest)

  densLen = Npts.(proposals)

  calclambdas = zeros(Ndim, Ndens)
  calcmu = zeros(Ndim, Ndens)
  destMu = zeros(Ndim)
  destCov = 0.0

  @inbounds @fastmath @simd for dim in 1:Ndim
    for j in 1:Ndens
      calclambdas[dim,j] = 1.0/getBW(proposals[j])[dim,labels[j]]
      calcmu[dim,j] = getPoints(proposals[j])[dim,labels[j]]
    end
    destCov = getLambda(calclambdas)
    destCov = 1.0/destCov
    # μ = 1/Λ * Λμ  ## i.e. already scaled to mean only
    destMu[dim] = getMu(calcmu[dim, :], calclambdas[dim, :], destCov)
  end

  # previous points
  pts = getPoints(dest)
  pts[:,smplIdx] = destMu

  manikde!(pts, manifolds)
end

# Returns the covariance (square), not deviation
function calcCovarianceBasic(M::AbstractManifold, ptsArr::Vector{P}) where P
  #TODO double check the maths,. it looks like its working at least for groups
  μ = mean(M, ptsArr)
  Xcs = vee.(Ref(M), Ref(μ), log.(Ref(M), Ref(μ), ptsArr))
  Σ = mean(Xcs .* transpose.(Xcs))
  @debug "calcCovarianceBasic" μ
  @debug "calcCovarianceBasic" Σ
  # TODO don't know what to do here so keeping as before, #FIXME it will break
  # a change between this and previous is that a full covariance matrix is returned
  msst = Σ
  msst_ = 0 < sum(1e-10 .< msst) ? maximum(msst) : 1.0
  return msst_
end

"""
    $SIGNATURES

Calculate covariance weighted mean as product of incoming Gaussian points ``μ_`` and coordinate covariances ``Σ_``.

Notes
- Return both weighted mean and new covariance (teh congruent product)
- More efficient helper function allows passing keyword inverse covariances `Λ_` instead.
- Assume `size(Σ_[1],1) == manifold_dimension(M)`.
- calc lambdas first and use to calculate mean product second.
- https://ccrma.stanford.edu/~jos/sasp/Product_Two_Gaussian_PDFs.html
- Pennec, X. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements, HAL Archive, 2011, Inria, France.
"""
function calcProductGaussians(M::AbstractManifold,
                              μ_::Union{<:AbstractVector{P},<:NTuple{N,P}}, # point type commonly known as P
                              Σ_::Union{Nothing,<:AbstractVector{S},<:NTuple{N,S}};
                                dim::Integer=manifold_dimension(M),
                              Λ_ = inv.(Σ_),
                              ) where {N,P,S<:AbstractMatrix{<:Real}}
  #
  # calc sum of covariances
  Λ = zeros(MMatrix{dim,dim})
  Λ .= sum(Λ_)

  # Tangent space reference around the evenly weighted mean of incoming points
  u0 = mean(M, μ_)

  # calc the covariance weighted delta means of incoming points and covariances
  ΛΔμ = zeros(MVector{dim})
  for (s,u) in zip(Λ_, μ_)
    # require vee as per Pennec, Caesar Ref [3.6]
    Δuvee = vee(M, u0, log(M, u0, u))
    ΛΔμ += s*Δuvee
  end

  # calculate the delta mean
  Δμ = Λ \ ΛΔμ

  # return new mean and covariance
  return exp(M, u0, hat(M, u0, Δμ)), inv(Λ)
end

# additional support case where covariances are passed as diagonal-only vectors
# still pass nothing, to avoid stack overflow.  Only Λ_ is needed further
calcProductGaussians( M::AbstractManifold,
                      μ_::Union{<:AbstractVector{P},<:NTuple{N,P}},
                      Σ_::Union{<:AbstractVector{S},<:NTuple{N,S}};
                        dim::Integer=manifold_dimension(M),
                      Λ_ = map(s->diagm( 1.0 ./ s), Σ_),
                      ) where {N,P,S<:AbstractVector} = calcProductGaussians(M, μ_, nothing; dim=dim, Λ_=Λ_ )
#

calcProductGaussians( M::AbstractManifold,
                      μ_::Union{<:AbstractVector{P},<:NTuple{N,P}};
                        dim::Integer=manifold_dimension(M),
                      Λ_ = diagm.( (1.0 ./ μ_) ),
                      ) where {N,P} = calcProductGaussians(M, μ_, nothing; dim=dim, Λ_=Λ_ )
#

# """
#     $SIGNATURES

# Once a Gibbs product is available, this function can be used to update the product assuming some change to the input
# to some or some or all of the input density kernels.

# Notes
# - This function does not resample a new posterior sample pairing of inputs, only updates with existing
# """
# function _updateMetricTreeDensityProduct( npd0::BallTreeDensity,
#                                           trees::Array{BallTreeDensity,1},
#                                           anFcns,
#                                           anParams;
#                                           Niter::Int=3,
#                                           addop::Tuple=(+,),
#                                           diffop::Tuple=(-,),
#                                           getMu::Tuple=(getEuclidMu,),
#                                           getLambda::T4=(getEuclidLambda,),
#                                           glbs = makeEmptyGbGlb(),
#                                           addEntropy::Bool=true )
#   #

# end

#


text:
using Pkg

Pkg.activate(".")

import Accounts: Domain, API, Infrastructure

using .Domain, .API, .Infrastructure

const FILE_ACCOUNTS = "./test_accounts.txt"

donald_email = create(EMAIL, "donald@duckcity.com")
donald_work = create(WORK,
  """
  Donalds Hardware Store
  attn. Donald Duck
  1190 Seven Seas Dr
  FL 32830 Lake Buena Vista
  USA
  """
)

addresses = [donald_email, donald_work]

donald = create("Donald Duck", addresses)

add_to_file(FILE_ACCOUNTS, [donald])

result = read_from_file(FILE_ACCOUNTS)

println(result)


text:
struct SubProblem{F, T <: Real, A <: AbstractVector, S <: LQSolver}
    id::Int
    π::T

    solver::S
    feasibility_solver::S

    h::Tuple{A,A}
    x::A
    y::A
    masterterms::Vector{Tuple{Int,Int,T}}

    function (::Type{SubProblem})(model::JuMP.Model, id::Integer, π::Real, x::AbstractVector, y₀::AbstractVector, masterterms::Vector{Tuple{Int,Int,R}}, optimsolver::MPB.AbstractMathProgSolver, F::Bool) where R <: Real
        T = promote_type(eltype(x), eltype(y₀), R, Float32)
        x_ = convert(AbstractVector{T}, x)
        y₀_ = convert(AbstractVector{T}, y₀)
        masterterms_ = convert(Vector{Tuple{Int,Int,T}}, masterterms)
        A = typeof(x_)
        solver = LQSolver(model,optimsolver)
        subproblem = new{F,T,A,typeof(solver)}(id,
                                               π,
                                               solver,
                                               LQSolver(model, optimsolver,Val{F}()),
                                               (convert(A, MPB.getconstrLB(solver.lqmodel)),
                                                convert(A, MPB.getconstrUB(solver.lqmodel))),
                                               x_,
                                               y₀_,
                                               masterterms_
                                               )
        return subproblem
    end

    function (::Type{SubProblem})(model::JuMP.Model, parent::JuMP.Model, id::Integer, π::Real, x::AbstractVector, y₀::AbstractVector, optimsolver::MPB.AbstractMathProgSolver, F::Bool)
        T = promote_type(eltype(x), eltype(y₀), Float32)
        x_ = convert(AbstractVector{T}, x)
        y₀_ = convert(AbstractVector{T}, y₀)
        A = typeof(x_)

        solver = LQSolver(model, optimsolver)

        subproblem = new{F,T,A,typeof(solver)}(id,
                                               π,
                                               solver,
                                               LQSolver(model,optimsolver,Val{F}()),
                                               (convert(A, MPB.getconstrLB(solver.lqmodel)),
                                                convert(A, MPB.getconstrUB(solver.lqmodel))),
                                               x_,
                                               y₀_,
                                               Vector{Tuple{Int,Int,T}}()
                                               )
        parse_subproblem!(subproblem, model, parent)
        return subproblem
    end
end

function parse_subproblem!(subproblem::SubProblem, model::JuMP.Model, parent::JuMP.Model)
    for (i, constr) in enumerate(model.linconstr)
        for (j, var) in enumerate(constr.terms.vars)
            if var.m == parent
                # var is a first stage variable
                push!(subproblem.masterterms, (i,var.col,-constr.terms.coeffs[j]))
            end
        end
    end
end

function update_subproblem!(subproblem::SubProblem{true}, x::AbstractVector)
    lb = MPB.getconstrLB(subproblem.solver.lqmodel)
    ub = MPB.getconstrUB(subproblem.solver.lqmodel)
    for i in [term[1] for term in unique(term -> term[1], subproblem.masterterms)]
        lb[i] = subproblem.h[1][i]
        ub[i] = subproblem.h[2][i]
    end
    for (i,j,coeff) in subproblem.masterterms
        lb[i] += coeff*x[j]
        ub[i] += coeff*x[j]
    end
    MPB.setconstrLB!(subproblem.solver.lqmodel, lb)
    MPB.setconstrUB!(subproblem.solver.lqmodel, ub)
    MPB.setconstrLB!(subproblem.feasibility_solver.lqmodel, lb)
    MPB.setconstrUB!(subproblem.feasibility_solver.lqmodel, ub)
    subproblem.x[:] = x
    return nothing
end
function update_subproblem!(subproblem::SubProblem{false}, x::AbstractVector)
    lb = MPB.getconstrLB(subproblem.solver.lqmodel)
    ub = MPB.getconstrUB(subproblem.solver.lqmodel)
    for i in [term[1] for term in unique(term -> term[1], subproblem.masterterms)]
        lb[i] = subproblem.h[1][i]
        ub[i] = subproblem.h[2][i]
    end
    for (i,j,coeff) in subproblem.masterterms
        lb[i] += coeff*x[j]
        ub[i] += coeff*x[j]
    end
    MPB.setconstrLB!(subproblem.solver.lqmodel, lb)
    MPB.setconstrUB!(subproblem.solver.lqmodel, ub)
    subproblem.x[:] = x
    return nothing
end
update_subproblems!(subproblems::Vector{<:SubProblem}, x::AbstractVector) = map(prob -> update_subproblem!(prob,x), subproblems)

function get_solution(subproblem::SubProblem)
    return copy(subproblem.y), getredcosts(subproblem.solver), getduals(subproblem.solver), getobjval(subproblem.solver)
end

function solve(subproblem::SubProblem)
    subproblem.solver(subproblem.y)
    solvestatus = status(subproblem.solver)
    if solvestatus == :Optimal
        subproblem.y[:] = getsolution(subproblem.solver)
        return OptimalityCut(subproblem)
    elseif solvestatus == :Infeasible
        @warn @sprintf("Subproblem %d is infeasible at the given first-stage variable. Rerun procedure with checkfeas=true to use feasibility cuts.", subproblem.id)
        return Infeasible(subproblem)
    elseif solvestatus == :Unbounded
        @warn @sprintf("Subproblem %d is unbounded at the given first-stage variable", subproblem.id)
        return Unbounded(subproblem)
    else
        error(@sprintf("Subproblem %d was not solved properly, returned status code: %s", subproblem.id,string(solvestatus)))
    end
end

function (subproblem::SubProblem{true})()
    subproblem.feasibility_solver(subproblem.y)
    w = getobjval(subproblem.feasibility_solver)
    if w > 0
        return FeasibilityCut(subproblem)
    end
    return solve(subproblem)
end
function (subproblem::SubProblem{false})()
    return solve(subproblem)
end

function (subproblem::SubProblem)(x::AbstractVector)
    update_subproblem!(subproblem, x)
    subproblem.solver(subproblem.y)
    solvestatus = status(subproblem.solver)
    if solvestatus == :Optimal
        subproblem.y[:] = getsolution(subproblem.solver)
        return getobjval(subproblem.solver)
    elseif solvestatus == :Infeasible
        @warn @sprintf("Subproblem %d is infeasible at the given first-stage variable", subproblem.id)
        return Inf
    elseif solvestatus == :Unbounded
        @warn @sprintf("Subproblem %d is unbounded at the given first-stage variable", subproblem.id)
        return -Inf
    else
        error(@sprintf("Subproblem %d was not solved properly, returned status code: %s", subproblem.id, string(solvestatus)))
    end
end


text:
# ===== VSB Verifications: Vorticity Coupling of VSB and VPM =====

println(" === BEGINNING Verifications7.jl === ")

# === Imports ===
println("   IMPORTING PACKAGES...")
using VSB
using SimpleVPM
using LinearAlgebra
using Printf

cd("/Users/Damyn/Documents/BYU/FLOW Lab/VSB/verifications")

# === Boundary Geometry ===
println("   GENERATING GEOMETRY...")
NPTS = 20                              # Number of Points on Boundary
del_theta = (2*pi)/(NPTS - 1)           # Step in theta around body
R = 1.0                                 # Cylinder Radius
k_hat = [0.0, 0.0, 1.0]                 # Z direction unit vector
X_coor = zeros(NPTS)                    # X-coordinates
Y_coor = zeros(NPTS)                    # Y-coordinates
X_tHat = zeros(NPTS)                    # X Tangent Vector
Y_tHat = zeros(NPTS)                    # X Tangent Vector
X_nHat = zeros(NPTS)                    # Y Normal Vector
Y_nHat = zeros(NPTS)                    # Y Normal Vector
phase = 0.0 * del_theta                 # Phase shift
for i=1:NPTS
    X_coor[i] = R*cos((i-1)*del_theta + phase)
    Y_coor[i] = R*sin((i-1)*del_theta + phase)

    X_tHat[i] = sin((i-1)*del_theta + phase)
    Y_tHat[i] = -cos((i-1)*del_theta + phase)

    X_nHat[i] = cos((i-1)*del_theta + phase)
    Y_nHat[i] = sin((i-1)*del_theta + phase)
end
body_pts = [[X_coor[i], Y_coor[i], 0.0] for i in 1:NPTS]
t_hats = [[X_tHat[i], Y_tHat[i], 0.0] for i in 1:NPTS]
n_hats = [[X_nHat[i], Y_nHat[i], 0.0] for i in 1:NPTS]
body = VSB.Boundary(body_pts,t_hats,n_hats)

# === Set up SimpleVPM ===
println("   SETTING UP SVPM...")
nu = 1e0                                                      # Kinematic Viscosity
dt = 1e0                                                      # Time step (place holder here)
overlap = 2.0                                                 # Particle overlap
sigma = overlap * del_theta                                   # Guassian Spreading
magU_inf = 1.0                                                # Magnitude of Free-Stream Velocity
U_inf(t) = magU_inf * [1.0, 0.0, 0.0]                         # Free-Stream velocity

println("   CREATING PARTICLE FIELD...")
pfield = SimpleVPM.ParticleField(nu,
                                 zeta_gauserf,                # Basis function
                                 g_gauserf;                   # Regularizing function of zeta
                                 transposed=true,             # Transposed scheme
                                 relax=!true,                 # Pedrizzeti relaxation scheme
                                 rlxf=0.3,                    # Relaxation factor
                                 integration="euler",         # Time integration scheme
                                 Uinf=U_inf)                  # Free-Stream velocity function

# === Velocity function ===
function U_slip(X)
    if length(pfield.particles) == 0
        return U_inf(0)
    else
        return SimpleVPM.U(pfield, X) + U_inf(0)
    end
end

# === Calculate parameters ===
println("   CALCULATING PARAMS...")
params = VSB.Parameters(body, U_slip, verbose = false, sigma=sigma)
alphas = params.alphas
betas = VSB.CalcVSDiffusionCoefs(body, alphas, dt, pfield.nu, sigma=sigma)
betas[1] = 0.0
betas[NPTS] = 0.0

println("betas: ",betas)

# === Add particles on surface ===
println("   ADDING PARTICLES TO FIELD...")
offset = 0             # Normal direction offset for new particle locations
for i in 1:NPTS
    X = body_pts[i] + offset .* n_hats[i]
    Gamma = betas[i] .* k_hat
    #Gamma = 1.0 .* k_hat
    vol = 1.0
    particle = SimpleVPM.Particle(X, Gamma, sigma, vol)
    SimpleVPM.addparticle(pfield, particle)
end

SimpleVPM.save(pfield,"pfield_testing")

# === Vorticity calculations ===
println("   BEGIN VORTICITY CALCULATIONS:")
println()
println()
# Method 1: VSB from diffused vortex sheet
omega_VSB(X) = (1/sigma^3)*VSB.CalcDiffusion(body, betas, X, sigma=sigma) .* k_hat

# Method 2: VPM RBF
omega_RBF(X) = SimpleVPM.omega_approx(pfield, X)

# Method 3: Definition of vorticity via VPM
function omega_VPM(i)
    junk, vorticity = SimpleVPM.vortexstretching(pfield, pfield)
    return vorticity[i]
end

# === Results ===
for i in 1:NPTS
    X = body_pts[i]
    XP = pfield.particles[i].X
    VSB = omega_VSB(X)
    RBF = omega_RBF(X)
    VPM = omega_VPM(i)

    println("=== i: ",@sprintf("%3.0i",i)," ======================================")
    println("            X = [",@sprintf("%3.2f",X[1]),", ",@sprintf("%3.2f",X[2]),"]")
    println("           XP = [",@sprintf("%3.2f",XP[1]),", ",@sprintf("%3.2f",XP[2]),"]")

    # println("     gamma[i]")
    println("      beta[i] = ",betas[i] .* k_hat)
    println("     Gamma[i] = ",pfield.particles[i].Gamma)

    println("    omega_VSB = [",@sprintf("%6.4f",VSB[1]),", ",@sprintf("%6.4f",VSB[2]),", ",@sprintf("%6.4f",VSB[3]),"]")
    println("    omega_RBF = [",@sprintf("%6.4f",RBF[1]),", ",@sprintf("%6.4f",RBF[2]),", ",@sprintf("%6.4f",RBF[3]),"]")
    println("    omega_VPM = [",@sprintf("%6.4f",VPM[1]),", ",@sprintf("%6.4f",VPM[2]),", ",@sprintf("%6.4f",VPM[3]),"]")

    println("  |omega_VSB| = ",@sprintf("%6.2f",norm(omega_VSB(X))))
    println("  |omega_RBF| = ",@sprintf("%6.2f",norm(omega_RBF(X))))
    println("  |omega_VPM| = ",@sprintf("%6.2f",norm(omega_VPM(i))))
    println()
end


text:
module Delayed

import Base
import Distributions
using DocStringExtensions

import ..IndexObjects
import ..SpacePrint
import ..Types

"""DistDelayed objects represent functions that draw output from distributions"""
abstract type AbstractDistDelayed <: Types.AbstractDelayed end
abstract type LogNormalQuantDist <: AbstractDistDelayed end
abstract type LogUniformQuantDist <: AbstractDistDelayed end
const NestedFloat = Union{Types.AbstractDelayed, Float64}

"""Switch objects represent a choice between options"""
abstract type AbstractSwitch <: Types.AbstractDelayed end

"""Allow Delayed objects to contain other Delayed objects"""
const NestedFloat = Union{Types.AbstractDelayed, Float64}
const NestedInt = Union{Types.AbstractDelayed, Int}
const NestedReal = Union{Types.AbstractDelayed, Real}

include("operators.jl")
include("categoricalindex.jl")
include("float.jl")
include("lognormal.jl")
include("logquantnormal.jl")
include("quantlognormal.jl")
include("normal.jl")
include("quantnormal.jl")
include("params.jl")
include("randindex.jl")
include("uniform.jl")
include("loguniform.jl")
include("quantuniform.jl")
include("logquantuniform.jl")
include("quantloguniform.jl")

function SpacePrint.spaceprint(
    item::Types.AbstractDelayed; index::Int = 1, tab::String = "", corner::String = "",
    final::Bool = true
)::Nothing
    println(tab, corner, index, ": ", typeof(item))
    for (i, propertyname) in enumerate(propertynames(item))
        finalitem = i == length(propertynames(item))
        println(tab, SpacePrint.fillerchar(final), "├", propertyname)
        SpacePrint.spaceprint(
            getproperty(item, propertyname); index = 1,
            tab = string(tab, SpacePrint.fillerchar(final)),
            corner = SpacePrint.cornerchar(finalitem), final = finalitem
        )
    end

    return nothing
end

end # module Delayed


text:
@testset "Params delete!" begin
  w = rand(2,3)
  b = rand(2)
  ps = Params([w,b])
  delete!(ps, b)
  @test length(ps.order) == length(ps.params) == 1
  @test first(ps.order) == first(ps.params) == w
end


text:
using Genie.Configuration, Logging
server_sett = TOML.parsefile("config/env/settings.toml")
port = server_sett["port"]
host = server_sett["host"]
admin = server_sett["admin"]
const config = Settings(
  server_port                     = Int64(port),
  server_host                     = host,
  log_level                       = Logging.Debug,
  log_to_file                     = true,
  server_handle_static_files      = true
)

ENV["JULIA_REVISE"] = "off"


text:
# broadcasting

using Base.Broadcast: BroadcastStyle, Broadcasted

struct CuArrayStyle{N} <: AbstractGPUArrayStyle{N} end
CuArrayStyle(::Val{N}) where N = CuArrayStyle{N}()
CuArrayStyle{M}(::Val{N}) where {N,M} = CuArrayStyle{N}()

BroadcastStyle(::Type{<:CuArray{T,N}}) where {T,N} = CuArrayStyle{N}()

Base.similar(bc::Broadcasted{CuArrayStyle{N}}, ::Type{T}) where {N,T} =
    similar(CuArray{T}, axes(bc))

Base.similar(bc::Broadcasted{CuArrayStyle{N}}, ::Type{T}, dims) where {N,T} =
    CuArray{T}(undef, dims)

# broadcasting type ctors isn't GPU compatible
Broadcast.broadcasted(::CuArrayStyle{N}, f::Type{T}, args...) where {N, T} =
    Broadcasted{CuArrayStyle{N}}((x...) -> T(x...), args, nothing)


text:
function gradient!(out::AbstractArray, f_at_face, f_at_bface, bf2c,bfn,bcv,f2c,fn,cv,floop::Face2CellLoop)
    fill!(out,zero(eltype(out)))
    NBF = nbfaces(floop)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        flux = bfn[i] * f_at_bface[i] * bcv[i]
        out[j] += flux
    end

    NF = nfaces(floop)
    @inbounds for i=1:NF
        j1,j2 = f2c[i]
        flux = fn[i] * f_at_face[i]
        v = cv[i]
        out[j1] += flux*v[1]
        out[j2] -= flux*v[2]
    end
end

@inline gradient!(out, ff, fbf, floop) = gradient!(out, ff, fbf, floop.bf2c,floop.bfn,floop.bcv,floop.f2c,floop.fn,floop.cv,floop)

@inline function gradient(ff,fbf,mesh)
    out = zeros(eltype(mesh.nodes),length(mesh.cells))
    gradient!(out,ff,fbf,mesh.f2cloops)
    return out
end

function div!(out::Vector, f_at_face, f_at_bface, bf2c, bfn, bcv, f2c, fn, cv, floops::Face2CellLoop)
    fill!(out,zero(eltype(out)))
    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        el = bf2c[i]
        j = el[1]
        n = bfn[i]
        flux = (f_at_bface[i] ⋅ n) * bcv[i]
        out[j] += flux
    end

    NF = nfaces(floops)
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        n = fn[i]
        flux = f_at_face[i] ⋅ n
        v = cv[i]
        out[j1] += flux * v[1]
        out[j2] -= flux * v[2]
    end
end

@inline div!(out, ff, fbf, floop) = div!(out, ff, fbf, floop.bf2c,floop.bfn,floop.bcv,floop.f2c,floop.fn,floop.cv,floop)


text:
module WikiSort

const DEFAULT_CACHE_SIZE = 512

immutable WikiSorter{T}
    cache::Vector{T}
    cache_size::Int
end

function call{T}(::Type{WikiSorter{T}}, cache_size::Integer=DEFAULT_CACHE_SIZE)
    return WikiSorter{T}(fill(zero(T), cache_size), cache_size)
end

function reverse!(arr::AbstractVector, r::UnitRange)
    for i = (div(length(r), 2) - 1):-1:0
        arr[r.start + i], arr[r.stop - i] = arr[r.stop - i], arr[r.start + i]
    end
end

function block_swap!(arr::AbstractVector, start_idx1::Integer, start_idx2::Integer, count::Integer)
    for i = 1:count
        # swap
        arr[start_idx1 + i], arr[start_idx2 + i] = arr[start_idx2 + i], arr[start_idx1 + i]
    end
end

# rotate the values in an array ([0 1 2 3] becomes [1 2 3 0] if we rotate by 1)
# this assumes that 0 <= amount <= length(r)
function rotate!{T}(s::WikiSorter{T}, arr::AbstractVector{T}, r::UnitRange, amount::Integer, use_cache=true)
    if length(r) == 0
        return
    end

    if amount >= 0
        split = r.start + amount
    else
        split = r.stop + amount + 1
    end

    range1 = r.start:(split - 1)
    range2 = split:r.stop

    if use_cache
        if length(range1) <= length(range2)
            if length(range1) <= s.cache_size
                # if cache not null
                # unsafe_copy!
                copy!(s.cache, 1, arr, range1.start, length(range1))
                copy!(arr, range1.start, arr, range2.start, length(range2))
                copy!(arr, range1.start + length(range2), s.cache, 1, length(range1))
                return
            end
        else
            if length(range2) <= s.cache_size
                # ditto
                copy!(s.cache, 1, arr, range2.start, length(range2))
                copy!(arr, range2.stop - length(range1) + 1, arr, range1.start, length(range1))
                copy!(arr, range1.start, s.cache, 1, length(range2))
                return
            end
        end
    end

    reverse!(arr, range1)
    reverse!(arr, range2)
    reverse!(arr, r)
    return
end

function rotate!(arr::AbstractVector, r::UnitRange, amount::Integer)
    Base.reverse!(arr, r.start, r.start + amount - 1)
    Base.reverse!(arr, r.start + amount, r.stop)
    Base.reverse!(arr, r.start, r.stop)
end

function binaryfirst(arr::AbstractVector, r::UnitRange, value)
    start = r.start
    stop = r.stop
    while start < stop
        mid = start + div(stop - start, 2)
        if arr[mid] < value
            start = mid + 1
        else
            stop = mid
        end
    end

    if start == r.stop && arr[start] < value
        start += 1
    end

    return start
end

function binarylast(arr::AbstractVector, r::UnitRange, value)
    start = r.start
    stop = r.stop
    while start < stop
        mid = start + div(stop - start, 2)
        if arr[mid] <= value
            start = mid + 1
        else
            stop = mid
        end
    end

    if start == r.stop && arr[start] <= value
        start += 1
    end

    return start
end

function insertion_sort!(arr::AbstractVector, r::UnitRange)
    for i = (r.start + 1):r.stop
        temp = arr[i]
        j = i
        while j > r.start && arr[j - 1] > temp
            arr[j] = arr[j - 1]
            j -= 1
        end
        arr[j] = temp
    end
end

end # module


text:
@testset "Embedding operators" begin

@testset "RandomBound" begin
    @testset "does nothing if within bounds" begin
        @test apply!(RandomBound(RectSearchSpace(1, (0.0, 1.0))), [0.0], [0.0]) == [0.0]

        @test apply!(RandomBound(RectSearchSpace([(0.0, 1.0), (10.0, 15.0)])),
                     [0.0, 11.4], [0.1, 12.3] ) == [0.0, 11.4]
    end

    @testset "bounds if lower than min bound" begin
        @test apply!(RandomBound(RectSearchSpace([(0.0, 1.0)])),
                     [-0.1], [0.0]) == [0.0]

        res = apply!(RandomBound(RectSearchSpace([(0.0, 1.0)])),
                     [-0.1], [0.5])
        @test 0.0 <= res[1] <= 0.5

        res = apply!(RandomBound(RectSearchSpace([(-11.0, 1.0), (0.0, 1.0)])),
                     [-11.1, 0.5], [-10.8, 0.5])
        @test -10.8 >= res[1] >= -11.0
        @test res[2] == 0.5

        res = apply!(RandomBound(RectSearchSpace([(30.0, 60.0), (-102.0, -1.0)])),
                     [50.4, -103.1], [49.6, -101.4])
        @test res[1] == 50.4
        @test -101.4 >= res[2] >= -102.0
    end

    @testset "bounds if higher than max bound" begin
        @test apply!(RandomBound(RectSearchSpace([(0.0, 1.0)])),
                     [1.1], [1.0]) == [1.0]

        res = apply!(RandomBound(RectSearchSpace([(-10.0, 96.0)])), [97.0], [95.0])
        @test 95.0 <= res[1] <= 96.0
    end

    @testset "MixedPrecisionRectSearchSpace" begin
        ss = RectSearchSpace([(0.0, 1.0), (1.0, 2.0), (2.0, 3.0)], dimdigits=[-1, 0, 2])

        @test apply!(RandomBound(ss), [0.0, 1.4, 2.1], [0.0, 1.45, 2.0]) == [0.0, 1.0, 2.1]
        @test apply!(RandomBound(ss), [0.0, 5.0, 0.123], [0.0, 1.5, 2.0]) == [0.0, 2.0, 2.0]
        @test apply!(RandomBound(ss), [1.12345, 1.12345, 2.1234], [1.0, 2.0, 3.0]) == [1.0, 1.0, 2.12]
    end
end

end


text:
using BinaryProvider

if arch(platform_key_abi()) != :x86_64
    @warn "FastTransforms has only been tested on x86_64 architectures."
end

ft_build_from_source = get(ENV, "FT_BUILD_FROM_SOURCE", "false")

if ft_build_from_source == "true"
    println("Building from source.")
    if Sys.isapple()
        script = raw"""
        brew update
        brew install gcc@8 fftw mpfr
        rm -rf FastTransforms
        git clone -b v0.2.7 https://github.com/MikaelSlevinsky/FastTransforms.git FastTransforms
        cd FastTransforms
        make lib CC=gcc-8 FT_USE_APPLEBLAS=1
        cd ..
        ln -sf FastTransforms/libfasttransforms.dylib libfasttransforms.dylib
        """
        run(`/bin/bash -c $(script)`)
    elseif Sys.islinux()
        script = raw"""
        sudo apt-get update
        sudo add-apt-repository ppa:ubuntu-toolchain-r/test
        sudo apt-get install gcc-8 libblas-dev libopenblas-base libfftw3-dev libmpfr-dev
        rm -rf FastTransforms
        git clone -b v0.2.7 https://github.com/MikaelSlevinsky/FastTransforms.git FastTransforms
        cd FastTransforms
        make lib CC=gcc-8
        cd ..
        ln -sf FastTransforms/libfasttransforms.so libfasttransforms.so
        """
        run(`/bin/bash -c $(script)`)
    else
        @warn "FastTransforms could not be built from source with the current build.jl script. Have you considered filing an issue? https://github.com/JuliaApproximation/FastTransforms.jl/issues"
    end
else
    println("Installing by downloading binaries.")
    if Sys.isapple()
        run(`brew update`)
        run(`brew install gcc@8 fftw mpfr`)
        const libfasttransforms = joinpath(dirname(@__DIR__), "deps", "libfasttransforms.dylib")
        GCC = BinaryProvider.detect_compiler_abi().gcc_version
        println("Building with ", GCC, ".")
        const release = "https://github.com/MikaelSlevinsky/FastTransforms/releases/download/v0.2.7/libfasttransforms.v0.2.7"
        if GCC == :gcc4
            download(release*".gcc-4.9.dylib", libfasttransforms)
        elseif GCC == :gcc5
            download(release*".gcc-5.dylib", libfasttransforms)
        elseif GCC == :gcc6
            download(release*".gcc-6.dylib", libfasttransforms)
        elseif GCC == :gcc7
            download(release*".gcc-7.dylib", libfasttransforms)
        elseif GCC == :gcc8
            download(release*".gcc-8.dylib", libfasttransforms)
        elseif GCC == :gcc9
            download(release*".gcc-9.dylib", libfasttransforms)
        else
            @warn "Please ensure you have a version of gcc from gcc-4.9 to gcc-9."
        end
    elseif Sys.islinux()
        run(`sudo add-apt-repository ppa:ubuntu-toolchain-r/test`)
        run(`sudo apt-get update`)
        run(`sudo apt-get install gcc-8 libblas-dev libopenblas-base libfftw3-dev libmpfr-dev`)
        const libfasttransforms = joinpath(dirname(@__DIR__), "deps", "libfasttransforms.so")
        GCC = BinaryProvider.detect_compiler_abi().gcc_version
        println("Building with ", GCC, ".")
        const release = "https://github.com/MikaelSlevinsky/FastTransforms/releases/download/v0.2.7/libfasttransforms.v0.2.7"
        if GCC == :gcc4
            download(release*".gcc-4.9.so", libfasttransforms)
        elseif GCC == :gcc5
            download(release*".gcc-5.so", libfasttransforms)
        elseif GCC == :gcc6
            download(release*".gcc-6.so", libfasttransforms)
        elseif GCC == :gcc7
            download(release*".gcc-7.so", libfasttransforms)
        elseif GCC == :gcc8
            download(release*".gcc-8.so", libfasttransforms)
        elseif GCC == :gcc9
            download(release*".gcc-9.so", libfasttransforms)
        else
            @warn "Please ensure you have a version of gcc from gcc-4.9 to gcc-9."
        end
    else
        @warn "FastTransforms could not be installed by downloading binaries. Have you tried building from source?"
    end
end


text:
using KLU
using KLU: increment!, KLUITypes, decrement, klu!, KLUFactorization,
klu_analyze!, klu_factor!
using Test
using SparseArrays: SparseMatrixCSC, sparse, nnz
using LinearAlgebra

@testset "KLU Wrappers" begin
    Ap = increment!([0,4,1,1,2,2,0,1,2,3,4,4])
    Ai = increment!([0,4,0,2,1,2,1,4,3,2,1,2])
    Ax = [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.]
    A0 = sparse(Ap, Ai, Ax)
    A1 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
                increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
                [2.,1.,3.,4.,-1.,-3.,3.,9.,2.,1.,4.,2.], 5, 5)
    @testset "Core functionality for $Tv elements" for Tv in (Float64, ComplexF64)
        @testset "Core functionality for $Ti indices" for Ti ∈ Base.uniontypes(KLUITypes)
            A = convert(SparseMatrixCSC{Tv, Ti}, A0)
            # test the raw vector construction method.
            klua = klu(size(A, 1), decrement(A.colptr), decrement(A.rowval), A.nzval)
            @test klua.F == klu(A).F
            @test nnz(klua) == 18
            R = Diagonal(Tv == ComplexF64 ? complex.(klua.Rs) : klua.Rs)
            @test R \ A[klua.p, klua.q] ≈ (klua.L * klua.U + klua.F)

            b = [8., 45., -3., 3., 19.]
            x = klua \ b
            @test x ≈ float([1:5;])
            @test A*x ≈ b

            z = complex.(b)
            x = ldiv!(klua, z)
            @test x ≈ float([1:5;])
            @test z === x
            # Can't match UMFPACK's ldiv!(<OUTPUT>, <KLU>, <INPUT>)
            # since klu_solve(<KLU>, <INPUT>) modifies <INPUT>, and has no field for <OUTPUT>.
            @test A*x ≈ b

            b = [8., 20., 13., 6., 17.]
            x = klua'\b
            @test x ≈ float([1:5;])

            @test A'*x ≈ b
            z = complex.(b)
            x = ldiv!(adjoint(klua), z)
            @test x ≈ float([1:5;])
            @test x === z

            @test A'*x ≈ b
            x = transpose(klua) \ b
            @test x ≈ float([1:5;])
            @test transpose(A) * x ≈ b

            x = ldiv!(transpose(klua), complex.(b))
            @test x ≈ float([1:5;])
            @test transpose(A) * x ≈ b

            @inferred klua\fill(1, size(A, 2))

            @testset "Permutation vectors" begin
                #Just to test this works, we'll use the existing permutation vectors.
                @test klu_analyze!(klua, klua.p, klua.q).common.status == 0
                klu_factor!(klua)
                x = klua \ b
                @test A*x ≈ b
            end
            @testset "Utility functions" begin
                K = KLUFactorization(A);
                @test size(K) == (5, 5)
                @test size(K, 3) == 1
                @test_throws ArgumentError K.symbolic
                @test_throws ArgumentError K.numeric
                klu_analyze!(K);
                @test K.symbolic.nz == 12
                @test K.symbolic.nzoff == 4
                klu_factor!(K)
                @test K.numeric.lnz == 7 == K.numeric.unz
                @test K.nblocks == 3
                @test propertynames(K) == (:lnz, :unz, :nzoff, :L, :U, :F, :q, :p, :Rs, :symbolic, :numeric,)
                @test_throws ArgumentError size(K, -1)
                @test size(K, 3) == 1
            end
            @testset "Refactorization" begin
                B = convert(SparseMatrixCSC{Tv, Ti}, A1)
                b = Tv[8., 45., -3., 3., 19.]
                F = klu(A)
                klu!(F, B)
                @test F\b ≈ B\b ≈ Matrix(B)\b
                @test klu!(F, B.nzval)\b ≈ B\b ≈ Matrix(B)\b #test just supplying nzval for recompute
            end

            @testset "Singular matrix" begin
                S = sparse([1 2; 0 0])
                S = convert(SparseMatrixCSC{Tv, Ti}, S)
                @test_throws SingularException klu(S)
            end
        end
    end
    @testset "REPL printing of KLU" begin
        A = sparse([1, 2], [1, 2], Float64[1.0, 1.0])
        F = klu(A)
        facstring = sprint((t, s) -> show(t, "text/plain", s), F)
        lstring = sprint((t, s) -> show(t, "text/plain", s), F.L)
        ustring = sprint((t, s) -> show(t, "text/plain", s), F.U)
        fstring = sprint((t, s) -> show(t, "text/plain", s), F.F)
        @test facstring == "$(summary(F))\nL factor:\n$lstring\nU factor:\n$ustring\nF factor:\n$fstring"
    end
end

@testset "Issue #4" begin
    A = SparseMatrixCSC(15, 15,
    [1, 8, 12, 16, 21, 22, 27, 34, 37, 39, 41, 44, 47, 50, 55, 58],
    [3, 4, 6, 8, 11, 13, 15, 2, 6, 8, 15, 1, 2, 5, 12, 3, 5, 6, 7, 9, 14, 7, 10, 11, 13, 15, 4, 5, 6, 8, 10, 12, 13, 8, 12, 15, 2, 13, 1, 7, 3, 4, 14, 2, 4, 5, 1, 5, 9, 4, 7, 9, 11, 14, 13, 14, 15],
    [0.2775474841561938, 0.19549953706849155, 0.7221976371086005, 0.4339373082200655, 0.983079431343046, 0.10918778088879799, 0.3728676112188065, 0.9134045061777432, 0.14560891622463457, 0.7715210431553383, 0.2945501295372417, 0.6134722502122134, 0.8777181195348973, 0.6995382425541914, 0.9562490972786235, 0.27001502642215325, 0.8573661029146233, 0.13020432565448115, 0.9221068910751316, 0.17087414970038983, 0.7062975193151109, 0.7668596005709167, 0.46967704631299334, 0.31764226298560083, 0.39054386892157833, 0.36610203401046015, 0.16689896372140534, 0.9624322297755521, 0.1478381603984824, 0.45423514524961806, 0.5564610482579242, 0.1844671322175948, 0.0823893170743808, 0.25409993152712307, 0.10475245199943273, 0.5863595004922162, 0.14733912690513562, 0.6504152422320895, 0.4339054908933866, 0.27614384058497166, 0.4019619228414033, 0.8631491210976487, 0.20159747073826084, 0.3273367915690062, 0.23866880928640288, 0.9557759456784265, 0.016351125161178537, 0.5320355909884844, 0.9010930260468242, 0.3780686420068593, 0.6375477164214856, 0.9850645305956225, 0.5366242762582065, 0.08835652070698918, 0.9877090693717305, 0.9775298646022268, 0.9759511830494418]
    )
    @test klu(A).F.nzval ≈ [0.20623152093700403, 0.09038754099133628, 1.0]
end


text:
#=
This file is auto-generated. Do not edit.
=#

"""As implemented in Milano's Book, Page 397"""
mutable struct VSCDCLine <: DCBranch
    name::String
    available::Bool
    activepower_flow::Float64
    arc::Arc
    rectifier_taplimits::Min_Max
    rectifier_xrc::Float64
    rectifier_firing_angle::Min_Max
    inverter_taplimits::Min_Max
    inverter_xrc::Float64
    inverter_firing_angle::Min_Max
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function VSCDCLine(name, available, activepower_flow, arc, rectifier_taplimits, rectifier_xrc, rectifier_firing_angle, inverter_taplimits, inverter_xrc, inverter_firing_angle, _forecasts=InfrastructureSystems.Forecasts(), )
    VSCDCLine(name, available, activepower_flow, arc, rectifier_taplimits, rectifier_xrc, rectifier_firing_angle, inverter_taplimits, inverter_xrc, inverter_firing_angle, _forecasts, InfrastructureSystemsInternal())
end

function VSCDCLine(; name, available, activepower_flow, arc, rectifier_taplimits, rectifier_xrc, rectifier_firing_angle, inverter_taplimits, inverter_xrc, inverter_firing_angle, _forecasts=InfrastructureSystems.Forecasts(), )
    VSCDCLine(name, available, activepower_flow, arc, rectifier_taplimits, rectifier_xrc, rectifier_firing_angle, inverter_taplimits, inverter_xrc, inverter_firing_angle, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function VSCDCLine(::Nothing)
    VSCDCLine(;
        name="init",
        available=false,
        activepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        rectifier_taplimits=(min=0.0, max=0.0),
        rectifier_xrc=0.0,
        rectifier_firing_angle=(min=0.0, max=0.0),
        inverter_taplimits=(min=0.0, max=0.0),
        inverter_xrc=0.0,
        inverter_firing_angle=(min=0.0, max=0.0),
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get VSCDCLine name."""
get_name(value::VSCDCLine) = value.name
"""Get VSCDCLine available."""
get_available(value::VSCDCLine) = value.available
"""Get VSCDCLine activepower_flow."""
get_activepower_flow(value::VSCDCLine) = value.activepower_flow
"""Get VSCDCLine arc."""
get_arc(value::VSCDCLine) = value.arc
"""Get VSCDCLine rectifier_taplimits."""
get_rectifier_taplimits(value::VSCDCLine) = value.rectifier_taplimits
"""Get VSCDCLine rectifier_xrc."""
get_rectifier_xrc(value::VSCDCLine) = value.rectifier_xrc
"""Get VSCDCLine rectifier_firing_angle."""
get_rectifier_firing_angle(value::VSCDCLine) = value.rectifier_firing_angle
"""Get VSCDCLine inverter_taplimits."""
get_inverter_taplimits(value::VSCDCLine) = value.inverter_taplimits
"""Get VSCDCLine inverter_xrc."""
get_inverter_xrc(value::VSCDCLine) = value.inverter_xrc
"""Get VSCDCLine inverter_firing_angle."""
get_inverter_firing_angle(value::VSCDCLine) = value.inverter_firing_angle
"""Get VSCDCLine _forecasts."""
get__forecasts(value::VSCDCLine) = value._forecasts
"""Get VSCDCLine internal."""
get_internal(value::VSCDCLine) = value.internal


text:
# Load the circulation and grid
grd, T_Circulation_unit = Circulation.load()
T_Circulation = ustrip.(T_Circulation_unit) # strip units for now

# Define useful constants and arrays
iwet = indices_of_wet_boxes(grd)
nb = number_of_wet_boxes(grd)
v = volumevec(grd)
z = depthvec(grd)
ztop = topdepthvec(grd) # strip units for now
# And matrices
DIV = DIVO(grd)
Iabove = buildIabove(grd)

@testset "Circulation and grid" begin
    @testset "wet3D" begin
        wet3D = iswet(grd)
        @test wet3D isa BitArray{3}
        @test wet3D == grd.wet3D
    end
    @testset "grid" begin
        @test grd isa OceanGrid
    end
    @testset "Transport matrix" begin
        T = T_Circulation
        @test T isa SparseMatrixCSC
        println("Divergence timescale: ", norm(ones(nb)) / norm(T * ones(nb)) * u"s" |> u"Myr")
        println("Mass conservation timescale: ", norm(v) / norm(T' * v) * u"s" |> u"Myr")
    end
    # TODO add tests for consistency of wet3D, grid, and T
end

@testset "Constants, vectors, and matrices types" begin
    # Test their types
    @testset "iwet (indices of wet boxes)" begin
        @test iwet isa Vector{<:Int}
        @test length(iwet) == nb
    end
    @testset "nb (number of wet boxes)" begin
        @test nb isa Int
        @test nb == sum(vec(grd.wet3D))
    end
    @testset "v (vector of volumes)" begin
        @test v isa Vector{Float64}
        @test length(v) == nb
        @test all(v .≥ 0)
        @test v * u"m^3" == array_of_volumes(grd)[iwet]
    end
    @testset "z (vector of depths of center of boxes)" begin
        @test z isa Vector{<:Real}
        @test length(z) == nb
    end
    @test ztop isa Vector{<:Real}
    @test DIV isa SparseMatrixCSC
    @test Iabove isa SparseMatrixCSC
end


text:
using FastPow, Test

@testset "integer powers" begin
    let x = 3.7, z = 5.2, u = 0.988
        y = sin(x)^6 + z^8.3 + u^9
        @fastpow y2 = sin(x)^6 + z^8.3 + u^9
        @test y ≈ y2
    end

    @test @fastpow(0.73^25) ≈ big"0.73"^25 rtol=3e-15
    @test @fastpow(0.6^-7) ≈ big"0.6"^-7 rtol=1e-15
    @test @fastpow(0.6^0) === 1.0
    @test @fastpow(0.6^1) === 0.6
end

struct DummyType; end
Base.:^(x::Int, ::DummyType) = "$x^dummy"
const dummy = DummyType()

@testset "1^p" begin
    for p in (-10,-9,0,9,10,0.2,2//3,0.3+0.4im)
        @test @fastpow(1^p) === one(p)
    end
    @test @fastpow(1^dummy) == "1^dummy"
end

@testset "(-1)^p" begin
    for n in (-10,-9,0,9,10)
        correct = iseven(n) ? 1 : -1
        @test @fastpow((-1)^n) === correct
        @test @fastpow((-1)^float(n)) == correct
        @test @fastpow((-1)^Rational(n)) == correct
    end
    for p in (0.2, 2//3)
        @test_throws DomainError @fastpow((-1)^p)
    end
    z = 0.3 + 0.4im
    @test @fastpow((-1)^z) ≈ (-1)^z ≈ 0.16728929223461353 + 0.23025395732014076im
    @test @fastpow((-1)^dummy) == "-1^dummy"
end

@testset "2^p" begin
    for p in Any[Int64.(0:65); big(100)]
        @test @fastpow(2^p)::typeof(2^p) == 2^p
    end
    p = -2
    @test_throws DomainError @fastpow(2^p)
    @test @fastpow(2^-2) === 0.25
    @test @fastpow(2^-100) === 7.888609052210118e-31
    for p in (0.3, -0.5, 4//3, 0.3+4im)
        @test @fastpow(2^p) == exp2(p)
    end
    @test @fastpow(2^dummy) == "2^dummy"
end

@testset "10^p" begin
    for p in Int64.(0:30)
        @test @fastpow(10^p) === 10^p
    end
    p = -2
    @test_throws DomainError @fastpow(10^p)
    @test @fastpow(10^-2) == 0.01
    for p in (0.3, -0.5, 4//3, 0.3+4im)
        @test @fastpow(10^p) == exp10(p)
    end
    @test @fastpow(10^dummy) == "10^dummy"
end


text:
# Misc data processing functionality

"""
    zeropad(fid, axis, pad)

Zero pad an FID for Fourier interpolation
"""
function zeropad(signal::AxisArray, axis=Axis{:time}, pad=2)
    if pad == 1
        return signal
    else
        dim = axisdim(signal, axis)
        padsize = [size(signal)...]
        padsize[dim] *= (pad-1)
        padded = cat(signal, zeros(padsize...), dims=dim)
        t = AxisArrays.axes(signal, axis).val
        tnew = first(t) .+ (0:size(padded,dim)-1)*step(t)
        newaxes = [AxisArrays.axes(signal)...]
        newaxes[dim] = axis(tnew)
        AxisArray(padded, newaxes...)
    end
end

"""
    simple_averaging(fids)

Do simple averaging from FIDs
"""
function simple_averaging(fids::AxisArray)
    d = axisdim(fids, Axis{:phase_cycle})
    averaged_signal = dropdims(mean(fids, dims=d), dims=d)
end

"""
    spectrum(signal::AxisArray)

Compute the spectrum from time domain signal via Fourier Transform.
"""
function spectrum(signal::AxisArray)
    names = axisnames(signal)
    if names == (:time,)
        spec = fftshift(fft(signal.data))
        f = frequency_axis(signal)
        AxisArray(spec, Axis{:freq}(f))
    elseif names == (:time2, :time1)
        spec = fftshift(fft(signal.data))
        f1 = frequency_axis(signal, Axis{:time1})
        f2 = frequency_axis(signal, Axis{:time2})
        AxisArray(spec, Axis{:freq2}(f2), Axis{:freq1}(f1))
    end
end

# Internal function for downsampling and truncating of acqusition data
function downsample_and_truncate(t, z, cutpre, cutpost, downsample)
    if downsample > 1
        # Do downsampling if requested, via naive Fourier filtering with square
        # window.  This appears to be the way Siemens implement this as well.
        # It's a very spectrum-focused way to do things and is a bit of an
        # abuse in the time domain.
        #
        # TODO: Do something useful if it's not a whole number of samples,
        # instead of an error?
        s1 = fft(z, 1)
        N = size(s1,1)
        n = Int(N/downsample/2)
        s2 = [s1[1:n,:]; s1[end-n+1:end,:]] ./ downsample
        z = ifft(s2, 1)
        t = t[1:downsample:end]
        cutpre  = Int(cutpre/downsample)
        cutpost = Int(cutpost/downsample)
    end
    # Add descriptive axes to the array and remove
    t[cutpre+1:end-cutpost], z[cutpre+1:end-cutpost,:]
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using Multifractals
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@test 1 == 2
@test hello("world") == "Hello, world" #just to show Henri made a change