text:
#=

# 406: 1D Weird Surface Reaction
 ([source code](SOURCE_URL))

Species $A$ and $B$ exist in the interior of the domain, species $C$
lives a the boundary $\Gamma_1$.  We assume a heterogeneous reaction scheme
where $A$ reacts to $B$ with a rate depending on $\nabla A$ near the surface

```math
\begin{aligned}
      A &\leftrightarrow B\\
\end{aligned}
```

In $\Omega$, both $A$ and $B$ are transported through diffusion:

```math
\begin{aligned}
\partial_t u_B - \nabla\cdot D_A \nabla u_A & = f_A\\
\partial_t u_B - \nabla\cdot D_B \nabla u_B & = 0\\
\end{aligned}
```
Here, $f(x)$ is a source term creating $A$.
On $\Gamma_2$, we set boundary conditions
```math
\begin{aligned}
D_A \nabla u_A & = 0\\
u_B&=0
\end{aligned}
```
describing no normal flux for $A$ and zero concentration of $B$.
On $\Gamma_1$, we use the mass action law to describe the boundary reaction and
the evolution of the boundary concentration $C$. We assume that there is a limited
amount of surface sites $S$ for species C, so in fact A has to react with a free
surface site in order to become $C$ which reflected by the factor $1-u_C$. The same
is true for $B$.
```math
\begin{aligned}
R_{AB}(u_A, u_B)&=k_{AB}^+exp(u_A'(0))u_A - k_{AB}^-exp(-u_A'(0))u_B\\
- D_A \nabla u_A  +  R_{AB}(u_A, u_B)& =0 \\
- D_B \nabla u_B  -  R_{AB}(u_A, u_B)& =0 \\
\end{aligned}
```

=#
module Example406_WeirdReaction
using Printf
using VoronoiFVM
using SparseArrays
using ExtendableGrids
using GridVisualize

function main(;n=10,
              Plotter=nothing,
              verbose=false,
              tend=1,
              unknown_storage=:sparse,
              autodetect_sparsity=true)

    h=1.0/convert(Float64,n)
    X=collect(0.0:h:1.0)
    N=length(X)

    grid=VoronoiFVM.Grid(X)
    ## By default, \Gamma_1 at X[1] and \Gamma_2 is at X[end]

    ## Species numbers
    iA=1
    iB=2
    iC=3

    ## Diffusion flux for species A and B
    D_A=1.0
    D_B=1.0e-2
    function flux!(f,u,edge)
        f[iA]=D_A*(u[iA,1]-u[iA,2])
        f[iB]=D_B*(u[iB,1]-u[iB,2])
    end

    ## Storage term of species A and B
    function storage!(f,u,node)
        f[iA]=u[iA]
        f[iB]=u[iB]
    end

    ## Source term for species a around 0.5
    function source!(f,node)
        x1=node[1]-0.5
        f[iA]=exp(-100*x1^2)
    end

    ## Reaction constants (p = + , m = -)
    ## Choosen to prefer path A-> B
    kp_AB=1.0
    km_AB=0.1

    function breaction!(f,u,node)
        if  node.region==1
            R=kp_AB*exp(u[iC])*u[iA] - exp(-u[iC])*km_AB*u[iB]
            f[iA]+=R
            f[iB]-=R
        end
    end

    ## This generic operator works on the full solution seen as linear vector, and indexing
    ## into it needs to be performed with the help of idx (defined below for a solution vector)
    ## Its sparsity is detected automatically using SparsityDetection.jl
    ## Here, we calculate the gradient of u_A at the boundary and store the value in u_C which
    ## is then used as a parameter in the boundary reaction
    function generic_operator!(f,u,sys)
        f.=0
        f[idx[iC,1]]=u[idx[iC,1]]  + 0.1*(u[idx[iA,1]]-u[idx[iA,2]])/(X[2]-X[1])
    end

    # If we know the sparsity pattern, we can here create a
    # sparse matrix with values set to 1 in the nonzero
    # slots. This allows to circumvent the
    # autodetection which may takes some time.
    function generic_operator_sparsity(sys)
        idx=unknown_indices(unknowns(sys))
        sparsity=spzeros(num_dof(sys),num_dof(sys))
        sparsity[idx[iC,1],idx[iC,1]]=1
        sparsity[idx[iC,1],idx[iA,1]]=1
        sparsity[idx[iC,1],idx[iA,2]]=1
        sparsity
    end

    if autodetect_sparsity
        physics=VoronoiFVM.Physics(
            breaction=breaction!,
            generic=generic_operator!,
            flux=flux!,
            storage=storage!,
            source=source!
        )
    else
        physics=VoronoiFVM.Physics(
            breaction=breaction!,
            generic=generic_operator!,
            generic_sparsity=generic_operator_sparsity,
            flux=flux!,
            storage=storage!,
            source=source!
        )
    end

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    ## Enable species in bulk resp
    enable_species!(sys,iA,[1])
    enable_species!(sys,iB,[1])

    ## Enable surface species
    enable_boundary_species!(sys,iC,[1])

    ## Set Dirichlet bc for species B on \Gamma_2
    boundary_dirichlet!(sys,iB,2,0.0)

    ## Initial values
    inival=unknowns(sys)
    inival.=0.0
    U=unknowns(sys)
    idx=unknown_indices(U)

    tstep=0.01
    time=0.0
    T=Float64[]
    u_C=Float64[]

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    p=GridVisualizer(Plotter=Plotter,layout=(2,1))
    while time<tend
        time=time+tstep
        solve!(U,inival,sys,tstep=tstep,control=control)
        inival.=U
        if verbose
            @printf("time=%g\n",time)
        end
        ## Record  boundary pecies
        push!(T,time)
        push!(u_C,U[iC,1])

        scalarplot!(p[1,1],grid,U[iA,:],label="[A]",title=@sprintf("max_A=%.5f max_B=%.5f u_C=%.5f",maximum(U[iA,:]),maximum(U[iB,:]),u_C[end]),color=:red)
        scalarplot!(p[1,1],grid,U[iB,:], label="[B]",clear=false,color=:blue)
        scalarplot!(p[2,1],copy(T),copy(u_C),label="[C]",clear=true,show=true)
    end
    return U[iC,1]
end

function test()
    testval=0.007027597470502758
    main(unknown_storage=:sparse) ≈ testval &&
        main(unknown_storage=:dense) ≈ testval &&
        main(unknown_storage=:sparse,autodetect_sparsity=false) ≈ testval &&
        main(unknown_storage=:dense,autodetect_sparsity=false) ≈ testval
end

end


text:
using Compat

################################################################################
# Export definitions supported by our export script                            #
################################################################################
export RTP_PAYLOAD_SIZE

# typealias
export ToxAv
export ToxAVCallback
export ToxAvAudioCallback
export ToxAvVideoCallback
export ToxAvCallbackID
export ToxAvCallType
export ToxAvCallState
export ToxAvError
export ToxAvCapabilities

# const
export av_OnInvite
export av_OnRinging
export av_OnStart
export av_OnCancel
export av_OnReject
export av_OnEnd
export av_OnRequestTimeout
export av_OnPeerTimeout
export av_OnPeerCSChange
export av_OnSelfCSChange
export av_TypeAudio
export av_TypeVideo
export av_CallNonExistent
export av_CallInviting
export av_CallStarting
export av_CallActive
export av_CallHold
export av_CallHungUp
export av_ErrorNone
export av_ErrorUnknown
export av_ErrorNoCall
export av_ErrorInvalidState
export av_ErrorAlreadyInCallWithPeer
export av_ErrorReachedCallLimit
export av_ErrorInitializingCodecs
export av_ErrorSettingVideoResolution
export av_ErrorSettingVideoBitrate
export av_ErrorSplittingVideoPayload
export av_ErrorEncodingVideo
export av_ErrorEncodingAudio
export av_ErrorSendingPayload
export av_ErrorCreatingRtpSessions
export av_ErrorNoRtpSession
export av_ErrorInvalidCodecState
export av_ErrorPacketTooLarge
export av_AudioEncoding
export av_AudioDecoding
export av_VideoEncoding
export av_VideoDecoding

export ToxAvCSettings

################################################################################
# This is the code generated by Clang (Anonymous enums removed)                #
################################################################################
const RTP_PAYLOAD_SIZE = 65535

typealias ToxAv Void
typealias ToxAVCallback Ptr{Void}
typealias ToxAvAudioCallback Ptr{Void}
typealias ToxAvVideoCallback Ptr{Void}

# begin enum ToxAvCallbackID
typealias ToxAvCallbackID Uint32
const av_OnInvite = (UInt32)(0)
const av_OnRinging = (UInt32)(1)
const av_OnStart = (UInt32)(2)
const av_OnCancel = (UInt32)(3)
const av_OnReject = (UInt32)(4)
const av_OnEnd = (UInt32)(5)
const av_OnRequestTimeout = (UInt32)(6)
const av_OnPeerTimeout = (UInt32)(7)
const av_OnPeerCSChange = (UInt32)(8)
const av_OnSelfCSChange = (UInt32)(9)
# end enum ToxAvCallbackID

# begin enum ToxAvCallType
typealias ToxAvCallType Uint32
const av_TypeAudio = (UInt32)(192)
const av_TypeVideo = (UInt32)(193)
# end enum ToxAvCallType

# begin enum ToxAvCallState
typealias ToxAvCallState Cint
const av_CallNonExistent = (Int32)(-1)
const av_CallInviting = (Int32)(0)
const av_CallStarting = (Int32)(1)
const av_CallActive = (Int32)(2)
const av_CallHold = (Int32)(3)
const av_CallHungUp = (Int32)(4)
# end enum ToxAvCallState

# begin enum ToxAvError
typealias ToxAvError Cint
const av_ErrorNone = (Int32)(0)
const av_ErrorUnknown = (Int32)(-1)
const av_ErrorNoCall = (Int32)(-20)
const av_ErrorInvalidState = (Int32)(-21)
const av_ErrorAlreadyInCallWithPeer = (Int32)(-22)
const av_ErrorReachedCallLimit = (Int32)(-23)
const av_ErrorInitializingCodecs = (Int32)(-30)
const av_ErrorSettingVideoResolution = (Int32)(-31)
const av_ErrorSettingVideoBitrate = (Int32)(-32)
const av_ErrorSplittingVideoPayload = (Int32)(-33)
const av_ErrorEncodingVideo = (Int32)(-34)
const av_ErrorEncodingAudio = (Int32)(-35)
const av_ErrorSendingPayload = (Int32)(-40)
const av_ErrorCreatingRtpSessions = (Int32)(-41)
const av_ErrorNoRtpSession = (Int32)(-50)
const av_ErrorInvalidCodecState = (Int32)(-51)
const av_ErrorPacketTooLarge = (Int32)(-52)
# end enum ToxAvError

# begin enum ToxAvCapabilities
typealias ToxAvCapabilities Uint32
const av_AudioEncoding = (UInt32)(1)
const av_AudioDecoding = (UInt32)(2)
const av_VideoEncoding = (UInt32)(4)
const av_VideoDecoding = (UInt32)(8)
# end enum ToxAvCapabilities

################################################################################
# Added by hand                                                                #
################################################################################
immutable ToxAvCSettings
	call_type::ToxAvCallType
	video_bitrate::UInt32
	max_video_width::UInt16
	max_video_height::UInt16
	audio_bitrate::UInt32
	audio_frame_duration::UInt16
	audio_sample_rate::UInt32
	audio_channels::UInt32
end

function ToxAvCSettings()
	# some stupid default values
    ToxAvCSettings(av_TypeAudio, 0, 0, 0, 0, 0, 0, 0 )
end

ext: jl


text:
export getManifoldPartial

# forcing ProductManifold to use ProductRepr as accompanying representation
const _PartiableRepresentationProduct = Union{Nothing,<:ProductRepr}
# forcing ProductManifold to use ProductRepr as accompanying representation
const _PartiableRepresentationFlat{T} = Union{Nothing,<:AbstractVector{T}}
# More general representation for Manifold Factors or Groups
const _PartiableRepresentation = Union{<:_PartiableRepresentationProduct,<:_PartiableRepresentationFlat,<:AbstractMatrix}

## COMMON UTILS FOR PARTIAL MANIFOLDS

function _checkManifoldPartialDims( M::AbstractManifold,
                                    partial::AbstractVector{Int},
                                    offset::Base.RefValue{Int},
                                    doError::Bool=true)
  #
  mask = 0 .< (partial .- offset[]) .<= manifold_dimension(M)
  doError && !any(mask) && error("Unknown manifold partial=$(partial .- offset[]) in dimension=$(manifold_dimension(M)) of $M")
  return mask
