text:
export Bar, NotchedBar

"""
    Bar
"""
struct Bar <: PDBGeom
    obj
end

function Bar(;bounds=[0.0 10; 0 2;0 2])
    Bar(Cuboid(bounds))
end

function PDMesh.create(geom::Bar; resolution=0.1, rand_=0.05, type=1)
    block = geom.obj
    return PDMesh.create(block; resolution=resolution, rand_=rand_, type=type)
end

"""
    NotchedBar
"""
struct NotchedBar <: PDBGeom
    obj
    notch
end

function NotchedBar(;bounds=[0.0 7.0; 0 1;0 1], notch=nothing)
    NotchedBar(Cuboid(bounds), notch)
end

function PDMesh.create(geom::NotchedBar; resolution=0.1, rand_=0.05, type=1)
    obj = geom.obj
    if geom.notch !== nothing
        obj = PDMesh.delete(obj, out -> begin x=out[1]; prod((x .<= geom.notch[:, 2]) .& (x .>= geom.notch[:, 1] ), dims=1) end)
    end
    return PDMesh.create(obj; resolution=resolution, rand_=rand_, type=type)
end

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark


text:
"""
    linreg(X, y)

    Objective function of Linear Regression on data matrix X and label vector y.
"""
linreg(X::Array{Float64,2}, y::Array{Float64,1}) = Objective(
    w::Array{Float64,1} -> _linreg_f(X, y, w),  # f
    w::Array{Float64,1} -> _linreg_g(X, y, w),  # g
    f_b = (w::Array{Float64,1}, b::Array{Int,1}) -> _linreg_f(X, y, w, b),
    g_b = (w::Array{Float64,1}, b::Array{Int,1}) -> _linreg_g(X, y, w, b),
    dim = length(y)
)

function _linreg_f(X::Array{Float64,2}, y::Array{Float64,1}, w::Array{Float64,1},
                   b::Array{Int,1}=Int[])
    if length(b) == 0
        return vecnorm(y - X*w)^2 / (2*length(y))
    else
        return vecnorm(y[b] - X[b,:]*w)^2 / (2*length(b))
    end
end

