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

mutable struct IoK8sApiCoreV1PodSpec <: SwaggerModel
    activeDeadlineSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: activeDeadlineSeconds
    affinity::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Affinity } # spec name: affinity
    automountServiceAccountToken::Any # spec type: Union{ Nothing, Bool } # spec name: automountServiceAccountToken
    containers::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Container} } # spec name: containers
    dnsConfig::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodDNSConfig } # spec name: dnsConfig
    dnsPolicy::Any # spec type: Union{ Nothing, String } # spec name: dnsPolicy
    enableServiceLinks::Any # spec type: Union{ Nothing, Bool } # spec name: enableServiceLinks
    hostAliases::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1HostAlias} } # spec name: hostAliases
    hostIPC::Any # spec type: Union{ Nothing, Bool } # spec name: hostIPC
    hostNetwork::Any # spec type: Union{ Nothing, Bool } # spec name: hostNetwork
    hostPID::Any # spec type: Union{ Nothing, Bool } # spec name: hostPID
    hostname::Any # spec type: Union{ Nothing, String } # spec name: hostname
    imagePullSecrets::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1LocalObjectReference} } # spec name: imagePullSecrets
    initContainers::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Container} } # spec name: initContainers
    nodeName::Any # spec type: Union{ Nothing, String } # spec name: nodeName
    nodeSelector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: nodeSelector
    priority::Any # spec type: Union{ Nothing, Int32 } # spec name: priority
    priorityClassName::Any # spec type: Union{ Nothing, String } # spec name: priorityClassName
    readinessGates::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodReadinessGate} } # spec name: readinessGates
    restartPolicy::Any # spec type: Union{ Nothing, String } # spec name: restartPolicy
    runtimeClassName::Any # spec type: Union{ Nothing, String } # spec name: runtimeClassName
    schedulerName::Any # spec type: Union{ Nothing, String } # spec name: schedulerName
    securityContext::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodSecurityContext } # spec name: securityContext
    serviceAccount::Any # spec type: Union{ Nothing, String } # spec name: serviceAccount
    serviceAccountName::Any # spec type: Union{ Nothing, String } # spec name: serviceAccountName
    shareProcessNamespace::Any # spec type: Union{ Nothing, Bool } # spec name: shareProcessNamespace
    subdomain::Any # spec type: Union{ Nothing, String } # spec name: subdomain
    terminationGracePeriodSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: terminationGracePeriodSeconds
    tolerations::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Toleration} } # spec name: tolerations
    volumes::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Volume} } # spec name: volumes

    function IoK8sApiCoreV1PodSpec(;activeDeadlineSeconds=nothing, affinity=nothing, automountServiceAccountToken=nothing, containers=nothing, dnsConfig=nothing, dnsPolicy=nothing, enableServiceLinks=nothing, hostAliases=nothing, hostIPC=nothing, hostNetwork=nothing, hostPID=nothing, hostname=nothing, imagePullSecrets=nothing, initContainers=nothing, nodeName=nothing, nodeSelector=nothing, priority=nothing, priorityClassName=nothing, readinessGates=nothing, restartPolicy=nothing, runtimeClassName=nothing, schedulerName=nothing, securityContext=nothing, serviceAccount=nothing, serviceAccountName=nothing, shareProcessNamespace=nothing, subdomain=nothing, terminationGracePeriodSeconds=nothing, tolerations=nothing, volumes=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("activeDeadlineSeconds"), activeDeadlineSeconds)
        setfield!(o, Symbol("activeDeadlineSeconds"), activeDeadlineSeconds)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("affinity"), affinity)
        setfield!(o, Symbol("affinity"), affinity)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("automountServiceAccountToken"), automountServiceAccountToken)
        setfield!(o, Symbol("automountServiceAccountToken"), automountServiceAccountToken)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("containers"), containers)
        setfield!(o, Symbol("containers"), containers)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("dnsConfig"), dnsConfig)
        setfield!(o, Symbol("dnsConfig"), dnsConfig)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("dnsPolicy"), dnsPolicy)
        setfield!(o, Symbol("dnsPolicy"), dnsPolicy)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("enableServiceLinks"), enableServiceLinks)
        setfield!(o, Symbol("enableServiceLinks"), enableServiceLinks)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("hostAliases"), hostAliases)
        setfield!(o, Symbol("hostAliases"), hostAliases)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("hostIPC"), hostIPC)
        setfield!(o, Symbol("hostIPC"), hostIPC)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("hostNetwork"), hostNetwork)
        setfield!(o, Symbol("hostNetwork"), hostNetwork)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("hostPID"), hostPID)
        setfield!(o, Symbol("hostPID"), hostPID)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("hostname"), hostname)
        setfield!(o, Symbol("hostname"), hostname)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("imagePullSecrets"), imagePullSecrets)
        setfield!(o, Symbol("imagePullSecrets"), imagePullSecrets)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("initContainers"), initContainers)
        setfield!(o, Symbol("initContainers"), initContainers)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("nodeName"), nodeName)
        setfield!(o, Symbol("nodeName"), nodeName)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("nodeSelector"), nodeSelector)
        setfield!(o, Symbol("nodeSelector"), nodeSelector)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("priority"), priority)
        setfield!(o, Symbol("priority"), priority)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("priorityClassName"), priorityClassName)
        setfield!(o, Symbol("priorityClassName"), priorityClassName)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("readinessGates"), readinessGates)
        setfield!(o, Symbol("readinessGates"), readinessGates)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("restartPolicy"), restartPolicy)
        setfield!(o, Symbol("restartPolicy"), restartPolicy)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("runtimeClassName"), runtimeClassName)
        setfield!(o, Symbol("runtimeClassName"), runtimeClassName)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("schedulerName"), schedulerName)
        setfield!(o, Symbol("schedulerName"), schedulerName)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("securityContext"), securityContext)
        setfield!(o, Symbol("securityContext"), securityContext)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("serviceAccount"), serviceAccount)
        setfield!(o, Symbol("serviceAccount"), serviceAccount)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("serviceAccountName"), serviceAccountName)
        setfield!(o, Symbol("serviceAccountName"), serviceAccountName)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("shareProcessNamespace"), shareProcessNamespace)
        setfield!(o, Symbol("shareProcessNamespace"), shareProcessNamespace)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("subdomain"), subdomain)
        setfield!(o, Symbol("subdomain"), subdomain)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("terminationGracePeriodSeconds"), terminationGracePeriodSeconds)
        setfield!(o, Symbol("terminationGracePeriodSeconds"), terminationGracePeriodSeconds)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("tolerations"), tolerations)
        setfield!(o, Symbol("tolerations"), tolerations)
        validate_property(IoK8sApiCoreV1PodSpec, Symbol("volumes"), volumes)
        setfield!(o, Symbol("volumes"), volumes)
        o
    end
end # type IoK8sApiCoreV1PodSpec

const _property_map_IoK8sApiCoreV1PodSpec = Dict{Symbol,Symbol}(Symbol("activeDeadlineSeconds")=>Symbol("activeDeadlineSeconds"), Symbol("affinity")=>Symbol("affinity"), Symbol("automountServiceAccountToken")=>Symbol("automountServiceAccountToken"), Symbol("containers")=>Symbol("containers"), Symbol("dnsConfig")=>Symbol("dnsConfig"), Symbol("dnsPolicy")=>Symbol("dnsPolicy"), Symbol("enableServiceLinks")=>Symbol("enableServiceLinks"), Symbol("hostAliases")=>Symbol("hostAliases"), Symbol("hostIPC")=>Symbol("hostIPC"), Symbol("hostNetwork")=>Symbol("hostNetwork"), Symbol("hostPID")=>Symbol("hostPID"), Symbol("hostname")=>Symbol("hostname"), Symbol("imagePullSecrets")=>Symbol("imagePullSecrets"), Symbol("initContainers")=>Symbol("initContainers"), Symbol("nodeName")=>Symbol("nodeName"), Symbol("nodeSelector")=>Symbol("nodeSelector"), Symbol("priority")=>Symbol("priority"), Symbol("priorityClassName")=>Symbol("priorityClassName"), Symbol("readinessGates")=>Symbol("readinessGates"), Symbol("restartPolicy")=>Symbol("restartPolicy"), Symbol("runtimeClassName")=>Symbol("runtimeClassName"), Symbol("schedulerName")=>Symbol("schedulerName"), Symbol("securityContext")=>Symbol("securityContext"), Symbol("serviceAccount")=>Symbol("serviceAccount"), Symbol("serviceAccountName")=>Symbol("serviceAccountName"), Symbol("shareProcessNamespace")=>Symbol("shareProcessNamespace"), Symbol("subdomain")=>Symbol("subdomain"), Symbol("terminationGracePeriodSeconds")=>Symbol("terminationGracePeriodSeconds"), Symbol("tolerations")=>Symbol("tolerations"), Symbol("volumes")=>Symbol("volumes"))
const _property_types_IoK8sApiCoreV1PodSpec = Dict{Symbol,String}(Symbol("activeDeadlineSeconds")=>"Int64", Symbol("affinity")=>"IoK8sApiCoreV1Affinity", Symbol("automountServiceAccountToken")=>"Bool", Symbol("containers")=>"Vector{IoK8sApiCoreV1Container}", Symbol("dnsConfig")=>"IoK8sApiCoreV1PodDNSConfig", Symbol("dnsPolicy")=>"String", Symbol("enableServiceLinks")=>"Bool", Symbol("hostAliases")=>"Vector{IoK8sApiCoreV1HostAlias}", Symbol("hostIPC")=>"Bool", Symbol("hostNetwork")=>"Bool", Symbol("hostPID")=>"Bool", Symbol("hostname")=>"String", Symbol("imagePullSecrets")=>"Vector{IoK8sApiCoreV1LocalObjectReference}", Symbol("initContainers")=>"Vector{IoK8sApiCoreV1Container}", Symbol("nodeName")=>"String", Symbol("nodeSelector")=>"Dict{String, String}", Symbol("priority")=>"Int32", Symbol("priorityClassName")=>"String", Symbol("readinessGates")=>"Vector{IoK8sApiCoreV1PodReadinessGate}", Symbol("restartPolicy")=>"String", Symbol("runtimeClassName")=>"String", Symbol("schedulerName")=>"String", Symbol("securityContext")=>"IoK8sApiCoreV1PodSecurityContext", Symbol("serviceAccount")=>"String", Symbol("serviceAccountName")=>"String", Symbol("shareProcessNamespace")=>"Bool", Symbol("subdomain")=>"String", Symbol("terminationGracePeriodSeconds")=>"Int64", Symbol("tolerations")=>"Vector{IoK8sApiCoreV1Toleration}", Symbol("volumes")=>"Vector{IoK8sApiCoreV1Volume}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodSpec }) = collect(keys(_property_map_IoK8sApiCoreV1PodSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1PodSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodSpec[property_name]

function check_required(o::IoK8sApiCoreV1PodSpec)
    (getproperty(o, Symbol("containers")) === nothing) && (return false)
    true
end

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


text:
using MultivariateSeries

X = @ring x1 x2
n = length(X)
r = 5

L = monoms(X,5)

Xi0 = rand(2,r)
w0  = rand(r)

sigma = series(w0, Xi0, L)

w, Xi = decompose(sigma)

sigma1 = series(w, Xi, L)

@assert norm(sigma-sigma1)<1.e-6

ext: jl


text:
using PdfFraud
using Test

include("data.jl")


text:
module ApproXD

import Base.show, Base.convert

	using StatsBase: sample
	using SparseArrays
	using LinearAlgebra

	# load files
	include("bspline.jl")
	include("approx.jl")
	include("fspacexd.jl")
	include("lininterp.jl")
	# couldn't get PyPlot to install properly on appveyor.
	# so windows users can't do the plots. there not that important anyway.
	# include("plotting.jl")

	export BSpline,
		   Lininterp,
	       show,
	       getBasis,
	       getNumCoefs,
	       getNumKnots,
	       getTensorCoef,
	       FSpaceXD,
	       getValue,
	       setValue,
	       setindex!,
	       getValue!,
	       getValue,
	       hitmiss,
	       resetCache!,
	       setGrid!

end # module


text:
using SpinMC

#set up MPI
using MPI
MPI.Initialized() || MPI.Init()
commSize = MPI.Comm_size(MPI.COMM_WORLD)
commRank = MPI.Comm_rank(MPI.COMM_WORLD)
commSize > 1 || error("This example should be run on an MPI communicator of two ranks or more.")

#build lattice
a1 = (3/2, sqrt(3)/2)
a2 = (3/2, -sqrt(3)/2)
uc = UnitCell(a1,a2)

b1 = addBasisSite!(uc, (0.0, 0.0))
b2 = addBasisSite!(uc, (1.0, 0.0))

M = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
addInteraction!(uc, b1, b2, M, (0, 0))
addInteraction!(uc, b1, b2, M, (0, -1))
addInteraction!(uc, b1, b2, M, (-1, 0))

L = (16,16)
lattice = Lattice(uc, L)

#init MC simulation
thermalizationSweeps = 50000
measurementSweeps = 50000

tmin = 5.0
tmax = 10.0
beta = (commSize == 1) ? 1.0/tmin : 1.0 / (reverse([ tmax * (tmin / tmax)^(n/(commSize-1)) for n in 0:commSize-1 ])[commRank+1])

m = MonteCarlo(lattice, beta, thermalizationSweeps, measurementSweeps)

#run simulation
run!(m, outfile="simulation.h5")


text:
"""
Define a dimension of an array in terms of a sequence of sub-vectors of specified
dimensions. When used with a `ViewArray`, it allows views into an array to be generated on
demand.
"""
struct ViewAxis
    dims::Vector{Int64}
    indices::Vector{UnitRange{Int64}}
    dirty::Base.RefValue{Bool}
end
function ViewAxis(n::Integer)
    return ViewAxis(zeros(Int64, n), Vector{UnitRange{Int64}}(undef, n), Ref(true))
end
function ViewAxis(dims::Vector)
    return ViewAxis(dims, Vector{UnitRange{Int64}}(undef, length(dims)), Ref(true))
end
Base.getindex(axis::ViewAxis, idx) = axis.dims[idx]

function Base.setindex!(axis::ViewAxis, dim, idx)
    axis.dirty[] = true
    axis.dims[idx] = dim
    return dim
end

"""
Update the internal indicies of the `ViewAxis` if the dimensions have been changed. This
function is automatically called when a `ViewArray` is created from the ViewAxis.
"""
function update_indices!(axis::ViewAxis)
    if axis.dirty[]
        idx0 = 0
        for i in eachindex(axis.dims)
            idx1 = idx0 + axis.dims[i]
            axis.indices[i] = (idx0 + 1):idx1
            idx0 = idx1
        end
    end
    return nothing
end

function Base.show(io::IO, mime::MIME"text/plain", axis::ViewAxis)
    return print(io, "$ViewAxis($(axis.dims))")
end

"""
A `ViewArray` generates views into an underlying array on demand when indexed. The
size/shape of the views is determined by the `ViewAxis` used to create the `ViewArray`.
"""
struct ViewArray{N,A}
    axes::NTuple{N,ViewAxis}
    array::A

    function ViewArray(axes::NTuple{N,ViewAxis}, array::A) where {N,A}
        if length(axes) > 2
            throw(ErrorException("Not implemented"))
        end
        for (i, viewaxis) in enumerate(axes)
            update_indices!(viewaxis)
            if viewaxis.indices[end][end] != size(array, i)
                throw(ErrorException("Size mismatch in dimension $i"))
            end
        end
        return new{N,A}(axes, array)
    end
end
ViewArray(axes::ViewAxis, array) = ViewArray((axes,), array)

Base.getindex(va::ViewArray{1}, ::Colon) = va.array
Base.getindex(va::ViewArray{2}, ::Colon) = va.array
Base.getindex(va::ViewArray{1}, i1) = view(va.array, va.axes[1].indices[i1])
function Base.getindex(va::ViewArray{2}, i1, i2)
    return view(va.array, va.axes[1].indices[i1], va.axes[2].indices[i2])
end

Base.setindex!(va::ViewArray{1}, value, i1) = va.array[va.axes[1].indices[i1]] = value
function Base.setindex!(va::ViewArray{2}, value, i1, i2)
    return va.array[va.axes[1].indices[i1], va.axes[2].indices[i2]] = value
end


text:
using Documenter, TravelingSalesmanExact

makedocs(;
    modules=[TravelingSalesmanExact],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/ericphanson/TravelingSalesmanExact.jl/blob/{commit}{path}#L{line}",
    sitename="TravelingSalesmanExact.jl",
    authors="Eric",
    assets=nothing,
)

deploydocs(;
    repo="github.com/ericphanson/TravelingSalesmanExact.jl",
)


text:
"""
    overapproximate(X::S, ::Type{S}) where {S<:LazySet}

Overapproximating a set of type `S` with type `S` is a no-op.

### Input

- `X` -- set
- `Type{S}` -- set type

### Output

The input set.
"""
function overapproximate(X::S, ::Type{S}) where {S<:LazySet}
    return X
end

"""
    overapproximate(S::LazySet{N},
                    ::Type{<:HPolygon},
                    [ε]::Real=Inf)::HPolygon where {N<:Real}

Return an approximation of a given 2D convex set.
If no error tolerance is given, or is `Inf`, the result is a box-shaped polygon.
Otherwise the result is an ε-close approximation as a polygon.

### Input

- `S`           -- convex set, assumed to be two-dimensional
- `HPolygon`    -- type for dispatch
- `ε`           -- (optional, default: `Inf`) error bound

### Output

A polygon in constraint representation.
"""
function overapproximate(S::LazySet{N},
                         ::Type{<:HPolygon},
                         ε::Real=Inf
                        )::HPolygon where {N<:Real}
    @assert dim(S) == 2
    if ε == Inf
        constraints = Vector{LinearConstraint{N}}(undef, 4)
        constraints[1] = LinearConstraint(DIR_EAST(N), ρ(DIR_EAST(N), S))
        constraints[2] = LinearConstraint(DIR_NORTH(N), ρ(DIR_NORTH(N), S))
        constraints[3] = LinearConstraint(DIR_WEST(N), ρ(DIR_WEST(N), S))
        constraints[4] = LinearConstraint(DIR_SOUTH(N), ρ(DIR_SOUTH(N), S))
        return HPolygon(constraints, sort_constraints=false)
    else
        return tohrep(approximate(S, ε))
    end
end

"""
    overapproximate(S::LazySet, ε::Real)::HPolygon

Alias for `overapproximate(S, HPolygon, ε)`.
"""
function overapproximate(S::LazySet,
                         ε::Real;
                        )::HPolygon
    return overapproximate(S, HPolygon, ε)
end

"""
    overapproximate(S::LazySet,
                    Type{<:Hyperrectangle})::Union{Hyperrectangle, EmptySet}

Return an approximation of a given set as a hyperrectangle.

### Input

- `S`              -- set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.
"""
function overapproximate(S::LazySet,
                         ::Type{<:Hyperrectangle};
                        )::Union{Hyperrectangle, EmptySet}
    return box_approximation(S)
end

"""
    overapproximate(S::CartesianProductArray{N, <:AbstractHyperrectangle{N}},
                    ::Type{<:Hyperrectangle}) where {N<:Real}

Return a tight overapproximation of the cartesian product array of a finite
number of convex sets with and hyperrectangle.

### Input

- `S`              -- cartesian product array of a finite number of convex set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This method falls back to the corresponding `convert` method. Since the sets wrapped
by the cartesian product array are hyperrectangles, it can be done efficiently
without overapproximation.
"""
function overapproximate(S::CartesianProductArray{N, <:AbstractHyperrectangle{N}},
                          ::Type{<:Hyperrectangle}) where {N<:Real}
    return convert(Hyperrectangle, S)
end

"""
    overapproximate(S::CartesianProduct{N, <:AbstractHyperrectangle{N}, <:AbstractHyperrectangle{N}},
                    ::Type{<:Hyperrectangle}) where {N<:Real}

Return a tight overapproximation of the cartesian product of two
hyperrectangles by a new hyperrectangle.

### Input

- `S`              -- cartesian product of two hyperrectangular sets
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This method falls back to the corresponding `convert` method. Since the sets wrapped
by the cartesian product are hyperrectangles, it can be done efficiently
without overapproximation.
"""
function overapproximate(S::CartesianProduct{N, <:AbstractHyperrectangle{N}, <:AbstractHyperrectangle{N}},
                          ::Type{<:Hyperrectangle}) where {N<:Real}
    return convert(Hyperrectangle, S)
end

"""
    overapproximate(lm::LinearMap{N, <:AbstractHyperrectangle{N}},
                    ::Type{Hyperrectangle}) where {N}

Return a tight overapproximation of the linear map of a hyperrectangular set
using a hyperrectangle.

### Input

- `S`              -- linear map of a hyperrectangular set
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

If `c` and `r` denote the center and vector radius of a hyperrectangle `H`,
a tight hyperrectangular overapproximation of `M * H` is obtained by transforming
`c ↦ M*c` and `r ↦ abs.(M) * c`, where `abs.(⋅)` denotes the element-wise absolute
value operator.
"""
function overapproximate(lm::LinearMap{N, <:AbstractHyperrectangle{N}},
                         ::Type{Hyperrectangle}) where {N<:Real}
    M, X = lm.M, lm.X
    center_MX = M * center(X)
    radius_MX = abs.(M) * radius_hyperrectangle(X)
    return Hyperrectangle(center_MX, radius_MX)
end

"""
    overapproximate(S::LazySet)::Union{Hyperrectangle, EmptySet}

Alias for `overapproximate(S, Hyperrectangle)`.
"""
overapproximate(S::LazySet;
               )::Union{Hyperrectangle, EmptySet} =
    overapproximate(S, Hyperrectangle)

"""
    overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
                    ::Type{<:Zonotope})::Zonotope where {N<:Real}

Overapproximate the convex hull of two zonotopes.

### Input

- `S`        -- convex hull of two zonotopes
- `Zonotope` -- type for dispatch

### Algorithm

This function implements the method proposed in [1].
The convex hull of two zonotopes ``Z₁`` and ``Z₂`` of the same order,
that we write

```math
Z_j = ⟨c^{(j)}, g^{(j)}_1, …, g^{(j)}_p⟩
```
for ``j = 1, 2``, can be overapproximated as follows:

```math
CH(Z_1, Z_2) ⊆ \\frac{1}{2}⟨c^{(1)}+c^{(2)}, g^{(1)}_1+g^{(2)}_1, …, g^{(1)}_p+g^{(2)}_p, c^{(1)}-c^{(2)}, g^{(1)}_1-g^{(2)}_1, …, g^{(1)}_p-g^{(2)}_p⟩.
```

If the zonotope order is not the same, this algorithm calls
`reduce_order` to reduce the order to the minimum of the arguments.

It should be noted that the output zonotope is not necessarily the minimal
enclosing zonotope, which is in general expensive in high dimensions. This is
further investigated in [2].

[1] Reachability of Uncertain Linear Systems Using Zonotopes, A. Girard.
    HSCC 2005.

[2] Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on
    Discrete Algorithms, pp. 803-812.
"""
function overapproximate(S::ConvexHull{N, Zonotope{N}, Zonotope{N}},
                         ::Type{<:Zonotope})::Zonotope where {N<:Real}
    Z1, Z2 = S.X, S.Y

    # reduce to the same order if possible
    if order(Z1) != order(Z2)
        min_order = min(order(Z1), order(Z2))
        Z1 = reduce_order(Z1, min_order)
        Z2 = reduce_order(Z2, min_order)
    end

    if order(Z1) >= order(Z2)
        c, G = _overapproximate_convex_hull_zonotope(Z1, Z2)
    else
        c, G = _overapproximate_convex_hull_zonotope(Z2, Z1)
    end
    return Zonotope(c, G)
end

# assumes that dim(Z1) == dim(Z2) and order(Z1) >= order(Z2)
function _overapproximate_convex_hull_zonotope(Z1::Zonotope{N}, Z2::Zonotope{N}) where {N}
    c = (Z1.center + Z2.center)/N(2)

    # the case of equal order is treated separately to avoid a slicing (this creates a copy)
    if order(Z1) == order(Z2)
        G = hcat(Z1.generators .+ Z2.generators,
                 Z1.center - Z2.center,
                 Z1.generators .- Z2.generators)/N(2)
    else
        G = hcat((Z1.generators[:, 1:ngens(Z2)] .+ Z2.generators)/N(2),
                 (Z1.center - Z2.center)/N(2),
                 (Z1.generators[:, 1:ngens(Z2)] .- Z2.generators)/N(2),
                 Z1.generators[:, ngens(Z2)+1:end])
    end
    return c, G
end

"""
    overapproximate(Z::Zonotope, ::Type{<:Hyperrectangle})::Hyperrectangle

Return a tight overapproximation of a zonotope with an axis-aligned box.

### Input

- `Z`              -- zonotope
- `Hyperrectangle` -- type for dispatch

### Output

A hyperrectangle.

### Algorithm

This function implements the method in [Section 5.1.2, 1]. A zonotope
``Z = ⟨c, G⟩`` can be overapproximated tightly by an axis-aligned box
(i.e. a `Hyperrectangle`) such that its center is ``c`` and the radius along
dimension ``i`` is the column-sum of the absolute values of the ``i``-th row
of ``G`` for ``i = 1,…, p``, where ``p`` is the number of generators of ``Z``.

[1] *Althoff, M., Stursberg, O., & Buss, M. (2010). Computing reachable sets of
hybrid systems using a combination of zonotopes and polytopes. Nonlinear analysis:
hybrid systems, 4(2), 233-249.*
"""
function overapproximate(Z::Zonotope, ::Type{<:Hyperrectangle})::Hyperrectangle
    r = sum(abs.(Z.generators), dims=2)[:]
    return Hyperrectangle(Z.center, r)
end

"""
    overapproximate(X::LazySet{N}, dir::AbstractDirections{N})::HPolytope{N}
        where {N}

Overapproximating a set with template directions.

### Input

- `X`   -- set
- `dir` -- (concrete) direction representation

### Output

An `HPolytope` overapproximating the set `X` with the directions from `dir`.
"""
function overapproximate(X::LazySet{N},
                         dir::AbstractDirections{N}
                        )::HPolytope{N} where {N}
    halfspaces = Vector{LinearConstraint{N}}()
    sizehint!(halfspaces, length(dir))
    H = HPolytope(halfspaces)
    for d in dir
        addconstraint!(H, LinearConstraint(d, ρ(d, X)))
    end
    return H
end

"""
    overapproximate(X::LazySet{N},
                    dir::Type{<:AbstractDirections})::HPolytope{N} where {N}

Overapproximating a set with template directions.

### Input

- `X`   -- set
- `dir` -- type of direction representation

### Output

A `HPolytope` overapproximating the set `X` with the directions from `dir`.
"""
function overapproximate(X::LazySet{N},
                         dir::Type{<:AbstractDirections}
                        )::HPolytope{N} where {N}
    return overapproximate(X, dir{N}(dim(X)))
end

"""
    overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real}

Return the overapproximation of a real unidimensional set with an interval.

### Input

- `S`        -- one-dimensional set
- `Interval` -- type for dispatch

### Output

An interval.

### Algorithm

The method relies on the exact conversion to `Interval`. Two support
function evaluations are needed in general.
"""
function overapproximate(S::LazySet{N}, ::Type{Interval}) where {N<:Real}
    @assert dim(S) == 1 "cannot overapproximate a $(dim(S))-dimensional set with an `Interval`"
    return convert(Interval, S)
end

function overapproximate_cap_helper(X::LazySet{N},             # compact set
                                    P::AbstractPolyhedron{N},  # polyhedron
                                    dir::AbstractDirections{N};
                                    kwargs...
                                   ) where {N<:Real}
    Hi = constraints_list(P)
    m = length(Hi)
    Q = HPolytope{N}()

    for di in dir
        ρ_X_Hi_min = ρ(di, X ∩ Hi[1], kwargs...)
        for i in 2:m
            ρ_X_Hi = ρ(di, X ∩ Hi[i], kwargs...)
            if ρ_X_Hi < ρ_X_Hi_min
                ρ_X_Hi_min = ρ_X_Hi
            end
        end
        addconstraint!(Q, HalfSpace(di, ρ_X_Hi_min))
    end
    return Q
end

"""
    overapproximate(cap::Intersection{N, <:LazySet, <:AbstractPolyhedron{N}},
                    dir::AbstractDirections{N};
                    kwargs...
                   ) where {N<:Real}

Return the overapproximation of the intersection between a compact set and a
polytope given a set of template directions.

### Input

- `cap`         -- intersection of a compact set and a polytope
- `dir`         -- template directions
- `kwargs`      -- additional arguments that are passed to the support function
                   algorithm

### Output

A polytope in H-representation such that the normal direction of each half-space
is given by an element of `dir`.

### Algorithm

Let `di` be a direction drawn from the set of template directions `dir`.
Let `X` be the compact set and let `P` be the polytope. We overapproximate the
set `X ∩ H` with a polytope in constraint representation using a given set of
template directions `dir`.

The idea is to solve the univariate optimization problem `ρ(di, X ∩ Hi)` for
each half-space in the set `P` and then take the minimum.
This gives an overapproximation of the exact support function.

This algorithm is inspired from [G. Frehse, R. Ray. Flowpipe-Guard Intersection
for Reachability Computations with Support
Functions](https://www.sciencedirect.com/science/article/pii/S1474667015371809).

### Notes

This method relies on having available the `constraints_list` of the polytope
`P`.

This method of overapproximations can return a non-empty set even if the original
intersection is empty.
"""
function overapproximate(cap::Intersection{N,
                                           <:LazySet,
                                           <:AbstractPolyhedron{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.X, cap.Y, dir; kwargs...)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolyhedron{N},
                                           <:LazySet},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.Y, cap.X, dir; kwargs...)
end

# disambiguation
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolytope{N},
                                           <:AbstractPolyhedron{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.X, cap.Y, dir; kwargs...)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolyhedron{N},
                                           <:AbstractPolytope{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.Y, cap.X, dir; kwargs...)
end

# disambiguation
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolytope{N},
                                           <:AbstractPolytope{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate_cap_helper(cap.X, cap.Y, dir; kwargs...)
end

"""
    overapproximate(cap::Intersection{N, <:HalfSpace{N}, <:AbstractPolytope{N}},
                    dir::AbstractDirections{N};
                    [kwargs]...
                   ) where {N<:Real}

Return the overapproximation of the intersection between a half-space and a
polytope given a set of template directions.

### Input

- `cap`         -- intersection of a half-space and a polytope
- `dir`         -- template directions
- `kwargs`      -- additional arguments that are passed to the support function
                   algorithm

### Output

A polytope in H-representation such that the normal direction of each half-space
is given by an element of `dir`.
"""
function overapproximate(cap::Intersection{N,
                                           <:HalfSpace{N},
                                           <:AbstractPolytope{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    H = HPolytope{N}()
    c = constraints_list(H)
    append!(c, constraints_list(cap.X))
    append!(c, constraints_list(cap.Y))
    return overapproximate(H, dir; kwargs...)
end

# symmetric method
function overapproximate(cap::Intersection{N,
                                           <:AbstractPolytope{N},
                                           <:HalfSpace{N}},
                         dir::AbstractDirections{N};
                         kwargs...
                        ) where {N<:Real}
    return overapproximate(swap(cap), dir; kwargs...)
end


text:
function ldlt!(a::Matrix{Float64}, d::Vector{Float64}; small=1.0e-10)
  n = size(a,1)
  for k in 1:n-1
    d[1] = a[1,1]
    if abs(a[k,k]) > small
      for i in k+1:n
        x=a[i,k]/a[k,k]
        for j=k+1:n
          a[i,j] -= a[k,j]*x
        end
        d[i]=a[i,i]
      end
    else
      println("Zero pivot found in row $k")
    end
  end
end


text:
"""
	grdgradient(cmd0::String="", arg1=[], kwargs...)

Compute the directional derivative in a given direction, or to find the direction [and the magnitude]
of the vector gradient of the data.

Full option list at [`grdgradient`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html)

Parameters
----------

- **A** : **azim** : -- Str or Number --    Flags = azim[/azim2]

    Azimuthal direction for a directional derivative.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html#a)
- **D** : **find_dir** : -- Str --      Flags = [a][c][o][n]

    Find the direction of the positive (up-slope) gradient of the data.
	[`-D`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html#d)
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdgradient(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html#g)
- **E** : **lambertian** : -- Str --    Flags = [m|s|p]azim/elev[+aambient][+ddiffuse][+pspecular][+sshine]

    Compute Lambertian radiance appropriate to use with grdimage and grdview.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/surface.html#e)
- **L** : **bc** : **boundary** : -- Str --

    Boundary condition flag.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html#l)
- **N** : **normalize** : -- Str --     Flags = [e|t][amp][+ssigma][+ooffset]

    Normalization. [Default is no normalization.] The actual gradients g are offset and scaled
    to produce normalized gradients.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html#n)
- $(GMT.opt_R)
- **S** : **slopegrid** : -- Str --

    Name of output grid file with scalar magnitudes of gradient vectors. Requires D but makes G optional.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grdgradient.html#s)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_f)
"""
function grdgradient(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("grdgradient", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_n(cmd, d)

    cmd = add_opt(cmd, 'A', d, [:A :azim])
    cmd = add_opt(cmd, 'D', d, [:D :find_dir])
    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
    cmd = add_opt(cmd, 'E', d, [:E :lambertian])
	cmd = add_opt(cmd, 'L', d, [:L :bc :boundary])
	cmd = add_opt(cmd, 'N', d, [:N :normalize])
    cmd = add_opt(cmd, 'S', d, [:S :slopegrid])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end
	return common_grd(d, cmd, got_fname, 1, "grdgradient", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdgradient(arg1=[], cmd0::String=""; kw...) = grdgradient(cmd0, arg1; kw...)


text:
function NodesSet()
    return NodesSet(())
end

function add(obj::NodesSet, arg0::BSPTree)
    return jcall(obj, "add", void, (BSPTree,), arg0)
end

function add_all(obj::NodesSet, arg0::Iterable)
    return jcall(obj, "addAll", void, (Iterable,), arg0)
end

function for_each(obj::Iterable, arg0::Consumer)
    return jcall(obj, "forEach", void, (Consumer,), arg0)
end

function iterator(obj::NodesSet)
    return jcall(obj, "iterator", JIterator, ())
end

function spliterator(obj::Iterable)
    return jcall(obj, "spliterator", Spliterator, ())
end


text:
using BenchmarkTools
using JuLIP
using JuLIP: AbstractAtoms, AbstractCalculator

function perfbm(id::AbstractString, at::AbstractAtoms, calc::AbstractCalculator;
                e = true, elist = true, f = true, flist = true)
   println("--------------------------------------------------------------------------")
   println(id)
   if e
      print("Energy Assembly (without nlist): ")
      @btime energy($calc, $at)
   end
   if elist
      print("Energy Assembly (with nlist):    ")
      @btime energy($calc,  rattle!($at, 0.001))
   end
   if f
      print("Force Assembly (without nlist):  ")
      @btime forces($calc, $at)
   end
   if flist
      print("Force Assembly (with nlist):     ")
      @btime forces($calc,  rattle!($at, 0.001))
   end
end

println()
println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
println("   JuLIP Performance Regression Tests")
println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
println()

perfbm("LENNARD-JONES",
         bulk(:Al, cubic=true) * (10,10,8),
         lennardjones(r0=rnn(:Al)) )

data = joinpath(dirname(@__FILE__), "..", "data") * "/"
perfbm("EAM (Splines)",
         bulk(:Fe, cubic=true) * (12,12,8),
         EAM(data * "pfe.plt", data * "ffe.plt", data * "F_fe.plt") )

perfbm("STILLINGER-WEBER",
         bulk(:Si, cubic=true) * (12,15,12),
         StillingerWeber())

##
using BenchmarkTools
using JuLIP
using JuLIP: AbstractAtoms, AbstractCalculator

at = bulk(:Al, cubic=true) * (10,10,8)
calc = lennardjones(r0=rnn(:Al))
tmp = JuLIP.alloc_temp(calc, at)
JuLIP.energy!(tmp, calc, at)

# Performance prior to restructuring
# TODO: fix performance of PairPotentials
# --------------------------------------------------------------------------
# LENNARD-JONES
# Energy Assembly (without nlist):   66.695 ms (41653 allocations: 47.65 MiB)
# Energy Assembly (with nlist):      68.778 ms (41675 allocations: 41.52 MiB)
# Force Assembly (without nlist):    73.158 ms (41657 allocations: 43.40 MiB)
# Force Assembly (with nlist):       72.561 ms (41677 allocations: 41.59 MiB)
# --------------------------------------------------------------------------
# EAM (Splines)
# Energy Assembly (without nlist):   69.074 ms (71535 allocations: 43.47 MiB)
# Energy Assembly (with nlist):      74.438 ms (71558 allocations: 44.49 MiB)
# Force Assembly (without nlist):    83.091 ms (69237 allocations: 44.23 MiB)
# Force Assembly (with nlist):       87.327 ms (69257 allocations: 44.51 MiB)
# --------------------------------------------------------------------------
# STILLINGER-WEBER
# Energy Assembly (without nlist):   78.205 ms (402680 allocations: 65.56 MiB)
# Energy Assembly (with nlist):      86.708 ms (410401 allocations: 69.01 MiB)
# Force Assembly (without nlist):    87.153 ms (411005 allocations: 67.34 MiB)
# Force Assembly (with nlist):       91.166 ms (410672 allocations: 69.42 MiB)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# priority list, for newly freed ambulances only
# (a priority list which uses all free ambulances can be described as a nested compliance table)
# From thesis: "Simulation optimisation and Markov models for dynamic ambulance redeployment"

# initialise data relevant to move up
function initPriorityList!(sim::Simulation, priorityListFilename::String)
	# read priority list
	priorityList = readPriorityListFile(priorityListFilename)
	initPriorityList!(sim, priorityList)
end

# initialise data relevant to move up
function initPriorityList!(sim::Simulation, priorityList::PriorityList)
	# shorthand names:
	@unpack ambulances, numStations = sim
	pld = sim.moveUpData.priorityListData

	checkPriorityList(priorityList, sim)
	pld.priorityList = priorityList

	pld.stationNumFreeAmbs = Vector{Int}(undef, numStations)

	# check that ambulance to station assignments follow priority list
	pld.stationNumFreeAmbs[:] .= 0
	for i = 1:sim.numAmbs
		pld.stationNumFreeAmbs[ambulances[i].stationIndex] += 1
	end
	if !all(pld.stationNumFreeAmbs[i] == sum(pld.priorityList .== i) for i = 1:numStations)
		@warn("Number of ambulances at each station does not match priority list")
	end
end

function priorityListMoveUp(sim::Simulation, newlyFreedAmb::Ambulance)
	@assert(sim.moveUpData.useMoveUp)

	# shorthand:
	pld = sim.moveUpData.priorityListData
	@unpack priorityList, stationNumFreeAmbs = pld
	@unpack ambulances, stations, numAmbs = sim

	# calculate the number of free ambulances at (or travelling to) each station
	stationNumFreeAmbs[:] .= 0
	for i = 1:numAmbs
		# do not count newly freed ambulance, it has not been assigned a station
		if isAmbMovable(ambulances[i]) && i != newlyFreedAmb.index
			stationNumFreeAmbs[ambulances[i].stationIndex] += 1
		end
	end
	numFreeAmbs = sum(stationNumFreeAmbs) + 1

	# return first item in priority list that is not already filled
	stationIndex = nullIndex
	for i = 1:numFreeAmbs
		j = priorityList[i] # station index
		stationNumFreeAmbs[j] -= 1
		if stationNumFreeAmbs[j] < 0
			stationIndex = j
			break
		end
	end
	@assert(stationIndex != nullIndex)

	return [newlyFreedAmb], [stations[stationIndex]]
end

# check that priority list is valid
function checkPriorityList(priorityList::PriorityList;
	numAmbs::Int = nullIndex, numStations::Int = nullIndex,
	stationCapacities::Union{Vector{Int},Nothing} = nothing)

	m = length(priorityList) # number of ambulances
	@assert(numAmbs == nullIndex || numAmbs == m)
	@assert(all(v -> v >= 1, priorityList)) # need positive integers (station indices)
	@assert(numStations == nullIndex || maximum(priorityList) <= numStations)
	if stationCapacities != nothing
		n = length(stationCapacities) # number of stations
		@assert(numStations == nullIndex || numStations == n)
		# check that priority list does not violate station capacity constraints
		for j = 1:n
			@assert(count(isequal(j), priorityList) <= stationCapacities[j])
		end
	end
	return true
end
function checkPriorityList(priorityList::PriorityList, sim::Simulation)
	checkPriorityList(priorityList; numAmbs = sim.numAmbs, numStations = sim.numStations, stationCapacities = [s.capacity for s in sim.stations])
end

# returns a randomly generated priority list
function makeRandPriorityList(numAmbs::Int, numStations::Int;
	stationCapacities::Union{Vector{Int},Nothing} = nothing,
	rng::AbstractRNG = GLOBAL_RNG)::PriorityList

	@assert(numStations > 0)
	if stationCapacities == nothing
		return rand(rng, 1:numStations, numAmbs)
	else
		@assert(numStations == length(stationCapacities))
		@assert(numAmbs <= sum(stationCapacities))
		remainingCapacity = copy(stationCapacities)
		unfilledStations = Set(findall(!iszero, remainingCapacity))
		priorityList = PriorityList(undef,numAmbs)
		priorityList[:] .= 0
		for i = 1:numAmbs
			j = rand(rng, unfilledStations) # station index
			priorityList[i] = j
			remainingCapacity[j] -= 1
			if remainingCapacity[j] == 0
				delete!(unfilledStations, j)
			end
		end
		return priorityList
	end
end


text:
using DataStructures

input = joinpath(@__DIR__, "input")
lines = readlines(input)

M = hcat(map(line -> collect(line), lines)...)

# p2 adjustment
newlines = []
append!(newlines, lines[1:3])
append!(newlines, ["  #D#C#B#A#  "])
append!(newlines, ["  #D#B#A#C#  "])
append!(newlines, lines[4:end])
N = hcat(map(line -> collect(line), newlines)...)

function p(m)
    y, x = size(m)
    for i in 1:x
        for j in 1:y
            print(m[CartesianIndex(j, i)])
        end
        println()
    end
end

function path(amphipod, space)
    x, y = Tuple(amphipod)
    dx, dy = Tuple(space)
    # x translation is the same for both cases
    path = [CartesianIndex(i, 2) for i in min(x, dx):max(x, dx)]
    # amphipod moving to hallway
    if dy != 2
        append!(path, [CartesianIndex(dx, i) for i in min(y, dy):max(y, dy)])
    # amphipod moving home
    else
        append!(path, [CartesianIndex(x, i) for i in min(y, dy):max(y, dy)])
    end
    unique!(path)
    filter!(idx -> idx != amphipod, path)
    return path
end

function cost(amphipod, p, m)
    a = m[amphipod]
    multiplier = 1
    if a == 'A'
        multiplier = 1
    elseif a == 'B'
        multiplier = 10
    elseif a == 'C'
        multiplier = 100
    elseif a == 'D'
        multiplier = 1000
    end
    return length(p) * multiplier
end

function desirable(amphipod, destination, m)
    a = m[amphipod]
    home = HOMES[a]
    y, x = Tuple(amphipod)
    dy, dx = Tuple(destination)

    # already "home". won't move
    amphipod in home && all(idx -> m[idx] == a, filter(idx -> idx > amphipod, home)) && return false

    # in hallway, will only move if destination is "home" and is empty and either the last space in room
    # or next to a similar amphipod
    if x == 2
        # destination not in hallway
        dx == 2 && return false
        # destination not in home row
        !(destination in home) && return false
        # destination is occupied
        m[destination] != '.' && return false
        # destination is deepest unoccupied spot
        destination != maximum(filter(idx -> m[idx] == '.', home)) && return false
        # room has non similar amphipods
        any(idx -> m[idx] != a, filter(idx -> idx > destination, home)) && return false
        return true
    # in side room, can always try to move to hallway
    elseif dx == 2
        return true
    # trying to move from hallway to hallway or room to room
    else
        return false
    end
end

function moves(m)
    # all legal spaces
    legal = LEGAL
    # all spaces currently occupied
    amphipods = filter(idx -> m[idx] in ['A', 'B', 'C', 'D'], CartesianIndices(m))
    # unoccupied spaces
    unoccupied = setdiff(legal, amphipods)

    acc = []
    for amphipod in amphipods
        # remove any side rooms that the amphipod doesn't like
        available = filter(idx -> desirable(amphipod, idx, m), unoccupied)
        for space in available
            pa = path(amphipod, space)
            if all(idx -> m[idx] == '.', pa)
                c = cost(amphipod, pa, m)
                n = copy(m)
                n[amphipod] = '.'
                n[space] = m[amphipod]
                push!(acc, (n, c))
            end
        end
    end
    return acc
end

# terrible heuristic. couldn't find a good heuristic for p2
function estimate(m, goal)
    return 0
end

function astar(start, goal, heuristic)
    open = PriorityQueue()
    open[start] = heuristic(start, goal)

    g = DefaultDict(typemax(Int))
    g[start] = 0

    # solution path, because i have no idea what i'm doing
    cf = Dict()

    rounds = 0

    while !isempty(open)
        current = dequeue!(open)
        # fin
        current == goal && return (g[current], cf)

        for move in moves(current)
            next, cost = move
            tentative = g[current] + cost
            if tentative < g[next]
                cf[next] = current
                g[next] = tentative
                open[next] = tentative + heuristic(next, goal)
            end
        end
    end

    return (nothing, nothing)
end

function solution(cf, goal)
    current = goal
    while current in keys(cf)
        current = cf[current]
        p(current)
        println()
    end
end

LEGAL = [
    CartesianIndex(2,2),
    CartesianIndex(3,2),
    CartesianIndex(5,2),
    CartesianIndex(7,2),
    CartesianIndex(9,2),
    CartesianIndex(11,2),
    CartesianIndex(12,2),
    CartesianIndex(4,3),
    CartesianIndex(4,4),
    CartesianIndex(6,3),
    CartesianIndex(6,4),
    CartesianIndex(8,3),
    CartesianIndex(8,4),
    CartesianIndex(10,3),
    CartesianIndex(10,4)
]

A = [CartesianIndex(4,3), CartesianIndex(4,4)]
B = [CartesianIndex(6,3), CartesianIndex(6,4)]
C = [CartesianIndex(8,3), CartesianIndex(8,4)]
D = [CartesianIndex(10,3), CartesianIndex(10,4)]
HOMES = Dict([('A', A), ('B', B), ('C', C), ('D', D)])

# construct goal
goal = copy(M)
goal[CartesianIndex(4,3)] = 'A'
goal[CartesianIndex(4,4)] = 'A'
goal[CartesianIndex(6,3)] = 'B'
goal[CartesianIndex(6,4)] = 'B'
goal[CartesianIndex(8,3)] = 'C'
goal[CartesianIndex(8,4)] = 'C'
goal[CartesianIndex(10,3)] = 'D'
goal[CartesianIndex(10,4)] = 'D'

p1, cf = astar(copy(M), goal, estimate)
# solution(cf, goal)

append!(LEGAL, [
    CartesianIndex(4,5),
    CartesianIndex(4,6),
    CartesianIndex(6,5),
    CartesianIndex(6,6),
    CartesianIndex(8,5),
    CartesianIndex(8,6),
    CartesianIndex(10,5),
    CartesianIndex(10,6)
])

append!(A, [CartesianIndex(4,5), CartesianIndex(4,6)])
append!(B, [CartesianIndex(6,5), CartesianIndex(6,6)])
append!(C, [CartesianIndex(8,5), CartesianIndex(8,6)])
append!(D, [CartesianIndex(10,5), CartesianIndex(10,6)])
HOMES = Dict([('A', A), ('B', B), ('C', C), ('D', D)])

# construct goal
goal = copy(N)
goal[CartesianIndex(4,3)] = 'A'
goal[CartesianIndex(4,4)] = 'A'
goal[CartesianIndex(6,3)] = 'B'
goal[CartesianIndex(6,4)] = 'B'
goal[CartesianIndex(8,3)] = 'C'
goal[CartesianIndex(8,4)] = 'C'
goal[CartesianIndex(10,3)] = 'D'
goal[CartesianIndex(10,4)] = 'D'
goal[CartesianIndex(4,5)] = 'A'
goal[CartesianIndex(4,6)] = 'A'
goal[CartesianIndex(6,5)] = 'B'
goal[CartesianIndex(6,6)] = 'B'
goal[CartesianIndex(8,5)] = 'C'
goal[CartesianIndex(8,6)] = 'C'
goal[CartesianIndex(10,5)] = 'D'
goal[CartesianIndex(10,6)] = 'D'

p2, cf = astar(copy(N), goal, estimate)
# solution(cf, goal)

println("-----------------------------------------------------------------------")
println("amphipod -- part one :: $p1")
println("amphipod -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 10607)
@assert(p2 == 59071)


text:
function get_dim(X1::Array{T,3}, X2::Matrix{T}, D::Matrix{T}) where T
    N, J, K1 = size(X1)
    J2, K2 = size(X2)
    N2, L = size(D)
    @assert N==N2 && J==J2
    return (N=N, J=J, K1=K1, K2=K2, L=L)
end

function get_kdim(X1::Array{T,3}, X2::Matrix{T}, D::Matrix{T}) where T
    N, J, K1, K2, L = get_dim(X1, X2, D)
    kdim = K1+K2*L+J-1
    return kdim
end

"""
Get random choice set of size NC from total J products

For each individividual, the actual choice y_i must be included in the dataset.
Other choices in the subset are randomly picked from fullset with equal probability
"""
function get_random_set(J::Integer, NC::Integer, Y::Vector{T}) where T<:Integer
    N = length(Y)
    C = Vector{Vector{T}}(undef, N)
    for (i, y) in enumerate(Y)
        subset = sample(1:J, NC; replace=false)
        if !(y in subset)
            subset[1] = y
        end
        C[i] = subset
    end
    return C
end

function initial_guess(m::DCModel)
    if !m.hetero_preference
        init_guess = randn(get_kdim(m.X1, m.X2, m.D))
    end
    return init_guess
end


text:
abstract type Perturbation end

mutable struct ThreePhaseFault <: Perturbation
    time::Float64
    Ybus::SparseMatrixCSC{Complex{Float64}, Int64}
end

function get_affect(::PSY.System, pert::ThreePhaseFault)
    return (integrator) -> begin
        integrator.p.Ybus = pert.Ybus
    end
end

mutable struct ControlReferenceChange <: Perturbation
    time::Float64
    device::PSY.DynamicInjection
    signal_index::Int64
    ref_value::Float64
end

function get_affect(system::PSY.System, pert::ControlReferenceChange)
    device = PSY.get_component(typeof(pert.device), system, PSY.get_name(pert.device))
    pert.device = device
    return (integrator) -> begin
        control_ref = PSY.get_ext(device)[CONTROL_REFS]
        return control_ref[pert.signal_index] = pert.ref_value
    end
end


text:
module TestLearningCompositesCore

# using Revise
using Test
using MLJBase
using ..Models
using ..TestUtilities
using CategoricalArrays
import Random.seed!
seed!(1234)

N =100
X = (x1=rand(N), x2=rand(N), x3=rand(N));
y = 2X.x1  - X.x2 + 0.05*rand(N);

@testset "error messages for invalid learning networks" begin

    Xs = source(X)
    ys = source(y)
    mach1 = machine(Standardizer(), Xs)
    W = transform(mach1, Xs)
    fit!(W, verbosity=0)
    @test_logs((:error, r"Failed"), @test_throws Exception W(34))

    mach2 = machine(DecisionTreeClassifier(), W, ys)
    yhat = predict(mach2, W)
    @test_logs((:error, r"Problem fitting"),
               (:info, r"Running type checks"),
               (:warn, MLJBase.alert_generic_scitype_mismatch(
                   scitype((X, y)),
                   MLJBase.fit_data_scitype(mach2.model),
                   typeof(mach2.model)
               )),
               (:info, r"^It seems an upstream"),
               (:error, r"^Problem fitting"),
               @test_throws Exception fit!(yhat, verbosity=-1))
    Xs = source()
    mach1 = machine(Standardizer(), Xs)
    W = transform(mach1, Xs)
    @test_logs((:error, r"Problem fitting"),
               (:warn, r"^Some learning network source"),
               (:info, r"Running type checks"),
               (:warn, MLJBase.alert_generic_scitype_mismatch(
                   Tuple{Nothing},
                   MLJBase.fit_data_scitype(mach1.model),
                   typeof(mach1.model)
               )),
               (:info, r"^It seems an upstream"),
               (:error, r"^Problem fitting"),
               @test_throws Exception fit!(W, verbosity=-1))
end

@testset "network #1" begin

    knn_ = KNNRegressor(K=7)

    # split the rows:
    allrows = eachindex(y);
    train, valid, test = partition(allrows, 0.7, 0.15);

    Xtrain = selectrows(X, train);
    ytrain = y[train];

    Xs = source(Xtrain);
    ys = source(ytrain);

    knn1 = machine(knn_, Xs, ys)
    @test_mach_sequence fit_only!(knn1) [(:train, knn1),]
    knn_.K = 5
    @test_logs((:info, r"Training"),
               fit_only!(knn1, rows=train[1:end-10], verbosity=2))
    @test_logs (:info, r"Training") fit_only!(knn1, verbosity=2)
    yhat = predict(knn1, Xs);
    yhat(selectrows(X, test))
    @test rms(yhat(selectrows(X, test)), y[test]) < 0.3

    @test_mach_sequence fit!(yhat) [(:skip, knn1),]
    pred = yhat();
end

@testset "network #2" begin

    N =100
    X = (x1=rand(N),
         x2=rand(N),
         x3=categorical(rand("yn",N)),
         x4=categorical(rand("yn",N)));

    y = 2X.x1  - X.x2 + 0.05*rand(N);
    X = source(X)
    y = source(y)

    hot = OneHotEncoder()
    hotM = machine(hot, X)
    W = transform(hotM, X)
    @test scitype(W) == CallableReturning{Table}

    knn = KNNRegressor()
    knnM = machine(knn, W, y)
    yhat = predict(knnM, W)
    @test scitype(yhat) ==
        MLJBase.CallableReturning{AbstractVector{Continuous}}

    @test_mach_sequence fit!(yhat) [(:train, hotM), (:train, knnM)]
    @test_mach_sequence fit!(yhat) [(:skip, hotM), (:skip, knnM)]

    hot.drop_last = true
    @test_mach_sequence fit!(yhat) [(:update, hotM), (:train, knnM)]
    @test_mach_sequence fit!(yhat) [(:skip, hotM), (:skip, knnM)]

    knn.K = 17
    @test_mach_sequence fit!(yhat) [(:skip, hotM), (:update, knnM)]

    # change rows:
    @test_mach_sequence fit!(yhat, rows=1:100) [(:train, hotM), (:train, knnM)]

    # change rows again:
    @test_mach_sequence fit!(yhat) [(:train, hotM), (:train, knnM)]

    # force:
    @test_mach_sequence fit!(yhat, force=true) [(:train, hotM), (:train, knnM)]

end

@testset "with parallel regressors and pre-processors" begin
    N =100
    X = (x1=rand(N),
         x2=rand(N),
         x3=categorical(rand("yn",N)),
         x4=categorical(rand("yn",N)));

    y = abs.(2X.x1  - X.x2 + 0.05*rand(N));
    X = source(X)
    y = source(y)

    cox = UnivariateBoxCoxTransformer()
    coxM = machine(cox, y, cache=true)
    z = transform(coxM, y)

    hot = OneHotEncoder()
    hotM = machine(hot, X, cache=true)
    W = transform(hotM, X)

    knn = KNNRegressor()
    knnM = machine(knn, W, z, cache=true)
    zhat1 = predict(knnM, W)

    @test_mach_sequence fit!(W) [(:train, hotM), ]
    @test_mach_sequence fit!(z) [(:train, coxM), ]
    evaluate!(knnM, measure=l2, verbosity=0);
    fit!(knnM, verbosity=0)

    cox.shift=true
    @test_mach_sequence(fit!(zhat1),
                   [(:update, coxM), (:skip, hotM), (:train, knnM)],
                   [(:skip, hotM), (:update, coxM), (:train, knnM)])

    # no training:
    @test_mach_sequence(fit!(zhat1),
                   [(:skip, coxM), (:skip, hotM), (:skip, knnM)],
                   [(:skip, hotM), (:skip, coxM), (:skip, knnM)])

    tree = DecisionTreeRegressor()
    treeM = machine(tree, W, z, cache=true)
    zhat2 = predict(treeM, W)

    zhat = 0.5*zhat1 + 0.5*zhat2
    @test elscitype(zhat) == Unknown

    @test_mach_sequence(fit!(zhat),
                   [(:skip, coxM), (:skip, hotM),
                     (:skip, knnM), (:train, treeM)],
                   [(:skip, hotM), (:skip, coxM),
                    (:skip, knnM), (:train, treeM)],
                   [(:skip, coxM), (:skip, hotM),
                    (:train, treeM), (:skip, knnM)],
                   [(:skip, hotM), (:skip, coxM),
                    (:train, treeM), (:skip, knnM)])

    yhat = inverse_transform(coxM, zhat)

    @test_mach_sequence(fit!(yhat),
                   [(:skip, coxM), (:skip, hotM),
                     (:skip, knnM), (:skip, treeM)],
                   [(:skip, hotM), (:skip, coxM),
                    (:skip, knnM), (:skip, treeM)],
                   [(:skip, coxM), (:skip, hotM),
                    (:skip, treeM), (:skip, knnM)],
                   [(:skip, hotM), (:skip, coxM),
                    (:skip, treeM), (:skip, knnM)])

    # error handling:
    MLJBase.rebind!(X, "junk5")
    @test_logs((:error, r""),
               (:error, r""),
               (:error, r""),
               (:error, r""),
               @test_throws Exception fit!(yhat, verbosity=-1, force=true))

end

@testset "network #3" begin

    N =100
    X = (x1=rand(N), x2=rand(N), x3=rand(N))
    y = 2X.x1  - X.x2 + 0.05*rand(N)

    XX = source(X)
    yy = source(y)

    # construct a transformer to standardize the target:
    uscale_ = UnivariateStandardizer()
    uscale = machine(uscale_, yy)

    # get the transformed inputs, as if `uscale` were already fit:
    z = transform(uscale, yy)

    # construct a transformer to standardize the inputs:
    scale_ = Standardizer()
    scale = machine(scale_, XX) # no need to fit

    # get the transformed inputs, as if `scale` were already fit:
    Xt = transform(scale, XX)

    # do nothing to the DataFrame
    Xa = node(identity, Xt)

    # choose a learner and make it machine:
    knn_ = KNNRegressor(K=7) # just a container for hyperparameters
    knn = machine(knn_, Xa, z) # no need to fit

    # get the predictions, as if `knn` already fit:
    zhat = predict(knn, Xa)

    # inverse transform the target:
    yhat = inverse_transform(uscale, zhat)

    # fit-through training:
    @test_mach_sequence(fit!(yhat, rows=1:50),
                   [(:train, uscale), (:train, scale), (:train, knn)],
                   [(:train, scale), (:train, uscale), (:train, knn)])

    @test_mach_sequence(fit!(yhat, rows=1:50),
                   [(:skip, uscale), (:skip, scale), (:skip, knn)],
                   [(:skip, scale), (:skip, uscale), (:skip, knn)])

    # change rows:
    @test_mach_sequence(fit!(yhat),
                   [(:train, uscale), (:train, scale), (:train, knn)],
                   [(:train, scale), (:train, uscale), (:train, knn)])

    knn_.K =67
    @test_mach_sequence(fit!(yhat),
                   [(:skip, uscale), (:skip, scale), (:update, knn)],
                   [(:skip, scale), (:skip, uscale), (:update, knn)])
end

@testset "network with machines sharing one model" begin

    N =100
    X = (x1=rand(N), x2=rand(N), x3=rand(N))
    y = 2X.x1  - X.x2 + 0.05*rand(N)

    XX = source(X)
    yy = source(y)

    # construct a transformer to standardize the target:
    uscale_ = UnivariateStandardizer()
    uscale = machine(uscale_, yy)

    # get the transformed inputs, as if `uscale` were already fit:
    z = transform(uscale, yy)

    # construct a transformer to standardize the inputs:
    xscale_ = Standardizer()
    xscale = machine(xscale_, XX) # no need to fit

    # get the transformed inputs, as if `scale` were already fit:
    Xt = transform(xscale, XX)

    # choose a learner and make two machines from it:
    knn_ = KNNRegressor(K=7) # just a container for hyperparameters
    knn1 = machine(knn_, Xt, z) # no need to fit
    knn2 = machine(knn_, Xt, z) # no need to fit

    # get the predictions, as if `knn` already fit:
    zhat1 = predict(knn1, Xt)
    zhat2 = predict(knn2, Xt)
    zhat = zhat1 + zhat2

    # inverse transform the target:
    yhat = inverse_transform(uscale, zhat)

    fit!(yhat, verbosity=0)

    θ = fitted_params(yhat)

end

mutable struct PermuteArgs <: MLJBase.Static
    permutation::NTuple{N,Int} where N
end

MLJBase.transform(p::PermuteArgs, _, args...) =
    Tuple([args[i] for i in p.permutation])
MLJBase.inverse_transform(p::PermuteArgs, _, args...) =
    Tuple([args[i] for i in sortperm(p.permutation |> collect)])

@testset "nodes wrapping Static transformers can be called with rows" begin

    y1 = source(10*ones(Int, 3))
    y2 = source(20*ones(Int, 3))
    y3 = source(30*ones(Int, 3))

    permuter = PermuteArgs((2, 3, 1))
    mach = machine(permuter)
    y = transform(mach, y1, y2, y3)

    @test y(rows=1:2) == (20*ones(Int, 2), 30*ones(Int, 2), 10*ones(Int, 2))

end

@testset "overloading methods for AbstractNode" begin

    A  = rand(3,7)
    As = source(A)
    @test MLJBase.matrix(MLJBase.table(As))() == A

    X = (x1 = [1,2,3], x2=[10, 20, 30], x3=[100, 200, 300])
    Xs = source(X)
    @test selectrows(Xs, 1)() == selectrows(X, 1)
    @test selectrows(Xs, 2:3)() == selectrows(X, 2:3)
    @test selectcols(Xs, 1)() == selectcols(X, 1)
    @test selectcols(Xs, 2:3)() == selectcols(X, 2:3)
    @test selectcols(Xs, :x1)() == selectcols(X, :x1)
    @test selectcols(Xs, [:x1, :x3])() == selectcols(X, [:x1, :x3])
    @test nrows(Xs)() == 3

    y = rand(4)
    ys = source(y)
    @test vcat(ys, ys)() == vcat(y, y)
    @test hcat(ys, ys)() == hcat(y, y)
    @test log(ys)() == log.(y)
    @test exp(ys)() == exp.(y)

    Z = (rand(4), rand(4), rand(4))
    Zs = source(Z)
    @test mean(Zs)() == mean.(Z)
    @test mode(Zs)() == mode.(Z)
    @test median(Zs)() == median.(Z)

    a, b, λ = rand(4), rand(4), rand()
    as, bs = source(a), source(b)
    @test (as + bs)() == a + b
    @test (λ * bs)() == λ * b

    X = source(1)
    y = source(7)
    @test glb(X, y)() == (1, 7)
    glb_node = @tuple X y
    @test glb_node() == (1, 7)

    X = source(1:10)
    Y = node(X->selectrows(X, 3:4), X)
    @test nrows_at_source(Y) == 10
    @test nrows(Y)() == 2

end

@testset "reformat logic - with upstream transformer learning from data" begin
    X = (x1=rand(5), x2=rand(5))
    y = categorical(collect("abaaa"))

    Xs = source(X)
    ys = source(y)

    std = Standardizer()
    mach1 = machine(std, Xs)
    W = transform(mach1, Xs)

    # a classifier with reformat front-end:
    clf = ConstantClassifier(testing=true)
    mach2 = machine(clf, W, ys; cache=true)
    yhat = predict(mach2, W)
    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(yhat, verbosity=0, rows=1:3))

    # training network with new rows changes upstream state of
    # classifier and hence retriggers reformatting of data:
    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(yhat, verbosity=0, rows=1:2))

    # however just changing classifier hyper-parameter avoids
    # reformatting and resampling:
    clf.bogus = 123
    @test_logs fit!(yhat, verbosity=0, rows=1:2)
end

mutable struct Scale <: MLJBase.Static
    scaling::Float64
end

function MLJBase.transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X * s.scaling
    MLJBase.table(s.scaling * MLJBase.matrix(X), prototype=X)
end

function MLJBase.inverse_transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X / s.scaling
    MLJBase.table(MLJBase.matrix(X) / s.scaling, prototype=X)
end

@testset "reformat logic - with upstream Static transformer" begin
    X = (x1=ones(5), x2=ones(5))
    y = categorical(collect("abaaa"))

    Xs = source(X)
    ys = source(y)

    scaler = Scale(2.0)
    mach1 = machine(scaler)
    W = transform(mach1, Xs)

    # a classifier with reformat front-end:
    clf = ConstantClassifier(testing=true)
    mach2 = machine(clf, W, ys, cache=true)
    yhat = predict(mach2, W)
    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(yhat, verbosity=0, rows=1:3))

    # training network with new rows does not change upstream state of
    # classifier, because `scaler isa Static` and no reformatting of
    # data:
    @test_logs((:info, "resampling X, y"),
               fit!(yhat, verbosity=0, rows=1:2)
               )

    # however changing an upstream hyperparameter forces reforamatting:
    scaler.scaling = 3.0
    @test_logs((:info, "reformatting X, y"),
               (:info, "resampling X, y"),
               fit!(yhat, verbosity=0, rows=1:2))

    # however just changing classifier hyper-parameter avoids
    # reformatting and resampling:
    clf.bogus = 123
    @test_logs fit!(yhat, verbosity=0, rows=1:2)
end

@testset "@node" begin
    X1 = source(4)
    X2 = source(5)
    Xp = source(1:10)

    add(a, b, c) = a + b + c
    N = @node add(X1, 1, X2)
    @test N() == 10

    N = @node tuple(X1, 5, X1)
    @test N() == (4, 5, 4)

    Y = @node selectrows(Xp, 3:4)
    @test Y() == 3:4
    @test Y([:one, :two, :three, :four]) == [:three, :four]
end

end

true


text:
using Test
using Distributions
using Random

using EnsembleKalmanProcesses.DataContainers

@testset "DataContainers" begin
    rng = Random.MersenneTwister(2021)

    parameter_samples = rand(rng, MvNormal(2, 0.1), 10) #10 samples of 4D params
    data_samples = rand(rng, MvNormal(12, 2), 10) #10 samples of 12D data
    data_samples_short = data_samples[:, 1:(end - 1)]

    new_parameter_samples = rand(rng, MvNormal(2, 0.1), 10) #10 samples of 4D params
    new_data_samples = rand(rng, MvNormal(12, 2), 10) #10 samples of 12D data
    new_data_samples_short = new_data_samples[:, 1:(end - 1)]

    #test DataContainer
    parameter_samples_T = permutedims(parameter_samples, (2, 1))
    idata = DataContainer(parameter_samples)
    idata_T = DataContainer(parameter_samples_T, data_are_columns = false)

    @test get_data(idata) == parameter_samples
    @test get_data(idata_T) == parameter_samples
    @test size(idata) == size(parameter_samples)

    #test Paired Container
    idata = DataContainer(parameter_samples)
    odata = DataContainer(data_samples)
    odata_short = DataContainer(data_samples_short)

    @test_throws DimensionMismatch PairedDataContainer(parameter_samples, data_samples_short)
    iopairs = PairedDataContainer(parameter_samples, data_samples)

    @test_throws DimensionMismatch PairedDataContainer(idata, odata_short)
    iopairs2 = PairedDataContainer(idata, odata)

    @test get_data(iopairs) == get_data(iopairs2)
    @test get_inputs(iopairs) == parameter_samples
    @test get_outputs(iopairs) == data_samples

    # test deep-copying
    parameter_samples[1, 1] += 1
    @test !isequal(get_data(idata), parameter_samples)
    @test !isequal(get_data(idata_T), parameter_samples)
    retrieved_samples = get_data(idata)
    retrieved_samples[1, 1] += 1
    @test !isequal(retrieved_samples, get_data(idata))
    retrieved_samples_T = get_data(idata_T)
    retrieved_samples_T[1, 1] += 1
    @test !isequal(retrieved_samples_T, get_data(idata_T))

end


text:
using Random
import StatsBase
export Daisy, Land, DaisyWorld

mutable struct Daisy <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    breed::Symbol
    age::Int
    albedo::Float64 # 0-1 fraction
end

mutable struct Land <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    temperature::Float64
end

const DaisyWorld = ABM{Union{Daisy, Land}};

"""
``` julia
daisyworld(;
    griddims = (30, 30),
    max_age = 25,
    init_white = 0.2,
    init_black = 0.2,
    albedo_white = 0.75,
    albedo_black = 0.25,
    surface_albedo = 0.4,
    solar_change = 0.005,
    solar_luminosity = 1.0,
    scenario = :default,
    seed = 165
)
```
Same as in [Daisyworld](@ref).

To access the `Daisy` and `Land` types, simply call
``` julia
using Agents.Models: Daisy, Land
```
"""
function daisyworld(;
        griddims = (30, 30),
        max_age = 25,
        init_white = 0.2, # % cover of the world surface of white breed
        init_black = 0.2, # % cover of the world surface of black breed
        albedo_white = 0.75,
        albedo_black = 0.25,
        surface_albedo = 0.4,
        solar_change = 0.005,
        solar_luminosity = 1.0, # initial luminosity
        scenario = :default,
        seed = 165
    )

    Random.seed!(seed)
    space = GridSpace(griddims, moore = true, periodic = true)
    properties = Dict(
        :max_age => max_age,
        :surface_albedo => surface_albedo,
        :solar_luminosity => solar_luminosity,
        :solar_change => solar_change,
        :scenario => scenario,
        :tick => 0
    )
    ## create a scheduler that only schedules Daisies
    daisysched(model) = [a.id for a in allagents(model) if a isa Daisy]
    model = ABM(Union{Daisy, Land}, space;
        scheduler = daisysched, properties = properties, warn = false
    )

    ## fill model with `Land`: every grid cell has 1 land instance
    fill_space!(Land, model, 0.0) # zero starting temperature

    ## Populate with daisies: each cell has only one daisy (black or white)
    white_nodes = StatsBase.sample(1:nv(space), Int(init_white * nv(space)); replace = false)
    for n in white_nodes
        wd = Daisy(nextid(model), vertex2coord(n, space), :white, rand(0:max_age), albedo_white)
        add_agent_pos!(wd, model)
    end
    allowed = setdiff(1:nv(space), white_nodes)
    black_nodes = StatsBase.sample(allowed, Int(init_black * nv(space)); replace = false)
    for n in black_nodes
        wd = Daisy(nextid(model), vertex2coord(n, space), :black, rand(0:max_age), albedo_black)
        add_agent_pos!(wd, model)
    end

    return model, daisyworld_agent_step!, daisyworld_model_step!
end

function update_surface_temperature!(node::Int, model::DaisyWorld)
    ids = get_node_contents(node, model)
    ## All grid points have at least one agent (the land)
    absorbed_luminosity = if length(ids) == 1
        ## Set luminosity via surface albedo
        (1 - model.surface_albedo) * model.solar_luminosity
    else
        ## more than 1 agents: daisy exists
        ## Set luminosity via daisy albedo
        (1 - model[ids[2]].albedo) * model.solar_luminosity
    end
    ## We expect local heating to be 80 ᵒC for an absorbed luminosity of 1,
    ## approximately 30 for 0.5 and approximately -273 for 0.01.
    local_heating = absorbed_luminosity > 0 ? 72 * log(absorbed_luminosity) + 80 : 80
    ## Surface temperature is the average of the current temperature and local heating.
    T0 = model[ids[1]].temperature
    model[ids[1]].temperature = (T0 + local_heating) / 2
end

function diffuse_temperature!(node::Int, model::DaisyWorld)
    ratio = get(model.properties, :ratio, 0.5) # diffusion ratio
    ids = space_neighbors(node, model)
    meantemp = sum(model[i].temperature for i in ids if model[i] isa Land)/8
    land = model[get_node_contents(node, model)[1]] # land at current node
    ## Each neighbor land patch is giving up 1/8 of the diffused
    ## amount to each of *its* neighbors
    land.temperature = (1 - ratio)*land.temperature + ratio*meantemp
end

function propagate!(node::Int, model::DaisyWorld)
    ids = get_node_contents(node, model)
    if length(ids) > 1
        daisy = model[ids[2]]
        temperature = model[ids[1]].temperature
        ## Set optimum growth rate to 22.5 ᵒC, with bounds of [5, 40]
        seed_threshold = (0.1457 * temperature - 0.0032 * temperature^2) - 0.6443
        if rand() < seed_threshold
            ## Collect all adjacent cells that have no daisies
            empty_neighbors = Int[]
            neighbors = node_neighbors(node, model)
            for n in neighbors
                if length(get_node_contents(n, model)) == 1
                    push!(empty_neighbors, n)
                end
            end
            if !isempty(empty_neighbors)
                ## Seed a new daisy in one of those cells
                seeding_place = vertex2coord(rand(empty_neighbors), model)
                a = Daisy(nextid(model), seeding_place, daisy.breed, 0, daisy.albedo)
                add_agent_pos!(a, model)
            end
        end
    end
end

function daisyworld_agent_step!(agent::Daisy, model::DaisyWorld)
    agent.age += 1
    agent.age >= model.max_age && kill_agent!(agent, model)
end

daisyworld_agent_step!(agent::Land, model::DaisyWorld) = nothing

function daisyworld_model_step!(model)
    for n in nodes(model)
        update_surface_temperature!(n, model)
        diffuse_temperature!(n, model)
        propagate!(n, model)
    end
    model.tick += 1
    solar_activity!(model)
end

function solar_activity!(model::DaisyWorld)
    if model.scenario == :ramp
        if model.tick > 200 && model.tick <= 400
            model.solar_luminosity += model.solar_change
        end
        if model.tick > 500 && model.tick <= 750
            model.solar_luminosity -= model.solar_change/2
        end
    elseif model.scenario == :change
        model.solar_luminosity += model.solar_change
    end
end


text:
module OneDTmazeERExperiment

import Flux: Descent, ADAM
import Random

using GVFHordes
using Curiosity
using MinimalRLCore
using SparseArrays
using ProgressMeter

const BU = Curiosity.BaselineUtils
const ODTMU = Curiosity.OneDTMazeUtils
const SRCU = Curiosity.SRCreationUtils

default_args() =
    Dict(
        "logger_interval" => 100,

        # Behaviour Items
        # "behaviour_eta" => 0.1/8,
        "behaviour_gamma" => 0.9,
        "behaviour_learner" => "RoundRobin",
        "behaviour_update" => "TB",
        "behaviour_reward_projector" => "maze",
        "behaviour_rp_tilings" => 1,
        "behaviour_rp_tiles" => 16,
        "behaviour_trace" => "ReplacingTraces",
        "behaviour_opt" => "Auto",
        "behaviour_lambda" => 0.95,
        "behaviour_alpha_init" => 0.1,
        "exploration_param" => 1.0,
        "exploration_strategy" => "epsilon_greedy",
        "ϵ_range" => (0.4,0.1),
        "decay_period" => 5000,
        "warmup_steps" => 12000,
        "behaviour_w_init" => 4,

        # Demon Attributes
        # "demon_eta" => 0.1/8,
        "demon_discounts" => 0.9,
        "demon_learner" => "SR",
        "demon_update" => "TB",
        "demon_interest_set" => "oned_tmaze",
        "demon_policy_type" => "greedy_to_cumulant",
        "demon_opt" => "Auto+ADAM",
        "demon_lambda" => 0.9,
        "demon_trace"=> "AccumulatingTraces",

        "demon_eta" => 0.001,
        "demon_alpha_init"=>0.1,
        "demon_beta" => (0.99, 0.999),

        #shared
        # "num_tiles" => 2,
        # "num_tilings" =>8,
        "tiling_structure" => [2,8],
        "demon_rep" => "ideal",
        # "demon_rep" => "tilecoding",
        "demon_num_tiles" => 8,
        "demon_num_tilings" => 1,
        "eta" => 0.05,
        "alpha_init" => 0.1,
        "replay_size" => 10000,
        "batch_size" => 4,

        # Environment Config
        "constant_target"=> [-10.0,10.0],
        "cumulant_schedule" => "DrifterDistractor",
        "cumulant" => 1,
        "distractor" => (1.0, 5.0),
        "drifter" => (1.0, sqrt(0.01)),
        "exploring_starts"=>"beg",
        "env_step_penalty" => -0.01,

        # Agent and Logger
        "horde_type" => "regular",
        "intrinsic_reward" => "weight_change",
        # "logger_keys" => [LoggerKey.TTMAZE_ERROR],
        "save_dir" => "OneDTMazeExperiment",
        "seed" => 1,
        "steps" => 20000,
        "use_external_reward" => true,
        "random_first_action" => false,
        "logger_keys" => [LoggerKey.ONEDTMAZEERROR,LoggerKey.ONEDTMAZEERROR_DPI,LoggerKey.ONEDTMAZEERROR_UNIFORM, LoggerKey.ONED_GOAL_VISITATION, LoggerKey.EPISODE_LENGTH, LoggerKey.INTRINSIC_REWARD]
    )

function construct_agent(parsed)
    action_space = 4
    obs_size = 6

    intrinsic_reward_type = parsed["intrinsic_reward"]
    use_external_reward = parsed["use_external_reward"]

    # Unpack Tiling Config
    if "tiling_structure" ∈ keys(parsed)
        parsed["num_tilings"] = parsed["tiling_structure"][1]
        parsed["num_tiles"] = parsed["tiling_structure"][2]
    end

    # Unpack joint eta param
    if "eta" in keys(parsed)
        prefixes = ["behaviour","demon"]
        for prefix in prefixes
            parsed[join([prefix, "eta"], "_")] = parsed["eta"]
        end
    end

    if "alpha_init" in keys(parsed)
        prefixes = ["behaviour","demon"]
        for prefix in prefixes
            parsed[join([prefix, "alpha_init"], "_")] = parsed["alpha_init"]
        end
    end

    if parsed["demon_opt"] == "Auto+Descent"
        parsed["demon_auto_alpha_init"] = parsed["demon_auto_alpha_init"] / parsed["num_tilings"]
    end

    # Unpack Auto Init
    if parsed["demon_opt"] == "Auto"
        parsed["demon_alpha_init"] =
            parsed["demon_alpha_init"] / parsed["num_tilings"]
    end

    if "behaviour_opt" ∈ keys(parsed) && parsed["behaviour_opt"] == "Auto"
        parsed["behaviour_alpha_init"] =
            parsed["behaviour_alpha_init"] / parsed["num_tilings"]
    end
    if parsed["behaviour_opt"] == "Descent"
        if "behaviour_eta" in keys(parsed)
            parsed["behaviour_eta"] = parsed["behaviour_eta"] / parsed["num_tilings"]
        end
    end

    if parsed["behaviour_w_init"] != 0.0
        parsed["behaviour_w_init"] = parsed["behaviour_w_init"] / parsed["num_tilings"]
    end

    fc = Curiosity.FeatureSubset(
        Curiosity.SparseTileCoder(parsed["num_tilings"], parsed["num_tiles"], 2),
        1:2)
    feat_size = size(fc)

    demon_projected_fc = if "demon_rep" ∉ keys(parsed)
        nothing
    elseif parsed["demon_rep"] == "tilecoding"
        Curiosity.FeatureProjector(Curiosity.FeatureSubset(
                Curiosity.SparseTileCoder(parsed["demon_num_tilings"], parsed["demon_num_tiles"], 2),
            1:2), false)
    elseif parsed["demon_rep"] == "ideal"
        Curiosity.FeatureSubset(ODTMU.IdealDemonFeatures(), 1:2)
    elseif parsed["demon_rep"] == "ideal_state_action"
        t = ODTMU.IdealStateActionDemonFeatures(4)
        Curiosity.FeatureSubset(t, 1:2)
    elseif parsed["demon_rep"] == "ideal_martha"
        Curiosity.FeatureProjector(Curiosity.FeatureSubset(ODTMU.MarthaIdealDemonFeatures(4), 1:2), false)
    else
        throw(ArgumentError("Not a valid demon projection rep for SR"))
    end

    demons = ODTMU.create_demons(parsed, demon_projected_fc)

    demon_learner = Curiosity.get_linear_learner(parsed,
                                feat_size,
                                action_space,
                                demons,
                                "demon",
                                demon_projected_fc)

    exploration_strategy = Curiosity.get_exploration_strategy(parsed, 1:action_space)

    behaviour_reward_projector = if "behaviour_reward_projector" ∉ keys(parsed)
        nothing
    elseif parsed["behaviour_reward_projector"] == "tilecoding"
        Curiosity.FeatureProjector(Curiosity.FeatureSubset(
            Curiosity.SparseTileCoder(parsed["behaviour_rp_tilings"], parsed["behaviour_rp_tiles"], 2),
            1:2), false)
    elseif parsed["behaviour_reward_projector"] == "base"
        Curiosity.FeatureProjector(fc, false)
    elseif parsed["behaviour_reward_projector"] == "identity"
        Curiosity.FeatureProjector(Curiosity.FeatureSubset(identity, 1:2), false)
    elseif parsed["behaviour_reward_projector"] == "ideal"
        Curiosity.FeatureProjector(Curiosity.FeatureSubset(
            ODTMU.IdealDemonFeatures(), 1:2), true)
    elseif parsed["behaviour_reward_projector"] == "ideal_martha"
        Curiosity.FeatureProjector(Curiosity.FeatureSubset(
            ODTMU.MarthaIdealDemonFeatures(), 1:2), false)
    elseif parsed["behaviour_reward_projector"] == "maze"
        Curiosity.FeatureSubset(ODTMU.TMazeEncoding(),1:2)
    else
        throw(ArgumentError("Not a valid demon projection rep for SR"))
    end

    behaviour_learner, behaviour_demons, behaviour_discount = if parsed["behaviour_learner"] == "RoundRobin"
        ODTMU.RoundRobinPolicy(), nothing, 0.0
    elseif parsed["behaviour_learner"] == "RandomDemons"
        action_set = 1:action_space
        blearner = if parsed["demon_learner"] == "SR"
            #We want the demons specifying the cumulants and not predicting the SF
            BU.RandomDemons(demons.PredHorde, action_set)
        else
            BU.RandomDemons(demons, action_set)
        end
        blearner, nothing, 0.0
    else
        behaviour_num_tasks = 1
        num_SFs = 4
        num_demons = if parsed["behaviour_learner"] ∈ ["GPI"]
            # num_SFs * size(behaviour_reward_projector) * action_space + behaviour_num_tasks
            num_SFs * size(behaviour_reward_projector) + behaviour_num_tasks

        elseif parsed["behaviour_learner"] ∈ ["Q"]
            behaviour_num_tasks
        elseif parsed["behaviour_learner"] == "RoundRobin"
            0
        else
            throw(ArgumentError("Hacky thing not working, yay!"))
        end
        behaviour_learner = Curiosity.get_linear_learner(parsed,
                                                         size(fc),
                                                         action_space,
                                                         num_demons,
                                                         behaviour_num_tasks,
                                                         "behaviour",
                                                         behaviour_reward_projector)

        behaviour_demons = if behaviour_learner isa GPI
            @assert !(behaviour_reward_projector isa Nothing)
            bh_gvf = ODTMU.make_behaviour_gvf(behaviour_learner,
                                              0.0,
                                              fc,
                                              exploration_strategy)
            pred_horde = GVFHordes.Horde([bh_gvf])
            SF_policies = [ODTMU.GoalPolicy(i) for i in 1:4]
            SF_discounts = [ODTMU.GoalTermination(parsed["behaviour_gamma"]) for i in 1:4]
            num_SFs = length(SF_policies)
            SF_horde = SRCU.create_SF_horde(SF_policies, SF_discounts, behaviour_reward_projector, 1:action_space)
            Curiosity.GVFSRHordes.SRHorde(pred_horde, SF_horde, num_SFs, behaviour_reward_projector)
        elseif behaviour_learner isa QLearner
            bh_gvf = ODTMU.make_behaviour_gvf(behaviour_learner,
                                              parsed["behaviour_gamma"],
                                              fc,
                                              exploration_strategy)
            GVFHordes.Horde([bh_gvf])
        elseif behaviour_learner isa ODTMU.RoundRobinPolicy
            nothing
        else
            throw(ArgumentError("goes with which horde? " ))
        end

        behaviour_learner, behaviour_demons, parsed["behaviour_gamma"]
    end #end behaviour_learner, behaviour_demons, behaviour_discount = if parsed["behaviour_learner"] == "RoundRobin"

    random_first_action = parsed["random_first_action"]
    if behaviour_learner isa ODTMU.RoundRobinPolicy
        if random_first_action == false
            @warn "Round Robin with random first action is recommended"
        end
    elseif random_first_action == true
        @warn "Random first action is enabled"
    end

    AgentER(demons,
            feat_size,
            behaviour_learner,
            behaviour_demons,
            behaviour_discount,
            demon_learner,
            obs_size,
            action_space,
            intrinsic_reward_type,
            fc,
            use_external_reward,
            exploration_strategy,
            random_first_action,
            parsed["batch_size"],
            parsed["replay_size"])
end

function main_experiment(parsed=default_args(); progress=false, working=false)

    num_steps = parsed["steps"]
    Random.seed!(parsed["seed"])

    logger_init_dict = Dict(
        LoggerInitKey.TOTAL_STEPS => num_steps,
        LoggerInitKey.INTERVAL => parsed["logger_interval"],
        # LoggerInitKey.ENV => "tabular_tmaze"
    )

    Curiosity.experiment_wrapper(parsed, logger_init_dict, working) do parsed, logger

        cumulant_schedule = ODTMU.get_cumulant_schedule(parsed)
        env = OneDTMaze(cumulant_schedule, parsed["exploring_starts"], parsed["env_step_penalty"])

        agent = construct_agent(parsed)
        goal_visitations = zeros(4)

        eps = 1
        max_num_steps = num_steps
        steps = Int[]

        logger_start!(logger, env, agent)

        prg_bar = ProgressMeter.Progress(num_steps, "Step: ")
        while sum(steps) < max_num_steps
            cur_step = 0
            max_episode_steps = max_num_steps - sum(steps)
            tr, stp =
                run_episode!(env, agent, max_episode_steps) do (s, a, s_next, r, t)
                    #This is a callback for every timestep where logger can go
                    # agent is accesible in this scope
                    if t == true && working==true
                        goals = s_next[3:end]
                        f = findfirst(!iszero, goals)
                        goal_visitations[f] += 1
                    end
                    if progress
                        next!(prg_bar)
                    end

                    logger_step!(logger, env, agent, s, a, s_next, r, t)
                    cur_step+=1
                end
                logger_episode_end!(logger)
            push!(steps, stp)
            eps += 1
        end
        if working == true
            println("Goal Visitation: ", goal_visitations)
        end
        # agent
    end

end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function get_field(::Type{GeoMagneticFieldFactory}, arg0::GeoMagneticFieldFactory_FieldModel, arg1::jdouble)
    return jcall(GeoMagneticFieldFactory, "getField", GeoMagneticField, (GeoMagneticFieldFactory_FieldModel, jdouble), arg0, arg1)
end

function get_geo_magnetic_fields(::Type{GeoMagneticFieldFactory})
    return jcall(GeoMagneticFieldFactory, "getGeoMagneticFields", LazyLoadedGeoMagneticFields, ())
end

function get_igrf(::Type{GeoMagneticFieldFactory}, arg0::jdouble)
    return jcall(GeoMagneticFieldFactory, "getIGRF", GeoMagneticField, (jdouble,), arg0)
end

function get_wmm(::Type{GeoMagneticFieldFactory}, arg0::jdouble)
    return jcall(GeoMagneticFieldFactory, "getWMM", GeoMagneticField, (jdouble,), arg0)
end


text:
struct Summary
    μ
    σ
    n

    Summary(data::Matrix{Float64}) = new(
        mean(data,dims=2),
        std(data,corrected=true,dims=2) / √(size(data)[2]-1),
        size(data)[2]
    )

    Summary(data::Vector{Float64}) = new(
        data,
        zeros(length(data)),
        1
    )
end

# Saving/Loading data
save!(self::ConfigManager, data::Dict) = save!(self.dataManager, get_datafile(self), data)
save(self::ConfigManager, data::Dict) = save(self.dataManager, get_datafile(self), data)
load(self::ConfigManager) = load(self.dataManager, get_datafile(self))
load(self::ConfigManager, idx::Int) = load(parse!(self, idx))

# Parsing saved data

function get_run_data(self::ConfigManager, idx, run)
    parse!(self, idx, run)
    return load(self)
end

function get_data(self::ConfigManager, idx)
    runs = get_max_runs(self,idx)
    all_data = Dict()
    for r in 1:runs
        data = get_run_data(self, idx, r)
        for k in keys(data)
            if haskey(all_data, k)
                all_data[k] = hcat(all_data[k], data[k])
            else
                all_data[k] = data[k]
            end
        end
    end
    return all_data
end

function get_fields(self::ConfigManager)
    return keys(get_run_data(copy(self),1,1))
end

function get_summary(self::ConfigManager, idx::Int, key::String)
    data = get_data(self,idx)
    return Summary(data[key])
end

function get_best(self::ConfigManager, key; metric=mean, selector=argmin)
    return get_best(self, key, 1:self.total_combinations, metric=metric, selector=selector)
end

function get_best(self::ConfigManager, key, indices; metric=mean, selector=argmin)
    data = []
    parsed_indices = []
    for idx in indices
        d = get_summary(self, idx, key).μ
        m = metric(d)
        if isfinite(m)
            push!(data, metric(d))
            push!(parsed_indices, idx)
        end
    end
    if length(data) == 0
        return nothing # all configs diverged
    end
    best = selector(data)
    return parsed_indices[best]
end

function sensitivity(self::ConfigManager, key; metric=mean)
    best_idx = get_best(self,key)
    cfg = parse!(self, best_idx)

    keys = get_sweep_params(self)
    d = Dict( k=> self[k] for k in keys )

    # delete the key we're looking for sensitivity over
    delete!(d,key)

    indices = configs_with(self, d)
    sens_data = []
    for idx in indices
        value = metric(get_summary(self,idx).μ)
        push!(sens_data, value)
    end
    return sens_data
end

function get_max_runs(self::ConfigManager, idx::Int)
    cfg = parse!(self,idx,1)
    runs = readdir(get_expdir(self))
    return length(runs)
end


text:
# this generated the error
#
# ERROR: LoadError: MethodError: no method matching init_node_sr_belief(::POMCPOW.POWNodeFilter, ::SimplePOMDP, ::Int64, ::Int64, ::Int64, ::Float64, ::Int64)

struct SimplePOMDP <: POMDP{Int,Int,Float64}
    p_success::Float64
end

function POMDPs.transition(m::SimplePOMDP, s::Int, a::Int)
    return SparseCat(clamp.([s + a, s - a], 1, 7), [m.p_success, 1 - m.p_success])
end

function POMDPs.observation(m::SimplePOMDP, a::Int, sp::Int)
    return SparseCat([sp + 0.5, sp - 0.5], [0.5, 0.5])
end

function POMDPs.reward(m::SimplePOMDP, s::Int, a::Int, sp::Int)
    if s == 1 || s == 7
        return 10
    elseif s == 4
        return -10
    else
        return 0
    end
end

function POMDPs.initialstate_distribution(m::SimplePOMDP)
    return SparseCat(1:7, ones(7) ./ 7)
end

function POMDPs.discount(m::SimplePOMDP)
    return 0.9
end

POMDPs.actions(m::SimplePOMDP) = [-1, 1]
POMDPs.states(m::SimplePOMDP) = 1:7
POMDPs.actionindex(m::SimplePOMDP, a::Int) = a == 1 ? 1 : 2
POMDPs.n_states(m::SimplePOMDP) = 7
POMDPs.n_actions(m::SimplePOMDP) = 2

pomdp = SimplePOMDP(0.7)

s = 5
b = SparseCat([s], [1.0])
solver = POMCPOWSolver()

planner = solve(solver, pomdp)
a, info = action_info(planner, b)


text:
export genH

function genH(n::Integer, c::Char)
	if c == 'X'
		return [chain(n, put(n, i=>X)) for i = 1:n]
	elseif c == 'Y'
		return [chain(n, put(n, i=>Y)) for i = 1:n]
	elseif c == 'Z'
		return [chain(n, put(n, i=>Z)) for i = 1:n]
	end
end


text:
function interpolator_error(T::Type, f::Function, xa, xb, n::Int)
    ϕ = T(f, xa, xb, n)
    x = collect(LinRange(xa, xb, 100))
    fy = f.(x)
    ϕy = ϕ.(x)
    maximum(abs.(ϕy .- fy))
end

function interpolator_error(T::Type, f::Function, xa, xb, nx::Int, ya, yb, ny::Int)
    Φ = T(f, xa, xb, nx, ya, yb, ny)
    x = collect(LinRange(xa, xb, 50))
    y = collect(LinRange(ya, yb, 50))
    X = x .* ones(length(y))'
    Y = y' .* ones(length(x))
    maximum(@. abs(Φ(X, Y) - f(X, Y)) )
end

T1 = [LinearInterpolator,
      CubicInterpolator,
      CubicSplineInterpolator,
      ChebyshevInterpolator]
T2 = [BilinearInterpolator,
      BicubicInterpolator,
      BichebyshevInterpolator]

for T ∈ T1
    @test interpolator_error(T, x->sin(x^2), -π, π, 128) < 1e-2
end

for T ∈ T2
    @test interpolator_error(T, (x,y)->sin(x)+cos(y), -π, π, 128, -π, π, 128) < 1e-2
end


text:
module GeneratedTables

using GeneratedTypes

import Base: getindex, unsafe_getindex, setindex!, unsafe_setindex!

export Cell, Column, Row, Table
export name, eltypes

include("cell.jl")
include("column.jl")
include("row.jl")
include("table.jl")
include("extra_constructors.jl")
include("show.jl")

end # module


text:
"""
clip_bbox(x::GeoArray, b::bbox)
clip_bbox(x::GeoArray, y::GeoArray)
"""
function st_clip(ga::GeoArray, b::bbox)
    x, y = st_dim(ga)
    I_x = findall((x .<= b.xmax) .& (x .>= b.xmin))
    I_y = findall((y .<= b.ymax) .& (y .>= b.ymin))
    ga[I_x, I_y]
end

function st_clip(x::GeoArray, y::GeoArray)
    st_clip(x, st_bbox(y))
end

export st_clip