end

## EXTRACT REPRESENTATION PARTIALS

# do nothing case
_getReprPartial(M::MB.AbstractManifold, ::Nothing, w...; kw...) = nothing

# drop back to regular vector
function _getReprPartial( M::MB.AbstractManifold,
                          repr::AbstractVector{T},
                          partial::AbstractVector{Int}, # total partial from user over all Factors
                          offset::Base.RefValue{Int}=Ref(0),
                          mask::BitVector=_checkManifoldPartialDims(M,partial,offset, doError) ) where {T <: Number}
  #
  ret = zeros(T, sum(mask))
  for (i,p) in enumerate(partial .- offset)
    ret[i] = repr[p]
  end
  return ret
end

function _getReprPartial( M::Union{<:typeof(SpecialOrthogonal(2)), <:Rotations{2}},
                          repr::AbstractMatrix{T},
                          partial::AbstractVector{Int}, # total partial from user over all Factors
                          offset::Base.RefValue{Int}=Ref(0),
                          mask::BitVector=_checkManifoldPartialDims(M,partial,offset, false) ) where {T <: Number}
  #
  @assert sum(mask) == 1 "Can only return the same point represenation matrix for SpecialOrthogonal(2) / Rotations(2)"
  return repr
end

## EXTRACT PARTIAL MANIFOLD

function getManifoldPartial(M::Union{<:Euclidean{Tuple{N}},<:TranslationGroup},
                            partial::AbstractVector{Int},
                            repr::_PartiableRepresentationFlat{T}=nothing,
                            offset::Base.RefValue{Int}=Ref(0);
                            doError::Bool=true) where {N, T <: Number}
  #
  mask = _checkManifoldPartialDims(M,partial,offset, doError)
  offset[] += manifold_dimension(M)
  len = sum(mask)
  repr_p = repr === nothing ? nothing : zeros(T,len)
  return (TranslationGroup(len),repr_p)
end

function getManifoldPartial(M::Circle,
                            partial::AbstractVector{Int},
                            repr::_PartiableRepresentation=nothing,
                            offset::Base.RefValue{Int}=Ref(0);
                            doError::Bool=true )
  #
  mask = _checkManifoldPartialDims(M,partial,offset,doError)
  offset[] += manifold_dimension(M)
  return (M,repr)
end

function getManifoldPartial(M::Rotations{2},
                            partial::AbstractVector{Int},
                            repr::_PartiableRepresentation=nothing,
                            offset::Base.RefValue{Int}=Ref(0);
                            doError::Bool=true )
  #
  mask = _checkManifoldPartialDims(M,partial,offset,doError)
  offset[] += manifold_dimension(M)
  return (M,repr)
end

function getManifoldPartial(M::TranslationGroup{Tuple{N}},
                            partial::AbstractVector{Int},
                            repr::_PartiableRepresentationFlat{T}=nothing,
                            offset::Base.RefValue{Int}=Ref(0);
                            doError::Bool=true ) where {N, T <: Number}
  #
  mask = _checkManifoldPartialDims(M,partial,offset,doError)
  offset[] += manifold_dimension(M)
  len = sum(mask)
  repr_p = repr === nothing ? nothing : zeros(T,len)
  return (TranslationGroup(len),repr_p)
end

function getManifoldPartial(M::typeof(SpecialOrthogonal(2)),
                            partial::AbstractVector{Int},
                            repr::_PartiableRepresentation=nothing,
                            offset::Base.RefValue{Int}=Ref(0);
                            doError::Bool=true )
  #
  mask = _checkManifoldPartialDims(M,partial,offset,doError)
  offset[] += manifold_dimension(M)
  return (M,repr)
end

"""
    $SIGNATURES

A so-called full dimension manifold can possibly be reduced to smaller partial manifolds over
some of the dimensions, returning a new programatically generated `<:AbstractManifold`.
This funciton can optionally also reduce a point representation for the desired
partial dimensions too.

Example
```julia
using Manifolds
using ApproxManifoldProducts

# a familiar manifold of translation and rotation in 2D
M = SpecialEuclidean(2)

# make a new partial of only the translation components
M_x, _ = getManifoldPartial(M,[1;])
# returns a new TranslationGroup(1) corresponding to just x dimension

# representation is semidirect product of translation and rotation matrix
u0 = ProductRepr([0.0;0],[1 0; 0 1.0])

# known coordinates are [x,y,θ], eg
#   vee(M,identity_element(M,u0),log(M,identity_element(M,u0),u0))
#   [0;0;0] in this example

# make another new partial of only the rotation component
M_rot, u_rot = getManifoldPartial(M,[3],u0)
# returns SpecialOrthogonal(2) information

# make another new partial of only  y and θ
M_yθ, u_yθ = getManifoldPartial(M,[2;3],u0)
# returns new manifold information as ProductArray(TranslationGroup(1),SpecialOrthogonal(2))
```

Notes
- Partial dimensions of interest are defined via a `AbstractVector{Int}` of coordinate dimensions.
- assumed to follow coordinates as transition step towards more general solution
- assume ProductManifold is the only way to stitch multiple manifolds together
- Still experimental.

DevNotes
- FIXME any semidirect product or action information is lost and naive Product manifold assumptions are currently made.

Related

[`getManifold`](@ref), [`AbstractManifold`], [`ProductManifold`], [`GroupManifold`], [`ProductRepr`]
"""
function getManifoldPartial(M::ProductManifold,
                            partial::AbstractVector{Int},
                            repr::_PartiableRepresentationProduct=nothing,
                            offset::Base.RefValue{Int}=Ref(0);
                            doError::Bool=true )
  #
  _checkManifoldPartialDims(M,partial,offset,doError)

  # loop through the ProductManifold components
  ManiArr = []
  ReprArr = []
  for (i,m) in enumerate(M.manifolds)
    mask = _checkManifoldPartialDims(m,partial,offset, false)
    if any(mask)
      Mp = if repr === nothing
        # decide if representation should also be updated or left as nothing
        Mp, = getManifoldPartial(m, partial, nothing, offset, doError=false)
        Mp
      else
        # hard assumption that repr::ProductRepr to go along with M::ProductManifold
        Mp, Rp = getManifoldPartial(m, partial, repr.parts[i], offset, doError=false)
        push!(ReprArr, Rp)
        Mp
      end
      push!(ManiArr, Mp)
    else
      offset[] += manifold_dimension(m)
    end
  end

  # trivial case, drop the ProductManifold for single element
  if length(ManiArr) == 1
    repr_p = repr === nothing ? nothing : ReprArr[1]
    return (ManiArr[1],repr_p)
  elseif 1 < length(ManiArr)
    repr_p = repr === nothing ? nothing : ProductRepr(ReprArr...)
    return (ProductManifold(ManiArr...), repr_p)
  end
  error("partial manifold calculations should not reach here")
end

function getManifoldPartial(M::GroupManifold,
                            partial::AbstractVector{<:Integer},
                            repr::_PartiableRepresentation=nothing,
                            offset::Base.RefValue{<:Integer}=Ref(0);
                            doError::Bool=true )
  #
  # mask the desired coordinate dimensions
  mask = _checkManifoldPartialDims(M,partial,offset, doError)

  if sum(mask) == manifold_dimension(M)
    # asking for all coordinate dimensions as offered by M
    return (M,repr)
  end
  # recursion may need to branch for ProductManifold
  # Note loss of the Group operation information at this time
  getManifoldPartial(M.manifold, partial, repr, offset, doError=doError)
end

#


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "SDL2_ttf"
version = v"2.0.15"

# Collection of sources required to build SDL2_ttf
sources = [
    ArchiveSource("https://www.libsdl.org/projects/SDL_ttf/release/SDL2_ttf-$(version).tar.gz",
                  "a9eceb1ad88c1f1545cd7bd28e7cbc0b2c14191d40238f531a15b01b1b22cd33"),
    DirectorySource("./bundled"),
]

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

atomic_patch -p1 ../patches/configure_in-v2.0.15.patch
atomic_patch -p1 ../patches/Makefile_in_dont_build_programs.patch
touch NEWS README AUTHORS ChangeLog
# For some reasons, the first time `autoreconf` may fail,
# but with some encouragement it can do it
autoreconf -vi || autoreconf -vi
export CPPFLAGS="-I${prefix}/include/SDL2"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --enable-shared \
    --disable-static \
    "${FLAGS[@]}"
make -j${nproc}
make install

if [[ "${target}" == *-mingw* ]]; then
    # As usual, build system for Windows is wrecked
    # and the shared library is not built at all
    cd ${prefix}/lib
    ar x libSDL2_ttf.a
    cc -shared -o ${libdir}/SDL2_ttf.dll SDL_ttf.o ${libdir}/libfreetype-6.dll ${libdir}/SDL2.dll
    rm SDL_ttf.o
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libSDL2_ttf", "SDL2_ttf"], :libsdl2_ttf)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("SDL2_jll"),
    Dependency("FreeType2_jll"),
    # The following libraries aren't needed for the build, but libSDL2_ttf is
    # dynamically linked to them regardless.
    Dependency("libpng_jll"),
    # TODO: Next time you build a new version of this package, make sure
    # HarfBuzz is actually needed, it doesn't look like it's the case.
    Dependency("HarfBuzz_jll"),
    Dependency("Graphite2_jll"),
    Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
    Dependency("PCRE_jll"),
]

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


text:
using Statistics
using Base.Threads: @threads

# prepare a dataset
X = rand(Int(1.25e6), 100)

function get_edges(X::AbstractMatrix{T}, nbins=250) where {T}
    edges = Vector{Vector{T}}(undef, size(X,2))
    @threads for i in 1:size(X, 2)
    # for i in 1:size(X, 2)
        edges[i] = quantile(view(X, :,i), (1:nbins)/nbins)
        if length(edges[i]) == 0
            edges[i] = [minimum(view(X, :,i))]
        end
    end
    return edges
end

println("num threads: ", Threads.nthreads())

println("trial 1: ")
edges = get_edges(X, 128);

println("trial 2: ")
edges = get_edges(X, 128);

println("trial 3: ")
edges = get_edges(X, 128);

println("trial 4: ")
edges = get_edges(X, 128);

println("trial 5: ")
edges = get_edges(X, 128);

println("trial 6: ")
edges = get_edges(X, 128);

println("trial 7: ")
edges = get_edges(X, 128);

println("trial 8: ")
edges = get_edges(X, 128);


text:
# This file contains methods that are dispatched to when using neural networks
# with a preallocated cache (called CachedNetworks). Those are used to greatly
# improve performance, especially in

export cached, vectorize_gradient, weights, grad_cache

"""
    NNCache{N}

The base abstract type holding a cache for the network type `N<:NeuralNetwork`.
"""
abstract type NNCache{N}
    #valid::Bool
end

"""
    invalidate_cache!(net)

Invalidates the cache of a cached Neural Network. Usually called when weights
are changed to signal the fact that lookup tables must be cleared.
"""
@inline invalidate_cache!(net::NNCache) = net.valid = false
@inline invalidate_cache!(net::Nothing) = nothing # if not implemented

"""
    CachedNet{N, NC}

A Cached Version of Neural Network `N` using Cache `NC`. It behaves the same as
the network `N`, but uses a preallocated cache to speed up computation of the
gradient and the function.
"""
struct CachedNet{N, NC} <: NeuralNetwork
    net::N
    cache::NC
end

"""
    cached(net)

Constructs a cached version of the neural network `net`, preallocating several
arrays for intermediate results for improved performance. This object will
behave identically to a standard network `net`.
"""
cached(net::NeuralNetwork) = CachedNet(net, cache(net))
cached(net::CachedNet) = throw("Can't build cached net of cached network!")

# When copying shallow copy the net but deepcopy the der_vec
"""
    copy(cnet::CachedNet)

Copy a cached network, building a shallow copy of the network and a deep-copy
of the cache.
"""
copy(cnet::CachedNet) = CachedNet(cnet.net, deepcopy(cnet.cache))

"""
    cache(net)

Constructs the `NNCache{typeof(net)}` object that holds the cache for this network.
If it has not been implemented returns nothing.
"""
cache(net) = nothing

"""
    weights(net)

Access the underlying structure holding the weights of the network. This is usefull
when dealing with Symmetrized networks or cached networks, as it gives access to
the wrapped data structure.
"""
weights(net) = net
weights(cnet::CachedNet) = cnet.net
grad_cache(net::CachedNet) = grad_cache(net.net)

@inline (cnet::CachedNet)(σ...) = logψ(cnet, σ...)
# When you call logψ on a cached net use the cache to compute the net
@inline logψ(cnet::CachedNet, σ...) = cnet.net(cnet.cache, config(σ)...)

#@inline ∇logψ(n::CachedNet, σ) = logψ_and_∇logψ(n, σ)[2]
function logψ_and_∇logψ(n::CachedNet, σ::Vararg{N,V}) where {N,V}
    #@warn "Inefficient calling logψ_and_∇logψ for cachedNet"
    ∇lnψ   = grad_cache(n)
    lψ, _g = logψ_and_∇logψ!(∇lnψ, n, σ...);
    return (lψ, ∇lnψ)
end

# Declare the two functions, even if config(blabla)=blabla, because of a shitty
# Julia's performance bug #32761
# see https://github.com/JuliaLang/julia/issues/32761
@inline logψ_and_∇logψ!(der, n::CachedNet, σ::State) = logψ_and_∇logψ!(der, n, config(σ))
@inline function logψ_and_∇logψ!(der, n::CachedNet, σ::NTuple{N,AbstractArray}) where N
    lψ = logψ_and_∇logψ!(der, n.net, n.cache, config(σ)...)
    return (lψ, der)
end
@inline function logψ_and_∇logψ!(der, n::CachedNet, σ::Vararg{AbstractArray,N}) where N
    lψ = logψ_and_∇logψ!(der, n.net, n.cache, σ...);
    return (lψ, der)
end

## Optimisation of cachednet
update!(opt, cnet::CachedNet, Δ, state=nothing) = (update!(opt, weights(cnet), weights(Δ), state); invalidate_cache!(cnet.cache))
apply!(opt, val1::Union{NeuralNetwork, CachedNet}, val2::Union{NeuralNetwork, CachedNet}, args...) =
    apply!(weights(val1), weights(val2), args...)

# The common operations are forwarded to the underlying network.
input_type(cnet::CachedNet) = input_type(cnet.net)
out_type(cnet::CachedNet) = out_type(cnet.net)
input_shape(net::CachedNet) = input_shape(cnet.net)
random_input_state(cnet::CachedNet) = random_input_state(cnet.net)
is_analytic(cnet::CachedNet) = is_analytic(cnet.net)
num_params(cnet::CachedNet) = num_params(cnet.net)

# Visualization
function Base.show(io::IO, m::CachedNet)
    print(io, "CachedNet{$(m.net)}")
end

Base.show(io::IO, ::MIME"text/plain", m::CachedNet) = print(io, "CachedNet{$(m.net)}")

## Particoular kinds of networks
# Thing for Matrix Neural Network
# A MatrixNeuralNetwork is a network for a state that has a bi-partite state,
# representing row and column of a density matrix in the same space.
abstract type MatrixNeuralNetwork <: NeuralNetwork end
const MatrixNet   = Union{MatrixNeuralNetwork, CachedNet{<:MatrixNeuralNetwork}}
#random_input_state(net) = sample_input_type(net).([rand(0:1) for i=1:input_length(model)])

# Thing for Pure Neural Network
# A MatrixNeuralNetwork is a network for a state that has a bi-partite state,
# representing row and column of a density matrix in the same space.
abstract type PureNeuralNetwork <: NeuralNetwork end
const PureNet   = Union{PureNeuralNetwork, CachedNet{<:PureNeuralNetwork}}
#random_input_state(net) = sample_input_type(net).([rand(0:1) for i=1:input_length(model)])

# Thing for Pure-state (Closed system) Neural Network
# A MatrixNeuralNetwork is a network for a state that has a bi-partite state,
# representing row and column of a density matrix in the same space.
abstract type KetNeuralNetwork <: NeuralNetwork end
const KetNet   = Union{KetNeuralNetwork, CachedNet{<:KetNeuralNetwork}}
#random_input_state(net) = sample_input_type(net).([rand(0:1) for i=1:input_length(model)])

# This overrides the standard behaviour of net(σ...) because Vector unpacking
# should not happen
@inline logψ(cnet::CachedNet{<:KetNeuralNetwork}, σ) = cnet.net(cnet.cache, config(σ))
function logψ_and_∇logψ(net::KetNeuralNetwork, σ)
    σ = config(σ)
    y, back = forward(net -> net(σ), net)

    # This computes the gradient, which is the conjugate of the derivative
    der = back(Int8(1))[1]

    # We take the conjugate of every element
    for key=keys(der)
        conj!(der[key])
    end
    der = add_vector_field(net, der)

    return y, der
end

function logψ_and_∇logψ(n::CachedNet{<:KetNeuralNetwork}, σ)
    ∇lnψ = grad_cache(n)
    lnψ = logψ_and_∇logψ!(∇lnψ, n.net, n.cache, config(σ));
    return (lnψ, ∇lnψ)
end


text:
"""
Stockmayer potential
T --- temperature [K]
n --- densty [mol/l]
del --- charge  []
sigma --- [nm]
eps --- epsilon/kB [K]
Mm -- molar mass  []

"""
function StECH(T::Number, n::Number, del::Number, sigma::Number, eps::Number, Mm::Number)
omdel11=[0, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5];
omtemp11=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3, 3.5,
4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 50, 75, 100 ]

oms11=[ 4.0079  4.002  4.655  5.521  6.454  8.214  9.824  11.31;
 3.13  3.164  3.355  3.721  4.198  5.23  6.225  7.16;
 2.6494  2.657  2.77  3.002  3.319  4.054  4.785  5.483;
 2.3144  2.32  2.402  2.572  2.812  3.386  3.972  4.539;
 2.0661  2.073  2.14  2.278  2.472  2.946  3.437  3.918;
 1.8767  1.885  1.944  2.06  2.225  2.628  3.054  3.474;
 1.7293  1.738  1.791  1.893  2.036  2.388  2.763  3.137;
 1.6122  1.622  1.67  1.76  1.886  2.198  2.535  2.872;
 1.5175  1.527  1.572  1.653  1.765  2.044  2.349  2.657;
 1.4398  1.45  1.49  1.564  1.665  1.917  2.196  2.478;
1.3204  1.33  1.364  1.425  1.509  1.72  1.956  2.199;
 1.2336  1.242  1.272  1.324  1.394  1.573  1.777  1.99;
 1.1679  1.176  1.202  1.246  1.306  1.461  1.639  1.827;
1.1166  1.124  1.146  1.185  1.237  1.372  1.53  1.698;
1.0753  1.082  1.102  1.135  1.181  1.3  1.441  1.592;
 1.0006  1.005  1.02  1.044  1.08  1.17  1.278  1.397;
0.95003  0.9538  0.9656  0.9852  1.012  1.082  1.168  1.265;
0.91311  0.9162  0.9256  0.9413  0.9626  1.019  1.09  1.17;
 0.88453  0.8871  0.8948  0.9076  0.9252  0.9721  1.031  1.098;
0.84277  0.8446  0.8501  0.8592  0.8716  0.9053  0.9483  0.9984;
0.81287  0.8142  0.8183  0.8251  0.8344  0.8598  0.8927  0.9316;
0.78976  0.7908  0.794  0.7993  0.8066  0.8265  0.8526  0.8836;
0.77111  0.772  0.7745  0.7788  0.7846  0.8007  0.8219  0.8474;
0.75553  0.7562  0.7584  0.7619  0.7667  0.78  0.7976  0.8189;
0.7422  0.7428  0.7446  0.7475  0.7515  0.7627  0.7776  0.7957;
0.72022  0.7206  0.722  0.7241  0.7271  0.7354  0.7464  0.76;
0.70254  0.7029  0.7039  0.7055  0.7078  0.7142  0.7228  0.7334;
0.68776  0.688  0.6888  0.6901  0.6919  0.697  0.704  0.7125;
0.6751  0.6753  0.676  0.677  0.6785  0.6827  0.6884  0.6955;
0.66405  0.6642  0.6648  0.6657  0.6669  0.6704  0.6752  0.6811;
0.64136  0.6415  0.6418  0.6425  0.6433  0.6457  0.649  0.6531;
0.6235  0.6236  0.6239  0.6243  0.6249  0.6267  0.6291  0.6321;
0.60882  0.6089  0.6091  0.6094  0.6099  0.6112  0.6131  0.6154;
0.5964  0.5964  0.5966  0.5969  0.5972  0.5983  0.5998  0.6017;
0.57626  0.5763  0.5764  0.5766  0.5768  0.5775  0.5785  0.5798;
0.54146  0.5415  0.5416  0.5416  0.5418  0.5421  0.5424  0.5429;
0.51803  0.5181  0.5182  0.5184  0.5184  0.5185  0.5186  0.5187]

omdel22=[0, 0.25, 0.5, 0.75, 1.0, 1.5, 2, 2.5]
omtemp22=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4,
5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 50, 75, 100 ]

oms22=[4.1005  4.266  4.833  5.742  6.729  8.624  10.34  11.89;
3.2626  3.305  3.516  3.914  4.433  5.57  6.637  7.618;
2.8399  2.836  2.936  3.168  3.511  4.329  5.126  5.874;
2.531  2.522  2.586  2.749  3.004  3.64  4.282  4.895;
2.2837  2.277  2.329  2.46  2.665  3.187  3.727  4.249;
2.0838  2.081  2.13  2.243  2.417  2.862  3.329  3.786;
1.922  1.924  1.97  2.072  2.225  2.614  3.028  3.435;
1.7902  1.795  1.84  1.934  2.07  2.417  2.788  3.156;
1.6823  1.689  1.733  1.82  1.944  2.258  2.596  2.933;
1.5929  1.601  1.644  1.725  1.838  2.124  2.435  2.746;
1.4551  1.465  1.504  1.574  1.67  1.913  2.181  2.451;
1.3551  1.365  1.4  1.461  1.544  1.754  1.989  2.228;
1.28  1.289  1.321  1.374  1.447  1.63  1.838  2.053;
1.2219  1.231  1.259  1.306  1.37  1.532  1.718  1.912;
1.1757  1.184  1.209  1.251  1.307  1.451  1.618  1.795;
1.0933  1.1  1.119  1.15  1.193  1.304  1.435  1.578;
1.0388  1.044  1.059  1.083  1.117  1.204  1.31  1.428;
0.99963  1.004  1.016  1.035  1.062  1.133  1.22  1.319;
0.96988  0.9732  0.983  0.9991  1.021  1.079  1.153  1.236;
0.92676  0.9291  0.936  0.9473  0.9628  1.005  1.058  1.121;
0.89616  0.8979  0.903  0.9114  0.923  0.9545  0.9955  1.044;
0.87272  0.8741  0.878  0.8845  0.8935  0.9181  0.9505  0.9893;
0.85379  0.8549  0.858  0.8632  0.8703  0.8901  0.9164  0.9482;
0.83795  0.8388  0.8414  0.8456  0.8515  0.8678  0.8895  0.916;
0.82435  0.8251  0.8273  0.8308  0.8356  0.8493  0.8676  0.8901;
0.80184  0.8024  0.8039  0.8065  0.8101  0.8201  0.8337  0.8504;
0.78363  0.784  0.7852  0.7872  0.7899  0.7976  0.8081  0.8212;
0.76834  0.7687  0.7696  0.7712  0.7733  0.7794  0.7878  0.7983;
0.75518  0.7554  0.7562  0.7575  0.7592  0.7642  0.7711  0.7797;
0.74364  0.7438  0.7445  0.7455  0.747  0.7512  0.7569  0.7642;
0.71982  0.72  0.7204  0.1211  0.7221  0.725  0.7289  0.7339;
0.70097  0.7011  0.7014  0.7019  0.7026  0.7047  0.7076  0.7112;
0.68545  0.6855  0.6858  0.6861  0.6867  0.6883  0.6905  0.6932;
0.67232  0.6724  0.6726  0.6728  0.6733  0.6745  0.6762  0.6784;
0.65099  0.651  0.6512  0.6513  0.6516  0.6524  0.6534  0.6546;
0.61397  0.6141  0.6143  0.6145  0.6147  0.6148  0.6148  0.6147;
0.5887  0.5889  0.5894  0.59  0.5903  0.5901  0.5895  0.5885]

spline11 = Spline2D(omtemp11, omdel11, oms11)
spline22 = Spline2D(omtemp22, omdel22, oms22)

sigw::Float64=sigma*10^(-9)	#to [m]

#condstants
Na::Float64 = 6.02 * 10.0^23
kb::Float64 = 1.38*10.0^-23;

#dimensionless units
tb::Float64 = T/eps
nb::Float64 = n*Na*1000.0*sigw^3

#integrals
omega11::Float64 = evaluate(spline11, tb, del)
omega22::Float64 = evaluate(spline22, tb, del)

#println("$(omega11) Na $(Na) $(pi) sigma $(sigma) tb $(tb) nb $(nb)")

D0z::Float64=3.0/(8.0*sqrt(pi)*nb)*sqrt(tb)/omega11
#println("D0 $(D0)")
D0::Float64=D0z*sigw/sqrt(Mm/1000/Na/(eps*kb))*10^9
#println("D $(D)")

nu0z::Float64=5.0/(16.0*sqrt(pi))*sqrt(tb)/omega22
#out.nu0=nu0*sqrt(mw/Na*(epsw*kb))/sigw^2*10^6;
nu0::Float64=nu0z*sqrt(Mm/1000/Na*(eps*kb))/sigw^2*10^6

la0z=75/(64*sqrt(pi))*sqrt(tb)/omega22
la0=la0z*kb/sigw^2/sqrt(Mm/1000/Na/(eps*kb))

return [D0z, D0, nu0z, nu0, la0z, la0, omega11, omega22]
end


text:
# Buggy and incomplete LISP-ish interpreter based mostly on the
# SICP Chapter 4 metacircular
# evaluator with a few changes to better suit Julia implementation.
#
# I also added a simple parser since the host language is not S-expression
# based. The parser generates arrays, which can cause some bugs... maybe
# at some point I'll make the interpreter properly cons-list based
# but I have to go shower now.

module Debug
  export debug
  function debug(msg)
    print_with_color(:blue, "DEBUG>> $(msg)\n")
  end
end

module Parser
  # Simple parser based on Peter Novig's "Lispy" parser.
  # Takes S-expression input
  #
  #   (+ (- 3 1) 4)
  #
  # And returns an AST array:
  #
  #   [:+, [:-, 3, 1], 4]
  #

  using Debug
  export parse

  function parse(program)
    debug("parse($(program))")
    program |> tokenize |> to_ast
  end

  # PRIVATE

  function tokenize(program)
    debug("tokenize($(program))")
    program                   |>
    x->replace(x, "(", " ( ") |>
    x->replace(x, ")", " ) ") |>
    x->replace(x, "\n", " ")  |>
    split
  end

  function to_ast(tokens::Array)
    debug("to_ast($(tokens))")

    if length(tokens) == 0
      error("SYNTAXERR: empty program given")
    end

    token = shift!(tokens)

    if "(" == token
      ast = []

      while tokens[1] != ")"
        push!(ast, to_ast(tokens))
      end

      shift!(tokens) # pop off ")"
      debug("... to_ast() -> $(ast)")
      return ast
    elseif ")" == token
      error("SYNTAXERR unexpected `)`")
    else
      atom(token)
    end
  end

  # Numbers become numbers. Everything else is a symbol.
  function atom(token)
    debug("atom($(token))")
    token                                                                |>
    x -> tryparse(Int, x)                                                |>
    maybeint -> (isnull(maybeint) ? tryparse(Float64, token) : maybeint) |>
    maybenum -> isnull(maybenum) ? Symbol(token) : maybenum.value
  end
end

module Environments
  using Debug
  export Environment, AbstractEnvironment, EmptyEnvironment
  export theglobalenvironment,
         setvariable!,
         definevariable!,
         lookup_variable_value

  # CONSTS

  const PRIMITIVE_PROCEDURES =
    [(:+, +),
     (:-, -),
     (:*, *),
     (:car, x->x[1]),
     (:cdr, x->x[2]),
     (:cons, (x,y)->(debug("push($(x), $(y))"); push!(Array{Any}([x]),y)))] |>
     tups -> map(tup -> (tup[1], [:primitive, tup[2]]), tups) |> # tag procs
     Dict{Symbol, Any}

  # TYPES

  abstract AbstractEnvironment
  type EmptyEnvironment <: AbstractEnvironment end

  type Environment <: AbstractEnvironment
    values::Dict
    outer::AbstractEnvironment
  end

  # FUNCTIONS

  function theglobalenvironment()
    debug("theglobalenvironment()")
    setupenvironment()
  end

  function definevariable!(env::Environment, var::Symbol, val)
    debug("definevariable!(env, $(var), $(val))")
    env.values[var] = val
  end

  function setvariable!(env::Environment, var::Symbol, val)
    debug("setvariable!(env, $(var), $(val))")
    if haskey(env.values, var)
      env.values[var] = val
    else
      setvariable!(env.outer, var, val)
    end
  end

  function setvariable!(env::EmptyEnvironment, var, val)
    debug("setvariable!(env, $(var), $(val))")
    error("ERROR: unbound variable $(var), could not set $(val)")
  end

  function lookup_variable_value(env::Environment, var::Symbol)
    debug("lookup_variable_value($(env), $(var))")
    if haskey(env.values, var)
      env.values[var]
    else
      lookup_variable_value(env.outer, var)
    end
  end

  function lookup_variable_value(env::EmptyEnvironment, var)
    debug("lookup_variable_value(env::EmptyEnvironment, $(var))")
    error("ERROR: unbound variable $(var)")
  end

  # PRIVATE

  function setupenvironment()
    debug("setupenvironment()")
    initialenvironment = Environment(PRIMITIVE_PROCEDURES, EmptyEnvironment())
    definevariable!(initialenvironment, :t, true)
    definevariable!(initialenvironment, :f, false)
    initialenvironment
  end
end

module Interpreter
  using Debug
  using Environments
  import Parser
  export eval

  # FUNCTIONS

  function eval(expression)
    eval(expression, theglobalenvironment())
  end

  function eval(expression::String, environment)
    Parser.parse(expression) |>
    parsed -> eval(parsed, environment)
  end

  function eval(expression::Array, environment)
    analyze(expression)(environment)
  end

  # PRIVATE

  function analyze(exp)
    debug("analyze($(exp))")
    if is_self_evaluating(exp)
      debug("analyze_self_evaluating($(exp))")
      env -> exp
    elseif is_variable(exp)
      debug("analyze_variable($(exp))")
      env -> lookup_variable_value(env, exp)
    # '(quote text)
    elseif is_tagged_list(exp, :quote)
      debug("analyze_quote($(exp))")
      env -> exp[2]                      # TODO: should be a cons list
    elseif is_tagged_list(exp, :set!)
      analyze_assignment(exp)
    elseif is_tagged_list(exp, :define)
      analyze_definition(exp)
    elseif is_tagged_list(exp, :if)
      analyze_if(exp)
    elseif is_tagged_list(exp, :lambda)
      analyze_lambda(exp)
    elseif is_tagged_list(exp, :begin)
      analyze_sequence(exp[2:end])
    elseif isa(exp, Array)
      analyze_application(exp)
    else
      error("SYNTAXERR unknown expression type $(exp)")
    end
  end

  function is_self_evaluating(exp)
    debug("is_self_evaluating($(exp))")
    (isa(exp, Number) || isa(exp, String)) ? true : false
  end

  function is_tagged_list(exp, tag::Symbol)
    debug("is_tagged_list($(exp), $(tag))")
    exp[1] == tag
  end

  function is_variable(exp)
    debug("is_variable($(exp))")
    isa(exp, Symbol)
  end

  function analyze_sequence(exps)
    debug("analyze_sequence($(exps))")
    sequentially = (proc1, proc2) -> (env) -> (proc1(env); proc2(env))
    procs = map(analyze, exps)

    function loop(procs)
      debug("... loop($(procs))")
      if length(procs) == 1
        procs[1]
      else
        loop(sequentially(procs[1], procs[2]), procs[3:end])
      end
    end

    if length(procs) == 0
      error("ANALYZEERR empty sequence")
    else
      loop(procs)
    end
  end

  # (set! varname value)
  function analyze_assignment(exp)
    debug("analyze_assignment($(exp))")
    varname = exp[2]
    # If the value is also a compound expression, return lambda that will
    # evaluate it.
    valueproc = analyze(exp[3])
    env -> (setvariable!(env, varname, valueproc(env)); :ok)
  end

  function analyze_definition(exp)
    debug("analyze_definition($(exp))")
    defname = definitionvariable(exp)
    # Similar to analyze_assignment
    valueproc = definitionvalue(exp) |> analyze
    env -> (definevariable!(env, defname, valueproc(env)); :ok)
  end

  # (if predicate consequent alternative)
  function analyze_if(exp)
    debug("analyze_if($(exp))")
    predicateproc = analyze(exp[2])
    consequentproc = analyze(exp[3])
    alternativeproc = analyze(exp[4])
    env -> predicateproc(env) ? consequentproc(env) : alternativeproc(env)
  end

  # (lambda parameters body)
  function analyze_lambda(exp)
    debug("analyze_lambda($(exp))")
    args = exp[2]
    bodyproc = analyze_sequence(exp[3:end])
    # PROC
    env -> [:procedure, args, bodyproc, env]
  end

  # (operator operands)
  function analyze_application(exp)
    debug("analyze_application($(exp))")
    fproc = analyze(exp[1])
    aprocs = map(analyze, exp[2:end])
    env -> execute_application(fproc(env), map(aproc->aproc(env), aprocs))
  end

  function execute_application(proc, args)
    debug("execute_application($(proc), $(args))")
    if is_tagged_list(proc, :primitive)
      debug("... $(proc[1])($(args))")
      proc[2](args...) # proc = [:primitive, lambda]
    elseif is_tagged_list(proc, :procedure)
      # See PROC
      args_dict = Dict(zip(proc[2], args))
      new_env = Environment(args_dict, proc[4]) # Extend env and pass to proc.
      proc[3](new_env)
    else
      error("ERR unknown procedure $(proc)")
    end
  end

  function definitionvariable(exp)
    debug("definitionvariable($(exp))")
    if isa(exp[2], Symbol)
      exp[2]               # (define VAR value)
    else
      exp[2][1]            # (define (VAR args) lambda)
    end
  end

  function definitionvalue(exp)
    debug("definitionvalue($(exp))")
    if isa(exp[2], Symbol)
      exp[3]
    else
      makelambda(exp[2][2:end], exp[3])
    end
  end

  function makelambda(args, body)
    debug("makelambda($(args), $(body))")
    [:lambda, args, body]
  end
end

module REPL
  using Debug
  using Interpreter
  using Environments
  using Parser
  export run

  # CONSTS

  const  INPUT_PROMPT = "--- INPUT ---"
  const  OUTPUT_PROMPT = "--- RESULT ---"

  # FNUCTIONS

  function run()
    printforio(INPUT_PROMPT)
    input = readline(STDIN) |> Parser.parse
    output = Interpreter.eval(input, theglobalenvironment())
    printforio(OUTPUT_PROMPT)
    printresult(output)
    run()
  end

  # PRIVATE

  function printforio(text)
    print_with_color(:yellow, "\n$(text)\n")
  end

  function printresult(text)
    print_with_color(:red, "$(text)")
  end
end

REPL.run()


text:
function inner_rows!(C,A,B)
    for i in 1:100, j in 1:100
        C[i,j] = A[i,j] + B[i,j]
    end
end

function inner_cols!(C,A,B)
  for j in 1:100, i in 1:100
    C[i,j] = A[i,j] + B[i,j]
  end
end

max_stars_repo_name: jvaverka/BookTest

max_issues_repo_name: jvaverka/BookTest

max_forks_repo_name: jvaverka/BookTest


text:
#================================================================

MonteCarloHash for hashing function spaces.

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

#========================
Global constants
========================#

const _MONTECARLOHASH_DEFAULT_N_SAMPLES = 1024

#========================
Typedefs
========================#

struct MonteCarloHash{F, H <: LSHFunction, D, T, S} <: LSHFunction
    similarity :: F

    discrete_hashfn :: H
    μ :: D

    # TODO: make typeof(volume) and typeof(p) match the element type
    # returned by the sampler μ.

    # The volume of the L^p space we're embedding, defined as
    #
    #   volume = ∫_Ω dx
    #
    volume :: T

    # The order of the L^p space that we're embedding
    p :: T

    n_samples :: Int64
    sample_points :: Vector{S}

    ### Internal constructors
    function MonteCarloHash(similarity::F, discrete_hashfn::H, μ, volume,
                            p, n_samples) where {F, H<:LSHFunction}
        sample_points = [μ() for ii = 1:n_samples]

        T = eltype(μ())
        volume = T(volume)
        p = T(p)

        new{F,H,typeof(μ),T,eltype(sample_points)}(
            similarity, discrete_hashfn, μ, volume,
            p, n_samples, sample_points
        )
    end
end

### External MonteCarloHash constructors
const _valid_MonteCarloHash_similarities = (
    # Function space similarities
    (L1, L2, cossim),
    # Discrete-space similarities corresponding to function space similarities
    (ℓ1, ℓ2, cossim),
    # Order of L^p space that the similarity applies to
    (1, 2, 2),
)

# TODO: restrict similarities. E.g. Jaccard should not be an available similarity
@doc """
    MonteCarloHash(sim, ω, args...; volume=1.0, n_samples=$(_MONTECARLOHASH_DEFAULT_N_SAMPLES), kws...)

Samples a hash function from an LSH family for the similarity `sim` defined over the function space ``L^p_{\\mu}(\\Omega)``. `sim` may be one of the following:
$(
join(
    ["- `" * sim * "`" for sim in (_valid_MonteCarloHash_similarities[1] .|>
                                   string |>
                                   collect |>
                                   sort!)
    ],
    "\n"
)
)

Given an input function ``f\\in L^p_{\\mu}(\\Omega)``, `MonteCarloHash` works by sampling ``f`` at some randomly-selected points in ``\\Omega``, and then hashing those samples.

# Arguments
- `sim`: the similarity statistic you want to hash on.
- `ω`: a function that takes no inputs and samples a single point from ``\\Omega``. Alternatively, it can be viewed as a random variable with probability measure

```math
\\frac{\\mu}{\\text{vol}_{\\mu}(\\Omega)} = \\frac{\\mu}{\\int_{\\Omega} d\\mu}
```

- `args...`: arguments to pass on when building the `LSHFunction` instance underlying the returned `MonteCarloHash` struct.
- `volume::Real` (default: `1.0`): the volume of the space ``\\Omega``, defined as

```math
\\text{vol}_{\\mu}(\\Omega) = \\int_{\\Omega} d\\mu
```

- `n_samples::Integer` (default: `$(_MONTECARLOHASH_DEFAULT_N_SAMPLES)`): the number of points to sample from each function that is hashed by the `MonteCarloHash`. Larger values of `n_samples` tend to capture the input function better and will thus be more likely to achieve desirable collision probabilities.
- `kws...`: keyword arguments to pass on when building the `LSHFunction` instance underlying the returned `MonteCarloHash` struct.

# Examples
Create a hash function for cosine similarity for functions in ``L^2([-1,1])``:

```jldoctest; setup = :(using LSHFunctions)
julia> μ() = 2*rand()-1;   # μ samples a random point from [-1,1]

julia> hashfn = MonteCarloHash(cossim, μ, 50; volume=2.0);

julia> n_hashes(hashfn)
50

julia> similarity(hashfn) == cossim
true

julia> hashtype(hashfn)
$(cossim |> LSHFunction |> hashtype)
```

Create a hash function for ``L^2`` distance in the function space ``L^2([0,2\\pi])``. Hash the functions `f(x) = cos(x)` and `f(x) = x/(2π)` using the returned `MonteCarloHash`.

```jldoctest; setup = :(using LSHFunctions, Random; Random.seed!(0))
julia> μ() = 2π * rand(); # μ samples a random point from [0,2π]

julia> hashfn = MonteCarloHash(L2, μ, 3; volume=2π);

julia> hashfn(cos)
3-element Array{Int32,1}:
 -1
  3
  0

julia> hashfn(x -> x/(2π))
3-element Array{Int32,1}:
 -1
 -2
 -1
```

Create a hash function with a different number of sample points.

```jldoctest; setup = :(using LSHFunctions; μ() = rand())
julia> μ() = rand();  # Samples a random point from [0,1]

julia> hashfn = MonteCarloHash(cossim, μ; volume=1.0, n_samples=512);

julia> length(hashfn.sample_points)
512
```
"""
MonteCarloHash(similarity, args...; kws...) =
    MonteCarloHash(SimilarityFunction(similarity), args...; kws...)

for (fn_space_simfn, simfn, p) in zip(_valid_MonteCarloHash_similarities...)
    quote
        # Add dispatch for case in which we specify the similarity function
        # to be $fn_space_simfn
        function MonteCarloHash(
                    sim::SimilarityFunction{$fn_space_simfn},
                    μ,
                    args...;
                    n_samples::Integer=_MONTECARLOHASH_DEFAULT_N_SAMPLES,
                    volume=1.0,
                    kws...)

            discrete_hashfn = LSHFunction($simfn, args...; kws...)
            MonteCarloHash($fn_space_simfn, discrete_hashfn, μ, volume,
                           $p, n_samples)
        end
    end |> eval
end

# Implementation of MonteCarloHash for invalid similarity functions. Just throws
# an error.
# Necessary because otherwise the first MonteCarloHash constructor will go into
# infinite recursion if it receives an invalid similarity function.
function MonteCarloHash(sim::SimilarityFunction, args...; kws...)
    valid_sims = _valid_MonteCarloHash_similarities[1] .|>
                 string  |>
                 collect |>
                 sort!
    valid_sims = join(valid_sims, ", ", " or ")
    ErrorException("similarity must be $(valid_sims)") |> throw
end

#========================
MonteCarloHash helper functions
========================#

# Embed an input function f in the discrete space
function get_samples(hashfn::MonteCarloHash, f)
    α = (hashfn.volume / hashfn.n_samples)^(1/hashfn.p)
    α * f.(hashfn.sample_points)
end

#========================
LSHFunction API compliance
========================#

hashtype(hashfn::MonteCarloHash) =
    hashtype(hashfn.discrete_hashfn)

similarity(hashfn::MonteCarloHash) = hashfn.similarity

n_hashes(hashfn::MonteCarloHash) =
    n_hashes(hashfn.discrete_hashfn)

# TODO: this may not be true
collision_probability(hashfn::MonteCarloHash, args...; kws...) =
    collision_probability(hashfn.discrete_hashfn, args...; kws...)

#========================
SymmetricLSHFunction API compliance
========================#

(hashfn::MonteCarloHash)(f) = index_hash(hashfn, f)

#========================
AsymmetricLSHFunction API compliance
========================#

function index_hash(hashfn::MonteCarloHash, f)
    samples = get_samples(hashfn, f)
    index_hash(hashfn.discrete_hashfn, samples)
end

function query_hash(hashfn::MonteCarloHash, f)
    samples = get_samples(hashfn, f)
    query_hash(hashfn.discrete_hashfn, samples)
end

#========================
MonteCarloHash API
========================#

# Compute the similarity between f and g in the embedded space
function embedded_similarity(hashfn::MonteCarloHash, f, g)
    samples_f = get_samples(hashfn, f)
    samples_g = get_samples(hashfn, g)
    simfun = similarity(hashfn)
    simfun(samples_f, samples_g)
end


text:
import RobotDynamics: jacobian!

""""
Specifies whether the constraint is an equality or inequality constraint.
Valid subtypes are `Equality`, `Inequality` ⟺ `NegativeOrthant`, and `SecondOrderCone`.

The sense of a constraint can be queried using `sense(::AbstractConstraint)`

If `sense(con) <: Conic` (i.e. not `Equality`), then the following operations are supported:
* `Base.in(::Conic, x::StaticVector)`. i.e. `x ∈ cone`
* `projection(::Conic, x::StaticVector)`
* `∇projection(::Conic, J, x::StaticVector)`
* `∇²projection(::Conic, J, x::StaticVector, b::StaticVector)`
"""
abstract type ConstraintSense end
abstract type Conic <: ConstraintSense end

"""
Equality constraints of the form ``g(x) = 0`.
Type singleton, so it is created with `Equality()`.
"""
struct Equality <: ConstraintSense end
"""
Inequality constraints of the form ``h(x) \\leq 0``.
Type singleton, so it is created with `Inequality()`.
Equivalent to `NegativeOrthant`.
"""
struct NegativeOrthant <: Conic end
const Inequality = NegativeOrthant

struct PositiveOrthant <: Conic end

"""
The second-order cone is defined as
``\\|x\\| \\leq t``
where ``x`` and ``t`` are both part of the cone.
TrajectoryOptimization assumes the scalar part ``t`` is
the last element in the vector.
"""
struct SecondOrderCone <: Conic end

dualcone(::NegativeOrthant) = NegativeOrthant()
dualcone(::PositiveOrthant) = PositiveOrthant()
dualcone(::SecondOrderCone) = SecondOrderCone()

projection(::NegativeOrthant, x) = min.(0, x)
projection(::PositiveOrthant, x) = max.(0, x)

function projection(::SecondOrderCone, x::StaticVector)
    # assumes x is stacked [v; s] such that ||v||₂ ≤ s
    s = x[end]
    v = pop(x)
    a = norm(v)
    if a <= -s          # below the cone
        return zero(x)
    elseif a <= s       # in the cone
        return x
    elseif a >= abs(s)  # outside the cone
        return 0.5 * (1 + s/a) * push(v, a)
    else
        throw(ErrorException("Invalid second-order cone projection"))
    end
end

function ∇projection!(::SecondOrderCone, J, x::StaticVector{n}) where n
    s = x[end]
    v = pop(x)
    a = norm(v)
    if a <= -s                               # below cone
        J .*= 0
    elseif a <= s                            # in cone
        J .*= 0
        for i = 1:n
            J[i,i] = 1.0
        end
    elseif a >= abs(s)                       # outside cone
        # scalar
        b = 0.5 * (1 + s/a)
        dbdv = -0.5*s/a^3 * v
        dbds = 0.5 / a

        # dvdv = dbdv * v' + b * oneunit(SMatrix{n-1,n-1,T})
        for i = 1:n-1, j = 1:n-1
            J[i,j] = dbdv[i] * v[j]
            if i == j
                J[i,j] += b
            end
        end

        # dvds
        J[1:n-1,n] .= dbds * v

        # ds
        dsdv = dbdv * a + b * v / a
        dsds = dbds * a
        ds = push(dsdv, dsds)
        J[n,:] .= ds
    else
        throw(ErrorException("Invalid second-order cone projection"))
    end
    return J
end

function Base.in(x, ::SecondOrderCone)
    s = x[end]
    v = pop(x)
    a = norm(v)
    return a <= s
end

function cone_status(::SecondOrderCone, x)
    s = x[end]
    v = pop(x)
    a = norm(v)
    if a <= -s
        return :below
    elseif a <= s
        return :in
    elseif a > abs(s)
        return :outside
    else
        return :invalid
    end
end

function ∇²projection!(::SecondOrderCone, hess, x::StaticVector, b::StaticVector)
    n = length(x)
    s = x[end]
    v = pop(x)
    bs = b[end]
    bv = pop(b)
    a =  norm(v)

    if a <= -s
        return hess .= 0
    elseif a <= s
        return hess .= 0
    elseif a > abs(s)
        dvdv = -s/norm(v)^2/norm(v)*(I - (v*v')/(v'v))*bv*v' +
            s/norm(v)*((v*(v'bv))/(v'v)^2 * 2v' - (I*(v'bv) + v*bv')/(v'v)) +
            bs/norm(v)*(I - (v*v')/(v'v))
        dvds = 1/norm(v)*(I - (v*v')/(v'v))*bv;
        dsdv = bv'/norm(v) - v'bv/norm(v)^3*v'
        dsds = 0
        hess[1:n-1,1:n-1] .= dvdv*0.5
        hess[1:n-1,n] .= dvds*0.5
        hess[n:n,1:n-1] .= dsdv*0.5
        hess[n,n] = 0
        return hess
        # return 0.5*[dvdv dvds; dsdv dsds]
    else
        throw(ErrorException("Invalid second-order cone projection"))
    end
end

function ∇projection!(::NegativeOrthant, J, x::StaticVector{n}) where n
    for i = 1:n
        J[i,i] = x <= 0 ? 1 : 0
    end
end

function ∇²projection!(::NegativeOrthant, hess, x::StaticVector, b::StaticVector)
    hess .= 0
end

Base.in(::NegativeOrthant, x::StaticVector) = all(x->x<=0, x)

"""
    AbstractConstraint

Abstract vector-valued constraint for a trajectory optimization problem.
May be either inequality or equality (specified by `sense(::AbstractConstraint)::ConstraintSense`),
and be function of single or adjacent knotpoints.

Interface:
Any constraint type must implement the following interface:
```julia
n = state_dim(::MyCon)
m = control_dim(::MyCon)
p = Base.length(::MyCon)
sense(::MyCon)::ConstraintSense
c = evaluate(::MyCon, args...)
jacobian!(∇c, ::MyCon, args...)
```

All constraints are categorized into the following type tree:
```text
                        AbstractConstraint
                        ↙                ↘
           StageConstraint               CoupledConstraint
            ↙        ↘                       ↙           ↘
StageConstraint ControlConstraint CoupledStateConstraint CoupledControlConstraint
```

The state and control dimensions (where applicable) can be queried using
`state_dim(::AbstractConstraint)` and `control_dim(::AbstractConstraint)`.
The dimensions of a constraint can be verified using [`check_dims`](@ref).
The width of the constraint Jacobian is given by [`get_inds`](@ref) or [`widths`](@ref).

The number of constraint values associated with the constraint (length of the constraint vector)
is given with `length(::AbstractConstraint)`.

# Evaluation methods
Refer to the doc strings for the following methods for more information on the required
signatures.
* [`evaluate`](@ref)
* [`jacobian!`](@ref)
* [`∇jacobian!`](@ref)
"""
abstract type AbstractConstraint end

"Only a function of states and controls at a single knotpoint"
abstract type StageConstraint <: AbstractConstraint end
"Only a function of states at a single knotpoint"
abstract type StateConstraint <: StageConstraint end
"Only a function of controls at a single knotpoint"
abstract type ControlConstraint <: StageConstraint end
"Only a function of states and controls at two adjacent knotpoints"
abstract type CoupledConstraint <: AbstractConstraint end
"Only a function of states at adjacent knotpoints"
abstract type CoupledStateConstraint <: CoupledConstraint end
"Only a function of controls at adjacent knotpoints"
abstract type CoupledControlConstraint <: CoupledConstraint end

const StateConstraints =
    Union{StageConstraint,StateConstraint,CoupledConstraint,CoupledStateConstraint}
const ControlConstraints =
    Union{StageConstraint,ControlConstraint,CoupledConstraint,CoupledControlConstraint}

"Get constraint sense (Inequality vs Equality)"
sense(::C) where {C<:AbstractConstraint} = throw(NotImplemented(:sense, Symbol(C)))

"Dimension of the state vector"
RobotDynamics.state_dim(::C) where {C<:StateConstraint} =
    throw(NotImplemented(:state_dim, Symbol(C)))

"Dimension of the control vector"
RobotDynamics.control_dim(::C) where {C<:ControlConstraint} =
    throw(NotImplemented(:control_dim, Symbol(C)))

"Return the constraint value"
evaluate(::C) where {C<:AbstractConstraint} = throw(NotImplemented(:evaluate, Symbol(C)))

"Length of constraint vector"
Base.length(::C) where {C<:AbstractConstraint} = throw(NotImplemented(:length, Symbol(C)))

# widths(con::StageConstraint, n=state_dim(con), m=control_dim(con)) = (n+m,)
# widths(con::StateConstraint, n=state_dim(con), m=0) = (n,)
# widths(con::ControlConstraint, n=0, m=control_dim(con)) = (m,)
# widths(con::CoupledConstraint, n=state_dim(con), m=control_dim(con)) = (n+m, n+m)
# widths(con::CoupledStateConstraint, n=state_dim(con), m=0) = (n,n)
# widths(con::CoupledControlConstraint, n=0, m=control_dim(con)) = (m,m)

"Upper bound of the constraint, as a vector, which is 0 for all constraints
(except bound constraints)"
@inline upper_bound(con::AbstractConstraint) =
    upper_bound(sense(con)) * @SVector ones(length(con))
@inline upper_bound(::Inequality) = 0.0
@inline upper_bound(::Equality) = 0.0

"Upper bound of the constraint, as a vector, which is 0 equality and -Inf for inequality
(except bound constraints)"
@inline lower_bound(con::AbstractConstraint) =
    lower_bound(sense(con)) * @SVector ones(length(con))
@inline lower_bound(::Inequality) = -Inf
@inline lower_bound(::Equality) = 0.0

"""
    primal_bounds!(zL, zU, con::AbstractConstraint)

Set the lower `zL` and upper `zU` bounds on the primal variables imposed by the constraint
`con`. Return whether or not the vectors `zL` or `zU` could be modified by `con`
(i.e. if the constraint `con` is a bound constraint).
"""
primal_bounds!(zL, zU, con::AbstractConstraint) = false

"Is the constraint a bound constraint or not"
@inline is_bound(con::AbstractConstraint) = false

"Check whether the constraint is consistent with the specified state and control dimensions"
@inline check_dims(con::StateConstraint, n, m) = state_dim(con) == n
@inline check_dims(con::ControlConstraint, n, m) = control_dim(con) == m
@inline check_dims(con::AbstractConstraint, n, m) =
    state_dim(con) == n && control_dim(con) == m

get_dims(con::Union{StateConstraint,CoupledStateConstraint}, nm::Int) =
    state_dim(con), nm - state_dim(con)
get_dims(con::Union{ControlConstraint,CoupledControlConstraint}, nm::Int) =
    nm - control_dim(con), control_dim(con)
get_dims(con::AbstractConstraint, nm::Int) = state_dim(con), control_dim(con)

con_label(::AbstractConstraint, i::Int) = "index $i"

"""
    get_inds(con::AbstractConstraint)

Get the indices of the joint state-control vector that are used to calculate the constraint.
If the constraint depends on more than one time step, the indices start from the beginning
of the first one.
"""
get_inds(con::StateConstraint, n, m) = (1:n,)
get_inds(con::ControlConstraint, n, m) = (n .+ (1:m),)
get_inds(con::StageConstraint, n, m) = (1:n+m,)
get_inds(con::CoupledConstraint, n, m) = (1:n+m, n+m+1:2n+2m)

"""
    widths(::AbstractConstraint)
    widths(::AbstractConstraint, n, m)

Return a tuple of the widths of the Jacobians for a constraint. If `n` and `m` are not passed
in, they are assumed to be consistent with those returned by `state_dim` and `control_dim`.
"""
@inline widths(con::AbstractConstraint, n, m) = length.(get_inds(con, n, m))
@inline widths(con::StageConstraint) = (state_dim(con) + control_dim(con),)
@inline widths(con::StateConstraint) = (state_dim(con),)
@inline widths(con::ControlConstraint) = (control_dim(con),)
@inline widths(con::CoupledConstraint) =
    (state_dim(con) + control_dim(con), state_dim(con) + control_dim(con))
@inline widths(con::CoupledStateConstraint) = (state_dim(con), state_dim(con))
@inline widths(con::CoupledControlConstraint) = (control_dim(con), control_dim(con))

############################################################################################
# 								EVALUATION METHODS 										   #
############################################################################################
"""
    evaluate!(vals, con::AbstractConstraint, Z, [inds])

Evaluate constraints for entire trajectory. This is the most general method used to evaluate
constraints along the trajectory `Z`, and should be the one used in other functions.
The `inds` argument determines at which knot points the constraint is evaluated.

If `con` is a `StageConstraint`, this will call `evaluate(con, z)` by default, or
`evaluate(con, z1, z2)` if `con` is a `CoupledConstraint`.
"""
function evaluate!(
    vals::Vector{<:AbstractVector},
    con::StageConstraint,
    Z::AbstractTrajectory,
    inds = 1:length(Z),
)
    for (i, k) in enumerate(inds)
        vals[i] .= evaluate(con, Z[k])
    end
end

function evaluate!(
    vals::Vector{<:AbstractVector},
    con::CoupledConstraint,
    Z::AbstractTrajectory,
    inds = 1:length(Z)-1,
)
    for (i, k) in enumerate(inds)
        vals[i] .= evaluate(con, Z[k], Z[k+1])
    end
end

"""
    jacobian!(∇c, con::AbstractConstraint, Z, [inds, is_const])

Evaluate constraints for entire trajectory. This is the most general method used to evaluate
constraints along the trajectory `Z`, and should be the one used in other functions.
The `inds` argument determines at which knot points the constraint is evaluated.
The optional `is_const` argument is a `BitArray` of the same size as `∇c`, and captures
the output of `jacobian!`, which should return a Boolean specifying if the Jacobian is
constant or not.

The values are stored in `∇c`, which should be a matrix of matrices. If `con` is a
`StageConstraint`, `size(∇c,2) = 1`, and `size(∇c,2) = 2` if `con` is a `CoupledConstraint`.

If `con` is a `StageConstraint`, this will call `jacobian!(∇c, con, z)` by default, or
`jacobian!(∇c, con, z1, z2, i)` if `con` is a `CoupledConstraint`.
"""
function jacobian!(
    ∇c::VecOrMat{<:AbstractMatrix},
    con::StageConstraint,
    Z::AbstractTrajectory,
    inds = 1:length(Z),
    is_const = BitArray(undef, size(∇c))
)
    for (i, k) in enumerate(inds)
        is_const[i] = jacobian!(∇c[i], con, Z[k])
    end
end

function jacobian!(
    ∇c::VecOrMat{<:AbstractMatrix},
    con::CoupledConstraint,
    Z::AbstractTrajectory,
    inds = 1:size(∇c, 1),
    is_const = BitArray(undef, size(∇c))
)
    for (i, k) in enumerate(inds)
        is_const[i,1] = jacobian!(∇c[i, 1], con, Z[k], Z[k+1], 1)
        is_const[i,2] = jacobian!(∇c[i, 2], con, Z[k], Z[k+1], 2)
    end
end

"""
    ∇jacobian!(G, con::AbstractConstraint, Z, λ, inds, is_const, init)
    ∇jacobian!(G, con::AbstractConstraint, Z::AbstractKnotPoint, λ::AbstractVector)

Evaluate the second-order expansion of the constraint `con` along the trajectory `Z`
after multiplying by the lagrange multiplier `λ`.
The optional `is_const` argument is a `BitArray` of the same size as `∇c`, and captures
the output of `jacobian!`, which should return a Boolean specifying if the Jacobian is
constant or not. The `init` flag will force re-calculation of constant Jacobians when true.

The method for each constraint should calculate the Jacobian of the vector-Jacobian product,
    and therefore should be of size n × n if the input dimension is n.

Importantly, this method should ADD and not overwrite the contents of `G`, since this term
is dependent upon all the constraints acting at that time step.
"""
function ∇jacobian!(
    G::VecOrMat{<:AbstractMatrix},
    con::StageConstraint,
    Z::AbstractTrajectory,
    λ::Vector{<:AbstractVector},
    inds = 1:length(Z),
    is_const = ones(Bool, length(inds)),
    init::Bool = false,
)
    for (i, k) in enumerate(inds)
        if init || !is_const[i]
            is_const[i] = ∇jacobian!(G[i], con, Z[k], λ[i])
        end
    end
end

function ∇jacobian!(
    G::VecOrMat{<:AbstractMatrix},
    con::CoupledConstraint,
    Z::AbstractTrajectory,
    λ::Vector{<:AbstractVector},
    inds = 1:length(Z),
    is_const = ones(Bool, length(inds)),
    init::Bool = false,
)
    for (i, k) in enumerate(inds)
        if init || !is_const[i]
            is_const[i] = ∇jacobian!(G[i, 1], con, Z[k], Z[k+1], λ[i], 1)
            is_const[i] = ∇jacobian!(G[i, 2], con, Z[k], Z[k+1], λ[i], 2)
        end
    end
end

# Default methods for converting KnotPoints to states and controls for StageConstraints
"""
    evaluate(con::AbstractConstraint, z)       # stage constraint
    evaluate(con::AbstractConstraint, z1, z2)  # coupled constraint

Evaluate the constraint `con` at knot point `z`. By default, this method will attempt to call

    evaluate(con, x)

if `con` is a `StateConstraint`,

    evaluate(con, u)

if `con` is a `ControlConstraint`, or

    evaluate(con, x, u)

if `con` is a `StageConstraint`. If `con` is a `CoupledConstraint` the constraint should
define

    evaluate(con, z1, z2)

"""
@inline evaluate(con::StateConstraint, z::AbstractKnotPoint) = evaluate(con, state(z))
@inline evaluate(con::ControlConstraint, z::AbstractKnotPoint) = evaluate(con, control(z))
@inline evaluate(con::StageConstraint, z::AbstractKnotPoint) =
    evaluate(con, state(z), control(z))

"""
    jacobian!(∇c, con::AbstractConstraint, z, i=1)       # stage constraint
    jacobian!(∇c, con::AbstractConstraint, z1, z2, i=1)  # coupled constraint

Evaluate the constraint `con` at knot point `z`. By default, this method will attempt to call

    jacobian!(∇c, con, x)

if `con` is a `StateConstraint`,

    jacobian!(∇c, con, u)

if `con` is a `ControlConstraint`, or

    jacobian!(∇c, con, x, u)

if `con` is a `StageConstraint`. If `con` is a `CoupledConstraint` the constraint should
define

    jacobian!(∇c, con, z, i)

where `i` determines which Jacobian should be evaluated. E.g. if `i = 1`, the Jacobian
with respect to the first knot point's stage and controls is calculated.

# Automatic Differentiation
If `con` is a `StateConstraint` or `ControlConstraint` then this method is automatically
defined using ForwardDiff.
"""
jacobian!(∇c, con::StateConstraint, z::AbstractKnotPoint, i = 1) =
    jacobian!(∇c, con, state(z))
jacobian!(∇c, con::ControlConstraint, z::AbstractKnotPoint, i = 1) =
    jacobian!(∇c, con, control(z))
jacobian!(∇c, con::StageConstraint, z::AbstractKnotPoint, i = 1) =
    jacobian!(∇c, con, state(z), control(z))

# ForwardDiff jacobians that are of only state or control
function jacobian!(∇c, con::StageConstraint, x::StaticVector)
    eval_c(x) = evaluate(con, x)
    ∇c .= ForwardDiff.jacobian(eval_c, x)
    return false
end

@inline ∇jacobian!(G, con::StateConstraint, z::AbstractKnotPoint, λ, i = 1) =
    ∇jacobian!(G, con, state(z), λ)
@inline ∇jacobian!(G, con::ControlConstraint, z::AbstractKnotPoint, λ, i = 1) =
    ∇jacobian!(G, con, control(z), λ)
@inline ∇jacobian!(G, con::StageConstraint, z::AbstractKnotPoint, λ, i = 1) =
    ∇jacobian!(G, con, state(z), control(z), λ)

function ∇jacobian!(G, con::StageConstraint, x::StaticVector, λ)
    eval_c(x) = evaluate(con, x)'λ
    G_ = ForwardDiff.hessian(eval_c, x)
    G .+= G_
    return false
end

function ∇jacobian!(
    G,
    con::StageConstraint,
    x::StaticVector{n},
    u::StaticVector{m},
    λ,
) where {n,m}
    ix = SVector{n}(1:n)
    iu = SVector{m}(n .+ (1:m))
    eval_c(z) = evaluate(con, z[ix], z[iu])'λ
    G .+= ForwardDiff.hessian(eval_c, [x; u])
    return false
end

function gen_jacobian(con::AbstractConstraint, i = 1)
    ws = widths(con)
    p = length(con)
    C1 = SizedMatrix{p,ws[i]}(zeros(p, ws[i]))
end

function gen_views(∇c::AbstractMatrix, con::StateConstraint, n = state_dim(con), m = 0)
    view(∇c, :, 1:n), view(∇c, :, n:n-1)
end

function gen_views(∇c::AbstractMatrix, con::ControlConstraint, n = 0, m = control_dim(con))
    view(∇c, :, 1:0), view(∇c, :, 1:m)
end

function gen_views(
    ∇c::AbstractMatrix,
    con::AbstractConstraint,
    n = state_dim(con),
    m = control_dim(con),
)
    if size(∇c, 2) < n + m
        view(∇c, :, 1:n), view(∇c, :, n:n-1)
    else
        view(∇c, :, 1:n), view(∇c, :, n .+ (1:m))
    end
end


text:
using Tables
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

t0 = 0.0
δt = 0.1
tmax = 40.0
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25]; # β,c,γ

function euler(f, u0, p, δt, t0, tmax)
    t = t0 # Initialize time
    u = copy(u0) # Initialize struct parametric inherited
    du = zeros(length(u0)) # Initialize derivatives
    f(du,u,p,t)
    sol = [] # Store output
    times = [] # Store times
    push!(sol,copy(u))
    push!(times,t)
    # Main loop
    while t < tmax
        t = t + δt # Update time
        u .= u .+ du.*δt # Update state
        sir_ode!(du,u,p,t) # Update derivative
        push!(sol,copy(u)) # Store output
        push!(times,t) # Store time
    end
    sol = hcat(sol...) # Convert to matrix
    return times, sol
end;

times, sol = euler(sir_ode!, u0, p, δt, t0, tmax);

df = DataFrame(Tables.table(sol'))
rename!(df,["S","I","R"])
df[!,:t] = times;

@df df plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark euler(sir_ode!, u0, p, δt, t0, tmax)


text:
using Automa
using Automa.RegExp: @re_str
const re = Automa.RegExp

ab = re"ab*"
c = re"c"
pattern = re.cat(ab, c)

ab.actions[:enter] = [:enter_ab]
ab.actions[:exit]  = [:exit_ab]
ab.actions[:all]   = [:all_ab]
ab.actions[:final] = [:final_ab]
c.actions[:enter]  = [:enter_c]
c.actions[:exit]   = [:exit_c]
c.actions[:final]  = [:final_c]

write("actions.dot", Automa.machine2dot(Automa.compile(pattern)))
run(`dot -Tpng -o src/figure/actions.png actions.dot`)


text:
import Colors: XYZ, xyY
import Unitful.@u_str
using UnitfulColors
using Test

@testset "UnitfulColors.jl" begin
    @test colormatch(612.2u"nm") == colormatch(612.2)
    @test typeof(colormatch(600u"THz")) <:  XYZ
    @test typeof(colormatch(2.3u"eV"))  <:  XYZ
    @test typeof(colormatch(5000u"K"))  <:  XYZ
    @test typeof(colormatch(1.5e4u"K",:spectrum)) <:  XYZ
end


text:
function InitialNBody16(T=Float64)

GmMoon = parse(BigFloat,"0.109318945074237400e-10")
GmEarth = parse(BigFloat,"0.888769244512563400e-9")

Gm = [parse(BigFloat,"0.295912208285591100e-3"),
      parse(BigFloat,"0.491248045036476000e-10"),
      parse(BigFloat,"0.724345233264412000e-9"),
      GmEarth,
      parse(BigFloat,"0.954954869555077000e-10"),
      parse(BigFloat,"0.282534584083387000e-6"),
      parse(BigFloat,"0.845970607324503000e-7"),
      parse(BigFloat,"0.129202482578296000e-7"),
      parse(BigFloat,"0.152435734788511000e-7"),
      parse(BigFloat,"0.217844105197418000e-11"),
      GmMoon,
# Ceres, Pallas, Vesta, Iris, Bamberga
      parse(BigFloat,"0.140047655617234400e-12"),
      parse(BigFloat,"0.310444819893871300e-13"),
      parse(BigFloat,"0.385475018780881000e-13"),
      parse(BigFloat,"0.213643444257140700e-14"),
      parse(BigFloat,"0.138862658985619900e-14")
  ]

N=length(Gm)

q = [
parse(BigFloat,"0.00450250878464055477"),parse(BigFloat,"0.00076707642709100705"),parse(BigFloat,"0.00026605791776697764"),
parse(BigFloat,"0.36176271656028195477"),parse(BigFloat,"-0.09078197215676599295"),parse(BigFloat,"-0.08571497256275117236"),
parse(BigFloat,"0.61275194083507215477"),parse(BigFloat,"-0.34836536903362219295"),parse(BigFloat,"-0.19527828667594382236"),
parse(BigFloat,"0.12051741410138465477"),parse(BigFloat,"-0.92583847476914859295"),parse(BigFloat,"-0.40154022645315222236"),
parse(BigFloat,"-0.11018607714879824523"),parse(BigFloat,"-1.32759945030298299295"),parse(BigFloat,"-0.60588914048429142236"),
parse(BigFloat,"-5.37970676855393644523"),parse(BigFloat,"-0.83048132656339789295"),parse(BigFloat,"-0.22482887442656542236"),
parse(BigFloat,"7.89439068290953155477"),parse(BigFloat,"4.59647805517127300705"),parse(BigFloat,"1.55869584283189997764"),
parse(BigFloat,"-18.26540225387235944523"),parse(BigFloat,"-1.16195541867586999295"),parse(BigFloat,"-0.25010605772133802236"),
parse(BigFloat,"-16.05503578023336944523"),parse(BigFloat,"-23.94219155985470899295"),parse(BigFloat,"-9.40015796880239402236"),
parse(BigFloat,"-30.48331376718383944523"),parse(BigFloat,"-0.87240555684104999295"),parse(BigFloat,"8.91157617249954997764")
]

v = [
parse(BigFloat,"-0.00000035174953607552"),parse(BigFloat,"0.00000517762640983341"),parse(BigFloat,"0.00000222910217891203"),
parse(BigFloat,"0.00336749397200575848"),parse(BigFloat,"0.02489452055768343341"),parse(BigFloat,"0.01294630040970409203"),
parse(BigFloat,"0.01095206842352823448"),parse(BigFloat,"0.01561768426786768341"),parse(BigFloat,"0.00633110570297786403"),
parse(BigFloat,"0.01681126830978379448"),parse(BigFloat,"0.00174830923073434441"),parse(BigFloat,"0.00075820289738312913"),
parse(BigFloat,"0.01448165305704756448"),parse(BigFloat,"0.00024246307683646861"),parse(BigFloat,"-0.00028152072792433877"),
parse(BigFloat,"0.00109201259423733748"),parse(BigFloat,"-0.00651811661280738459"),parse(BigFloat,"-0.00282078276229867897"),
parse(BigFloat,"-0.00321755651650091552"),parse(BigFloat,"0.00433581034174662541"),parse(BigFloat,"0.00192864631686015503"),
parse(BigFloat,"0.00022119039101561468"), parse(BigFloat,"-0.00376247500810884459"),parse(BigFloat,"-0.00165101502742994997"),
parse(BigFloat,"0.00264276984798005548"),parse(BigFloat,"-0.00149831255054097759"),parse(BigFloat,"-0.00067904196080291327"),
parse(BigFloat,"0.00032220737349778078"),parse(BigFloat,"-0.00314357639364532859"),parse(BigFloat,"-0.00107794975959731297")
]

dqMoon = [
parse(BigFloat,"-0.00080817735147818490"), parse(BigFloat,"-0.0019946299854970130"),parse(BigFloat,"-0.00108726268307068900")
]
dvMoon = [
parse(BigFloat,"0.00060108481561422370"), parse(BigFloat,"-0.00016744546915764980"), parse(BigFloat,"-0.00008556214140094871")
];

qq = copy(q)
vv = copy(v)

qEM = q[10:12]
qEarth = qEM - GmMoon/GmEarth*dqMoon
qMoon = qEarth + dqMoon
qq[10:12] = qEarth

vEM = v[10:12]
vEarth = vEM - GmMoon/GmEarth*dvMoon
vMoon = vEarth + dvMoon
vv[10:12] = vEarth

# Ceres, Pallas, Vesta, Iris, Bamberga
qsun=qq[1:3]
qCeres=qsun+[parse(BigFloat,"1.438681809676469747"),parse(BigFloat,"-2.204373633189407045"),parse(BigFloat,"-1.326397853361325874")]
qPallas=qsun+[parse(BigFloat,"0.203832272462290465"),parse(BigFloat,"-3.209619436062307152"),parse(BigFloat,"0.623843179079393351")]
qVesta=qsun+[parse(BigFloat,"0.182371836377417107"),parse(BigFloat,"2.386628211277654010"),parse(BigFloat,"0.924596062836265498")]
qIris=qsun+[parse(BigFloat,"1.892475267790300286"),parse(BigFloat,"-0.848414748075139946"),parse(BigFloat,"-0.157159319044464590")]
qBamberga=qsun+[parse(BigFloat,"1.398759064223541682"), parse(BigFloat,"-1.287476729008325105"), parse(BigFloat,"-0.669098428660833799")]

vsun=vv[1:3]
vCeres=vsun+[parse(BigFloat,"0.008465406136316316"), parse(BigFloat,"0.004684247977335608"),parse(BigFloat,"0.000466157738595739")]
vPallas=vsun+[parse(BigFloat,"0.008534313855651248"), parse(BigFloat,"-0.000860659210123161"), parse(BigFloat,"-0.000392901992572746")]
vVesta=vsun+[parse(BigFloat,"-0.010174496747119257"), parse(BigFloat,"0.000041478190529952"), parse(BigFloat,"0.001344157634155624")]
vIris=vsun+[parse(BigFloat,"0.002786950314570632"), parse(BigFloat,"0.011314057384917047"), parse(BigFloat,"0.004975132577079665")]
vBamberga=vsun+[parse(BigFloat,"0.007164363244556328"), parse(BigFloat,"0.009219958777618218"), parse(BigFloat,"0.006857861727407507")]

qq = [qq; qMoon; qCeres; qPallas; qVesta; qIris; qBamberga]
vv = [vv; vMoon; vCeres; vPallas; vVesta; vIris; vBamberga]

q0 = reshape(qq,3,:)
v0 = reshape(vv,3,:)

q0bar = [sum(Gm .* q0[j,:])/sum(Gm) for j in 1:3]
v0bar = [sum(Gm .* v0[j,:])/sum(Gm) for j in 1:3]

q0 = q0 .- q0bar
v0 = v0 .- v0bar

u0 = Array{T}(undef,2,3,N)
u0[1,:,:] = v0
u0[2,:,:] = q0

return u0, convert.(T,Gm)
end


text:
import ..UncertainValues:
    getquantileindex,
    support,
    AbstractUncertainScalarKDE

import StatsBase:
    sample, sample!,
    Weights

import Distributions:
    Uniform

import ..SamplingConstraints:
	NoConstraint,
	TruncateLowerQuantile,
	TruncateUpperQuantile,
	TruncateQuantiles,
	TruncateMinimum,
	TruncateMaximum,
	TruncateRange,
    TruncateStd,
    fallback

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::NoConstraint)

Resample without contraints (use the full distribution representing the value)

## Example

```julia
some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

# Resample the uncertain value by resampling the full distribution once.
resample(uncertainval, NoConstraint())
```
"""
resample(uv::AbstractUncertainScalarKDE, constraint::NoConstraint) = resample(uv)

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::NoConstraint, n::Int)

Resample without contraints (use the full distribution representing the value)

## Example

```julia
some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

# Resample the uncertain value by resampling the full distribution n times
resample(uncertainval, NoConstraint(), n)
```
"""
resample(uv::AbstractUncertainScalarKDE, constraint::NoConstraint, n::Int) = resample(uv, n)

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateLowerQuantile)