function _linreg_g(X::Array{Float64,2}, y::Array{Float64,1}, w::Array{Float64,1},
                   b::Array{Int,1}=Int[])
    if length(b) == 0
        return vec(X' * (X*w - y)) ./ length(y)
    else
        return vec(X[b,:]' * (X[b,:]*w - y[b])) ./ length(b)
    end
end

"""
    linreg_predict(w, X)

    Predict the examples in X with Linear Regression under parameter vector w.
"""
linreg_predict(w::Array{Float64,1}, X::Array{Float64,2}) = sign(vec(X*w))

ext: jl

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
# ellipse to rotated and shifted ellipse, initial and final landmark positions observed
using Revise
using Polynomials.PolyCompat
using BridgeLandmarks
const BL = BridgeLandmarks
using Distributions
using LinearAlgebra
using Random
using FileIO
using JLD2

Random.seed!(9)
workdir = @__DIR__
cd(workdir)
outdir = workdir

# read data
dat = load("data_exp1_hard.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]
n = dat["n"]

# set pars
p_ms = Pars_ms(δmom=0.01/n, cinit=0.2, γinit=2.0/√n, σobs = 0.001, skip_saveITER=1)
#p_ahs = Pars_ahs(δmom=0.1, cinit=0.02, γinit=0.2, db=[2.5,2.5],stdev=0.75, σobs = 0.01, skip_saveITER=1)
ainit = mean(norm.(diff(xobs0)))/1.0
#ups = [:rmrw_mom, :innov, :parameter]
#ups = [:sgd_mom]#, :innov, :parameter]
#ups = [:rmrw_mom, :parameter, :innov]
ups = [:mala_mom, :parameter, :innov]

# run algorithm
@time landmarksmatching(xobs0,xobsT; ITER=10,pars=p_ms, outdir=outdir, ainit=ainit, updatescheme=ups)
#@time landmarksmatching(xobs0,xobsT; ITER=10,pars=p_ahs, outdir=outdir, ainit=ainit, updatescheme=ups)

plotlandmarksmatching(outdir)


text:
@testset "QuGAN" begin
    include("QuGAN.jl")
end

@testset "QCBM" begin
    include("QCBM.jl")
end


text:
#####################################################
#Train functions

#Train classifier in S(D)AE with classifier
function train!(sae,smallimgs,
	 			labels,
				smallimgstest,
				labelstest,
				conf::configClass)
  		errors = zeros(3,length(conf.hidden_sizes),conf.n_inits)
		nets = Array{AE, 1}()
	#get hidden repr. of last AE in stack (n_ae is global!):
	global n_ae = sae.nae+1
	for i in 1:length(conf.hidden_sizes)
		if conf.hidden_sizes[i] == 0
			print("-> shallow classifier...\n")
		else
			print(string("-> classifier with 1 hidden layer with",conf.hidden_sizes[i]," hidden neurons...\n"))
		end
		for j in 1:conf.n_inits
			if conf.hidden_sizes[i] == 0
				net = Classifier([size(sae.reps[end][1])[1], 10])
			else
				net = Classifier([size(sae.reps[end][1])[1], conf.hidden_sizes[i], 10])
			end
      		@time loss = learn!(net, sae, gethiddenreps, getlabel, conf.iterations, conf.learningrates, conf.weight_decay,
                            nonlinearity = conf.nonlinearity, nonlinearity_diff = conf.nonlinearity_diff,
                            nonlin_diff = conf.nonlin_diff)
			push!(nets,deepcopy(net))
			error_train, error_test = printerrors(net, sae.reps[end][1], labels, sae.reps[end][2], labelstest, conf)
      		errors[:,i,j] = [conf.hidden_sizes[i], error_train, error_test]
     	end
  	end
  	errors, nets
end

function trainSAEc_full_patchy!(smallimgs,smallimgstest,labels,labelstest,
									n_out_per_population,conf;
									patchsize = 10,
									n_of_populations = 100,
									stride = 3, im_size = 28,
									no_recurrence_after_training = false,
									Földiak_model = false,
									rand_patch_pos = true,
									colorimage = false)

	n_in = size(smallimgs)[1]
	AE_s_p = AE_sparse_patchy(n_in, n_out_per_population, n_of_populations;
							patchsize = patchsize, stride = stride,
							im_size = im_size, colorimage = colorimage)
	rand_patch_pos && (global global_patch_positions = getrandompatchpositions(AE_s_p; colorimage = colorimage))
	confae = configAE_sparse(conf.configSAE_sparse.task,conf.configSAE_sparse.hidden_sizes[1],
				conf.configSAE_sparse.n_inits,conf.configSAE_sparse.iterations[1],
				conf.configSAE_sparse.learningrates[1],conf.configSAE_sparse.lambdas[1],
				conf.configSAE_sparse.weight_decay,conf.configSAE_sparse.noise_type,
				conf.configSAE_sparse.noise[1],conf.configSAE_sparse.nonlinearity[1],
				conf.configSAE_sparse.nonlinearity_diff[1],true)

	print("train sparse coding layer...\n")
	learnAE_Brito_patchy!(AE_s_p, getsmallimg, confae.iterations, confae.learningrates,
						confae.weight_decay, confae.noise_type, confae.noise, confae.lambda;
						nonlinearity = confae.nonlinearity, nonlinearity_diff = confae.nonlinearity_diff,
						nonlin_diff = confae.nonlin_diff, BPuptolayer = AE_s_p.AE_sparse_patches[1].nl,
						Földiak_model = Földiak_model, rand_patch_pos = rand_patch_pos,
						colorimage = colorimage)

	#construct SAE with trained patchy SAE
	sae = SAE_sparse_patchy(AE_s_p.n_AE_sparse_patches,n_in,n_out_per_population;
							n_samples = [size(smallimgs)[2], size(smallimgstest)[2]])
	sae.nae = 1
	sae.aes = [deepcopy(AE_s_p)]

	print("generate sparse hidden representations...\n")
	generatehiddenreps_Brito!(sae,1, smallimgs, smallimgstest, confae;
		no_recurrence = no_recurrence_after_training,
		Földiak_model = Földiak_model, rand_patch_pos = rand_patch_pos,
		colorimage = colorimage)
	for reps in sae.reps[1]
		reps .*= 1/maximum(reps)
	end
	print("train top-end classifier...\n")
	errors, nets = train!(sae,smallimgs,labels,smallimgstest,labelstest,conf.configClass)

	return sae, errors
end


text:
abstract type Distribution end

function ==(d1::D, d2::D) where D<:Distribution
    ret = true
    for fld in fieldnames(D)
        ret = ret && getfield(d1, fld) == getfield(d2, fld)
    end
    return ret
end

# DISTRIBUTION AS TRAIT

# fallback
isdistribution(::Type{<:Any}) = false

# for distributions in Distributions.jl:
isdistribution(::Type{<:Distributions.Distribution}) = true

# for MLJ custom distibutions defined below:
isdistribution(::Type{Distribution}) = true

isdistribution(d) = isdistribution(typeof(d))

## UNIVARIATE NOMINAL PROBABILITY DISTRIBUTION

"""
    UnivariateNominal(prob_given_label)

A discrete univariate distribution whose finite support is the set of keys of the
provided dictionary, `prob_given_label`. The dictionary values specify
the corresponding probabilities, which must be nonnegative and sum to
one.

    UnivariateNominal(labels, p)

A discrete univariate distribution whose finite support is the
elements of the vector `labels`, and whose corresponding probabilities
are elements of the vector `p`.

````julia
d = UnivariateNominal(["yes", "no", "maybe"], [0.1, 0.2, 0.7])
pdf(d, "no") # 0.2
mode(d) # "maybe"
rand(d, 5) # ["maybe", "no", "maybe", "maybe", "no"]
d = fit(UnivariateNominal, ["maybe", "no", "maybe", "yes"])
pdf(d, "maybe") ≈ 0.5 # true
````

"""
struct UnivariateNominal{L,T<:Real} <: Distribution
    prob_given_label::Dict{L,T}
    function UnivariateNominal{L,T}(prob_given_label::Dict{L,T}) where {L,T<:Real}
        p = values(prob_given_label) |> collect
        Distributions.@check_args(UnivariateNominal, Distributions.isprobvec(p))
        return new{L,T}(prob_given_label)
    end
end
UnivariateNominal(prob_given_label::Dict{L,T}) where {L,T<:Real} =
    UnivariateNominal{L,T}(prob_given_label)

function UnivariateNominal(labels::Vector{L}, p::Vector{T}) where {L,T<:Real}
        Distributions.@check_args(UnivariateNominal, length(labels)==length(p))
        prob_given_label = Dict{L,T}()
        for i in eachindex(p)
            prob_given_label[labels[i]] = p[i]
        end
        return  UnivariateNominal(prob_given_label)
end

function average(dvec::Vector{UnivariateNominal{L,T}}; weights=nothing) where {L,T}

    n = length(dvec)

    Distributions.@check_args(UnivariateNominal, weights == nothing || n==length(weights))

    if weights == nothing
        weights = fill(1/n, n)
    else
        weights = weights/sum(weights)
    end

    # get all labels:
    labels = reduce(union, [keys(d.prob_given_label) for d in dvec])

    z = Dict{L,T}([x => zero(T) for x in labels]...)
    prob_given_label_vec = map(dvec) do d
        merge(z, d.prob_given_label)
    end

    # initialize the prob dictionary for the distribution sum:
    prob_given_label = Dict{L,T}()
    for x in labels
        prob_given_label[x] = zero(T)
    end

    # sum up:
    for x in labels
        for k in 1:n
            prob_given_label[x] += weights[k]*prob_given_label_vec[k][x]
        end
    end

    return UnivariateNominal(prob_given_label)

end

function Distributions.mode(d::UnivariateNominal)
    dic = d.prob_given_label
    p = values(dic)
    max_prob = maximum(p)
    m = first(first(dic)) # mode, just some label for now
    for (x, prob) in dic
        if prob == max_prob
            m = x
            break
        end
    end
    return m
end

Distributions.pdf(d::UnivariateNominal, x) = d.prob_given_label[x]

"""
    _cummulative(d::UnivariateNominal)

Return the cummulative probability vector `[0, ..., 1]` for the
distribution `d`, using whatever ordering is used in the dictionary
`d.prob_given_label`. Used only for to implement random sampling from
`d`.

"""
function _cummulative(d::UnivariateNominal{L,T}) where {L,T<:Real}
    p = collect(values(d.prob_given_label))
    K = length(p)
    p_cummulative = Array{T}(undef, K + 1)
    p_cummulative[1] = zero(T)
    p_cummulative[K + 1] = one(T)
    for i in 2:K
        p_cummulative[i] = p_cummulative[i-1] + p[i-1]
    end
    return p_cummulative
end

"""
_rand(p_cummulative)

Randomly sample the distribution with discrete support `1:n` which has
cummulative probability vector `p_cummulative=[0, ..., 1]` (of length
`n+1`). Does not check the first and last elements of `p_cummulative`
but does not use them either.

"""
function _rand(p_cummulative)
    real_sample = rand()
    K = length(p_cummulative)
    index = K
    for i in 2:K
        if real_sample < p_cummulative[i]
            index = i - 1
            break
        end
    end
    return index
end

function Base.rand(d::UnivariateNominal)
    p_cummulative = _cummulative(d)
    labels = collect(keys(d.prob_given_label))
    return labels[_rand(p_cummulative)]
end

function Base.rand(d::UnivariateNominal, n::Int)
    p_cummulative = _cummulative(d)
    labels = collect(keys(d.prob_given_label))
    return [labels[_rand(p_cummulative)] for i in 1:n]
end

function Distributions.fit(d::Type{<:UnivariateNominal}, v::AbstractVector{L}) where L
    N = length(v)
    prob_given_label = Dict{L,Float64}()
    count_given_label = Distributions.countmap(v)
    for (x, c) in count_given_label
        prob_given_label[x] = c/N
    end
    return UnivariateNominal(prob_given_label)
end

# if fitting to categorical array, must include missing labels with prob zero
function Distributions.fit(d::Type{<:UnivariateNominal}, v::CategoricalVector{L,R,V}) where {L,R,V}
    N = length(skipmissing(v) |> collect)
    prob_given_label = Dict{V,Float64}() # V is the type of levels
    for x in levels(v)
        prob_given_label[x] = 0
    end
    count_given_label = Distributions.countmap(skipmissing(v) |> collect)
    for (x, c) in count_given_label
        prob_given_label[x] = c/N
    end
    return UnivariateNominal(prob_given_label)
end

## NORMAL DISTRIBUTION ARITHMETIC

# *(lambda::Number, d::Distributions.Normal) = Distributions.Normal(lambda*d.μ, lambda*d.σ)
# function +(ds::Distributions.Normal...)
#     μ = sum([d.μ for d in ds])
#     σ = sum([d.σ^2 for d in ds]) |> sqrt
#     return Distributions.Normal(μ, σ)
# end

# hack for issue #809 in Distributions.jl:
Distributions.fit(::Type{Distributions.Normal{T}}, args...) where T = Distributions.fit(Distributions.Normal, args...)


text:
"""
    a(m,q)

Compute the characteristic value \$a_m(q)\$ corresponding to an even (cosine-elliptic) solution to the angular Mathieu equation. `m` is the order of the solution and must be an integer \$m=0,1,2,3…\$.

Pass a vector, e.g. `m=0:4`, to compute multiple characteristic values efficiently.
"""
function a(m::Order,q::Number,N::Integer=matsize2(m,q))
    ie = map(iseven,m)
    all(ie)  && return ap(div.(m,2),q,N) # all even
    !any(ie) && return aa(div.(m,2),q,N) # all odd
    a = zeros(length(m))
    a[ie]   = ap(div.(filter(iseven,m),2),q,N)
    a[.!ie] = aa(div.(filter(isodd, m),2),q,N)
    return a
end

"""
    b(m,q)

Compute the characteristic value \$b_m(q)\$ corresponding to an odd (sine-elliptic) solution to the angular Mathieu equation. `m` is the order of the solution and must be an integer, excluding 0, \$m=1,2,3,4…\$.

Pass a vector, e.g. `m=1:5`, to compute multiple characteristic values efficiently.
"""
function b(m::Order,q::Number,N::Integer=matsize2(m,q))
    ie = map(iseven,m)
    all(ie)  && return bp(div.(m.-2,2),q,N) # all even
    !any(ie) && return ba(div.(m.-1,2),q,N) # all odd
    b = zeros(length(m))
    b[ie]   = bp(div.(filter(iseven,m).-2,2),q,N)
    b[.!ie] = ba(div.(filter(isodd, m).-1,2),q,N)
    return b
end

"""
    ap(n,q)

Compute the characteristic value \$a_{2n}(q)\$ corresponding to an even (cosine-elliptic), π-periodic solution to the angular Mathieu equation, for \$n=0,1,2,3,…\$.

Pass a vector, e.g. `n=0:4`, to compute multiple characteristic values efficiently.
"""
function ap(n::Order,q::Number,N::Integer=matsize1(n,q))
    C = mat_C1(q,N)
    a = eigvals(C)[n.+1] # stev aka STERF
    #a = LAPACK.stebz!('I','E',0.0,0.0,1,maximum(n)+1,2*eps(q),C.dv,C.ev)[1][n+1]
    return a
end

"""
    aa(n,q)

Compute the characteristic value \$a_{2n+1}(q)\$ corresponding to an even (cosine-elliptic), π-antiperiodic (2π-periodic) solution to the angular Mathieu equation, for \$n=0,1,2,3,…\$.

Pass a vector, e.g. `n=0:4`, to compute multiple characteristic values efficiently.
"""
function aa(n::Order,q::Number,N::Integer=matsize1(n,q))
    C = mat_C2(q,N)
    a = eigvals(C)[n.+1]
    return a
end

"""
    ba(n,q)

Compute the characteristic value \$b_{2n+1}(q)\$ corresponding to an odd (sine-elliptic), π-antiperiodic (2π-periodic) solution to the angular Mathieu equation, for \$n=0,1,2,3,…\$.

Pass a vector, e.g. `n=0:4`, to compute multiple characteristic values efficiently.
"""
function ba(n::Order,q::Number,N::Integer=matsize1(n,q))
    C = mat_C3(q,N)
    b = eigvals(C)[n.+1]
    return b
end

"""
    bp(n,q)

Compute the characteristic value \$b_{2n+2}(q)\$ corresponding to an odd (sine-elliptic), π-periodic solution to the angular Mathieu equation, for \$n=0,1,2,3,…\$.

Pass a vector, e.g. `n=0:4`, to compute multiple characteristic values efficiently.
"""
function bp(n::Order,q::Number,N::Integer=matsize1(n,q))
    C = mat_C4(q,N)
    b = eigvals(C)[n.+1]
    return b
end

"""
    char_per(m,q)

Compute the characteristic value \$a_{m}(q)\$ for even `m`, or \$b_{m+1}(q)\$ for odd `m`. These characteristic values correspond respectively to even and odd π-periodic solutions to the angular Mathieu equation, for \$n=0,1,2,3,…\$.

|   m |  value |
| --- |--------|
|   0 |    a_0 |
|   1 |    b_2 |
|   2 |    a_2 |
|   3 |    b_4 |
|   4 |    a_4 |
|   ⋮ |       ⋮ |

Pass a vector, e.g. `m=0:4`, to compute multiple characteristic values efficiently.
"""
function char_per(m::Order,q::Number,N::Integer=matsize2(m,q))
    ie = map(iseven,m)
    all(ie)  && return ap(div.(m,2),q,N)
    !any(ie) && return bp(div.(m,2),q,N)
    a = zeros(length(m))
    a[ie]   = ap(div.(filter(iseven,m),2),q,N)
    a[.!ie] = bp(div.(filter(isodd, m),2),q,N)
    return a
end

"""
    char_aper(m,q)

Compute the characteristic value \$a_{m+1}(q)\$ for even `m`, or \$b_{m}(q)\$ for odd `m`. These characteristic values correspond respectively to even and odd π-antiperiodic (2π-periodic) solutions to the angular Mathieu equation, for \$n=0,1,2,3,…\$.

|   m |  value |
| --- |--------|
|   0 |    a_1 |
|   1 |    b_1 |
|   2 |    a_3 |
|   3 |    b_3 |
|   4 |    a_5 |
|   ⋮ |       ⋮ |

Pass a vector, e.g. `m=0:4`, to compute multiple characteristic values efficiently.
"""
function char_aper(m::Order,q::Number,N::Integer=matsize2(m,q))
    ie = map(iseven,m)
    all(ie)  && return aa(div.(m,2),q,N)
    !any(ie) && return ba(div.(m,2),q,N)
    a = zeros(length(m))
    a[ie]   = aa(div.(filter(iseven,m),2),q,N)
    a[.!ie] = ba(div.(filter(isodd, m),2),q,N)
    return a
end


text:
#=
Computes asset prices in an endowment economy when the endowment obeys
geometric growth driven by a finite state Markov chain.  The transition
matrix of the Markov chain is P, and the set of states is s.  The
discount factor is beta, and gamma is the coefficient of relative risk
aversion in the household's utility function.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-06-27

References
----------

http://quant-econ.net/jl/markov_asset.html
=#

"""
A class to compute asset prices when the endowment follows a finite Markov chain

##### Fields

- `bet::Float64` : Discount factor in (0, 1)
- `P::Matrix{Float64}` A valid stochastic matrix
- `s::Vector{Float64}` : Growth rate of consumption in each state
- `gamma::Float64` : Coefficient of risk aversion
- `n::Int(size(P, 1))`: The numberof states
- `P_tilde::Matrix{Float64}` : modified transition matrix used in computing the
price of the lucas tree
- `P_check::Matrix{Float64}` : modified transition matrix used in computing the
price of the consol

"""
type AssetPrices
    bet::Real
    P::Matrix
    s::Vector
    gamm::Real
    n::Int
    P_tilde::Matrix
    P_check::Matrix
end

"""
Construct an instance of `AssetPrices`, where `n`, `P_tilde`, and `P_check` are
computed automatically for you. See also the documentation for the type itself
"""
function AssetPrices(bet::Real, P::Matrix, s::Vector, gamm::Real)
    P_tilde = P .* s'.^(1-gamm)
    P_check = P .* s'.^(-gamm)
    return AssetPrices(bet, P, s, gamm, size(P, 1), P_tilde, P_check)
end

"""
Computes the function v such that the price of the lucas tree is v(lambda)C_t

##### Arguments

- `ap::AssetPrices` : An instance of the `AssetPrices` type

##### Returns

- `v::Vector{Float64}` : the pricing function for the lucas tree

"""
function tree_price(ap::AssetPrices)
    # Simplify names
    P, s, gamm, bet, P_tilde = ap.P, ap.s, ap.gamm, ap.bet, ap.P_tilde

    # Compute v
    I = eye(ap.n)
    O = ones(ap.n)
    v = bet .* ((I - bet .* P_tilde)\ (P_tilde * O))
    return v
end

"""
Computes price of a consol bond with payoff zeta

##### Arguments

- `ap::AssetPrices` : An instance of the `AssetPrices` type
- `zeta::Float64` : Per period payoff of the consol

##### Returns

- `pbar::Vector{Float64}` : the pricing function for the lucas tree

"""
function consol_price(ap::AssetPrices, zet::Real)
    # Simplify names
    P, s, gamm, bet, P_check = ap.P, ap.s, ap.gamm, ap.bet, ap.P_check

    # Compute v
    I = eye(ap.n)
    O = ones(ap.n)
    v = bet .* ((I - bet .* P_check)\ (P_check * (zet .*O)))
    return v
end

"""
Computes price of a call option on a consol bond, both finite and infinite
horizon

##### Arguments

- `zeta::Float64` : Coupon of the console
- `p_s::Float64` : Strike price
- `T::Vector{Int}(Int[])`: Time periods for which to store the price in the
finite horizon version
- `epsilon::Float64` : Tolerance for infinite horizon problem

##### Returns

- `w_bar::Vector{Float64}` Infinite horizon call option prices
- `w_bars::Dict{Int, Vector{Float64}}` A dictionary of key-value pairs {t: vec},
where t is one of the dates in the list T and vec is the option prices at that
date

"""
function call_option(ap::AssetPrices, zet::Real, p_s::Real,
                     T::Vector{Int}=Int[], epsilon=1e-8)
    # Simplify names, initialize variables
    P, s, gamm, bet, P_check = ap.P, ap.s, ap.gamm, ap.bet, ap.P_check

    # Compute consol price
    v_bar = consol_price(ap, zet)

    # Compute option price
    w_bar = zeros(ap.n)

    err = epsilon + 1.0
    t = 0
    w_bars = Dict{Int, Vector{eltype(P)}}()
    while err > epsilon
        if t in T w_bars[t] = w_bar end
        # Maximize across columns
        w_bar_new = max(bet .* (P_check * w_bar), v_bar.-p_s)
        # Find maximal difference of each component
        err = Base.maxabs(w_bar - w_bar_new)
        # Update
        w_bar = w_bar_new
        t += 1
    end
    return w_bar, w_bars
end

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
module runtests

using Dates
using DataFrames
using DataArrays

tests = ["abstractFuncs.jl",
         "preservingFuncs.jl",
         "constraints.jl",
         "constructors.jl",
         "convert.jl",
         "dataframe_extensions.jl",
         "doctests.jl",
         "getindex.jl",
         "setindex.jl",
         "getEntries.jl",
         "iterators.jl",
         "map.jl",
         "io.jl",
         "timematr.jl",
         "timenum.jl",
         "timedata.jl",
         "operators.jl"]

println("Running TimeData tests:")

for t in tests
    println("")
    println("----------------------")
    println(" * $(t)")
    println("----------------------")
    println("")
    include(string(Pkg.dir("TimeData"), "/test/", t))
end

end


text:
function _readgeom(wktstring::String, wktreader::WKTReader, context::GEOSContext = _context)
    result = GEOSWKTReader_read_r(context.ptr, wktreader.ptr, pointer(wktstring))
    if result == C_NULL
        error("LibGEOS: Error in GEOSWKTReader_read_r while reading $wktstring")
    end
    result
end
_readgeom(wktstring::String, context::GEOSContext = _context) =
    _readgeom(wktstring, WKTReader(context), context)

function _readgeom(wkbbuffer::Vector{Cuchar}, wkbreader::WKBReader, context::GEOSContext = _context)
    result = GEOSWKBReader_read_r(context.ptr, wkbreader.ptr, pointer(wkbbuffer), length(wkbbuffer))
    if result == C_NULL
        error("LibGEOS: Error in GEOSWKBReader_read_r")
    end
    result
end
_readgeom(wkbbuffer::Vector{Cuchar}, context::GEOSContext = _context) =
    _readgeom(wkbbuffer, WKBReader(context), context)

function _writegeom(geom::GEOSGeom, wktwriter::WKTWriter, context::GEOSContext = _context)
    GEOSWKTWriter_write_r(context.ptr, wktwriter.ptr, geom)
end

function _writegeom(geom::GEOSGeom, wkbwriter::WKBWriter, context::GEOSContext = _context)
    wkbsize = Ref{Csize_t}()
    result = GEOSWKBWriter_write_r(context.ptr, wkbwriter.ptr, geom, wkbsize)
    unsafe_wrap(Array, result, wkbsize[], own=true)
end

_writegeom(geom::GEOSGeom, context::GEOSContext = _context) =
    _writegeom(geom, WKTWriter(context), context)

# -----
# Coordinate Sequence functions
# -----
"""
    createCoordSeq(size::Integer; ndim::Integer=2) -> Ptr{Ptr{Void}}

Create a Coordinate sequence with ``size'' coordinates of ``dims'' dimensions (Return NULL on exception)
"""
function createCoordSeq(size::Integer, context::GEOSContext = _context; ndim::Integer=2)
    @assert ndim >= 2
    result = GEOSCoordSeq_create_r(context.ptr, size, ndim)
    if result == C_NULL
        error("LibGEOS: Error in GEOSCoordSeq_create")
    end
    result
end

# Clone a Coordinate Sequence (Return NULL on exception)
function cloneCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSCoordSeq_clone_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSCoordSeq_clone")
    end
    result
end

function destroyCoordSeq(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSCoordSeq_destroy_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSCoordSeq_destroy")
    end
    result
end

# Set ordinate values in a Coordinate Sequence (Return 0 on exception)
function setX!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
    result = GEOSCoordSeq_setX_r(context.ptr, ptr, i-1, value)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_setX")
    end
    result
end

function setY!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
    result = GEOSCoordSeq_setY_r(context.ptr, ptr, i-1, value)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_setY")
    end
    result
end

function setZ!(ptr::GEOSCoordSeq, i::Integer, value::Real, context::GEOSContext = _context)
    result = GEOSCoordSeq_setZ_r(context.ptr, ptr, i-1, value)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_setZ")
    end
    result
end

"Get size info from a Coordinate Sequence (Return 0 on exception)"
function getSize(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    out = Ref{UInt32}()
    result = GEOSCoordSeq_getSize_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getSize")
    end
    Int(out[])
end

"Get dimensions info from a Coordinate Sequence (Return 0 on exception)"
function getDimensions(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    out = Ref{UInt32}()
    result = GEOSCoordSeq_getDimensions_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getDimensions")
    end
    Int(out[])
end

# convenience functions
function setCoordSeq!(ptr::GEOSCoordSeq, i::Integer, coords::Vector{Float64}, context::GEOSContext = _context)
    ndim = length(coords)
    @assert ndim >= 2
    setX!(ptr, i, coords[1], context)
    setY!(ptr, i, coords[2], context)
    ndim >= 3 && setZ!(ptr, i, coords[3], context)
    ptr
end

"""
    createCoordSeq(x::Real, y::Real) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a single 2D coordinate
"""
function createCoordSeq(x::Real, y::Real, context::GEOSContext = _context)
    coordinates = createCoordSeq(1, context, ndim=2)
    setX!(coordinates, 1, x, context)
    setY!(coordinates, 1, y, context)
    coordinates
end

"""
    createCoordSeq(x::Real, y::Real, z::Real) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a single 3D coordinate
"""
function createCoordSeq(x::Real, y::Real, z::Real, context::GEOSContext = _context)
    coordinates = createCoordSeq(1, context, ndim=3)
    setX!(coordinates, 1, x, context)
    setY!(coordinates, 1, y, context)
    setZ!(coordinates, 1, z, context)
    coordinates
end

"""
    createCoordSeq(coords::Vector{Float64}) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a single N dimensional coordinate
"""
function createCoordSeq(coords::Vector{Float64}, context::GEOSContext = _context)
    ndim = length(coords)
    @assert ndim >= 2
    coordinates = createCoordSeq(1, context, ndim=ndim)
    setCoordSeq!(coordinates, 1, coords, context)
end

"""
    createCoordSeq(coords::Vector{Float64}) -> Ptr{Ptr{Void}}

Create a createCoordSeq of a multiple N dimensional coordinate
"""
function createCoordSeq(coords::Vector{Vector{Float64}}, context::GEOSContext = _context)
    ncoords = length(coords)
    @assert ncoords > 0
    ndim = length(coords[1])
    coordinates = createCoordSeq(ncoords, context, ndim=ndim)
    for (i,coord) in enumerate(coords)
        setCoordSeq!(coordinates, i, coord, context)
    end
    coordinates
end

function getX(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSCoordSeq_getX_r(context.ptr, ptr, i-1, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getX")
    end
    out[]
end

function getX(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ncoords = getSize(ptr, context)
    xcoords = Array{Float64}(ncoords)
    start = pointer(xcoords)
    floatsize = sizeof(Float64)
    for i=0:ncoords-1
        GEOSCoordSeq_getX_r(context.ptr, ptr, i, start + i*floatsize)
    end
    xcoords
end

function getY(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSCoordSeq_getY_r(context.ptr, ptr, i-1, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getY")
    end
    out[]
end

function getY(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ncoords = getSize(ptr, context)
    ycoords = Array{Float64}(ncoords)
    start = pointer(ycoords)
    floatsize = sizeof(Float64)
    for i=0:ncoords-1
        GEOSCoordSeq_getY_r(context.ptr, ptr, i, start + i*floatsize)
    end
    ycoords
end

function getZ(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSCoordSeq_getZ_r(context.ptr, ptr, i-1, out)
    if result == 0
        error("LibGEOS: Error in GEOSCoordSeq_getZ")
    end
    out[]
end

function getZ(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ncoords = getSize(ptr, context)
    zcoords = Array{Float64}(ncoords)
    start = pointer(zcoords)
    floatsize = sizeof(Float64)
    for i=0:ncoords-1
        GEOSCoordSeq_getZ_r(context.ptr, ptr, i, start + i*floatsize)
    end
    zcoords
end

function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = _context)
    ndim = getDimensions(ptr, context)
    coord = Array{Float64}(undef, ndim)
    start = pointer(coord)
    floatsize = sizeof(Float64)
    GEOSCoordSeq_getX_r(context.ptr, ptr, i-1, start)
    GEOSCoordSeq_getY_r(context.ptr, ptr, i-1, start+floatsize)
    if ndim == 3
        GEOSCoordSeq_getZ_r(context.ptr, ptr, i-1, start+2*floatsize)
    end
    coord
end

function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    ndim = getDimensions(ptr, context)
    ncoords = getSize(ptr, context)
    coordseq = Vector{Float64}[]
    sizehint!(coordseq, ncoords)
    for i=1:ncoords
        push!(coordseq, getCoordinates(ptr, i, context))
    end
    coordseq
end

# -----
# Linear referencing functions -- there are more, but these are probably sufficient for most purposes
# -----
# (GEOSGeometry ownership is retained by caller)

# Return distance of point 'p' projected on 'g' from origin of 'g'. Geometry 'g' must be a lineal geometry
project(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = _context) =
    GEOSProject_r(context.ptr, g, p)

# Return closest point to given distance within geometry (Geometry must be a LineString)
function interpolate(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
    result = GEOSInterpolate_r(context.ptr, ptr, d)
    if result == C_NULL
        error("LibGEOS: Error in GEOSInterpolate")
    end
    result
end

projectNormalized(g::GEOSGeom, p::GEOSGeom, context::GEOSContext = _context) =
    GEOSProjectNormalized_r(context.ptr, g, p)

function interpolateNormalized(ptr::GEOSGeom, d::Real, context::GEOSContext = _context)
    result = GEOSInterpolateNormalized_r(context.ptr, ptr, d)
    if result == C_NULL
        error("LibGEOS: Error in GEOSInterpolateNormalized")
    end
    result
end

# -----
# Buffer related functions
# -----
# Computes the buffer of a geometry, for both positive and negative buffer distances.
# Since true buffer curves may contain circular arcs, computed buffer polygons can only be approximations to the true geometry.
# The user can control the accuracy of the curve approximation by specifying the number of linear segments with which to approximate a curve.

# Always returns a polygon. The negative or zero-distance buffer of lines and points is always an empty Polygon.
buffer(ptr::GEOSGeom, width::Real, quadsegs::Integer=8, context::GEOSContext = _context) =
    GEOSBuffer_r(context.ptr, ptr, width, Int32(quadsegs))

bufferWithStyle(ptr::GEOSGeom, width::Real, quadsegs::Integer=8, endCapStyle::GEOSBufCapStyles=GEOSBUF_CAP_ROUND, joinStyle::GEOSBufJoinStyles=GEOSBUF_JOIN_ROUND, mitreLimit::Real=5.0, context::GEOSContext = _context) =
    GEOSBufferWithStyle_r(context.ptr, ptr, width, Int32(quadsegs), Int32(endCapStyle), Int32(joinStyle), mitreLimit)

# GEOSBufferParams_create
# GEOSBufferParams_destroy
# GEOSBufferParams_setEndCapStyle
# GEOSBufferParams_setJoinStyle
# GEOSBufferParams_setMitreLimit
# GEOSBufferParams_setQuadrantSegments
# GEOSBufferParams_setSingleSided
# GEOSBufferWithParams
# GEOSOffsetCurve

# -----
# Geometry Constructors -- All functions return NULL on exception
# -----
# (GEOSCoordSequence* arguments will become ownership of the returned object.)

function createPoint(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSGeom_createPoint_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createPoint")
    end
    result
end
createPoint(x::Real, y::Real, context::GEOSContext = _context) =
    createPoint(createCoordSeq(x,y,context), context)
createPoint(x::Real, y::Real, z::Real, context::GEOSContext = _context) =
    createPoint(createCoordSeq(x,y,z,context), context)
createPoint(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
    createPoint(createCoordSeq(coords,context), context)
createPoint(coords::Vector{Float64}, context::GEOSContext = _context) =
    createPoint(createCoordSeq(Vector{Float64}[coords],context), context)

function createLinearRing(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSGeom_createLinearRing_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createLinearRing")
    end
    result
end
createLinearRing(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
    GEOSGeom_createLinearRing_r(context.ptr, createCoordSeq(coords, context))

function createLineString(ptr::GEOSCoordSeq, context::GEOSContext = _context)
    result = GEOSGeom_createLineString_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createLineString")
    end
    result
end
createLineString(coords::Vector{Vector{Float64}}, context::GEOSContext = _context) =
    GEOSGeom_createLineString_r(context.ptr, createCoordSeq(coords, context))

# Second argument is an array of GEOSGeometry* objects.
# The caller remains owner of the array, but pointed-to
# objects become ownership of the returned GEOSGeometry.
function createPolygon(shell::GEOSGeom, holes::Vector{GEOSGeom}, context::GEOSContext = _context)
    result = GEOSGeom_createPolygon_r(context.ptr, shell, pointer(holes), length(holes))
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createPolygon")
    end
    result
end

function createCollection(geomtype::GEOSGeomTypes, geoms::Vector{GEOSGeom}, context::GEOSContext = _context)
    result = GEOSGeom_createCollection_r(context.ptr, geomtype, pointer(geoms), length(geoms))
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createCollection")
    end
    result
end

function createEmptyCollection(geomtype::GEOSGeomTypes, context::GEOSContext = _context)
    result = GEOSGeom_createEmptyCollection_r(context.ptr, geomtype)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_createEmptyCollection")
    end
    result
end

# Memory management
function cloneGeom(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeom_clone_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_clone")
    end
    result
end

destroyGeom(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom_destroy_r(context.ptr, ptr)

# -----
# Topology operations - return NULL on exception.
# -----

function envelope(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSEnvelope_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSEnvelope")
    end
    result
end

"""
    minimumRotatedRectangle(geom)

Returns the minimum rotated rectangular POLYGON which encloses the input geometry. The rectangle
has width equal to the minimum diameter, and a longer length. If the convex hill of the input is
degenerate (a line or point) a LINESTRING or POINT is returned. The minimum rotated rectangle can
be used as an extremely generalized representation for the given geometry.
"""
function minimumRotatedRectangle(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSMinimumRotatedRectangle_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSMinimumRotatedRectangle")
    end
    result
end

function intersection(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSIntersection_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSIntersection")
    end
    result
end

# Returns a Geometry that represents the convex hull of the input geometry. The returned geometry contains the minimal number of points needed to represent the convex hull. In particular, no more than two consecutive points will be collinear.

# Returns:
# if the convex hull contains 3 or more points, a Polygon; 2 points, a LineString; 1 point, a Point; 0 points, an empty GeometryCollection.
function convexhull(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSConvexHull_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSConvexHull")
    end
    result
end

function difference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSDifference_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSDifference")
    end
    result
end

function symmetricDifference(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSSymDifference_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSymDifference")
    end
    result
end

function boundary(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSBoundary_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSBoundary")
    end
    result
end

function union(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSUnion_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSUnion")
    end
    result
end

function unaryUnion(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSUnaryUnion_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSUnaryUnion")
    end
    result
end

function pointOnSurface(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPointOnSurface_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSPointOnSurface")
    end
    result
end

function centroid(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetCentroid_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetCentroid")
    end
    result
end

function node(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSNode_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSNode")
    end
    result
end

# all arguments remain ownership of the caller (both Geometries and pointers)
function polygonize(geoms::Vector{GEOSGeom}, context::GEOSContext = _context)
    result = GEOSPolygonize_r(context.ptr, pointer(geoms), length(geoms))
    if result == C_NULL
        error("LibGEOS: Error in GEOSPolygonize")
    end
    result
end
# GEOSPolygonizer_getCutEdges
# GEOSPolygonize_full

function lineMerge(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSLineMerge_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSLineMerge")
    end
    result
end

function simplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSSimplify_r(context.ptr, ptr, tol)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSimplify")
    end
    result
end

function topologyPreserveSimplify(ptr::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSTopologyPreserveSimplify_r(context.ptr, ptr, tol)
    if result == C_NULL
        error("LibGEOS: Error in GEOSTopologyPreserveSimplify")
    end
    result
end

# Return all distinct vertices of input geometry as a MULTIPOINT.
# (Note that only 2 dimensions of the vertices are considered when testing for equality)
function uniquePoints(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeom_extractUniquePoints_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_extractUniquePoints")
    end
    result
end

# Find paths shared between the two given lineal geometries.
# Returns a GEOMETRYCOLLECTION having two elements:
# - first element is a MULTILINESTRING containing shared paths
#    having the _same_ direction on both inputs
#  - second element is a MULTILINESTRING containing shared paths
#    having the _opposite_ direction on the two inputs
# (Returns NULL on exception)
function sharedPaths(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSSharedPaths_r(context.ptr, g1, g2)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSharedPaths")
    end
    result
end

# Snap first geometry on to second with given tolerance
# (Returns a newly allocated geometry, or NULL on exception)
function snap(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSSnap_r(context.ptr, g1, g2, tol)
    if result == C_NULL
        error("LibGEOS: Error in GEOSSnap")
    end
    result
end

# Return a Delaunay triangulation of the vertex of the given geometry
#
# @param g the input geometry whose vertex will be used as "sites"
# @param tolerance optional snapping tolerance to use for improved robustness
# @param onlyEdges if non-zero will return a MULTILINESTRING, otherwise it will
#                  return a GEOMETRYCOLLECTION containing triangular POLYGONs.
#
# @return  a newly allocated geometry, or NULL on exception
function delaunayTriangulation(ptr::GEOSGeom, tol::Real=0.0, onlyEdges::Bool=false, context::GEOSContext = _context)
    result = GEOSDelaunayTriangulation_r(context.ptr, ptr, tol, Int32(onlyEdges))
    if result == C_NULL
        error("LibGEOS: Error in GEOSDelaunayTriangulation")
    end
    result
end

# -----
# Binary predicates - return 2 on exception, 1 on true, 0 on false
# -----
function disjoint(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSDisjoint_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSDisjoint")
    end
    result != 0x00
end

function touches(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSTouches_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSTouches")
    end
    result != 0x00
end

function intersects(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSIntersects_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSIntersects")
    end
    result != 0x00
end

function crosses(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSCrosses_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSCrosses")
    end
    result != 0x00
end

function within(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSWithin_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSWithin")
    end
    result != 0x00
end

function contains(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSContains_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSContains")
    end
    result != 0x00
end

function overlaps(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSOverlaps_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSOverlaps")
    end
    result != 0x00
end

function equals(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSEquals_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSEquals")
    end
    result != 0x00
end

function equalsexact(g1::GEOSGeom, g2::GEOSGeom, tol::Real, context::GEOSContext = _context)
    result = GEOSEqualsExact_r(context.ptr, g1, g2, tol)
    if result == 0x02
        error("LibGEOS: Error in GEOSEqualsExact")
    end
    result != 0x00
end

function covers(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSCovers_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSCovers")
    end
    result != 0x00
end

function coveredby(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSCoveredBy_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSCoveredBy")
    end
    result != 0x00
end

# -----
# Prepared Geometry Binary predicates - return 2 on exception, 1 on true, 0 on false
# -----

# GEOSGeometry ownership is retained by caller
function prepareGeom(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPrepare_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSPrepare")
    end
    result
end

destroyPreparedGeom(ptr::Ptr{GEOSPreparedGeometry}, context::GEOSContext = _context) =
    GEOSPreparedGeom_destroy_r(context.ptr, ptr)

function prepcontains(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedContains_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedContains")
    end
    result != 0x00
end

function prepcontainsproperly(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedContainsProperly_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedContainsProperly")
    end
    result != 0x00
end

function prepcoveredby(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedCoveredBy_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedCoveredBy")
    end
    result != 0x00
end

function prepcovers(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedCovers_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedCovers")
    end
    result != 0x00
end

function prepcrosses(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedCrosses_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedCrosses")
    end
    result != 0x00
end

function prepdisjoint(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedDisjoint_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedDisjoint")
    end
    result != 0x00
end

function prepintersects(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedIntersects_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedIntersects")
    end
    result != 0x00
end

function prepoverlaps(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedOverlaps_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedOverlaps")
    end
    result != 0x00
end

function preptouches(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedTouches_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedTouches")
    end
    result != 0x00
end

function prepwithin(g1::Ptr{GEOSPreparedGeometry}, g2::GEOSGeom, context::GEOSContext = _context)
    result = GEOSPreparedWithin_r(context.ptr, g1, g2)
    if result == 0x02
        error("LibGEOS: Error in GEOSPreparedWithin")
    end
    result != 0x00
end

# -----
# STRtree functions
# -----
# GEOSSTRtree_create
# GEOSSTRtree_insert
# GEOSSTRtree_query
# GEOSSTRtree_iterate
# GEOSSTRtree_remove
# GEOSSTRtree_destroy

# -----
# Unary predicate - return 2 on exception, 1 on true, 0 on false
# -----
function isEmpty(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisEmpty_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisEmpty")
    end
    result != 0x00
end

function isSimple(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisSimple_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisSimple")
    end
    result != 0x00
end

function isRing(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisRing_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisRing")
    end
    result != 0x00
end

function hasZ(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSHasZ_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSHasZ")
    end
    result != 0x00
end

# Call only on LINESTRING (return 2 on exception, 1 on true, 0 on false)
function isClosed(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisClosed_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisClosed")
    end
    result != 0x00
end

# -----
# Dimensionally Extended 9 Intersection Model related
# -----

# GEOSRelatePattern (return 2 on exception, 1 on true, 0 on false)
# GEOSRelate (return NULL on exception, a string to GEOSFree otherwise)
# GEOSRelatePatternMatch (return 2 on exception, 1 on true, 0 on false)
# GEOSRelateBoundaryNodeRule (return NULL on exception, a string to GEOSFree otherwise)

# -----
# Validity checking -- return 2 on exception, 1 on true, 0 on false
# -----

# /* These are for use with GEOSisValidDetail (flags param) */
# enum GEOSValidFlags {
#     GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE=1
# };

function isValid(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSisValid_r(context.ptr, ptr)
    if result == 0x02
        error("LibGEOS: Error in GEOSisValid")
    end
    result != 0x00
end

# * return NULL on exception, a string to GEOSFree otherwise
# GEOSisValidReason

# * Caller has the responsibility to destroy 'reason' (GEOSFree)
# * and 'location' (GEOSGeom_destroy) params
# * return 2 on exception, 1 when valid, 0 when invalid
# GEOSisValidDetail

# -----
# Geometry info
# -----

# Return NULL on exception, result must be freed by caller
# function geomType(ptr::GEOSGeom)
#     result = GEOSGeomType(ptr)
#     if result == C_NULL
#         error("LibGEOS: Error in GEOSGeomType")
#     end
#     result
# end

# Return -1 on exception
function geomTypeId(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomTypeId_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGeomTypeId")
    end
    result
end

# Return 0 on exception
function getSRID(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetSRID_r(context.ptr, ptr)
    if result == 0
        error("LibGEOS: Error in GEOSGeomTypeId")
    end
    result
end

setSRID(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSSetSRID_r(context.ptr, ptr)

# May be called on all geometries in GEOS 3.x, returns -1 on error and 1
# for non-multi geometries. Older GEOS versions only accept
# GeometryCollections or Multi* geometries here, and are likely to crash
# when fed simple geometries, so beware if you need compatibility with
# old GEOS versions.
function numGeometries(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetNumGeometries_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGeomTypeId")
    end
    result
end

# Call only on GEOMETRYCOLLECTION or MULTI*
# (Return a pointer to the internal Geometry. Return NULL on exception.)
# Returned object is a pointer to internal storage:
# it must NOT be destroyed directly.
# Up to GEOS 3.2.0 the input geometry must be a Collection, in
# later version it doesn't matter (i.e. getGeometryN(0) for a single will return the input).
function getGeometry(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
    result = GEOSGetGeometryN_r(context.ptr, ptr, n-1)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetGeometryN")
    end
    cloneGeom(result, context)
end
getGeometries(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom[getGeometry(ptr, i, context) for i=1:numGeometries(ptr, context)]

# Converts Geometry to normal form (or canonical form).
# Return -1 on exception, 0 otherwise.
function normalize!(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSNormalize_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSNormalize")
    end
    result
end

# Return -1 on exception
function numInteriorRings(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetNumInteriorRings_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGetNumInteriorRings")
    end
    result
end

# Call only on LINESTRING (returns -1 on exception)
function numPoints(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomGetNumPoints_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGeomGetNumPoints")
    end
    result
end

# Return -1 on exception, Geometry must be a Point.
function getGeomX(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeomGetX_r(context.ptr, ptr, out)
    if result == -1
        error("LibGEOS: Error in GEOSGeomGetX")
    end
    out[]
end

function getGeomY(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    result = GEOSGeomGetY_r(context.ptr, ptr, out)
    if result == -1
        error("LibGEOS: Error in GEOSGeomGetY")
    end
    out[]
end

# Return NULL on exception, Geometry must be a Polygon.
# Returned object is a pointer to internal storage: it must NOT be destroyed directly.
function interiorRing(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
    result = GEOSGetInteriorRingN_r(context.ptr, ptr, n-1)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetInteriorRingN")
    end
    cloneGeom(result, context)
end

function interiorRings(ptr::GEOSGeom, context::GEOSContext = _context)
    n = numInteriorRings(ptr, context)
    if n == 0
        return GEOSGeom[]
    else
        return GEOSGeom[interiorRing(ptr, i, context) for i=1:n]
    end
end

# Return NULL on exception, Geometry must be a Polygon.
# Returned object is a pointer to internal storage: it must NOT be destroyed directly.
function exteriorRing(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetExteriorRing_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGetExteriorRing")
    end
    cloneGeom(result, context)
end

# Return -1 on exception
function numCoordinates(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGetNumCoordinates_r(context.ptr, ptr)
    if result == -1
        error("LibGEOS: Error in GEOSGetNumCoordinates")
    end
    result
end

# Geometry must be a LineString, LinearRing or Point (Return NULL on exception)
function getCoordSeq(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeom_getCoordSeq_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_getCoordSeq")
    end
    result
end
# getGeomCoordinates(ptr::GEOSGeom) = getCoordinates(getCoordSeq(ptr))

# Return 0 on exception (or empty geometry)
getGeomDimensions(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom_getDimensions_r(context.ptr, ptr)

# Return 2 or 3.
getCoordinateDimension(ptr::GEOSGeom, context::GEOSContext = _context) =
    GEOSGeom_getCoordinateDimension_r(context.ptr, ptr)

# Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
function getPoint(ptr::GEOSGeom, n::Integer, context::GEOSContext = _context)
    result = GEOSGeomGetPointN_r(context.ptr, ptr, n-1)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeomGetPointN")
    end
    result
end

# Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
function startPoint(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomGetStartPoint_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeomGetStartPoint")
    end
    result
end

# Call only on LINESTRING, and must be freed by caller (Returns NULL on exception)
function endPoint(ptr::GEOSGeom, context::GEOSContext = _context)
    result = GEOSGeomGetEndPoint_r(context.ptr, ptr)
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeomGetEndPoint")
    end
    result
end

# -----
# Misc functions
# -----
function geomArea(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSArea_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSArea")
    end
    out[]
end

function geomLength(ptr::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSLength_r(context.ptr, ptr, out)
    if result == 0
        error("LibGEOS: Error in GEOSLength")
    end
    out[]
end

function geomDistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSDistance_r(context.ptr, g1, g2, out)
    if result == 0
        error("LibGEOS: Error in GEOSDistance")
    end
    out[]
end

function hausdorffdistance(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSHausdorffDistance_r(context.ptr, g1, g2, out)
    if result == 0
        error("LibGEOS: Error in GEOSHausdorffDistance")
    end
    out[]
end

function hausdorffdistance(g1::GEOSGeom, g2::GEOSGeom, densifyFrac::Real, context::GEOSContext = _context)
    out = Ref{Float64}()
    # Return 0 on exception, 1 otherwise
    result = GEOSHausdorffDistanceDensify_r(context.ptr, g1, g2, densifyFrac, out)
    if result == 0
        error("LibGEOS: Error in GEOSHausdorffDistanceDensify")
    end
    out[]
end

# Return 0 on exception, the closest points of the two geometries otherwise.
# The first point comes from g1 geometry and the second point comes from g2.
nearestPoints(g1::GEOSGeom, g2::GEOSGeom, context::GEOSContext = _context) =
    GEOSNearestPoints_r(context.ptr, g1, g2)

# -----
# Precision functions
# -----

"""
    getPrecision(geom)

Return the size of the geometry's precision grid, 0 for FLOATING precision.
"""
function getPrecision(geom::GEOSGeom, context::GEOSContext = _context)
    GEOSGeom_getPrecision_r(context.ptr, geom)
end

"""
    setPrecision(geom, gridSize; flags = 0)

Set the geometry's precision, optionally rounding all its coordinates to the
precision grid (if it changes).

Note that operations will always be performed in the precision
of the geometry with higher precision (smaller `gridSize`).
That same precision will be attached to the operation outputs.

* `gridSize` size of the precision grid, or 0 for FLOATING precision.
* `flags` The bitwise OR of one of more of the "precision options"
  * `LibGEOS.GEOS_PREC_KEEP_COLLAPSED` retain collapsed elements
  * `LibGEOS.GEOS_PREC_NO_TOPO` do not attempt at preserving the topology
"""
function setPrecision(geom::GEOSGeom, gridSize::Real, flags::Int, context::GEOSContext = _context)
    result = geomFromGEOS(GEOSGeom_setPrecision_r(context.ptr, geom, gridSize, flags))
    if result == C_NULL
        error("LibGEOS: Error in GEOSGeom_setPrecision_r")
    end
    result
end


text:
using Plots
using Reproduce
using JLD2
using FileIO

p = plot()
# folder = "TabularTMazeExperiment/RP_0_0x278a9a4b9fa34c2b/"
folder = "TabularTMazeExperiment/RP_0_0xf91eb689cd1d55f6/"
save_file = "plotting/chunlok/generated_plots/test_plot_error_map"

results_file = folder * "results.jld2"
@load results_file results
# println(results)

settings_file = folder * "settings.jld2"
println(settings_file)
settings = FileIO.load(settings_file)["parsed_args"]

# print(settings)
# print(sett)

error_map_data = results[:ttmaze_error_map]
eval_set = results[:ttmaze_error_map_eval_set]
all_q_ests = results[:ttmaze_error_map_est]

function valid_state_mask()
    """
    return a mask of valid states that is 9x7
    """
    world = [["G1", "0", "0", "0", "0", "0", "G3"],
             ["1", "0", "0", "0", "0", "0", "1"],
             ["1", "1", "1", "1", "1", "1", "1"],
             ["1", "0", "0", "1", "0", "0", "1"],
             ["G2", "0", "0", "1", "0", "0", "G4"],
             ["0", "0", "0", "1", "0", "0", "0"],
             ["0", "0", "0", "1", "0", "0", "0"],
             ["0", "0", "0", "1", "0", "0", "0"],
             ["0", "0", "0", "1", "0", "0", "0"]]

    world = permutedims(hcat(world...))
    valid_states = findall(x-> x!="0", world)
    return valid_states
end

state_map = valid_state_mask()

# @gif for time_step in 2:40
    # Last error map of the 3rd GVF (drifter)
    time_step = 40
    errors = error_map_data[time_step]
    q_ests = all_q_ests[time_step]

    # println(keys(eval_set))

    true_values = eval_set["ests"]
    actions = eval_set["actions"]

    imgs = []
    for action in 1:4
        error_map = zeros(9, 7)
        for (ind, obs) in enumerate(eval_set["states"])

            index_state = convert(Int32, obs[1])
            state = state_map[index_state]

            # println(size(ests))
            # println(size(data))

            error = errors[:, ind]
            true_value = true_values[:, ind]
            q_est = q_ests[:, ind]

            # println(size(q_ests))

            # println(size(error))
            gvf = 4
            # println(q_est)

            if (actions[ind] == action)
                error_map[state] += error[gvf]
                # error_map[state] += q_est[gvf]
                # error_map[state] += true_value[gvf]
            end
        end

        action_names = ["up", "right", "down", "left"]
        push!(imgs, heatmap(error_map, yflip=true, title="$(action_names[action]) $(time_step)"))
    end
    plot(imgs..., layout=4)
# end

# heatmap!(hcat(imgs...), yflip=true, layout=4)

savefig(save_file)

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# GCM Initial Perturbation
# This file contains helpers and lists currely avaiable options

using Distributions
using Random
using Distributions: Uniform
using Random: rand

abstract type AbstractPerturbation end
struct NoPerturbation <: AbstractPerturbation end
struct DeterministicPerturbation <: AbstractPerturbation end
struct RandomPerturbation <: AbstractPerturbation end

# Helper for parsing `--init-perturbation` command line argument
function parse_perturbation_arg(arg)
    if arg === nothing
        perturbation = nothing
    elseif arg == "deterministic"
        perturbation = DeterministicPerturbation()
    elseif arg == "zero"
        perturbation = NoPerturbation()
    elseif arg == "random"
        perturbation = RandomPerturbation()
    else
        error("unknown perturbation: " * arg)
    end

    return perturbation
end

function init_perturbation(::NoPerturbation, bl, state, aux, coords, t)
    FT = eltype(state)

    u′, v′, w′ = (FT(0), FT(0), FT(0))
    rand_pert = FT(1)

    return u′, v′, w′, rand_pert
end

# Velocity perturbation following
# Ullrich et al. (2016) Dynamical Core Model Intercomparison Project (DCMIP2016) Test Case Document
function init_perturbation(
    ::DeterministicPerturbation,
    bl,
    state,
    aux,
    coords,
    t,
)
    FT = eltype(state)

    # get parameters
    param_set = parameter_set(bl)
    _a::FT = planet_radius(param_set)

    φ = latitude(bl, aux)
    λ = longitude(bl, aux)
    z = altitude(bl, aux)

    # perturbation specific parameters
    z_t::FT = 15e3
    λ_c::FT = π / 9
    φ_c::FT = 2 * π / 9
    d_0::FT = _a / 6
    V_p::FT = 10

    F_z::FT = 1 - 3 * (z / z_t)^2 + 2 * (z / z_t)^3
    if z > z_t
        F_z = FT(0)
    end
    d::FT = _a * acos(sin(φ) * sin(φ_c) + cos(φ) * cos(φ_c) * cos(λ - λ_c))
    c3::FT = cos(π * d / 2 / d_0)^3
    s1::FT = sin(π * d / 2 / d_0)
    if 0 < d < d_0 && d != FT(_a * π)
        u′::FT =
            -16 * V_p / 3 / sqrt(3) *
            F_z *
            c3 *
            s1 *
            (-sin(φ_c) * cos(φ) + cos(φ_c) * sin(φ) * cos(λ - λ_c)) /
            sin(d / _a)
        v′::FT =
            16 * V_p / 3 / sqrt(3) * F_z * c3 * s1 * cos(φ_c) * sin(λ - λ_c) /
            sin(d / _a)
    else
        u′ = FT(0)
        v′ = FT(0)
    end
    w′ = FT(0)
    rand_pert = FT(1)

    return u′, v′, w′, rand_pert
end

function init_perturbation(::RandomPerturbation, bl, state, aux, coords, t)
    FT = eltype(state)
    u′, v′, w′ = (FT(0), FT(0), FT(0))
    rand_pert = FT(1.0 + rand(Uniform(-1e-3, 1e-3)))

    return u′, v′, w′, rand_pert
end


text:
using PushGP: hasstack, stackhas, pop2!

@testset "PushInterpreter" begin

@testset "hasstack" begin
    interp = PushInterpreter()
    @test hasstack(interp, Int)
    @test hasstack(interp, Float64)
    @test hasstack(interp, String)
    @test hasstack(interp, Bool)
    @test !hasstack(interp, Vector{Int})
    @test !hasstack(interp, Regex)
end # @testset "hasstack" begin

@testset "stackhas, execute! and pop2!" begin
    interp = PushInterpreter()
    for T in [Int, Float64, String, Bool]
        @test stackhas(interp, T, 0) == true
        for numelems in 1:10
            @test stackhas(interp, T, numelems) == false
        end
    end

    execute!(interp, [1]) # Should add one to the Int stack
    @test stackhas(interp, Int, 0) == true
    @test stackhas(interp, Int, 1) == true
    @test stackhas(interp, Int, 2) == false

    execute!(interp, [4]) # Should add one to the Int stack, now 2 elements on it
    @test stackhas(interp, Int, 0) == true
    @test stackhas(interp, Int, 1) == true
    @test stackhas(interp, Int, 2) == true
    @test stackhas(interp, Int, 3) == false

    a, b = pop2!(interp, Int)
    @test a == 4
    @test b == 1
end # @testset "stackhas" begin

end # @testset "PushInterpreter" begin


text:
function ground(video, action, framenumber)
    background("white")
    sethue("blue")
    return framenumber
end

function ground_color(color_bg, color_pen, framenumber)
    background(color_bg)
    sethue(color_pen)
    return framenumber
end

function latex_title(args...)
    fontsize(20)
    latex(L"E=mc^2", 0, -200)
end

function circ(p = O, color = "black")
    sethue(color)
    circle(p, 25, :fill)
    return p
end

function circ_ret_trans(p = O, color = "black")
    sethue(color)
    circle(p, 25, :fill)
    return Transformation(p, 0.0)
end

function path!(points, pos, color)
    sethue(color)
    push!(points, pos)
    circle.(points, 2, :fill)
end

function rad(p1, p2, color)
    sethue(color)
    line(p1, p2, :stroke)
end

@testset "Dancing circles (gif)" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:25,
                ground,
                Rotation(0.0),
                Translation(Point(25, 25), Point(25, 25)),
            ),
            Action(latex_title),
            Action(
                Rel(-24:0),
                :red_ball,
                (args...) -> circ(p1, "red"),
                Rotation(from_rot, to_rot),
            ),
            Action(
                1:25,
                :blue_ball,
                (args...) -> circ(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(
                1:25,
                (video, args...) -> path!(path_of_red, get_position(:red_ball), "red"),
            ),
            Action(:same, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:25, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "dancing.gif",
    )

    @test_reference "refs/dancing_circles_16.png" load("images/0000000016.png")
    @test isfile("dancing.gif")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
    rm("images/palette.bmp")
    rm("dancing.gif")
end

@testset "Dancing circles (mp4)" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:30,
                ground,
                Rotation(0.0),
                Translation(Point(25, 25), Point(25, 25)),
            ),
            Action(latex_title),
            Action(:red_ball, (args...) -> circ(p1, "red"), Rotation(from_rot, to_rot)),
            Action(
                1:30,
                :blue_ball,
                (args...) -> circ(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(
                1:30,
                (video, args...) -> path!(path_of_red, get_position(:red_ball), "red"),
            ),
            Action(:same, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:30, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "dancing.mp4",
        framerate = 1,
    )

    # 30 frames with a framerate of 1 should take about 30 seconds ;)
    @test isapprox(VideoIO.get_duration("dancing.mp4"), 30.0, atol = 0.1)
    rm("dancing.mp4")
end

@testset "Dancing circles layered" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            Action(
                1:25,
                ground,
                Rotation(π / 2, π / 2, O),
                Translation(Point(25, 25), Point(25, 25));
                in_global_layer = true,
            ),
            Action(1:25, latex_title),
            Action(1:25, :red_ball, (args...) -> circ(p1, "red"), Rotation(to_rot)),
            Action(
                1:25,
                :blue_ball,
                (args...) -> circ(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(
                1:25,
                (video, args...) -> path!(path_of_red, get_position(:red_ball), "red"),
            ),
            Action(1:25, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:25, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Dancing circles layered return Transformation" begin
    p1 = Point(100, 0)
    p2 = Point(100, 80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(
        video,
        [
            Action(
                1:25,
                ground,
                Rotation(π / 2, π / 2, O),
                Translation(Point(25, 25), Point(25, 25));
                in_global_layer = true,
            ),
            Action(1:25, latex_title),
            Action(
                1:25,
                :red_ball,
                (args...) -> circ_ret_trans(p1, "red"),
                Rotation(to_rot),
            ),
            Action(
                1:25,
                :blue_ball,
                (args...) -> circ_ret_trans(p2, "blue"),
                Rotation(to_rot, from_rot, :red_ball),
            ),
            Action(1:25, (video, args...) -> path!(path_of_red, pos(:red_ball), "red")),
            Action(1:25, (video, args...) -> path!(path_of_blue, pos(:blue_ball), "blue")),
            Action(1:25, (args...) -> rad(pos(:red_ball), pos(:blue_ball), "black")),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i in 1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Drawing grid" begin
    video = Video(500, 500)
    javis(
        video,
        [
            Action(1:40, ground),
            Action(1:10, draw_grid(direction = "BL", line_gap = 25)),
            Action(zero_lines(direction = "BL", line_thickness = 10)),
            Action(Rel(10), draw_grid(direction = "BR", line_gap = 25)),
            Action(zero_lines(direction = "BR", line_thickness = 10)),
            Action(Rel(10), draw_grid(direction = "TL", line_gap = 25)),
            Action(zero_lines(direction = "TL", line_thickness = 10)),
            Action(Rel(10), draw_grid(direction = "TR", line_gap = 25)),
            Action(zero_lines(direction = "TR", line_thickness = 10)),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/grid_drawing_bl.png" load("images/0000000008.png")
    @test_reference "refs/grid_drawing_br.png" load("images/0000000018.png")
    @test_reference "refs/grid_drawing_tl.png" load("images/0000000028.png")
    @test_reference "refs/grid_drawing_tr.png" load("images/0000000038.png")
    for i in 1:40
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Use returned angle" begin
    video = Video(500, 500)
    p = Point(100, 0)
    javis(
        video,
        [
            Action(1:10, ground),
            Action(:circ, (args...) -> circ_ret_trans(p), Rotation(2π)),
            Action(
                (args...) -> line(Point(-200, 0), Point(-200, -10 * ang(:circ)), :stroke),
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/circle_angle.png" load("images/0000000008.png")
    for i in 1:10
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

astar(args...) = star(Point(-100, -100), 30)
acirc(args...) = circle(Point(100, 100), 30)

@testset "morphing star2circle and back" begin
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:20,
                :framenumber,
                (args...) -> ground_color("white", "black", args[3]),
            ),
            Action(1:10, (args...) -> circle(Point(-100, 0), val(:framenumber), :fill)),
            Action(1:10, morph(astar, acirc)),
            Action(11:20, morph(acirc, astar)),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/star2circle5.png" load("images/0000000005.png")
    @test_reference "refs/star2circle15.png" load("images/0000000015.png")
    for i in 1:20
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "morphing star2circle and back with fill" begin
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(
                1:20,
                :framenumber,
                (args...) -> ground_color("white", "black", args[3]),
            ),
            Action(1:10, (args...) -> circle(Point(-100, 0), val(:framenumber), :fill)),
            Action(1:10, morph(astar, acirc; action = :fill)),
            Action(11:20, morph(acirc, astar; action = :fill)),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/star2circle_fill5.png" load("images/0000000005.png")
    @test_reference "refs/star2circle_fill15.png" load("images/0000000015.png")
    for i in 1:20
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_nicholas(args...)
    background("white")
    sethue("blue")
    setline(3)
end

function house_of_nicholas(; p1 = O, width = 100, color = "black")
    # sethue(color)
    #     .p1
    # .p2   .p3
    #
    # .p4   .p5
    width2 = div(width, 2)
    p2 = p1 + Point(-width2, width2)
    p3 = p1 + Point(width2, width2)
    p4 = p2 + Point(0, width)
    p5 = p3 + Point(0, width)
    line(p4, p5, :stroke)
    line(p5, p2, :stroke)
    line(p2, p4, :stroke)
    line(p4, p3, :stroke)
    line(p3, p1, :stroke)
    line(p1, p2, :stroke)
    line(p2, p3, :stroke)
    setline(8)
    line(p3, p5, :stroke)
end

@testset "House of Nicholas line_width" begin
    demo = Video(500, 500)
    javis(
        demo,
        [
            BackgroundAction(1:50, ground_nicholas),
            Action(
                (args...) -> house_of_nicholas();
                subactions = [
                    SubAction(1:25, appear(:fade_line_width)),
                    SubAction(Rel(25), disappear(:fade_line_width)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/nicholas15.png" load("images/0000000015.png")
    @test_reference "refs/nicholas25.png" load("images/0000000025.png")
    @test_reference "refs/nicholas35.png" load("images/0000000035.png")
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Squeezing a circle using scale" begin
    demo = Video(500, 500)
    javis(
        demo,
        [
            BackgroundAction(1:100, ground),
            Action(:start_scale, (args...) -> (1.0, 1.0)),
            Action(
                (args...) -> circ();
                subactions = [
                    SubAction(1:25, Scaling((1.0, 1.5))),
                    SubAction(Rel(25), Scaling((2.0, 1.0))),
                    SubAction(Rel(25), Scaling(:start_scale)),
                    SubAction(Rel(25), Scaling(2.0)),
                ],
            ),
            Action((args...) -> circ(Point(-100, 0))),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/squeeze15.png" load("images/0000000015.png")
    @test_reference "refs/squeeze25.png" load("images/0000000025.png")
    @test_reference "refs/squeeze35.png" load("images/0000000035.png")
    @test_reference "refs/squeeze65.png" load("images/0000000065.png")
    @test_reference "refs/squeeze85.png" load("images/0000000085.png")
    for i in 1:100
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_opacity(args...)
    background("white")
    sethue("black")
    setopacity(0.5)
end

function square_opacity(p1, w)
    setopacity(1.0)
    sethue("red")
    rect(p1, w, w, :fill)
end

@testset "Circle/square appear opacity" begin
    demo = Video(500, 500)
    javis(
        demo,
        [
            BackgroundAction(1:50, ground_opacity),
            Action(
                1:42,
                (args...) -> circ();
                subactions = [
                    SubAction(1:25, appear(:fade)),
                    SubAction(26:42, disappear(:fade)),
                ],
            ),
            Action(
                5:50,
                (args...) -> square_opacity(Point(-100, 0), 60);
                subactions = [
                    SubAction(1:15, linear(), appear(:fade)),
                    SubAction(Rel(20), linear(), Translation(100, 50)),
                    SubAction(Rel(5), disappear(:fade)),
                    # for global frames 46-50 it should still be disappeared
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/circleSquare07opacity.png" load("images/0000000007.png")
    @test_reference "refs/circleSquare25opacity.png" load("images/0000000025.png")
    @test_reference "refs/circleSquare42opacity.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) ==
          RGB{Float64}(500 * 500, 500 * 500, 500 * 500)
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Animations.jl translate()" begin
    video = Video(500, 500)
    circle_anim = Animation(
        [0, 0.3, 0.6, 1], # must go from 0 to 1
        [O, Point(150, 0), Point(150, 150), O],
        [sineio(), polyin(5), expin(8)],
    )
    javis(
        video,
        [
            BackgroundAction(1:150, ground),
            Action(
                (args...) -> circle(O, 25, :fill);
                subactions = [SubAction(1:150, circle_anim, translate())],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )
    @test_reference "refs/anim_circle020.png" load("images/0000000020.png")
    @test_reference "refs/anim_circle075.png" load("images/0000000075.png")
    @test_reference "refs/anim_circle142.png" load("images/0000000142.png")
    for i in 1:150
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Animations.jl @warn" begin
    video = Video(500, 500)
    circle_anim = Animation(
        [0, 0.3, 0.6, 1.2], # must go from 0 to 1
        [O, Point(150, 0), Point(150, 150), O],
        [sineio(), polyin(5), expin(8)],
    )
    # warning as animation goes to 1.2 but should go to 1.0
    @test_logs (:warn,) (:warn,) javis(
        video,
        [
            BackgroundAction(1:2, ground),
            Action(
                (args...) -> circle(O, 25, :fill);
                subactions = [SubAction(1:2, circle_anim, translate())],
            ),
        ],
        pathname = "",
    )
end

@testset "Animations.jl rotate, scale, translate" begin
    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()],
    )

    color_anim = Animation(
        [0, 0.5, 1], # must go from 0 to 1
        [Lab(colorant"red"), Lab(colorant"cyan"), Lab(colorant"black")],
        [sineio(), sineio()],
    )

    javis(
        video,
        [
            BackgroundAction(1:150, ground),
            Action(
                (args...) -> circle(O, 25, :fill);
                subactions = [
                    SubAction(1:10, sineio(), scale()),
                    SubAction(11:50, translate_anim, translate()),
                    SubAction(51:100, rotate_anim, rotate_around(Point(-150, 0))),
                    SubAction(101:140, translate_back_anim, translate()),
                    SubAction(141:150, rev(sineio()), scale()),
                    SubAction(1:150, color_anim, sethue()),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/animations_all_05.png" load("images/0000000005.png")
    @test_reference "refs/animations_all_25.png" load("images/0000000025.png")
    @test_reference "refs/animations_all_65.png" load("images/0000000065.png")
    @test_reference "refs/animations_all_125.png" load("images/0000000125.png")
    @test_reference "refs/animations_all_145.png" load("images/0000000145.png")

    for i in 1:150
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Rotate around center animation" begin
    rotate_anim = Animation([0.0, 1.0], [0.0, 2π], [sineio()])

    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(1:50, ground),
            BackgroundAction(1:50, (args...) -> scaleto(2)),
            Action(
                (args...) -> circ(Point(75, 0)),
                subactions = [SubAction(1:50, rotate_anim, rotate())],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/rotate_center25.png" load("images/0000000025.png")
    @test_reference "refs/rotate_center45.png" load("images/0000000045.png")

    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Scaling circle" begin
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(1:50, ground),
            Action(:set_scale, (args...) -> 2),
            Action(
                (args...) -> circ(),
                subactions = [
                    SubAction(1:15, Scaling(0.0, :set_scale)),
                    SubAction(36:50, Scaling(0.0)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/scalingCircle07.png" load("images/0000000007.png")
    @test_reference "refs/scalingCircle25.png" load("images/0000000025.png")
    @test_reference "refs/scalingCircle42.png" load("images/0000000042.png")

    # test using appear and disappear
    video = Video(500, 500)
    javis(
        video,
        [
            BackgroundAction(1:50, ground),
            BackgroundAction(1:50, (args...) -> scale(2)),
            Action(
                (args...) -> circ(),
                subactions = [
                    SubAction(1:15, appear(:scale)),
                    SubAction(36:50, disappear(:scale)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/scalingCircle07.png" load("images/0000000007.png")
    @test_reference "refs/scalingCircle25.png" load("images/0000000025.png")
    @test_reference "refs/scalingCircle42.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) ==
          RGB{Float64}(500 * 500, 500 * 500, 500 * 500)
    for i in 1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "animating text" begin
    # does only test that it doesn't fail but I wasn't able to test this on all platforms
    # with using reference tests
    video = Video(400, 300)

    javis(
        video,
        [
            BackgroundAction(1:100, ground),
            BackgroundAction(1:100, (args...) -> fontsize(30)),
            Action(
                1:100,
                (args...) -> text("Hello Stream!", -50, 50; halign = :centre);
                subactions = [
                    SubAction(1:15, sineio(), appear(:draw_text)),
                    SubAction(76:100, sineio(), disappear(:draw_text)),
                ],
            ),
            Action(
                1:100,
                (args...) ->
                    text("Hello World!", -50, -100; halign = :wrong, valign = :wrong);
                subactions = [
                    SubAction(1:15, sineio(), appear(:draw_text)),
                    SubAction(76:100, sineio(), disappear(:draw_text)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    img07 = load("images/$(lpad(7, 10, "0")).png")
    img30 = load("images/$(lpad(30, 10, "0")).png")
    img82 = load("images/$(lpad(30, 10, "0")).png")

    # does only test that it doesn't fail but I wasn't able to test this on all platforms
    # with using reference tests
    video = Video(400, 300)

    javis(
        video,
        [
            BackgroundAction(1:100, ground),
            BackgroundAction(1:100, (args...) -> fontsize(30)),
            Action(
                1:100,
                (args...) -> text("Hello Stream!", -50, 50; halign = :center);
                subactions = [
                    SubAction(1:15, sineio(), appear(:draw_text)),
                    SubAction(76:100, sineio(), disappear(:draw_text)),
                ],
            ),
            Action(
                1:100,
                (args...) -> text("Hello World!", -50, -100);
                subactions = [
                    SubAction(1:15, sineio(), appear(:draw_text)),
                    SubAction(76:100, sineio(), disappear(:draw_text)),
                ],
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )

    img_other07 = load("images/$(lpad(7, 10, "0")).png")
    img_other30 = load("images/$(lpad(30, 10, "0")).png")
    img_other82 = load("images/$(lpad(30, 10, "0")).png")

    @test img07 == img_other07
    @test img30 == img_other30
    @test img82 == img_other82

    for i in 1:100
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Following a path" begin
    video = Video(800, 600)

    anim = Animation([0, 1], [0.0, 2.0], [sineio()])

    color_anim = Animation(
        [0, 0.5, 1], # must go from 0 to 1
        [Lab(colorant"red"), Lab(colorant"cyan"), Lab(colorant"red")],
        [sineio(), sineio()],
    )

    actions = [
        Action(
            frame_start:(frame_start + 149),
            (args...) -> star(O, 20, 5, 0.5, 0, :fill);
            subactions = [
                SubAction(1:150, anim, follow_path(star(O, 100))),
                SubAction(1:150, color_anim, sethue()),
            ],
        ) for frame_start in 1:7:22
    ]

    javis(
        video,
        [BackgroundAction(1:180, ground), actions...],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/followPath10.png" load("images/0000000010.png")
    @test_reference "refs/followPath30.png" load("images/0000000030.png")
    @test_reference "refs/followPath100.png" load("images/0000000100.png")
    @test_reference "refs/followPath160.png" load("images/0000000160.png")

    # Following along a bezier path
    function simple_bezier()
        P1 = Point(-300, 0)
        CP1 = Point(-200, -200)
        CP2 = Point(200, -200)
        P2 = Point(300, 0)

        beziersegment = BezierPathSegment(P1, CP1, CP2, P2)
        beziertopoly(beziersegment)
    end

    video = Video(800, 600)

    anim = Animation([0, 1], [0.0, 1.0], [sineio()])

    color_anim = Animation(
        [0, 0.5, 1], # must go from 0 to 1
        [Lab(colorant"red"), Lab(colorant"cyan"), Lab(colorant"red")],
        [sineio(), sineio()],
    )

    actions = [
        Action(
            frame_start:(frame_start + 149),
            (args...) -> star(O, 20, 5, 0.5, 0, :fill);
            subactions = [
                SubAction(1:1, Translation(simple_bezier()[1] + Point(0, 3 * frame_start))),
                SubAction(1:10, appear(:fade)),
                SubAction(
                    11:150,
                    anim,
                    follow_path(
                        simple_bezier() .- (simple_bezier()[1] + Point(0, 3 * frame_start));
                        closed = false,
                    ),
                ),
                SubAction(1:150, color_anim, sethue()),
            ],
        ) for frame_start in 1:7:22
    ]

    javis(
        video,
        [BackgroundAction(1:180, ground), actions...],
        tempdirectory = "images",
        pathname = "",
    )

    @test_reference "refs/followPathBezier1.png" load("images/0000000001.png")
    @test_reference "refs/followPathBezier5.png" load("images/0000000005.png")
    @test_reference "refs/followPathBezier10.png" load("images/0000000010.png")
    @test_reference "refs/followPathBezier30.png" load("images/0000000030.png")
    @test_reference "refs/followPathBezier100.png" load("images/0000000100.png")
    @test_reference "refs/followPathBezier160.png" load("images/0000000160.png")

    for i in 1:180
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "test default kwargs" begin
    video = Video(500, 500)
    pathname = javis(video, [Action(1:10, ground), Action(1:10, morph(astar, acirc))])
    path, ext = splitext(pathname)
    @test ext == ".gif"
    @test isfile(pathname)
    rm(pathname)
end

@testset "test @error .mp3" begin
    video = Video(500, 500)
    @test_logs (:error,) javis(
        video,
        [Action(1:10, ground), Action(1:10, morph(astar, acirc))];
        pathname = "test.mp3",
    )
end


text:
using Nemo, GibbsTypePriors, RCall, JLD
R"library(tidyverse)"

to_plot = Float64.(GibbsTypePriors.Pkn_NGG_robust(1000, 1.2, 0.6; verbose = true))

R"p = tibble(Pkn = $to_plot) %>%
    rowid_to_column(var = 'k') %>%
    ggplot(aes(x = k, y = Pkn)) +
    theme_bw() +
    geom_point() +
    geom_line()"

R"pdf('test/graphical_tests/figures_graphical_tests/dist_Pkn_1000.pdf', height = 5, width = 5)
    plot(p)
    dev.off()"

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors