text:
module RecursiveFiniteDifferences

# TODO: clean this up :)
# This is simply the notebook also shown in the docs.
# I just wrapped it in a module, fixed the import of Trees and commented out some examples

include(joinpath(@__DIR__, "Trees.jl"))
using .Trees

### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ e7d76f0e-eea5-11eb-1bdf-9196c513ff0a
begin
	using StaticArrays
	using Symbolics
	using Symbolics: Sym
	using Parameters
end

# ╔═╡ 7b914a4f-89b4-4284-b8c6-d629157b9a7c
md"# Finite Differences for Taylor Models."

# ╔═╡ 44f2f27b-2ac6-4df1-b222-4cfa37f41985
begin
	const RVec = AbstractVector{<:Real}
	const RVecOrR = Union{Real, RVec}
end

# ╔═╡ 10b9b062-450e-43a2-9c7a-2fe6c7121ab7
md"""
### Abstract

For calculating finite difference gradients and hessians in Morbit we want to separate the grid construction from the function evaluation and derivative calculations.
To my knowledge, this is not easily possible using `FiniteDiff` or `FiniteDifferences.jl`.

This notebook explores a recursive calculation procedure starting from first order finite difference rules (interfaced by `FiniteDiffStamp`).
To calculate partial derivatives of ``f:ℝ^n \to ℝ^k`` at a fixed point ``x_0 ∈ ℝ^n`` a `DiffWrapper` can be used.
It is initialized devoid of any function evaluation information but only has symbolic representations of required evaluation sites relative to ``x_0``.
A `DiffWrapper` can hence be used to first construct a grid for ``x_0`` and at a later stage to approximate gradients or Hessians.
"""

# ╔═╡ bd83fc2a-d7bd-489e-83f3-ef5b327d5cca
md"""
## Finite Difference Rules

A finite difference rule is defined by a grid of integer offsets, a decision space stepsize ``h``, and coefficients to combine evaluation results at the grid points.

For details, please see [this Wikipedia page](https://en.wikipedia.org/wiki/Finite_difference_coefficient).

We abstract the basic information with the `FiniteDiffStamp` interface:
"""

# ╔═╡ a7d81ac5-0124-441a-8eae-f83cee2b81c3
begin
	abstract type FiniteDiffStamp end

	# integer offsets to build a grid along some variable dimension
	_grid( :: FiniteDiffStamp ) = nothing
	# stepsize to build a grid along some variable dimension in decision space
	_stepsize( :: FiniteDiffStamp ) = nothing

	# coefficients of the rule
	_coeff( :: FiniteDiffStamp ) = nothing

	# order of the method
	_order( :: FiniteDiffStamp ) = nothing

	# accuracy (exponent of residual in Taylor expansion)
	_accuracy( :: FiniteDiffStamp ) = nothing

	# zero index of grid
	_zi( :: FiniteDiffStamp ) = nothing

	Base.broadcastable( fds :: FiniteDiffStamp ) = Ref(fds)

	# evaluation of rule at sites in `X` (a vector or an iterable)
	function (stamp :: FiniteDiffStamp)(X)
		coeff = _coeff(stamp)
		@assert length(X) == length(coeff)
		h = _stepsize(stamp)
		return sum( c * x for (c,x) in zip(coeff,X) ) ./ h^(_order(stamp))
	end
end

# ╔═╡ fdcf00a4-836e-4404-ba1c-6058ae85d9ce
md"We could already define a finite difference function using a `FiniteDiffStamp`:"

# ╔═╡ 07e698bb-fac0-4951-85dc-5adb208cf5dd
md"Of course, some concrete implementations are still missing.
We provide some entries from the table in [Wikipedia](https://en.wikipedia.org/wiki/Finite_difference_coefficient)."

# ╔═╡ 9462b310-0842-4373-ac2c-f2eaef5e6e7f
begin
	stepsize( x :: F ) where F <: AbstractFloat = 10 * sqrt(eps(F))
	stepsize( F :: Type{<:AbstractFloat} ) = 10 * sqrt(eps(F))
	stepsize( x :: AbstractVector{F} ) where F<:AbstractFloat = stepsize(F)
	stepsize( x ) = stepsize(Float64)
end

# ╔═╡ c19325b6-7516-4acf-a5f1-fe7084bc113d
begin
	struct CFDStamp{N, F <: AbstractFloat} <: FiniteDiffStamp
		grid :: SVector{N,Int8}
		coeff :: SVector{N,Rational{Int16}}
		order :: UInt8
		accuracy :: UInt8

		stepsize :: F

		zero_index :: UInt8
	end

	_grid( stamp :: CFDStamp ) = stamp.grid
	_coeff( stamp :: CFDStamp ) = stamp.coeff
	_order( stamp :: CFDStamp ) = stamp.order
	_accuracy( stamp :: CFDStamp ) = stamp.accuracy
	_stepsize( stamp :: CFDStamp ) = stamp.stepsize
	_zi( stamp :: CFDStamp ) = stamp.zero_index

	function CFDStamp( order :: Int, accuracy :: Int, _stepsize_ = stepsize(Float64) )
		return CFDStamp( Val(order), Val( floor(Int, accuracy / 2) * 2 ), _stepsize_ )
	end

	function CFDStamp( ::Val{1}, ::Val{2}, _stepsize_ :: AbstractFloat )
		CFDStamp(
			SVector{3,Int8}([-1, 0, 1]),
			SVector{3,Rational{Int16}}([-1//2, 0, 1//2]),
			UInt8(1), UInt8(2), _stepsize_, UInt8(2),
		)
	end

	function CFDStamp( ::Val{1}, ::Val{4}, _stepsize_ :: AbstractFloat )
		CFDStamp(
			SVector{5,Int8}([-2,-1,0,1,2]),
			SVector{5,Rational{Int16}}([1//12, -2//3, 0, 2//3, -1//12]),
			UInt8(1), UInt8(4), _stepsize_, UInt8(3)
		)
	end

	function CFDStamp( ::Val{1}, ::Val{6}, _stepsize_ :: AbstractFloat )
		CFDStamp(
			SVector{7,Int8}(collect(-3:1:3)),
			SVector{7,Rational{Int16}}([-1//60, 3//20, -3//4, 0, 3//4, -3//20, 1//60]),
			UInt8(1), UInt8(6), _stepsize_, UInt8(4)
		)
	end

end

# ╔═╡ 2426f92a-21f3-4455-aeb7-fe6674aafe02
begin
	struct FFDStamp{N, F <: AbstractFloat} <: FiniteDiffStamp
		grid :: SVector{N,Int8}
		coeff :: SVector{N,Rational{Int16}}
		order :: UInt8
		accuracy :: UInt8

		stepsize :: F
	end

	_grid( stamp :: FFDStamp ) = stamp.grid
	_coeff( stamp :: FFDStamp ) = stamp.coeff
	_order( stamp :: FFDStamp ) = stamp.order
	_accuracy( stamp :: FFDStamp ) = stamp.accuracy
	_stepsize( stamp :: FFDStamp ) = stamp.stepsize
	_zi( stamp :: FFDStamp ) = 1

	function FFDStamp( order :: Int, accuracy :: Int, _stepsize_ = stepsize(Float64) )
		return FFDStamp( Val(order), Val( ceil(Int, accuracy / 2) ), _stepsize_ )
	end

	function FFDStamp( ::Val{1}, ::Val{1}, _stepsize_ :: AbstractFloat )
		FFDStamp(
			SVector{2,Int8}([0, 1]),
			SVector{2,Rational{Int16}}([-1, 1]),
			UInt8(1), UInt8(1), _stepsize_
		)
	end

	function FFDStamp( ::Val{1}, ::Val{2}, _stepsize_ :: AbstractFloat )
		FFDStamp(
			SVector{3,Int8}([0, 1, 2]),
			SVector{3,Rational{Int16}}([-3//2, 2, -1//2]),
			UInt8(1), UInt8(2), _stepsize_
		)
	end

	function FFDStamp( ::Val{1}, ::Val{3}, _stepsize_ :: AbstractFloat )
		FFDStamp(
			SVector{4,Int8}([0, 1, 2, 3]),
			SVector{4,Rational{Int16}}([-11/6, 3, -3//2, 1//3]),
			UInt8(1), UInt8(3), _stepsize_
		)
	end

end

# ╔═╡ 4ebb80a7-539d-49a0-981f-166f0bbf3964
begin
	struct BFDStamp{N, F <: AbstractFloat} <: FiniteDiffStamp
		grid :: SVector{N,Int8}
		coeff :: SVector{N,Rational{Int16}}
		order :: UInt8
		accuracy :: UInt8

		stepsize :: F
	end

	_grid( stamp :: BFDStamp ) = stamp.grid
	_coeff( stamp :: BFDStamp ) = stamp.coeff
	_order( stamp :: BFDStamp ) = stamp.order
	_accuracy( stamp :: BFDStamp ) = stamp.accuracy
	_stepsize( stamp :: BFDStamp ) = stamp.stepsize
	_zi( stamp :: BFDStamp ) = 1

	function BFDStamp( order :: Int, accuracy :: Int, _stepsize_ = stepsize(Float64) )
		ffd_stamp = FFDStamp( Val(order), Val( ceil(Int, accuracy / 2) ), _stepsize_ )
		if isodd(order)
			return BFDStamp(
				- _grid( ffd_stamp ),
				- _coeff( ffd_stamp ),
				_order(ffd_stamp), _accuracy(ffd_stamp), _stepsize(ffd_stamp)
			)
		else
			return  BFDStamp(
				_grid( ffd_stamp ), _coeff( ffd_stamp ),
				_order(ffd_stamp), _accuaracy(ffd_stamp), _stepsize(ffd_stamp)
			)
		end
	end

end

# ╔═╡ aeac724f-b87a-4f37-aa57-893d1d131f6d
"Evaluate all outputs of `f` with respect to ``x_i``."
function dx(stamp :: FiniteDiffStamp, f :: Function, x0 :: RVec, i = 1)
	h = _stepsize(stamp)
	grid = [ [x0[1:i-1];x0[i] + h*g; x0[i+1:end]] for g in _grid(stamp) ]
	evals = f.(grid)
	return stamp( evals )
end

# ╔═╡ a5e02956-c0f9-44a6-80da-abfb497bcfb9
md"Does it work?"

# ╔═╡ 09d78bbc-f568-459a-8bd0-05de84b04af2
begin
	cfd = CFDStamp(1,2)
	Func = x -> x[1]*x[2]^2

	# ∂₁Func = x[2]^2
	# ∂₂Func = x[1]
	dx( cfd, Func, [1,2], 1)
end

# ╔═╡ 04861d1b-38d9-4ba7-9fde-b10150f39280
md"## Tree structure"

# ╔═╡ 20907685-2bb6-49d5-9ff0-a175435fa591
md"""
Our goal is to recursively approximate high order derivatives of a function ``f\colon ℝ^n \to ℝ^k`` using a first order finite difference formula ``G_f``.
We assume, that ``G_f`` requires ``m`` values of derivative order ``i-1`` to approximate the partial derivative of order ``i``.

For simplicity, first consider ``k=1`` (the single output case).
* For forming ``∂_1 f(x_0)`` we would need values ``f(ξ_1^1), …, f(ξ_m^1)``, where ``ξ^1`` varies along the first coordinate and results from applying ``G_f`` to ``x_0``.
* For forming ``∂_1 ∂_1 f(x_0)`` we would instead need values approximating ``∂_1 f(ξ_1^1),…, ∂_1 f(ξ_m^1)``.
  Using the same rule to approximate the first order partial derivatives requires function evaluations at sites ``ξ_{1,1}^1, …, ξ_{1,m}^1, ξ_{2,1}^1, …, ξ_{m,m}^1`` resulting from applying ``G_f`` to each of ``ξ_i^1``.

This process can be recursed infinitely.
To actually compute the desired derivative approximation we would first have to construct the evaluation sites ``ξ_I^J``.
We can think of a tree-like structure, with a root node containing ``x_0`` where we add -- for each derivative order -- child nodes by applying ``G_f`` to the previous leaves. \
The tree is hence build from top to bottom but for evaluation we start at the leaves and resolve the lowest order derivative approximations first.
In a sense, this is a dynamic programming approach.
"""

# ╔═╡ 65d0f18a-6da9-4836-bf4c-8ff10e6a265f
md"""
Our tree is made out of `FDiffNode`s.
An `FDiffNode` has a field `x_sym` storing the symbolic representation of some point relative to ``x_0`` that we need a value for. \
If `N :: FDiffNode{T,X,C}` is a leave node, then `T` should be a mutable vector type and `vals :: T` is meant to store the result of ``f`` evaluated at ``x``, where ``x`` results from substituting ``x_0`` and a stepsize ``h`` into `x_sym`. \
If `N` is *not* a leave node, then `T` should be a vector type to store ``n`` vectors of `FDiffNode`s and each of those vectors should have ``m`` elements, so that it can be used to recursively approximate a derivative.

Note, that the variables `x_vars::Vector{Symbolics.Num}` and the finite difference rule are stored outside of the tree to keep it as simple as possible.

We also allow to use a cache of type ``C`` which should be Dict-Like (if not `nothing`).
It can be used in the `val` method to retrieve values that have already been calculated before.
"""

# ╔═╡ 6cdc6c1e-790e-4d2a-984c-36f5f0c98f63
md"We inherit from `Trees.Node` and define a `children` method` so that some nice iterators are available."

# ╔═╡ dc4d3c61-f2ec-4fb1-b84b-2c1957da3a7e
md"For leave nodes (indicated by `T` being a real vector type), we simply retrieve the function values if they are present. Else, `NaN` is returned.
We have a helper function `missing` to get the right type of `NaN` and the important `val` retrieval function."

# ╔═╡ c2cb0dd1-e36d-44ff-aa59-6365a15e7ab6
begin
	missing(T::Type{<:AbstractFloat}) = T(NaN)
	missing(T) = Float16(NaN)
end

# ╔═╡ 3f67de97-2a12-4ca7-82bd-16851f958989
md"Otherwise, we *recurse* and apply the finite difference rule `stamp` to compute the return values."

# ╔═╡ 18e51761-98ae-4814-b575-3bb4a40a06cb
function val( subnode_list, indices, stamp, output_index )
	coeff = _coeff( stamp )
	h = _stepsize( stamp )
	return stamp( val(sub_node, indices, stamp; output_index) for sub_node in subnode_list )
end

# ╔═╡ a5565a82-2b6d-4c43-909c-9bfeca9890cb
md"""The tree is build recursively from the top down by calling the `build_tree` function with a decreasing `order` value.
* An `order` value of 0 indicates that we need need a leave node for `x_sym` and the correct value container is initialized (and the cache deactivated).
* For a higher `order` value, we collect the ``n`` vectors of ``m`` sub nodes by calling `buid_tree` again.
"""

# ╔═╡ 1db29942-3ff0-49e8-b5c6-2480eadd608f
md"`build_tree` is called in default the `DiffWrapper` constructor to store the `tree` (out of `FDiffNode`s) for the derivative `order` stored within the container.
A `DiffWrapper` also stores the base point `x0` and the finite difference rule `stamp`.
"

# ╔═╡ af650705-0636-4728-a643-ab658d362290
begin
	"Helper function to get Symbolic variables for `x_1,…,x_n` and the stepsize `h`."
	function _get_vars(n :: Int)
		return Tuple(Num.(Symbolics.variable.(:x, 1:n))), Num(Symbolics.variable(:h))
	end

	_get_vars( x :: AbstractVector ) = _get_vars(length(x))
end

# ╔═╡ 7c091ff3-7d34-4ad8-b959-fe601397f091
begin
	vec_typex( x :: AbstractVector{Y}, stamp ) where Y = Vector{ Base.promote_op( +, Y, typeof(_stepsize(stamp))) }
	vec_typef( fx :: AbstractVector{Y}, stamp ) where Y = Vector{ Base.promote_op( *, Y, typeof(_coeff(stamp)[1])) }
end

# ╔═╡ 59211507-8901-4f4c-b2ab-c5d4fb950bf1
md"""
!!! note
    A `DiffWrapper` should be initialized with the right (floating point) vector types for `x0` and `fx0`. If `fx0` is not known, but the precision, use something like `fx0 = Float32[]`.
"""

# ╔═╡ eb61f081-bd96-4a25-b2ba-9a3bc9e87934
md"We exploit the meta data stored to forward the `val` method.
`indices` should be an array of `Int`s of length `dw.order`, indicating the partial derivatives we want.
Suppose, `dw.order == 2`, then `val( dw, [1,2] )` will give ``∂_1∂_2f_1(x_0)``.
"

# ╔═╡ f77f8375-5b8f-4a3e-b571-28ebccb204c4
md"""
!!! note
    The `tree` of a `DiffWrapper` needs to be initialized befor calling any derivative method.
    That means, the `x` and `vals` fields have to be set for each leave node.
	This can be done in one go, by calling `prepare_tree!( dw, f )`, or sequentially by calling first `substitute_leaves!(dw)` and and `set_leave_values(dw,f)`.
"""

# ╔═╡ fe8171dc-7efd-4acb-b1f9-0da5f1c6dc9d
md"From this, it is easy to define convenience functions for hessians and gradients.
For second order `DiffWrapper`s we can also make use of the fact, that each derivative rule has a grid point with index `zi` where ``x_0`` is not varied.
That means, we can collect the gradient from the leaves with the parent node of index `zi` (relative to the root)."

# ╔═╡ 5dd33903-f4a4-455c-8348-c3c98e5f6404
md"Does it work?"

# ╔═╡ 8ad29cb3-9bb5-48d6-b0a6-ceb46e30293e
md"""
#### Helpers for filling the Tree

I won't go into detail concering these helper functions. They mostly leverage the `Trees.Leaves` iterator and should be comprehensible by themselves.
"""

# ╔═╡ fef9850c-3cbd-47f8-a786-e420f01d9141
function _substitute_symbols!(root_iterator, x0, h, vars)
	for node in root_iterator
		_substitute_symbols!(node, x0, h, vars)
	end
end

# ╔═╡ ac3d101b-c9d1-4fdb-80a4-645464298ba3
function _set_node_values!( node, f :: Function )
	empty!(node.vals)
	append!(node.vals, f( node.x ))
end

# ╔═╡ 5d586923-a32e-4d16-a0a6-56bbc73680f2
md"""
### Morbit Example
"""

# ╔═╡ 3feed8b1-1338-47c7-960f-892481f96a5c
function unique_with_indices( x :: AbstractVector{T} ) where T
	unique_elems = T[]
	indices = Int[]
	for elem in x
		i = findfirst( e -> all( isequal.(e,elem) ), unique_elems )
		if isnothing(i)
			push!(unique_elems, elem)
			push!(indices, length(unique_elems) )
		else
			push!(indices, i)
		end
	end
	return unique_elems, indices
end

# ╔═╡ d41308c0-b7f1-4ed9-b27e-7c6fe87080e6
md"---"

# ╔═╡ f294835d-ee9a-42b1-b2b7-5816d65e0f18
md"`ingredients` thanks to [fonsp](https://github.com/fonsp/Pluto.jl/issues/115)"

# ╔═╡ e9af48fd-ac05-4358-8a33-c14008ca1cb9
#=
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end
=#

# ╔═╡ 6c6e63ae-28ce-4c1c-bb60-7e7826ffefbb
# Trees = ingredients( joinpath(@__DIR__, ".." ,"src", "Trees.jl") )

# ╔═╡ c3fa417a-898f-4273-a75c-974a42e80a79
@with_kw struct FDiffNode{T,X,C} <: Trees.Node
	x_sym :: Vector{Symbolics.Num} = []
	x :: X = Float64[]
	vals :: T = nothing
	cache :: C = nothing
end

# ╔═╡ aabc4fed-e62f-4668-9984-c9b72cad15e9
function val( n :: FDiffNode{<:AbstractVector{Y},X,C}, args...; output_index = 1) where {Y<:Real,X,C}
	if isempty(n.vals)
		return missing(Y)
	else
		n.vals[output_index]
	end
end

# ╔═╡ d21fed2f-f61a-42d6-aed3-430bf73fd170
function val( n :: FDiffNode{T,X,C}, indices, stamp; output_index = 1 ) where{T,X,C}
	copied_indices = copy(indices)
	i = popfirst!( copied_indices )

	if C <: Nothing
		return val( n.vals[i], copied_indices, stamp, output_index )
	end

	cache_key = [ i; copied_indices; output_index ]
	if haskey( n.cache, cache_key )
		ret_val = n.cache[cache_key]
		!isnan(ret_val) && return ret_val
	end

	ret_val = val( n.vals[i], copied_indices, stamp, output_index )
	n.cache[cache_key] = ret_val
	return ret_val
end

# ╔═╡ 17347e0e-08ee-4a29-83b3-f337e397f49b
function build_tree( x_sym, stamp, vars, order = 1; x_type = Vector{Float64}, val_type = Vector{Float64}, cache_type = Nothing )

	if order <= 0
		# return a leave node
		return FDiffNode(; x_sym, x = x_type(), vals = val_type(), cache = nothing )
	else
		# collect ``n`` subnode vectors of length ``m``
		grid = _grid( stamp )	# ``m`` grid points to define the points for sub_nodes
		m = length( grid )
		n = length( x_sym )
		h = vars[2]
		#sub_nodes = SizedVector{n, SizedVector{m, FDiffNode}}(undef)
		sub_nodes = SizedVector{n, SizedVector{m, FDiffNode}}(undef)
		for i = 1 : n
			sub_nodes[i] = [
				build_tree(
					# vary `x_sym` along variable `i`
					[x_sym[1:i-1]; x_sym[i] + h * g; x_sym[i+1:n]],
					stamp, vars, order - 1; x_type, val_type, cache_type )
				for g in grid
			]
		end
		return FDiffNode(; x_sym, x = x_type(), vals = sub_nodes, cache = cache_type() )
	end
end

# ╔═╡ a09ab593-8efd-4799-8f0f-aa2ceb1e135d
@with_kw struct DiffWrapper{X <: RVec,Y <: RVec,N,S,T,C}
	x0 :: X
	fx0 :: Y
	#vars :: Tuple{Tuple{N,Symbolics.Num}, Symbolics.Num} = _get_vars( x0 )
	vars :: N = _get_vars(x0)
	x0_sym :: Vector{Symbolics.Num} = [vars[1]...]
	stamp :: S = CFDStamp(1,2, stepsize(x0))
	order :: Int = 1

	cache_type :: C = Dict{Vector{UInt8}, eltype(vec_typef(fx0, stamp))}
	tree :: T = build_tree( x0_sym, stamp, vars, order; x_type = vec_typex(x0,stamp), val_type = vec_typef(fx0, stamp), cache_type = cache_type )
end

# ╔═╡ 5253756b-1a03-4d8f-854a-b03c5f396dc9
val( dw :: DiffWrapper, indices; output_index = 1 ) = val( dw.tree, indices, dw.stamp ; output_index)

# ╔═╡ dfc93b7e-980d-4f35-9d72-16cee2c747c2
function gradient( dw :: DiffWrapper; output_index = 1 )
	@assert 1 <= dw.order <= 2 "Gradient retrieval only implemented for DiffWrapper of order 1 and 2."
	return gradient( dw :: DiffWrapper, Val(dw.order); output_index )
end

# ╔═╡ a728763c-8995-4d09-9878-1db83afc9194
function gradient( dw :: DiffWrapper, ::Val{1}; output_index = 1 )
	n = length( dw.x0 )
	g = Vector{eltype(dw.fx0)}(undef, n)
	for i = 1 : n
		g[i] = val( dw, [i,]; output_index )
	end
	return g
end

# ╔═╡ 2f50ce22-7a20-4b24-ac57-5f3c3918079e
function gradient( dw :: DiffWrapper, ::Val{2}; output_index = 2 )
	n = length( dw.x0 )
	g = Vector{eltype(dw.fx0)}(undef, n)
	zi = _zi( dw.stamp )	# index of stamp grid point is zero
	node = dw.tree.vals[1][zi]
	for i = 1 : n
		g[i] = val( node, [i,], dw.stamp; output_index )
	end
	return g
end

# ╔═╡ 634067d4-15db-4aa5-8b5a-b0794ee88a39
function hessian( dw :: DiffWrapper; output_index = 1)
	@assert dw.order == 2 "Hessian only implemented for DiffWrapper of 2."
	n = length(dw.x0)
	H = Matrix{eltype(dw.fx0)}(undef, n, n)
	# TODO make sure that dw.fx0 == vec_typef(fx0, stamp)
	for i = 1 : n
		for j = 1 : n
			H[i,j] = val( dw, [i, j]; output_index )
		end
	end
	return H
end

# ╔═╡ dbdceb0d-e4bb-406c-b3cf-defffcc1ccf9
function substitute_symbols!(dw :: DiffWrapper)
	substitute_symbols!( dw.tree, dw.x0, _stepsize( dw.stamp ), dw.vars )
end

# ╔═╡ 8c11b8e1-c0fb-48b8-b9f1-8679e3d0e596
function substitute_leaves!(dw :: DiffWrapper)
	substitute_leaves!( dw.tree, dw.x0, _stepsize( dw.stamp ), dw.vars )
end

# ╔═╡ 351d4199-28fa-4733-89d1-266ee16137c4
function _substitute_symbols!( node :: FDiffNode, x0, h, vars )
	x_vars, h_var = vars
	empty!(node.x)
	append!(node.x, Symbolics.value.(substitute.(node.x_sym, (
		Dict((x_vars[i] => x0[i] for i = eachindex(x0))...,
					h_var=>h),
	))))
end

# ╔═╡ 94d3947d-fb1c-4e8d-9b16-bdf14a411193
begin
	Trees.children( n :: FDiffNode{<:RVec} ) = nothing
	Trees.children( n :: FDiffNode ) = Iterators.flatten( n.vals )
end

# ╔═╡ 98051445-7ffd-4a92-af38-07507efab445
function jacobian( dw :: DiffWrapper )
	@assert 1 <= dw.order <= 2 "Gradient retrieval only implemented for DiffWrapper of order 1 and 2."
	k = isempty( dw.fx0 ) ? length(first(Trees.Leaves(dw.tree)).vals) : length(dw.fx0)
	return transpose( hcat( collect(gradient(dw; output_index = m) for m = 1 : k )...) )
end

# ╔═╡ 0be3107e-a3c9-46f5-837d-c7f6d2b2d595
function substitute_symbols!(root :: FDiffNode, x0, h, vars )
	_substitute_symbols!( Trees.PreOrderDFS( root ), x0, h, vars )
end

# ╔═╡ 42e6c63e-0525-4ec6-88b9-b077130603ff
function substitute_leaves!(root :: FDiffNode, x0, h, vars )
	_substitute_symbols!( Trees.Leaves(root), x0, h, vars )
end

# ╔═╡ df8bd829-5c28-46c1-b171-1b5dc0358c76
function collect_leave_sites( dw :: DiffWrapper )
	return [ node.x for node in Trees.Leaves( dw.tree ) ]
end

# ╔═╡ 7b725d45-f6fa-4bd8-99d7-fcad166320e5
function set_leave_values!(dw :: DiffWrapper, f :: Function )
	for node in Trees.Leaves(dw.tree)
		_set_node_values!(node, f)
	end
end

# ╔═╡ 9cb6592b-ba7d-4b46-9171-1816fd7dff69
function set_leave_values!(dw :: DiffWrapper, leave_vals :: AbstractVector )
	for (i,node) in enumerate(Trees.Leaves(dw.tree))
		empty!(node.vals)
		append!(node.vals, leave_vals[i] )
	end
end

# ╔═╡ 8a1443d4-0beb-4f78-8ab4-ed47adf45e8c
#=
begin
	func = x -> [ sum(x.^2); exp( sum( x ) ) ]
	x_t = ones(Float32, 3)
	fx_t = func( x_t )
	diff_wrapper = DiffWrapper(; x0 = x_t, fx0 = fx_t, stamp = CFDStamp(1,2,stepsize(fx_t)), order = 2 )

	#prepare_tree!(diff_wrapper, func)

	# phase I: get sites
	substitute_leaves!( diff_wrapper )
	all_sites = collect_leave_sites( diff_wrapper )

	# phase II: evaluation

	## find unique sites to avoid costly objective evaluations
	unique_sites, set_indices = unique_with_indices( [[x_t,]; all_sites] )
	## now unique_sites[set_indices] == [[x_t,]; all_sites]

	@assert unique_sites[set_indices][2:end] == all_sites

	## and x_t will be the first site in unique sites

	## call `func` on new sites
	unique_evals = [ [fx_t,]; func.( unique_sites[ 2 : end ] )]
	all_evals = unique_evals[set_indices[2:end]]

	# phase III: set the leave values and get jacobian
	set_leave_values!(diff_wrapper, all_evals)

	hessian(diff_wrapper; output_index = 1), jacobian(diff_wrapper)
end
=#

# ╔═╡ bf982874-9ec2-4aa2-87ed-27f3d26d5d51
function prepare_tree!( dw :: DiffWrapper, f :: Function )
	x0 = dw.x0
	vars = dw.vars
	h = _stepsize( dw.stamp )

	for node in Trees.Leaves( dw.tree )
		_substitute_symbols!(node, x0, h, vars)
		_set_node_values!( node, f )
	end
end

# ╔═╡ e0154be3-4c57-4ab9-9218-13ab3a4f44b4
#=
begin
	# I used this cell for fiddling around.

	f = x -> [ 2 * x[1]^2 + x[2]^3; sum( x ) ]
	# ∂1 f1 = 4 * x[1]
	# ∂2 f1 = 3 * x[2]^2

	#H1 = [
	#	4.0 0.0
	#	0.0 6x[2]
	# ]

	dw = DiffWrapper(; x0 = [1.0; 1.0], fx0 = rand(2) , stamp = BFDStamp(1,3,1e-3), order = 2)
	# to deactivate caching: `cache_type = Nothing`

	prepare_tree!(dw,f)
	# this would require 2 tree traversals:
	# substitute_leaves!(dw)
	# set_leave_values!(dw,f)

	H = hessian(dw)
	H, sum( ([ 4 0; 0 6 ] .- H).^2 )
end
=#

# ╔═╡ 90b565b5-0826-4674-836a-376430823cc5
# jacobian(dw)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

[compat]
BenchmarkTools = "~1.1.1"
Parameters = "~0.12.2"
StaticArrays = "~1.2.7"
Symbolics = "~1.4.2"
"""

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

[[AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "a71d224f61475b93c9e196e83c17c6ac4dedacfa"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.18"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
git-tree-sha1 = "c31ebabde28d102b602bada60ce8922c266d205b"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.1.1"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f53ca8d41e4753c41cdafa6ec5f7ce914b34be54"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "0.10.13"

[[Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

[[CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.31.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[CompositeTypes]]
git-tree-sha1 = "d5b014b216dc891e81fea299638e4c10c657b582"
uuid = "b152e2b5-7a66-4b01-a709-34e65c35f657"
version = "0.1.2"

[[ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "85d2d9e2524da988bffaf2a381864e20d2dae08d"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.2.1"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Distributions]]
deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "3889f646423ce91dd1055a76317e9a1d3a23fff1"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.11"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[DomainSets]]
deps = ["CompositeTypes", "IntervalSets", "LinearAlgebra", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "6cdd99d0b7b555f96f7cb05aa82067ee79e7aef4"
uuid = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
version = "0.5.2"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[DynamicPolynomials]]
deps = ["DataStructures", "Future", "LinearAlgebra", "MultivariatePolynomials", "MutableArithmetics", "Pkg", "Reexport", "Test"]
git-tree-sha1 = "5e47c4d652ea67652b7c5945c79c46472397d47f"
uuid = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
version = "0.3.18"

[[EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "8041575f021cba5a099a456b4163c9a08b566a02"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.1.0"

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8c8eac2af06ce35973c3eadb4ab3243076a408e7"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.1"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IntervalSets]]
deps = ["Dates", "EllipsisNotation", "Statistics"]
git-tree-sha1 = "3cc368af3f110a767ac786560045dceddfc16758"
uuid = "8197267c-284f-5f27-9208-e0e47529a953"
version = "0.5.3"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[LabelledArrays]]
deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
git-tree-sha1 = "5e38cfdd771c34821ade5515f782fe00865d60b3"
uuid = "2ee39098-c373-598a-b85f-a56591580800"
version = "1.6.2"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "LinearAlgebra"]
git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.2.5"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.6"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MultivariatePolynomials]]
deps = ["DataStructures", "LinearAlgebra", "MutableArithmetics"]
git-tree-sha1 = "45c9940cec79dedcdccc73cc6dd09ea8b8ab142c"
uuid = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
version = "0.3.18"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.1"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.1"

[[RecursiveArrayTools]]
deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "a4bd5d7c4bf7effc1e7ab75d503928082f63bd71"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.16.0"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[RuntimeGeneratedFunctions]]
deps = ["ExprTools", "SHA", "Serialization"]
git-tree-sha1 = "5975a4f824533fa4240f40d86f1060b9fc80d7cc"
uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
version = "0.5.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SciMLBase]]
deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
git-tree-sha1 = "f0bf114650476709dd04e690ab2e36d88368955e"
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
version = "1.18.2"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "d5640fc570fb1b6c54512f0bd3853866bd298b3e"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.0"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.5.1"

[[Static]]
deps = ["IfElse"]
git-tree-sha1 = "62701892d172a2fa41a1f829f66d2b0db94a9a63"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.3.0"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "1b9a0f17ee0adde9e538227de093467348992397"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.7"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.8"

[[StatsFuns]]
deps = ["LogExpFunctions", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "30cd8c360c54081f806b1ee14d2eecbef3c04c49"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.8"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[SymbolicUtils]]
deps = ["AbstractTrees", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DynamicPolynomials", "IfElse", "LabelledArrays", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "TimerOutputs"]
git-tree-sha1 = "17fecd52e9a82ca52bfc9d28a2d31b33458a6ce0"
uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
version = "0.13.1"

[[Symbolics]]
deps = ["ConstructionBase", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "IfElse", "Latexify", "Libdl", "LinearAlgebra", "MacroTools", "NaNMath", "RecipesBase", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils", "TreeViews"]
git-tree-sha1 = "dae26a27018d0cad7efd585a9a0012c6a2752a88"
uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
version = "1.4.2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.4.4"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TimerOutputs]]
deps = ["ExprTools", "Printf"]
git-tree-sha1 = "209a8326c4f955e2442c07b56029e88bb48299c7"
uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
version = "0.5.12"

[[TreeViews]]
deps = ["Test"]
git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
version = "0.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

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

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.1"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─7b914a4f-89b4-4284-b8c6-d629157b9a7c
# ╠═e7d76f0e-eea5-11eb-1bdf-9196c513ff0a
# ╠═44f2f27b-2ac6-4df1-b222-4cfa37f41985
# ╟─10b9b062-450e-43a2-9c7a-2fe6c7121ab7
# ╟─bd83fc2a-d7bd-489e-83f3-ef5b327d5cca
# ╠═a7d81ac5-0124-441a-8eae-f83cee2b81c3
# ╟─fdcf00a4-836e-4404-ba1c-6058ae85d9ce
# ╠═aeac724f-b87a-4f37-aa57-893d1d131f6d
# ╟─07e698bb-fac0-4951-85dc-5adb208cf5dd
# ╠═9462b310-0842-4373-ac2c-f2eaef5e6e7f
# ╠═c19325b6-7516-4acf-a5f1-fe7084bc113d
# ╠═2426f92a-21f3-4455-aeb7-fe6674aafe02
# ╠═4ebb80a7-539d-49a0-981f-166f0bbf3964
# ╟─a5e02956-c0f9-44a6-80da-abfb497bcfb9
# ╠═09d78bbc-f568-459a-8bd0-05de84b04af2
# ╟─04861d1b-38d9-4ba7-9fde-b10150f39280
# ╟─20907685-2bb6-49d5-9ff0-a175435fa591
# ╟─65d0f18a-6da9-4836-bf4c-8ff10e6a265f
# ╠═c3fa417a-898f-4273-a75c-974a42e80a79
# ╟─6cdc6c1e-790e-4d2a-984c-36f5f0c98f63
# ╠═94d3947d-fb1c-4e8d-9b16-bdf14a411193
# ╟─dc4d3c61-f2ec-4fb1-b84b-2c1957da3a7e
# ╠═c2cb0dd1-e36d-44ff-aa59-6365a15e7ab6
# ╠═aabc4fed-e62f-4668-9984-c9b72cad15e9
# ╟─3f67de97-2a12-4ca7-82bd-16851f958989
# ╠═18e51761-98ae-4814-b575-3bb4a40a06cb
# ╠═d21fed2f-f61a-42d6-aed3-430bf73fd170
# ╟─a5565a82-2b6d-4c43-909c-9bfeca9890cb
# ╠═17347e0e-08ee-4a29-83b3-f337e397f49b
# ╟─1db29942-3ff0-49e8-b5c6-2480eadd608f
# ╠═af650705-0636-4728-a643-ab658d362290
# ╠═7c091ff3-7d34-4ad8-b959-fe601397f091
# ╠═a09ab593-8efd-4799-8f0f-aa2ceb1e135d
# ╟─59211507-8901-4f4c-b2ab-c5d4fb950bf1
# ╟─eb61f081-bd96-4a25-b2ba-9a3bc9e87934
# ╠═5253756b-1a03-4d8f-854a-b03c5f396dc9
# ╟─f77f8375-5b8f-4a3e-b571-28ebccb204c4
# ╟─fe8171dc-7efd-4acb-b1f9-0da5f1c6dc9d
# ╠═dfc93b7e-980d-4f35-9d72-16cee2c747c2
# ╠═a728763c-8995-4d09-9878-1db83afc9194
# ╠═2f50ce22-7a20-4b24-ac57-5f3c3918079e
# ╠═98051445-7ffd-4a92-af38-07507efab445
# ╠═634067d4-15db-4aa5-8b5a-b0794ee88a39
# ╟─5dd33903-f4a4-455c-8348-c3c98e5f6404
# ╠═e0154be3-4c57-4ab9-9218-13ab3a4f44b4
# ╠═90b565b5-0826-4674-836a-376430823cc5
# ╟─8ad29cb3-9bb5-48d6-b0a6-ceb46e30293e
# ╠═351d4199-28fa-4733-89d1-266ee16137c4
# ╠═fef9850c-3cbd-47f8-a786-e420f01d9141
# ╠═0be3107e-a3c9-46f5-837d-c7f6d2b2d595
# ╠═42e6c63e-0525-4ec6-88b9-b077130603ff
# ╠═dbdceb0d-e4bb-406c-b3cf-defffcc1ccf9
# ╠═8c11b8e1-c0fb-48b8-b9f1-8679e3d0e596
# ╠═df8bd829-5c28-46c1-b171-1b5dc0358c76
# ╠═ac3d101b-c9d1-4fdb-80a4-645464298ba3
# ╠═7b725d45-f6fa-4bd8-99d7-fcad166320e5
# ╠═9cb6592b-ba7d-4b46-9171-1816fd7dff69
# ╠═bf982874-9ec2-4aa2-87ed-27f3d26d5d51
# ╟─5d586923-a32e-4d16-a0a6-56bbc73680f2
# ╠═3feed8b1-1338-47c7-960f-892481f96a5c
# ╠═8a1443d4-0beb-4f78-8ab4-ed47adf45e8c
# ╟─d41308c0-b7f1-4ed9-b27e-7c6fe87080e6
# ╠═6c6e63ae-28ce-4c1c-bb60-7e7826ffefbb
# ╟─f294835d-ee9a-42b1-b2b7-5816d65e0f18
# ╟─e9af48fd-ac05-4358-8a33-c14008ca1cb9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002
end


text:
module Constraints

# Imports
using Dictionaries
using CompositionalNetworks
using ConstraintDomains

# Exports
export Constraint
export concept, error_f
export args_length, params_length
export symmetries
export usual_constraints, usual_symmetries
export learn_from_icn

# Includes internals
include("utils.jl")
include("constraint.jl")

# Includes learned errors from ICN
foreach(include, readdir(joinpath(dirname(pathof(Constraints)),"compositions"), join = true))

# Includes constraints
include("constraints/all_different.jl")
include("constraints/dist_different.jl")
include("constraints/equals.jl")
include("constraints/less.jl")
include("constraints/ordered.jl")
include("constraints/sequential_tasks.jl")

"""
    usual_constraints::Dict
Dictionary that contains all the usual constraints defined in Constraint.jl.
"""
const usual_constraints = Dict(
    :all_different => all_different,
    :all_equal => all_equal,
    :all_equal_param => all_equal_param,
    :always_true => Constraint(),
    :dist_different => dist_different,
    :eq => eq,
    :less_than_param => less_than_param,
    :minus_equal_param => minus_equal_param,
    :ordered => ordered,
    :sequential_tasks => sequential_tasks,
    :sum_equal_param => sum_equal_param,
)

# include learn script
include("learn.jl")

end

ext: jl


text:
"""
    appear(s::Symbol)

Appear can be used inside an [`Action`](@ref) to make an [`Object`](@ref) or an entire [`Object`](@ref) (including it's objects) to appear.

# Example
```julia
house = Object(101:200, (args...)->house_of_nicholas())
act!(house, Action(1:20, appear(:fade)))
act!(house, Action(81:100, disappear(:fade)))
```
In this case the `house_of_nicholas` will fade in during the first 20 frames
of the [`Object`](@ref) so `101-120`.

# Arguments
- `s::Symbol`: the symbol defines the animation of appearance
    The only symbols that are currently supported are:
    - `:fade_line_width` which increases the line width up to the default value
       or the value specified by [`setline`](@ref)
    - `:fade` which increases the opcacity up to the default value
       or the value specified by [`setopacity`](@ref)
    - `:scale` which increases the scale up to the default value `1`
       or the value specified by [`scale`](@ref)
    - `:draw_text` which only works for [`text`](@ref) and lets it appear from left to right.

For a layer only `appear(:fade)` is supported
"""
function appear(s::Symbol)
    (video, object, action, rel_frame) -> _appear(video, object, action, rel_frame, Val(s))
end

function _appear(video, object, action, rel_frame, symbol::Val{:fade_line_width})
    t = get_interpolation(action, rel_frame)
    object.current_setting.mul_line_width = t
end

function _appear(video, object, action, rel_frame, symbol::Val{:fade})
    t = get_interpolation(action, rel_frame)
    object.current_setting.mul_opacity = t
end

function _appear(video, layer::Layer, action, rel_frame, symbol::Val{:fade})
    t = get_interpolation(action, rel_frame)
    layer.current_setting.opacity = t
end

function _appear(video, object, action, rel_frame, symbol::Val{:scale})
    t = get_interpolation(action, rel_frame)
    object.current_setting.mul_scale = t
end

function _appear(video, object, action, rel_frame, symbol::Val{:draw_text})
    t = get_interpolation(action, rel_frame)
    object.opts[:draw_text_t] = t
end

"""
    disappear(s::Symbol)

Disappear can be used inside an [`Action`](@ref) to make an [`Object`](@ref) or an entire [`Layer`](@ref) (including it's objects) to disappear.

# Example
```julia
house = Object(101:200, (args...)->house_of_nicholas())
act!(house, Action(1:20, appear(:fade)))
act!(house, Action(81:100, disappear(:fade)))
```
In this case the `house_of_nicholas` will fade out during the last 20 frames
of the [`Object`](@ref) so `181-200`.

# Arguments
- `s::Symbol`: the symbol defines the animation of disappearance
    The only symbols that are currently supported are:
    - `:fade_line_width` which decreases the line width down to `0`
    - `:fade` which decreases the opacity down to `0`
    - `:scale` which decreases the scale down to `0`
    - `:draw_text` which only works for text and let the text disappear from right to left.

For a layer only `disappear(:fade)` is supported
"""
function disappear(s::Symbol)
    (video, object, action, rel_frame) ->
        _disappear(video, object, action, rel_frame, Val(s))
end

function _disappear(video, object, action, rel_frame, symbol::Val{:fade_line_width})
    t = get_interpolation(action, rel_frame)
    object.current_setting.mul_line_width = 1 - t
end

function _disappear(video, object, action, rel_frame, symbol::Val{:fade})
    t = get_interpolation(action, rel_frame)
    object.current_setting.mul_opacity = 1 - t
end

# for layer only disappear(:fade) is possible
function _disappear(video, layer::Layer, action, rel_frame, symbol::Val{:fade})
    t = get_interpolation(action, rel_frame)
    layer.current_setting.opacity = 1 - t
end

function _disappear(video, object, action, rel_frame, symbol::Val{:scale})
    t = get_interpolation(action, rel_frame)
    object.current_setting.mul_scale = 1 - t
end

function _disappear(video, object, action, rel_frame, symbol::Val{:draw_text})
    t = get_interpolation(action, rel_frame)
    object.opts[:draw_text_t] = 1 - t
end

"""
    translate()

Translate an [`Object`](@ref) or a [`Layer`](@ref) using an [`Action`](@ref) and an Animation defined
with Animations.jl.

If you're used to working with Animations.jl this should feel quite natural.
Instead of defining each movement in its own action it's possible to define it in one
by using an Animation.

# Example
```julia
using Javis, Animations

function ground(args...)
    background("black")
    sethue("white")
end

video = Video(500, 500)
circle_anim = Animation(
    [0.0, 0.3, 0.6, 1.0], # must go from 0 to 1
    # the circle will move from the origin to `Point(150, 0)` then `Point(150, 150)`
    # and back to the origin `O`.
    [O, Point(150, 0), Point(150, 150), O],
    [sineio(), polyin(5), expin(8)],
)

Background(1:150, ground)
obj = Object((args...)->circle(O, 25, :fill))
act!(obj, Action(1:150, circle_anim, translate()))

render(video)
```

Here `circle_anim` defines the movement of the circle. The most important part is that the
time in animations has to be from `0.0` to `1.0`.

This notation uses the Animations.jl library very explicitly. It's also possible to do the
same with:

```julia
obj = Object((args...)->circle(O, 25, :fill))
act!(obj, Action(1:50, sineio(), anim_translate(150, 0)))
act!(obj, Action(51:100, polyin(2), anim_translate(0, 150)))
act!(obj, Action(101:150, expin(8), anim_translate(-150, -150)))
```

which uses the `Action` syntax three times and only uses easing functions instead of
specifying the `Animation` directly. Have a look at [`anim_translate`](@ref) for details.
"""
function Luxor.translate()
    (video, object, action, rel_frame) -> _translate(video, object, action, rel_frame)
end

function _translate(video, object, action, rel_frame)
    p = get_interpolation(action, rel_frame)
    Luxor.translate(p)
end

function _translate(video, layer::Layer, action, rel_frame)
    p = get_interpolation(action, rel_frame)
    layer.position = p
end

"""
    rotate()

Rotate an [`Object`](@ref) or a [`Layer`](@ref) using an [`Action`](@ref) and an Animation defined
with Animations.jl.

If you're used to working with Animations.jl this should feel quite natural.
Instead of defining each movement in its own action it's possible to define it in one
by using an Animation.

# Example
```julia
using Javis, Animations

# define ground function here

video = Video(500, 500)
translate_anim = Animation(
    [0, 1], # must go from 0 to 1
    [O, Point(150, 0)],
    [sineio()],
)

translate_back_anim = Animation(
    [0, 1], # must go from 0 to 1
    [O, Point(-150, 0)],
    [sineio()],
)

rotate_anim = Animation(
    [0, 1], # must go from 0 to 1
    [0, 2π],
    [linear()],
)

Background(1:150, ground)
ball = Object((args...) -> circle(O, 25, :fill))
act!(ball, Action(1:10, sineio(), scale()))
act!(ball, Action(11:50, translate_anim, translate()))
act!(ball, Action(51:100, rotate_anim, rotate_around(Point(-150, 0))))
act!(ball, Action(101:140, translate_back_anim, translate()))
act!(ball, Action(141:150, rev(sineio()), scale()))

render(video)
```

which uses the `Action` syntax five times with both easing functions directly and animation objects.
The `rev(sineio())` creates an `Animation` which goes from `1.0` to `0.0`.
"""
function Luxor.rotate()
    (video, object, action, rel_frame) -> _rotate(video, object, action, rel_frame)
end

function _rotate(video, object, action, rel_frame)
    a = get_interpolation(action, rel_frame)
    Luxor.rotate(a)
end

function _rotate(video, layer::Layer, action, rel_frame)
    a = get_interpolation(action, rel_frame)
    layer.current_setting.rotation_angle = a
end

"""
    rotate_around(p)

Rotate an [`Object`](@ref) or a [`Layer`](@ref) using an [`Action`](@ref) and an Animation defined
with Animations.jl around a point `p`. For [`rotate`](@ref) it rotates around the current origin.

An example can be seen in [`rotate`](@ref).

# Arguments
- `p`: the point to rotate around
"""
function rotate_around(p)
    (video, object, action, rel_frame) ->
        _rotate_around(video, object, action, rel_frame, get_position(p))
end

function _rotate_around(video, object, action, rel_frame, p)
    # p should be global so without the translated start_pos
    p = p - object.start_pos
    # save the radius and start angle
    r = get(action.defs, :rotate_radius, distance(p, O))
    pnormed = get(action.defs, :pnormed, p / r)
    if !haskey(action.defs, :rotate_radius)
        action.defs[:rotate_radius] = r
        action.defs[:pnormed] = pnormed
    end
    i = get_interpolation(action, rel_frame)
    Luxor.translate(p)
    Luxor.rotate(i)
    Luxor.translate(-r * pnormed)
end

function _rotate_around(video, layer::Layer, action, rel_frame, p)
    # p should be global so without the translated start_pos
    p = p - layer.position
    # save the radius and start angle
    r = get(action.defs, :rotate_radius, distance(p, O))
    pnormed = get(action.defs, :pnormed, p / r)
    if !haskey(action.defs, :rotate_radius)
        action.defs[:rotate_radius] = r
        action.defs[:pnormed] = pnormed
    end
    i = get_interpolation(action, rel_frame)
    p1 = -r * pnormed
    push!(
        layer.current_setting.misc,
        :rotate_around => true,
        :angle => i,
        :translate => p,
        :translate_back => p1,
    )
end

"""
    scale()

Scale a function defined inside an [`Action`](@ref) using an Animation defined
with Animations.jl.

An example can be seen in [`rotate`](@ref).
"""
function scale()
    (video, object, action, rel_frame) -> _scale(video, object, action, rel_frame)
end

function _scale(video, object, action, rel_frame)
    s = get_interpolation(action, rel_frame)
    scaleto(s)
end

function _scale(video, layer::Layer, action, rel_frame)
    s = get_interpolation(action, rel_frame)
    layer.current_setting.scale = s
end

"""
    sethue()

Set the color of an [`Object`](@ref) using an [`Action`](@ref) and an Animation defined
with Animations.jl.

# Example
A possible animation would look like this:
```julia
color_anim = Animation(
    [0, 0.5, 1], # must go from 0 to 1
    [
        Lab(colorant"red"),
        Lab(colorant"cyan"),
        Lab(colorant"black"),
    ],
    [sineio(), sineio()],
)
```

An example on how to integrate this into an [`Action`](@ref) can be seen in [`rotate`](@ref).
Where this would be a valid Action: `Action(1:150, color_anim, sethue())`.
"""
function Luxor.sethue()
    (video, object, action, rel_frame) -> _sethue(video, object, action, rel_frame)
end

function _sethue(video, object, action, rel_frame)
    color = get_interpolation(action, rel_frame)
    Luxor.sethue(color)
end

"""
    setopacity()

Set the color of an [`Object`](@ref) or a [`Layer`](@ref) using an [`Action`](@ref) and an Animation defined
with Animations.jl.

# Example
A possible animation would look like this:
```julia
opacity_anim = Animation(
    [0, 0.5, 1], # must go from 0 to 1
    [
        0.0,
        0.3,
        0.7,
    ],
    [sineio(), sineio()],
)
```

An example on how to integrate this into an [`Action`](@ref) can be seen in [`rotate`](@ref).
Where this would be a valid Action: `Action(1:150, opacity_anim, setopacity())`.
"""
function setopacity()
    (video, object, action, rel_frame) -> _setopacity(video, object, action, rel_frame)
end

function _setopacity(video, object, action, rel_frame)
    opacity = get_interpolation(action, rel_frame)
    setopacity(opacity)
end

function _setopacity(video, layer::Layer, action, rel_frame)
    opacity = get_interpolation(action, rel_frame)
    layer.current_setting.opacity = opacity
end

"""
    follow_path(points::Vector{Point}; closed=true)

Can be applied inside an action such that the parent object follows a path.
It takes a vector of points which can be created as an example by calling
`circle(O, 50)`. Notice that the object is set to `:none`, the default.

# Example
```julia
Action(1:150, follow_path(star(O, 300)))
```

# Arguments
- `points::Vector{Point}` - the vector of points the object should follow

# Keywords
- `closed::Bool` default: true, sets whether the path is a closed path as for example when
    using a circle, ellipse or any polygon. For a bezier path it should be set to false.
"""
function follow_path(points::Vector{Point}; closed = true)
    (video, object, action, rel_frame) ->
        _follow_path(video, object, action, rel_frame, points; closed = closed)
end

function _follow_path(video, object, action, rel_frame, points; closed = closed)
    isfirstframe = rel_frame == first(get_frames(action))
    t = get_interpolation(action, rel_frame)
    # if not closed it should be always between 0 and 1
    if !closed
        t = clamp(t, 0.0, 1.0)
    end
    # if t is discrete and not 0.0 take the last point or first if closed
    if !isfirstframe && isapprox_discrete(t)
        if closed
            translate(points[1])
        else
            translate(points[end])
        end
        return
    end
    # get only the frobjectal part to be between 0 and 1
    t -= floor(t)
    if isfirstframe
        # compute the distances only once for performance reasons
        action.defs[:p_dist] = polydistances(points, closed = closed)
    end
    if isapprox(t, 0.0, atol = 1e-4)
        translate(points[1])
        return
    end
    pdist = action.defs[:p_dist]

    overshootpoint = get_polypoint_at(points, t; pdist = pdist)
    translate(overshootpoint)
end

"""
    change(s::Symbol, [val(s)])

Changes the keyword `s` of the parent [`Object`](@ref) from `vals[1]` to `vals[2]`
in an animated way if vals is given as a `Pair` otherwise it sets the keyword `s` to `val`.

# Arguments
- `s::Symbol` Change the keyword with the name `s`
- `vals::Pair` If vals is given i.e `0 => 25` it will be animated from 0 to 25.
    - The default is to use `0 => 1` or use the value given by the animation
    defined in the [`Action`](@ref)

# Example
```julia
Background(1:100, ground)
obj = Object((args...; radius = 25, color="red") -> object(O, radius, color), Point(100, 0))
act!(obj, Action(1:50, change(:radius, 25 => 0)))
act!(Action(51:100, change(:radius, 0 => 25)))
act!(Action(51:100, change(:color, "blue")))
```
"""
function change(s::Symbol, vals::Pair)
    (video, object, action, rel_frame) -> _change(video, object, action, rel_frame, s, vals)
end

function change(s::Symbol, val)
    (video, object, action, rel_frame) -> _change(video, object, action, rel_frame, s, val)
end

function change(s::Symbol)
    (video, object, action, rel_frame) -> _change(video, object, action, rel_frame, s)
end

function _change(video, object, action, rel_frame, s, vals::Pair)
    t = get_interpolation(action, rel_frame)
    # use the linear animation power by the Animation package to interpolate colors as well
    lin_anim = Animation([0.0, 1.0], interpolateable([vals...]))
    object.change_keywords[s] = at(lin_anim, t)
end

function _change(video, object, action, rel_frame, s, val)
    object.change_keywords[s] = val
end

function _change(video, object, action, rel_frame, s)
    val = get_interpolation(action, rel_frame)
    object.change_keywords[s] = val
end


text:
module Interact

using Reactive, Compat

import Base: mimewritable, writemime
export signal, Widget, InputWidget

# A widget
abstract Widget

# A widget that gives out a signal of type T
abstract InputWidget{T}  <: Widget

signal(w::InputWidget) = w.signal
signal(x::Signal) = x

function statedict(w)
    msg = Dict()
    attrs = @compat fieldnames(w)
    for n in attrs
        if n in [:signal, :label]
            continue
        end
        msg[n] = getfield(w, n)
    end
    msg
end

# Convert e.g. JSON values into Julia values
parse_msg{T <: Number}(::InputWidget{T}, v::AbstractString) = parse(T, v)
parse_msg(::InputWidget{Bool}, v::Number) = v != 0
parse_msg{T}(::InputWidget{T}, v) = convert(T, v)

"""
This should be overrided by communication providers
"""
function update_view end

function error_handler(sig, value, err)
    Reactive.print_error(sig, value, err, open("/tmp/Interact.log", "a"))
end

function recv_msg{T}(widget ::InputWidget{T}, value)
    # Hand-off received value to the signal graph
    parsed = parse_msg(widget, value)
    #println(STDERR, signal(widget))
    push!(signal(widget), parsed)
    widget.value = parsed
    if value != parsed
        update_view(widget)
    end
end

uuid4() = string(Base.Random.uuid4())

const id_to_widget = Dict{AbstractString, InputWidget}()
const widget_to_id = Dict{InputWidget, AbstractString}()

function register_widget(w::InputWidget)
    if haskey(widget_to_id, w)
        return widget_to_id[w]
    else
        id = uuid4()
        widget_to_id[w] = id
        id_to_widget[id] = w
        return id
    end
end

function get_widget(id::AbstractString)
    if haskey(id_to_widget, id)
        return id_to_widget[id]
    else
        warn("Widget with id $(id) does not exist.")
    end
end

include("widgets.jl")
include("compose.jl")
include("manipulate.jl")

const ijulia_setup_path = joinpath(dirname(@__FILE__), "IJulia", "setup.jl")

function __init__()
    if isdefined(Main, :IJulia)
        include(ijulia_setup_path)
    end
end

"""
Interact.jl allows you to use interactive widgets such as sliders, dropdowns and checkboxes to play with your Julia code.

Basic widgets:

- `slider(1:10)` creates a slider widget with the specified range of values
- `checkbox(false)` creates a check box
- `togglebutton(false)` creates a toggle button
- `button()` creates a button

Option (selection) widgets:

- `dropdown(["one", "two", "three"])` creates a drop-down widget
- `togglebuttons(["one", "two", "three"])` creates a toggle button group
- `radiobuttons(["one", "two", "three"])` creates a radio button group

Text boxes:

- `textbox(value)` creates a box accepting string or numeric text input
- `textarea` creates a box for extended text input

Output widgets:

- `latex`
- `progress`
"""
Interact

end # module


text:
"""
Depth First Walk (Tree Postwalk) on expressions, mutates expression in-place.
"""
function df_walk!(f, e, f_args...; skip=Vector{Symbol}(), skip_call=false)
    if !(e isa Expr) || e.head ∈ skip
        return f(e, f_args...)
    end
    #println("walking on", e)
    start = 1
    # skip walking on function names
    if skip_call && isexpr(e, :call)
        start = 2
    end
    e.args[start:end] = e.args[start:end] .|> x ->
        df_walk!(f, x, f_args...; skip=skip, skip_call=skip_call)
    return f(e, f_args...)
end

"""
Depth First Walk (Tree Postwalk) on expressions. Does not mutate expressions.
"""
function df_walk(f, e, f_args...; skip=Vector{Symbol}(), skip_call=false)
    if !(e isa Expr) || e.head ∈ skip
        return f(e, f_args...)
    end
    start = 1
    # skip walking on function names
    if skip_call && isexpr(e, :call)
        start = 2
    end

    ne = copy(e)
    ne.args[start:end] = ne.args[start:end] .|> x ->
        df_walk(f, x, f_args...; skip=skip, skip_call=skip_call)
    return f(ne, f_args...)
end

## Breadth First Walk on expressions

"""
Breadth First Walk (Tree Prewalk) on expressions mutates expression in-place.
"""
function bf_walk!(f, e, f_args...; skip=Vector{Symbol}(), skip_call=false)
    if !(e isa Expr) || e.head ∈ skip
        return f(e, f_args...)
    end
    e = f(e, f_args...)
    if !(e isa Expr) return e end
    start = 1
    # skip walking on function names
    if skip_call && isexpr(e, :call)
        start = 2
    end
    e.args[start:end] = e.args[start:end] .|> x ->
        bf_walk!(f, x, f_args...; skip=skip, skip_call=skip_call)
    return e
end

"""
Breadth First Walk (Tree Prewalk) on expressions. Does not mutate expressions.
"""
function bf_walk(f, e, f_args...; skip=Vector{Symbol}(), skip_call=false)
    if !(e isa Expr) || e.head ∈ skip
        return f(e, f_args...)
    end
    ne = copy(e)
    ne = f(e, f_args...)
    if !(ne isa Expr) return ne end
    start = 1
    # skip walking on function names
    if skip_call && isexpr(ne, :call)
        start = 2
    end
    ne.args[start:end] = ne.args[start:end] .|> x ->
        bf_walk(f, x, f_args...; skip=skip, skip_call=skip_call)
    return ne
end


text:
module LJH
export LJHGroup, LJHFile, LJH3File, ljhopen
using Nullables

"    ljh_get_header_dict(io::IO)
Return a Dict{String,String} mapping entries in the header. "
function ljh_get_header_dict(io::IO)
    headerdict = Dict()
    while true
        line = readline(io)
        if startswith(line, "#End of Header")
            break
        elseif eof(io)
          error("eof while reading header")
        elseif startswith(line, "#")
            continue
        else
            splitline=split(line,":")
            if length(splitline) != 2
              continue
            end
            a,b = splitline
            headerdict[strip(a)]=strip(b)
        end
    end
    headerdict
end

mutable struct LJHFile{VersionInt, FrameTime, PretrigNSamples, NRow, T<:IO}
    filename         ::String           # filename
    io               ::IO               # IO to read from LJH file
    headerdict       ::Dict             # LJH file header data
    datastartpos     ::Int
    record_nsamples  ::Int64            # number of sample per record
    channum          ::Int16            # channel number
    column           ::Int16            # column number
    row              ::Int16            # row number
    num_columns      ::Int16            # number of rows
    inverted_data    ::Bool             # data values inverted (downward-going pulses)?
end
struct LJHRecord{FrameTime, PretrigNSamples, NRow}
    data::Vector{UInt16}
    rowcount::Int64
    timestamp_usec::Int64
end
data(r::LJHRecord) = r.data
rowcount(r::LJHRecord) = r.rowcount
timestamp_usec(r::LJHRecord) = r.timestamp_usec
num_rows(r::LJHRecord{FrameTime, PretrigNSamples, NRow}) where {FrameTime, PretrigNSamples, NRow} = NRow
frametime(r::LJHRecord{FrameTime, PretrigNSamples, NRow}) where {FrameTime, PretrigNSamples, NRow} = FrameTime
frame1index(r::LJHRecord{FrameTime, PretrigNSamples, NRow}) where {FrameTime, PretrigNSamples, NRow} = div(rowcount(r), num_rows(r))
first_rising_sample(r::LJHRecord{FrameTime, PretrigNSamples, NRow}) where {FrameTime, PretrigNSamples, NRow} = PretrigNSamples
frameperiod(r::LJHRecord) = frametime(r)

Base.length(r::LJHRecord) = length(r.data)
import Base: ==
==(a::LJHRecord, b::LJHRecord) = a.data==b.data && a.rowcount == b.rowcount && a.timestamp_usec == b.timestamp_usec

@enum LJHVERSION LJH_20 LJH_21 LJH_22
VERSIONS = Dict(v"2.0.0"=>LJH_20, v"2.1.0"=>LJH_21, v"2.2.0"=>LJH_22,
v"2.2.1"=>LJH_22)

LJHFile(fname::String; inverted=false) = LJHFile(fname, open(fname,"r"); inverted=inverted)
function LJHFile(fname::String, io::IO; inverted=false)
    headerdict = ljh_get_header_dict(seekstart(io))
    datastartpos = position(io)
    # ioend = position(seekend(io))
    ioend = stat(io).size
    # seek(io,datastartpos)
    version = VersionNumber(headerdict["Save File Format Version"])
    version in keys(VERSIONS) || error("$fname has version $version, which is not in VERSIONS $VERSIONS")
    versionint = VERSIONS[version]
    record_nsamples = parse(Int,headerdict["Total Samples"])
    NRow = parse(Int16, get(headerdict,"Number of rows","1"))
    num_columns = parse(Int16, get(headerdict,"Number of columns","0"))
    FrameTime = parse(Float64,headerdict["Timebase"]) # frametime
    PretrigNSamples = parse(Int16,headerdict["Presamples"]) # pretrig_nsamples
    LJHFile{versionint, FrameTime, PretrigNSamples, NRow, typeof(io)}(
        fname, # filename
        io, # IO
        headerdict, #headerdict
        datastartpos, # datastartpos
        parse(Int16,headerdict["Total Samples"]), # record_nsamples
        round(Int16,parse(Float64,headerdict["Channel"])), # channel number
        parse(Int16,get(headerdict,"Column number (from 0-$(num_columns-1) inclusive)","0")),   #column, defaults to zero if not in file
        parse(Int16,get(headerdict,"Row number (from 0-$(NRow-1) inclusive)","0")), #row, defaults to zero if not in file
        num_columns,
        inverted)
end
LJHFile(f::LJHFile) = f

"header_nbytes(f::LJHFile)
    Return number of bytes in the header of each record in LJH file, based on version number"
header_nbytes(f::LJHFile{LJH_20}) = 6
header_nbytes(f::LJHFile{LJH_21}) = 6
header_nbytes(f::LJHFile{LJH_22}) = 16

"record_nbytes(f::LJHFile)
    Return number of bytes per record in LJH file, based on version number"
record_nbytes(f::LJHFile) = header_nbytes(f)+2*f.record_nsamples

"    ljh_number_of_records(f::LJHFile)
Return the number of complete records currently available to read from `f`."
function ljh_number_of_records(f::LJHFile)
    nbytes = stat(f.io).size-f.datastartpos
    number_of_records = div(nbytes, record_nbytes(f))
end

"    progresssize(f::Union{LJHFile, LJH3File)
Return the value that `progressposition` will return when the file is at `seekend(f.io)`.
Used to allow `ProgressBar` to work for both `LJHFile` and `LJH3File`."
progresssize(f::LJHFile) = length(f)
"    progressposition(f::LJHFile)
Return a value to pass to `update!(p::ProgressBar,position)`. Used to allow
`ProgressBar` to work for both `LJHFile` and `LJH3File`."
progressposition(f::LJHFile) = div(position(f.io)-f.datastartpos,record_nbytes(f))

# support for ljhfile[1:7] syntax
seekto(f::LJHFile, i::Int) = seek(f.io,f.datastartpos+(i-1)*record_nbytes(f))
Base.getindex(f::LJHFile,indexes::AbstractVector)=LJHGroupSlice(LJHGroup(f), indexes)
function Base.getindex(f::LJHFile,index::Int)
    seekto(f, index)
    pop!(f)
end
Base.pop!(f::LJHFile) = get(tryread(f)) # use get to raise error if Nullable is null.
_readrecord(f::LJHFile,i) = pop!(f)

Base.size(f::LJHFile) = (ljh_number_of_records(f),)
Base.length(f::LJHFile) = ljh_number_of_records(f)
Base.lastindex(f::LJHFile) = ljh_number_of_records(f)
function Base.eltype(f::LJHFile{V, FrameTime, PretrigNSamples, NRow}) where {V, FrameTime, PretrigNSamples, NRow}
     LJHRecord{FrameTime, PretrigNSamples, NRow}
end
# access as iterator
function Base.iterate(f::LJHFile)
    seekto(f, 1)
    return Base.iterate(f, 0)
end
function Base.iterate(f::LJHFile, j)
    j≥length(f) && return
    return pop!(f), j+1
end
filename(f::LJHFile) = f.filename
record_nsamples(f::LJHFile) = f.record_nsamples
function pretrig_nsamples(f::LJHFile{VersionInt, FrameTime, PretrigNSamples, NRow}) where {VersionInt, FrameTime, PretrigNSamples, NRow}
    convert(Int,PretrigNSamples)
end
channel(f::LJHFile) = f.channum
row(f::LJHFile) = f.row
column(f::LJHFile) = f.column
function frameperiod(f::LJHFile{VersionInt, FrameTime, PretrigNSamples, NRow}) where {VersionInt, FrameTime, PretrigNSamples, NRow}
    convert(Float64, FrameTime)
end
frametime(f::LJHFile) = frameperiod(f)
function num_rows(f::LJHFile{VersionInt, FrameTime, PretrigNSamples, NRow}) where {VersionInt, FrameTime, PretrigNSamples, NRow}
    convert(Int, NRow)
end
function Base.show(io::IO, g::LJHFile)
    print(io, "LJHFile $(filename(g))\n")
    inverted = ""
    if g.inverted_data
        inverted = "inverted "
    end
    print(io, "$(length(g)) $(inverted)records\n")
    print(io, "record_nsamples $(record_nsamples(g)), pretrig_nsamples $(pretrig_nsamples(g)).\n")
    print(io, "Channel $(channel(g)), row $(row(g)), column $(column(g)), frametime $(frametime(g)) s.\n")
end

# open and close
Base.close(f::LJHFile) = close(f.io)

"""
    tryread{T}(f::LJHFile)

Attempt to read an `LJHRecord` from `f`. Return a `Nullable{LJHRecord}` containing
that record if succesful, or one that `isnull` if not. On success the file position
moves forward, on failure the file position does not change.
"""
function tryread(f::LJHFile{LJHv, FrameTime, PretrigNSamples, NRow}) where {LJHv, FrameTime, PretrigNSamples, NRow}
    nbytes = record_nbytes(f)
    record = read(f.io, nbytes)
    if length(record) == nbytes
        rowcount, timestamp_usec = parse_record_header(f, record)
        data = reinterpret(UInt16, record[1+header_nbytes(f):nbytes])
        if f.inverted_data
            data = .~data
        end
        return Nullable(LJHRecord{FrameTime, PretrigNSamples, NRow}(data, rowcount, timestamp_usec))
    else
        seek(f.io, position(f.io)-length(record)) # go back to the start of the record
        return Nullable{LJHRecord{FrameTime, PretrigNSamples, NRow}}()
    end
end
# this version with nb_available seems like it should work, and possibly be
# more performant, but nb_available doesn't seem to work as I expect
# function tryread{T}(f::LJHFile{LJH_22,T})
#   n = nb_available(f.io)
#   if n<record_nbytes(f)
#     return Nullable{LJHRecord}()
#   else
#     rowcount = read(Int,f.io)
#     timestamp_usec = read(Int,f.io)
#     data = read(UInt16, f.io, f.record_nsamples)
#     return Nullable(LJHRecord(data, rowcount, timestamp_usec))
#   end
# end

function parse_record_header(f::LJHFile{LJH_22}, record::Vector{UInt8})
    rowcount = reinterpret(Int, record[1:8])[1]
    timestamp_usec = reinterpret(Int, record[9:16])[1]
    rowcount, timestamp_usec
end

# Used only for reading LJH version 2.1.0 files. This parses the ugly
# "encoded" version of the frame counter, which is converted into an
# approximate time, rounded to 4 microseconds, and divided into the integer
# and fractional parts of the millisecond. The latter are stored in bytes
# 3-6, and the former is in byte 1. Ignore byte 2 (it has some need to be
# 0 for backward compatibility). Ugly! That's why LJH 2.2.0 does something
# totally different.
function parse_record_header(f::LJHFile{LJH_21}, record::Vector{UInt8})
    frac = Int64(record[1])
    ms = UInt64(record[3]) |
         (UInt64(record[4])<<8) |
         (UInt64(record[5])<<16) |
         (UInt64(record[6])<<24)
    count_4usec = Int64(250ms+frac)
    ns_per_frame = round(Int64, frametime(f)*1e9)
    ns_per_4usec = Int64(4000)
    count_nsec = count_4usec*ns_per_4usec
    count_frame = cld(count_nsec, ns_per_frame)
    rowcount = count_frame*num_rows(f)+f.row
    return rowcount, 4*count_4usec
end

# Used only for reading LJH version 2.0.0 files. This parses the ugly
# "encoded" version of the frame counter, which is converted into an
# approximate time, rounded to the nearest millisecond and stored in bytes
# 3-6. Ignore bytes 1-2 (they have some need to be
# 0 for backward compatibility).
function parse_record_header(f::LJHFile{LJH_20}, record::Vector{UInt8})
    ms = UInt64(record[3]) |
         (UInt64(record[4])<<8) |
         (UInt64(record[5])<<16) |
         (UInt64(record[6])<<24)
    ns_per_frame = round(Int64, frametime(f)*1e9)
    ns_per_msec = Int64(1000000)
    count_nsec = ms*ns_per_msec
    count_frame = cld(count_nsec, ns_per_frame)
    rowcount = count_frame*num_rows(f)+f.row
    return rowcount, 1000*ms
end

watch_file(f::LJHFile, timeout_s::Real) = watch_file(f.filename, timeout_s)

"""Represent one or more LJHFiles as a seamless sequence that can be addressed
by record number from 1 to the sum of all records in the group."""
mutable struct LJHGroup
    ljhfiles::Vector{LJHFile}
    lengths::Vector{Int}
end
function LJHGroup(x::Vector)
    ljhfiles = LJHFile[LJHFile(f) for f in x]
    LJHGroup(ljhfiles, [length(f) for f in ljhfiles])
end
LJHGroup(x::LJHFile) = LJHGroup(LJHFile[x])
LJHGroup(x::String) = LJHGroup(LJHFile(x))
Base.length(g::LJHGroup) = sum(g.lengths)
Base.close(g::LJHGroup) = map(close, g.ljhfiles)
Base.open(g::LJHGroup) = map(open, g.ljhfiles)
fieldvalue(g::LJHGroup, s::Symbol) = unique([getfield(f, s) for f in g.ljhfiles])
channel(g::LJHGroup) = (@assert length(fieldvalue(g, :channum))==1; g.ljhfiles[1].channum)
record_nsamples(g::LJHGroup) = (@assert length(fieldvalue(g, :record_nsamples))==1; g.ljhfiles[1].record_nsamples)
pretrig_nsamples(g::LJHGroup) = (@assert length(unique(pretrig_nsamples.(g.ljhfiles)))==1; pretrig_nsamples(g.ljhfiles[1]))
frametime(g::LJHGroup) = (@assert length(unique(frametime.(g.ljhfiles)))==1; frametime(g.ljhfiles[1]))
column(g::LJHGroup) = (@assert length(fieldvalue(g, :column))==1; g.ljhfiles[1].column)
row(g::LJHGroup) = (@assert length(fieldvalue(g, :row))==1; g.ljhfiles[1].row)
num_columns(g::LJHGroup) = (@assert length(fieldvalue(g, :num_columns))==1; g.ljhfiles[1].num_columns)
num_rows(g::LJHGroup) = (@assert length(unique(num_rows.(g.ljhfiles)))==1; num_rows(g.ljhfiles[1]))
filenames(g::LJHGroup) = [f.filename for f in g.ljhfiles]
lengths(g::LJHGroup) = g.lengths
"    watch(g::LJHGroup, timeout_s)
calls `watch_file` on the last LJH file in `g`, passes `timeout_s` through."
watch(g::LJHGroup, timeout_s) = watch_file(last(g.ljhfiles).filename, timeout_s)
"    update_num_records(g::LJHGroup)
examine the underlying LJHFiles to determine if any grew. If the last one grew, update `g.lengths`, otherwise throw an error"
function update_num_records(g::LJHGroup)
    old_lengths = copy(g.lengths)
    new_lengths = Int[length(f) for f in g.ljhfiles]
    for i = 1:length(g.lengths)-1
        if old_lengths[i]!=new_lengths[i]
            error("a ljh file other than the last file in grew in length $g it was $(g.ljhfiles[i])")
        end
    end
    g.lengths=new_lengths
end
"    filenum_recordnum(g::LJHGroup, j::Int)
The record `g[j]` is actually `g.ljhfiles[i][k]`, return `i,k`."
function filenum_recordnum(g::LJHGroup, j::Int)
    for (i,len) in enumerate(g.lengths)
        j <= len ? (return i,j) : (j-=len)
    end
    1,1 # default return value in case of empty range
end
function Base.getindex(g::LJHGroup, i::Int)
    filenum, recordnum = filenum_recordnum(g,i)
    g.ljhfiles[filenum][recordnum]
end
Base.getindex(g::LJHGroup, slice::AbstractArray) = LJHGroupSlice(g, slice)
Base.lastindex(g::LJHGroup) = length(g)
function Base.iterate(g::LJHGroup)
    for f in g.ljhfiles seekto(f,1) end
    filenum, recordnum = filenum_recordnum(g,1)
    donefilenum, donerecordnum = filenum_recordnum(g, length(g))
    Base.iterate(g, (filenum, recordnum, donefilenum, donerecordnum))
end
function Base.iterate(g::LJHGroup, state)
    filenum, recordnum, donefilenum, donerecordnum = state
    if filenum>donefilenum || filenum==donefilenum && recordnum>donerecordnum
        return nothing
    end
    ljhrecord = pop!(g.ljhfiles[filenum])
    recordnum+=1
    recordnum > g.lengths[filenum] && (recordnum-=g.lengths[filenum];filenum+=1)
    ljhrecord, (filenum, recordnum, donefilenum, donerecordnum)
end
function Base.show(io::IO, g::LJHGroup)
    print(io, "LJHGroup with $(length(g.ljhfiles)) files, $(length(g)) records, split as $(lengths(g)),")
    print(io, "record_nsamples $(record_nsamples(g)),\n")
    print(io, "pretrig_nsamples $(pretrig_nsamples(g)).")
    print(io, "channel $(channel(g)), row $(row(g)), column $(column(g)), frametime $(frametime(g)) s.\n")
    print(io, "First filename $(g.ljhfiles[1].filename)")
end

"`LJHGroupSlice` is used to allow acces to ranges of LJH records, eg `[r.data for r in ljh[1:100]]`."
struct LJHGroupSlice{T<:AbstractArray}
    g::LJHGroup
    slice::T
    function LJHGroupSlice{T}(ljhgroup, slice) where T
        isempty(slice) || maximum(slice)<=length(ljhgroup) || error("$(maximum(slice)) is greater than nrec=$(length(ljhgroup)) in $ljhgroup")
        new(ljhgroup, slice)
    end
end
Base.length(g::LJHGroupSlice) = length(g.slice)
Base.lastindex(g::LJHGroupSlice) = length(g.slice)
LJHGroupSlice(ljhfile::LJHGroup, slice::T) where T <: AbstractArray = LJHGroupSlice{T}(ljhfile, slice)
function Base.iterate(g::LJHGroupSlice{T}) where T <: UnitRange
    for f in g.g.ljhfiles seekto(f,1) end
    isempty(g.slice) && return (2,2,1,1) # ensure done condition is immediatley met on empty range
    filenum, recordnum = filenum_recordnum(g.g, first(g.slice))
    donefilenum, donerecordnum = filenum_recordnum(g.g, last(g.slice))
    seekto(g.g.ljhfiles[filenum], recordnum)
    Base.iterate(g, (filenum, recordnum, donefilenum, donerecordnum))
end
function Base.iterate(g::LJHGroupSlice{T}, state) where T <: UnitRange
    filenum, recordnum, donefilenum, donerecordnum = state
    if filenum>donefilenum || filenum==donefilenum && recordnum>donerecordnum
        return nothing
    end
    ljhrecord = pop!(g.g.ljhfiles[filenum])
    recordnum+=1
    recordnum > g.g.lengths[filenum] && (recordnum-=g.g.lengths[filenum];filenum+=1)
    ljhrecord, (filenum, recordnum, donefilenum, donerecordnum)
end
record_nsamples(s::LJHGroupSlice) = record_nsamples(s.g)

const LJHLike = Union{LJHFile, LJHGroup, LJHGroupSlice}
"    get_data_rowcount_timestamp(g::LJHGroup)
Get all data from an `LJHGroup`, returned as a tuple of Vectors `(data, rowcount, timestamp_usec)`."
function get_data_rowcount_timestamp(g::LJHLike)
    data = Matrix{UInt16}(undef, record_nsamples(g),length(g))
    rowcount = zeros(Int64, length(g))
    timestamp_usec = zeros(Int64, length(g))
    get_data_rowcount_timestamp!(g,data,rowcount,timestamp_usec)
end
"    get_data_rowcount_timestamp!(g::LJHGroupSlice),data::Matrix{UInt16},rowcount::Vector{Int64},timestamp_usec::Vector{Int64})
Get all data from an `LJHGroupSlice`, pass in vectors of length `length(g)` and correct type to be filled with the answers."
function get_data_rowcount_timestamp!(g::LJHLike,data::Matrix{UInt16},rowcount::Vector{Int64},timestamp_usec::Vector{Int64})
    @assert size(data,2)==length(rowcount)==length(timestamp_usec)==length(g) "data, rowcount, timestap_usec: length mismatch: lengths $((size(data,2), length(rowcount), length(timestamp_usec), length(g)))"
    for (i, record) in enumerate(g)
        data[:,i] = record.data
        rowcount[i] = record.rowcount
        timestamp_usec[i] = record.timestamp_usec
    end
    data,rowcount,timestamp_usec
end

function finalize_longrecord!(longrecords, records, nsamples)
    v = Vector{UInt16}(undef, nsamples)
    i=1
    lasti=1
    for record in records
        ilast = min(i+length(record)-1, length(v))
        v[i:ilast] = data(record)[1:(ilast-i+1)]
        i=ilast+1
    end
    longrecord = LJH3Record{frameperiod(first(records))}(
        v,
        first_rising_sample(first(records)),
        frame1index(first(records)),
        timestamp_usec(first(records))
    )
    empty!(records)
    push!(longrecords, longrecord)
end

"""
     read_longrecords(ljh::LJHLike, nsamples;maxrecords=1, allowdiscontinuity=true, unevenfinalrecord=!allowdiscontinuity)

Return a `Vector{LJH3Record}` of at most `maxrecords` length, containing longrecords each with
`nsamples` samples. If `allowdiscontinuity` is `true`, each longrecord will have exactly `nsamples`
samples, even if the records in `ljh` were not continuous. If `allowdiscontinuity` is `false`
longrecords will have fewer samples when a gap between records in `ljh` exists. If `unevenfinalrecord` is `true`
the last returned longrecord may have fewer samples than `nsamples`, because `ljh` did not have enough
data for a full length longrecord.
"""
function read_longrecords(ljh::LJHLike, nsamples;maxrecords=1, allowdiscontinuity=true, unevenfinalrecord=!allowdiscontinuity)
    records = eltype(ljh)[] # records to be collated into a longrecord
    longrecords = LJH3Record{frameperiod(ljh)}[]
    for (i,record) in enumerate(ljh)
        if !isempty(records) && !allowdiscontinuity
            lastrec = last(records)
            if frame1index(record) != frame1index(lastrec)+length(lastrec)
                finalize_longrecord!(longrecords, records,sum(length.(records)))
                if length(longrecords)>=maxrecords
                    break
                end
                continue
            end
        end
        push!(records,record)
        if sum(length.(records)) >= nsamples
            finalize_longrecord!(longrecords, records, nsamples)
        end
        if length(longrecords)>=maxrecords
            break
        end
    end
    if !isempty(records) && length(longrecords) < maxrecords
        avail_nsamples = sum(length.(records))
        if avail_nsamples > nsamples || unevenfinalrecord
            finalize_longrecord!(longrecords, records,min(nsamples, avail_nsamples))
        end
    end
    return longrecords
end

"""
    create(filename::AbstractString, dt, npre, nsamp; version="2.2.0",
    channel=1, num_rows=32, num_cols=8)

Create a new file, write an LJH Header to it, and return an `LJHFile` ready for
writing with `write`.  `dt` in seconds. `npre` is the number of presamples.
`nsamp` is the number of samples per record. Versions "2.2.0", "2.1.0", and "2.0.0" are supported.
"""
function create(filename::AbstractString, dt, npre, nsamp; version="2.2.0", channel=1, num_rows=32, num_cols=8)
    f = open(filename,"w+")
    writeljhheader(f, dt, npre, nsamp; version=version, channel=channel, num_rows=num_rows, num_cols=num_cols)
    LJHFile(filename,seekstart(f))
end

"""
    writeljhheader(filename::AbstractString, dt, npre, nsamp; version="2.2.0",channel=1, num_rows=32, num_cols=8)
    writeljhheader(io::IO, dt, npre, nsamp; version="2.2.0",channel=1, num_rows=32, num_cols=8)

Write a header for an LJH file. `dt` in seconds. `npre` is the number of presamples.
`nsamp` is the number of samples per record. Versions "2.2.0", "2.1.0", and "2.0.0" are supported.
`create_ljh` is easier to use.
"""
function writeljhheader(filename::AbstractString, dt, npre, nsamp; version="2.2.0", channel=1, num_rows=32, num_cols=8)
    open(filename, "w") do f
    writeljhheader(f, dt, npre, nsamp; version=version)
    end #do
end

# the header here seems absurd, yes
# but it is intended to allow compatilibity with older LJH readers,
# like that in IGOR
function writeljhheader(io::IO, dt, npre, nsamp; version="2.2.0", channel=1, num_rows=32, num_cols=8)
    write(io,
"#LJH Memorial File Format
Save File Format Version: $version
Software Version: LJH file generated by Julia
Software Driver Version: n/a
Date: %(asctime)s GMT
Acquisition Mode: 0
Digitized Word Size in bytes: 2
Location: LANL, presumably
Cryostat: Unknown
Thermometer: Unknown
Temperature (mK): 100.0000
Bridge range: 20000
Magnetic field (mGauss): 100.0000
Detector:
Sample:
Excitation/Source:
Operator: Unknown
SYSTEM DESCRIPTION OF THIS FILE:
USER DESCRIPTION OF THIS FILE:
#End of description
Number of Digitizers: 1
Number of Active Channels: 1
Timestamp offset (s): 1304449182.876200
Digitizer: 1
Description: CS1450-1 1M ver 1.16
Master: Yes
Bits: 16
Effective Bits: 0
Anti-alias low-pass cutoff frequency (Hz): 0.000
Timebase: $dt
Number of samples per point: 1
Number of rows: $num_rows
Number of columns: $num_cols
Presamples: $npre
Total Samples: $nsamp
Trigger (V): 250.000000
Tigger Hysteresis: 0
Trigger Slope: +
Trigger Coupling: DC
Trigger Impedance: 1 MOhm
Trigger Source: CH A
Trigger Mode: 0 Normal
Trigger Time out: 351321
Use discrimination: No
Channel: $channel
Description: A (Voltage)
Range: 0.500000
Offset: -0.000122
Coupling: DC
Impedance: 1 Ohms
Inverted: No
Preamp gain: 1.000000
Discrimination level (%%): 1.000000
#End of Header\n"
    )
flush(io)
end

"    write(ljh::LJHFile{LJH_22},traces::Array{UInt16,2}, rowcounts::Vector{Int64}, times::Vector{Int64})"
function Base.write(ljh::LJHFile{LJH_22},traces::Array{UInt16,2}, rowcounts::Vector{Int64}, times::Vector{Int64})
    for j = 1:length(times)
        write(ljh, traces[:,j], rowcounts[j], times[j])
    end
end
"    write(ljh::LJHFile{LJH_22}, trace::Vector{UInt16}, rowcount::Int64, time::Int64)"
function Base.write(ljh::LJHFile{LJH_22}, trace::Vector{UInt16}, rowcount::Int64, timestamp_usec::Int64)
    write(ljh.io, rowcount, timestamp_usec, trace)
end
"    write(ljh::LJHFile{LJH_22}, record::LJHRecord)"
function Base.write(ljh::LJHFile{LJH_22}, record::LJHRecord)
  write(ljh, record.data, record.rowcount, record.timestamp_usec)
end

"    write(ljh::LJHFile,traces::Array{UInt16,2}, rowcounts::Vector{Int64})"
function Base.write(ljh::LJHFile,traces::Array{UInt16,2}, rowcounts::Vector{Int64})
    for j = 1:length(rowcounts)
        write(ljh, traces[:,j], rowcounts[j])
    end
end
"    write(ljh::LJHFile{LJH_21}, trace::Vector{UInt16}, rowcount::Int64)"
function Base.write(ljh::LJHFile{LJH_21}, trace::Vector{UInt16}, rowcount::Int64)
  lsync_us = 1e6*frametime(ljh)/num_rows(ljh)
  timestamp_us = round(Int, rowcount*lsync_us)
  timestamp_ms = Int32(div(timestamp_us, 1000))
  subms_part = round(UInt8, mod(div(timestamp_us,4), 250))
  z = Int8(0)
  write(ljh.io, subms_part)
  write(ljh.io, z) # Required to be zero.
  write(ljh.io, timestamp_ms)
  write(ljh.io, trace)
end
"    write(ljh::LJHFile{LJH_20}, trace::Vector{UInt16}, rowcount::Int64)"
function Base.write(ljh::LJHFile{LJH_20}, trace::Vector{UInt16}, rowcount::Int64)
  lsync_us = 1e6*frametime(ljh)/num_rows(ljh)
  timestamp_us = round(Int, rowcount*lsync_us)
  timestamp_ms = Int32(div(timestamp_us, 1000))
  z = Int8(0)
  write(ljh.io, z)
  write(ljh.io, z) # Required to be zero.
  write(ljh.io, timestamp_ms)
  write(ljh.io, trace)
end
"    Base.write(ljh::LJHFile, record::LJHRecord)"
function Base.write(ljh::LJHFile, record::LJHRecord)
  write(ljh, record.data, record.rowcount)
end

include("ljhutil.jl")
include("ljh3.jl")
end #module


text:
using LinearAlgebra

import ..Embeddings: AbstractEmbedding, Embedding
using ..GroupSlices: groupslices, firstinds, groupinds

include("../plot_recipes/helperfunctions_rectgrid.jl")

include("bin_encoding.jl")
include("marginal_visitation_frequency.jl")

struct RectangularPartition{D, N}
    rectangles::Vector{Rectangle}
end

Base.length(r::RectangularPartition) = length(r.rectangles)
Base.getindex(r::RectangularPartition, i) = r.rectangles[i]
Base.firstindex(r::RectangularPartition, i) = 1
Base.lastindex(r::RectangularPartition, i) = Base.length(r.rectangles)
Base.iterate(r::RectangularPartition, state = 1) = iterate(r.rectangles, state)

export RectangularPartition


text:
# arb_hypgeom.jl

module arb_hypgeom
export one_f_two_fast, one_f_two

using ArbNumerics

"""
Implementation of the hypergeometric function 1F2 using the expansion:

    1F2(a1; {b1, b2}; z) = ∑_{k=0}^∞ (a1)_k z^k / ((b1)_k ⋅ (b2)_k ⋅ k!)

where (x)_k is the Pochhammer rising factorial = Γ(x+k)/Γ(x) with Γ(x) the Gamma function.
"""

function one_f_two(a, b, x; prec = 64) # z > 0 & n >= 0

    # Initialise precision of ArbReal to prec.
    setextrabits(0)
    setprecision(ArbReal, prec + 8)

    a = ArbReal("$a")
    b = (ArbReal("$i") for i in b)
    x = ArbReal("$x")

    k = ArbReal("0")
    result = ArbReal("0.0")
    err = eps(result)  # Machine accuracy of specified precision prec.

    risingfact(w, n) = ArbReal(ArbNumerics.gamma(w + n) / ArbNumerics.gamma(w))

    while true
        term = ArbReal(risingfact(a, k) * x^k / (prod(risingfact.(b, k)) * ArbNumerics.gamma(k + 1)))

        # Break loop if term smaller than accuracy of result.
        if abs(term) < err
            break
        end

        result += term
        k += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if ball(result)[2] > err
            setprecision(ArbReal, precision(result) * 2)
            k = ArbReal("0")
            a = ArbReal("$a")
            b = (ArbReal("$i") for i in b)
            x = ArbReal("$x")
            result = ArbReal("0.0")
        end
    end
    ArbReal(result, bits = prec + 8)
end

"""
Implementation of specific 1F2 hypergeometric function used for integrals:

    ∫_0^1 cosh(zx) x^{2m} dx = 1F2(m+1/2; {1/2, m+3/2}; z^2/4)/(2m+1)
                             = ∑_{t=0}^∞ z^{2t} / ((2t+2m+1)⋅(2t)!)

    ∫_0^1 sinh(zx) x^{2m} dx = z ⋅ 1F2(m+1; {3/2, m+2}; z^2/4)/(2m+2)
                             = ∑_{t=0}^∞ z^{2t+1} / ((2t+2m+2)⋅(2t+1)!)

which we can combine into a generic summation:

    ∑_{t=0}^∞ z^{2t+h} / ((2t+2m+1+h)⋅(2t+h)!)

which gives the cosh version for h=0 and the sinh version for h=1. This specialised 1F2 converges faster than generic 1F2 algorithm, one_f_two.
"""

function one_f_two_fast(x, m, h; prec = 64) # z > 0 & n >= 0

    # Initialise precision of ArbReal to prec.
    p = prec
    setextrabits(0)
    setprecision(ArbReal, p)

    x = ArbReal("$x")
    m = ArbReal("$m")
    h = ArbReal("$h")

    k = ArbReal("0")
    result = ArbReal("0.0")
    term = ArbReal("1.0")
    err = eps(result)  # Machine accuracy of specified precision prec.

    while abs(midpoint(term)) > err * abs(midpoint(result))
        term = ArbReal(x^(2 * k + h) / ((2 * k + 2 * m + 1 + h) * ArbNumerics.gamma(2 * k + 1 + h)))
        result += term
        # println("term: k = ", k, "\nterm value: ", ball(ArbReal(term, bits = prec)), "\ncumulant result: ", ball(ArbReal(result, bits = prec)), "\n")
        k += 1

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if radius(result) > err * abs(midpoint(result))
            p *= 2
            setprecision(ArbReal, p)
            # println("Not precise enough. Error = ", abs(radius(result)/midpoint(result)), " > ", err, ". Increasing precision to ", p, " bits.\n")
            x = ArbReal("$x")
            m = ArbReal("$m")
            h = ArbReal("$h")
            k = ArbReal("0")
            result = ArbReal("0.0")
            term = ArbReal("1.0")
        end
    end
    # println("x: ", ArbReal(x, bits = prec), ". Final result: ", ArbReal(result, bits = prec))
    ArbReal(result, bits = prec)
end

# β = ArbReal("2.0")
# α = ArbReal("7.0")
# v = ArbReal("5.8")
# w = ArbReal("1.6")
# R = ArbReal((v^2 - w^2) / (w^2 * v))
# a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
# z = ArbReal("90.61")
# m = ArbReal("1.0")
# @time c = one_f_two_fast(z, m, 0; prec = 64)
# @show(c)

end # end module

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

"""
Some plots and tests of the above functions to check that they produce the correct values to the specified precision.
"""

# using Plots
# using QuadGK
# plotly()
# Plots.PlotlyBackend()
#
# cosh_integral(z, m) = quadgk(t -> (m + 1) * cosh(z * t) * t^m, BigFloat(0), BigFloat(1))[1]
#
# # x_range = 0:200
# # m = 1
# # @time I = [cosh_integral(x, m) for x in x_range]
# # @time F = abs.([arb_hypgeom_1f2(m / 2 + 1 / 2, (1 / 2, m / 2 + 3/2), x^2 / 4.0; prec = 128) for x in x_range])
# # p = plot(x_range, F, yaxis=:log)
# # plot!(x_range, I)
# # display(p)
# # @show(I, F)
# # println(m / 2 + 1 / 2, (1 / 2, m / 2 + 3/2), x_range^2 / 4.0)
#
# """
# Conclusion: The integral is a bit faster than the hypergeometric function, but not as accurate. So, probably use 1F2 since the time difference doesn't seem to be major!
# """
#
# # @time F_test = arb_hypgeom_1f2(1/2, (1/2, 1/2), 1/2; prec = 512)
# # setprecision(BigFloat, 512)
# # t_test = parse(BigFloat, "2.17818355660857086398922206782012528343129403292165693281081574094992093930206091916837394310903919052930938811320102906774905335969974706202490049810447181898779259195640840870264051023127290510625691429186995752638681515944856002645204867940492548077249544")
# # @show(F_test, t_test)
# # @show(isequal("$F_test"[end-1], "$t_test"[end-1]))
#
# """
# The 1F2 function above can accurate produce high precision values as shown in the above test. But it can be simplified and optimized mathematically initially as shown below.
# """
#
# x_range = 0:200
# m = 1
# h = 0
# @time L = abs.([arb_hypergeom_1f2_fast(x, m, h; prec = 128) for x in x_range])
# # @time F = abs.([arb_hypgeom_1f2(m + 1 / 2, (1 / 2, m + 3/2), x^2 / 4.0; prec = 128) / (2 * m + 1) for x in x_range])
# p = plot(x_range, L, yaxis=:log)
# # plot!(x_range, F)
# display(p)
# @show(L)


text:
global ADDIR
global PROC

PROC="GPU"

ADDIR=pwd()
cd("src")
include("$ADDIR/src/setup.jl")
#usegpu(true)
useproc(PROC)
using AutoDiff
cd("../")

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
@everywhere begin

using Dagger
import Dagger: ThreadProc

struct FakeProc <: Dagger.Processor
    owner::Int
end
FakeProc() = FakeProc(myid())
Dagger.get_parent(proc::FakeProc) = OSProc(proc.owner)
Base.:(==)(proc1::FakeProc, proc2::FakeProc) = proc1.owner == proc2.owner

struct FakeVal
    x
end
fakesum(xs...) = FakeVal(sum(map(y->y.x, xs)))

Dagger.iscompatible_func(proc::FakeProc, opts, f) = true
Dagger.iscompatible_arg(proc::FakeProc, opts, ::Type{<:Integer}) = true
Dagger.iscompatible_arg(proc::FakeProc, opts, ::Type{<:FakeVal}) = true
Dagger.move(from_proc::ThreadProc, to_proc::FakeProc, x::Integer) = FakeVal(x)
Dagger.execute!(proc::FakeProc, func, args...) = FakeVal(42+func(args...).x)

end