Resample `uv` by first truncating below the kernel density estimate of the
distribution furnishing the value at some lower quantile, then resampling
it once.

## Example

```julia
using UncertainData

some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

constraint = TruncateLowerQuantile(0.16)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateLowerQuantile)
    # Find the index of the kernel density estimated distribution
    # corresponding to the lower quantile at which we want to truncate.
    idx_lower_quantile = getquantileindex(uv, constraint.lower_quantile)

    # Box width
    δ = step(uv.range)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_lower_quantile:end] |> collect
    wts = Weights(uv.pdf[idx_lower_quantile:end])

    # Sample a box
    sampled_val = sample(range, wts)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateLowerQuantile,
        n::Int)

Resample `uv` by first truncating below the kernel density estimate of the
distribution furnishing the value at some lower quantile, then resampling
it `n` times.

## Example

```julia
using UncertainData

some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

constraint = TruncateLowerQuantile(0.16)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 500 times.
resample(uncertainval, constraint, 500)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateLowerQuantile,
        n::Int)
    # Find the index of the kernel density estimated distribution
    # corresponding to the lower quantile at which we want to truncate.
    idx_lower_quantile = getquantileindex(uv, constraint.lower_quantile)

    # Box width
    δ = step(uv.range)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_lower_quantile:end] |> collect
    wts = Weights(uv.pdf[idx_lower_quantile:end])

    # Sample n boxes according to estimated weights (pdf)
    sampled_vals = Vector{Float64}(undef, n)
    sample!(range, wts, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile)

Resample `uv` by first truncating above the kernel density estimate of the
distribution furnishing the value at some upper quantile, then resampling
it once.

## Example

```julia
using UncertainData

some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

constraint = TruncateUpperQuantile(0.78)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile)
    # Find the index of the kernel density estimated distribution
    # corresponding to the lower quantile at which we want to truncate.
    idx_upper_quantile = getquantileindex(uv, constraint.upper_quantile)

    # Box width
    δ = step(uv.range)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[1:idx_upper_quantile] |> collect
    wts = Weights(uv.pdf[1:idx_upper_quantile])

    # Sample a box
    sampled_val = sample(range, wts)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile,
        n::Int)

Resample `uv` by first truncating above the kernel density estimate of the
distribution furnishing the value at some upper quantile, then resampling
it `n` times.

## Example

```julia
using UncertainData

some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

constraint = TruncateLowerQuantile(0.16)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 500 times.
resample(uncertainval, constraint, 500)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile,
        n::Int)
    # Find the index of the kernel density estimated distribution
    # corresponding to the lower quantile at which we want to truncate.
    idx_upper_quantile = getquantileindex(uv, constraint.upper_quantile)

    # Box width
    δ = step(uv.range)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[1:idx_upper_quantile] |> collect
    wts = Weights(uv.pdf[1:idx_upper_quantile])

    # Sample n boxes according to estimated weights (pdf)
    sampled_vals = Vector{Float64}(undef, n)
    sample!(range, wts, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateUpperQuantile)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value both above and below at some quantile range,
then resampling it once.

## Example

```julia
using UncertainData

some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

constraint = TruncateQuantiles(0.1, 0.9)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateQuantiles)
    # Find the index of the kernel density estimated distribution
    # corresponding to the lower quantile at which we want to truncate.
    idx_lower_quantile = getquantileindex(uv, constraint.lower_quantile)
    idx_upper_quantile = getquantileindex(uv, constraint.upper_quantile)

    # Box width
    δ = step(uv.range)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_lower_quantile:idx_upper_quantile] |> collect
    wts = Weights(uv.pdf[idx_lower_quantile:idx_upper_quantile])

    # Sample a box
    sampled_val = sample(range, wts)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateQuantiles,
        n::Int)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value both above and below at some quantile range,
then resampling it `n` times.

## Example

```julia
using UncertainData

some_sample = rand(Normal(), 1000)

# Calling UncertainValue with a single vector of numbers triggers KDE estimation
uncertainval = UncertainValue(some_sample) # -> UncertainScalarKDE

constraint = TruncateQuantiles(0.1, 0.9)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 500 times.
resample(uncertainval, constraint, 500)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateQuantiles,
        n::Int)
    # Find the index of the kernel density estimated distribution
    # corresponding to the lower quantile at which we want to truncate.
    idx_upper_quantile = getquantileindex(uv, constraint.upper_quantile)
    idx_lower_quantile = getquantileindex(uv, constraint.lower_quantile)

    # Box width
    δ = step(uv.range)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_lower_quantile:idx_upper_quantile] |> collect
    wts = Weights(uv.pdf[idx_lower_quantile:idx_upper_quantile])

    # Sample n boxes according to estimated weights (pdf)
    sampled_vals = Vector{Float64}(undef, n)
    sample!(range, wts, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMaximum)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value at some maximum value,
then resampling it once.

## Example

```julia
# Uncertain value represented by a normal distribution with mean = 0 and
# standard deviation = 1.
uncertainval = UncertainValue(rand(Normal(0, 1), 1000))

constraint = TruncateMaximum(0.8) # accept no values larger than 1.1

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 700 times.
resample(uncertainval, constraint, 700)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMaximum)
    # Box width
    δ = step(uv.range)

    upper_bound = constraint.max
    lower_bound = minimum(uv.range)

    idx_max = findlast(uv.range .<= upper_bound)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[1:idx_max] |> collect
    wts = Weights(uv.pdf[1:idx_max])

    # Sample a box
    sampled_val = sample(range, wts)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMaximum, n::Int)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value at some minimum value,
then resampling it `n` times.

## Example

```julia
# Uncertain value represented by a normal distribution with mean = 0 and
# standard deviation = 1.
uncertainval = UncertainValue(rand(Normal(0, 1), 1000))

constraint = TruncateMaximum(0.8) # accept no values larger than 1.1

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 700 times.
resample(uncertainval, constraint, 700)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMaximum, n::Int)
    # Box width
    δ = step(uv.range)

    upper_bound = constraint.max
    lower_bound = minimum(uv.range)

    idx_max = findlast(uv.range .<= upper_bound)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[1:idx_max] |> collect
    wts = Weights(uv.pdf[1:idx_max])

    # Sample n boxes according to estimated pdf
    sampled_vals = Vector{Float64}(undef, n)
    sample!(range, wts, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMinimum)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value at some minimum value,
then resampling it once.

## Example

```julia
# Uncertain value represented by a normal distribution with mean = 0 and
# standard deviation = 1.
uncertainval = UncertainValue(rand(Normal(0, 1), 1000))

constraint = TruncateMinimum(0.2) # accept no values smaller than 0.2

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 700 times.
resample(uncertainval, constraint, 700)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMinimum)
    # Box width
    δ = step(uv.range)

    lower_bound = constraint.min
    upper_bound = maximum(uv.range)

    idx_min = findfirst(uv.range .>= lower_bound)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_min:end] |> collect
    wts = Weights(uv.pdf[idx_min:end])

    # Sample a box
    sampled_val = sample(range, wts)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMinimum, n::Int)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value at some minimum value,
then resampling it `n` times.

## Example

```julia
# Uncertain value represented by a normal distribution with mean = 0 and
# standard deviation = 1.
uncertainval = UncertainValue(rand(Normal(0, 1), 1000))

constraint = TruncateMinimum(0.2) # accept no values smaller than 0.2

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 700 times.
resample(uncertainval, constraint, 700)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateMinimum, n::Int)
    # Box width
    δ = step(uv.range)

    lower_bound = constraint.min
    upper_bound = maximum(uv.range)

    idx_min = findfirst(uv.range .>= lower_bound)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_min:end] |> collect
    wts = Weights(uv.pdf[idx_min:end])

    # Sample n boxes according to estimated pdf
    sampled_vals = Vector{Float64}(undef, n)
    sample!(range, wts, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateRange)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value at some minimum and maximum values,
then resampling it once.

## Example

```julia
# Uncertain value represented by a normal distribution with mean = 0 and
# standard deviation = 1.
uncertainval = UncertainValue(rand(Normal(0, 1), 1000))

# Only accept values in the range [-0.9, 1.2]
constraint = TruncateRange(-0.9, 1.2)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 300 times.
resample(uncertainval, constraint, 300)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateRange)
    # Box width
    δ = step(uv.range)

    lower_bound = constraint.min
    upper_bound = constraint.max

    idx_min = findfirst(uv.range .>= lower_bound)
    idx_max = findlast(uv.range .<= upper_bound)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_min:idx_max] |> collect
    wts = Weights(uv.pdf[idx_min:idx_max])

    # Sample a box
    sampled_val = sample(range, wts)

    # Sample uniformly from within the box
    rand(Uniform(sampled_val, sampled_val + δ))
end

"""
    resample(uv::AbstractUncertainScalarKDE, constraint::TruncateRange, n::Int)

Resample `uv` by first truncating the kernel density estimate of the
distribution furnishing the value at some minimum and maximum values,
then resampling it `n` times.

## Example

```julia
# Uncertain value represented by a normal distribution with mean = 0 and
# standard deviation = 1.
uncertainval = UncertainValue(rand(Normal(0, 1), 1000))

# Only accept values in the range [-0.9, 1.2]
constraint = TruncateRange(-0.9, 1.2)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 300 times.
resample(uncertainval, constraint, 300)
```
"""
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateRange, n::Int)
    # Box width
    δ = step(uv.range)

    lower_bound = constraint.min
    upper_bound = constraint.max

    idx_min = findfirst(uv.range .>= lower_bound)
    idx_max = findlast(uv.range .<= upper_bound)

    # Subset the values and weights (values of the pdf at those values)
    range = uv.range[idx_min:idx_max] |> collect
    wts = Weights(uv.pdf[idx_min:idx_max])

    # Sample n boxes according to estimated pdf
    sampled_vals = Vector{Float64}(undef, n)
    sample!(range, wts, sampled_vals)

    # Sample uniformly from within each box
    [rand(Uniform(sampled_vals[i], sampled_vals[i] + δ)) for i = 1:n]
end

# Resampling UncertainScalarKDE with TruncateStd is undefined, so fall back to quantiles
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateStd)
    resample(uv, fallback(uv, constraint))
end

# Resampling UncertainScalarKDE with TruncateStd is undefined, so fall back to quantiles
function resample(uv::AbstractUncertainScalarKDE, constraint::TruncateStd, n::Int)
    resample(uv, fallback(uv, constraint), n)
end