text:
using Test
using PathWeightSampling
using DiffEqJump
using Catalyst
using StaticArrays

traj = PathWeightSampling.Trajectory([SA[1, 4], SA[2, 5], SA[3, 6]], [1.0, 2.0, 3.0], [1, 1, 0])
traj_mat = PathWeightSampling.Trajectory([1 2 3; 4 5 6], [1.0, 2.0, 3.0], [1, 1, 0])

@test traj == traj_mat
@test length(traj) == 3 == length(traj.t)
@test traj[1] == [1, 4] == traj[:, 1]
@test traj[1, :] == [1, 2, 3]
@test traj[2, :] == [4, 5, 6]
@test traj[1:2] == PathWeightSampling.Trajectory([SA[1, 4], SA[2, 5]], [1.0, 2.0], [1, 1])
@test traj[:, :] == [1 2 3; 4 5 6]

@test collect(tuples(traj)) == [([1, 4], 1.0, 1), ([2, 5], 2.0, 1), ([3, 6], 3.0, 0)]
@test traj(0.0) == [1, 4]
@test traj(1.5) == [2, 5]
@test traj(2.9) == [3, 6]

@test traj([-1.0, 0.0, 0.5, 1.5, 3.5]) == [1 1 1 2 3; 4 4 4 5 6]

sn = @reaction_network begin
    0.005, S --> ∅
    0.25, ∅ --> S
end

u0 = [50.0]
tspan = (0.0, 100.0)
discrete_prob = DiscreteProblem(sn, u0, tspan, [])
jump_prob = JumpProblem(sn, discrete_prob, Direct())

integrator = init(jump_prob, SSAStepper(), tstops=Float64[])
ssa_iter = PathWeightSampling.SSAIter(integrator)

traj_sol = PathWeightSampling.collect_trajectory(ssa_iter)

@test length(traj_sol.i) == length(traj_sol)
@test traj_sol.syms == [Symbol("S(t)")]

sol = integrator.sol
for i in eachindex(sol)
    if i > 1
        @test sol.t[i] == traj_sol.t[i-1]
        @test sol[i] == traj_sol.u[i]
    end
end

for i in 1:length(traj_sol)-1
    if traj_sol.i[i] == 1
        @test (traj_sol.u[i+1] - traj_sol.u[i]) == [-1]
    else
        @test (traj_sol.u[i+1] - traj_sol.u[i]) == [1]
    end
end

# test SSAIterator
iterator = PathWeightSampling.SSAIter(init(jump_prob, SSAStepper()))
collected_values = collect(iterator)
times = getindex.(collected_values, 2)
@test times[begin] > 0.0
@test times[end] == 100.0
@test issorted(times)
@test allunique(times)

traj2 = PathWeightSampling.Trajectory([[1, 4], [2, 5], [3, 6]], [0.5, 1.5, 2.5], [2, 4, 0])

joint = PathWeightSampling.merge_trajectories(traj, traj2)
@test collect(joint) == [
    ([1, 4, 1, 4], 0.5, 2),
    ([1, 4, 2, 5], 1.0, 1),
    ([2, 5, 2, 5], 1.5, 4),
    ([2, 5, 3, 6], 2.0, 1),
    ([3, 6, 3, 6], 2.5, 0),
]

rn = @reaction_network begin
    0.01, S --> X + S
    0.01, X --> ∅
end

network = ModelingToolkit.extend(rn, sn)

u0 = [50.0, 50.0]
tspan = (0.0, 100.0)
discrete_prob = DiscreteProblem(network, u0, tspan, [])
jump_prob = JumpProblem(network, discrete_prob, Direct())

iter = PathWeightSampling.SSAIter(init(jump_prob, SSAStepper(), seed=1))
full = PathWeightSampling.collect_trajectory(iter)
iter = PathWeightSampling.SSAIter(init(jump_prob, SSAStepper(), seed=1))
partial = PathWeightSampling.collect_sub_trajectory(iter, [1])

@test full.syms == [Symbol("S(t)"), Symbol("X(t)")]
@test partial.syms == [Symbol("S(t)")]

sol = iter.integrator.sol
for i in eachindex(partial.t)
    @test partial.t[i] ∈ vcat(sol.t, 100.0)
    @test partial.u[i][1] ∈ sol[[1], :]
    t = partial.t[i]
    @test partial(t) == partial.u[min(i + 1, length(partial))]
    @test partial(t) == full(t)[[1]]
end

# test driven jump problem
import PathWeightSampling
using PathWeightSampling: DrivenJumpProblem

cb_traj = PathWeightSampling.Trajectory([[50], [110], [120], [130]], [30.0, 60.0, 90.0, 100.0])
u0 = [50.0]
tspan = (0.0, 100.0)
discrete_prob = DiscreteProblem(sn, u0, tspan, [])
jump_prob = JumpProblem(sn, discrete_prob, Direct())
djp = DrivenJumpProblem(jump_prob, cb_traj)
integrator = init(djp)
iter = PathWeightSampling.SSAIter(integrator)
traj = iter |> PathWeightSampling.collect_trajectory

for t in [30, 60, 90]
    @test t ∈ traj.t
end
@test traj(30.0)[1] == 110
@test traj(60.0)[1] == 120
@test traj(90.0)[1] == 130


text:
mutable struct Saver
	dirname::AbstractString
    iter::Int64
	Saver(dirname::AbstractString = pwd()) = (mkpath(dirname); new(dirname, 1))
end

macro stash!(group, args...)
    exs = [:($(esc(group))[$(string(arg))] = $(esc(arg))) for arg in args]
    return Expr(:block, exs...)
end

function save_data(saver_o, P, Λ, static_fmin, m_opt, f_min, bs_opt, τs)
	jldopen("$(saver_o.dirname)/dset_$(saver_o.iter).jld2", "w") do file
        @stash!(file, P, Λ, static_fmin, m_opt, f_min, bs_opt, τs)
    end
    saver_o.iter += 1
end

function save_data(saver_o, P, Λ, μ, α, Nₚ, β, Δτ, a, opt, static_f_min, m_opt, f_min, bs_opt, τs)
    jldopen("$(saver_o.dirname)/dset_$(saver_o.iter).jld2", "w") do file
        @stash!(file, P, Λ, μ, α, Nₚ, β, Δτ, a, opt, static_f_min, m_opt, f_min, bs_opt, τs)
    end
    saver_o.iter += 1
end

function process_stacktrace(st)
	str = "\t[1]\t"*string(st[1])
	for i = 2:length(st)
		str = str * "\n\t[$i]\t$(string(st[i]))"
	end
	return str
end

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

const IRange{T} = Union{AbstractVector{T}, T}

function get_optimum(W::Float64, m::Union{Int64, Nothing}, N::Int64, a::Float64, reptyp::AbstractRepulsionType, rdisp::ReducedDispersion, opt::Optim.Options, int_rtol::Real=1e-6, limits::Int64 = 200)
	t′ = -time()
	psamples_raw = get_psamples(rdisp, int_rtol, limits)
	psamples = widen(psamples_raw, W, m)|>separate_psamples
	Nₚ = length(psamples_raw)
	@info "Number of ε⁻ points is: $Nₚ.\n\n"
    params = ParamsB1(N, m, W, a, psamples_raw)
    params0 = ParamsNoB1(params.N, params.m, params.W, params.u)
    ####
    f_nm = precompute_step(zeros(params.N), params0, psamples)
    @info "Free energy of normal metal: $(f_nm).\n\n"
	####
    τs = get_τs(params0)
    bs0 = get_bs(params0)
	d0 = construct_objective(params0, psamples, bs0, reptyp)
	results0 = optimize(d0, bs0, LBFGS(m=120, linesearch = LineSearches.MoreThuente()), opt)
	f_seed = Optim.minimum(results0)
	bs_seed = Optim.minimizer(results0)
    #bs_seed′ = zeros(get_length(params))
    #bs_seed′[1:params.N] .= bs_seed
	t′ += time()
	@info "Free energy of without b1 is $f_seed. Calculation took $t′ s.\n\n"
	####
	#t′ = -time()
	#d = construct_objective(params, psamples, bs_seed′, ω₀)
	#results = optimize(d, bs_seed′, LBFGS(m=120, linesearch = LineSearches.MoreThuente()), opt)
	#f_final = Optim.minimum(results)
	#bs_final = Optim.minimizer(results)
	#t′ += time()
	#@info "Free energy with b1 is $f_final. Calculation took $t′ s.\n\n"
    return f_nm, f_seed, bs_seed, τs
end

function km_walkthrough_repul(W_range::AbstractVector, m_range::AbstractVector, N::Int64, a::Float64, reptyp_range::IRange{T}, rdisp::ReducedDispersion, saver_o::Saver, opt::Optim.Options, int_rtol::Float64 = 1e-7, limits::Int64 = 200, i₀::Int64 = 1) where {T<:AbstractRepulsionType}
	saver_o.iter = i₀
	tups = vec([tup for tup in product(W_range, m_range, reptyp_range)])[i₀:end]
    N_tot = length(tups)
    t0 = time()
    for i in eachindex(tups)
        t1 = time()
        W, m, reptyp = tups[i]
		@info "I am currently dealing with $i-th tuple of (W, m, reptyp), which is (W, m, reptyp) = ($W, $(trm(m)), $reptyp).\n\n"
		f_nm, f_seed, bs_seed, τs = get_optimum(W, m, N, a, reptyp, rdisp, opt, int_rtol, limits)
		jldopen("$(saver_o.dirname)/dset_$(saver_o.iter).jld2", "w") do file
			@stash!(file, W, m, N, a, reptyp, rdisp, opt, int_rtol, limits, f_nm, f_seed, bs_seed, τs)
		end
		saver_o.iter += 1
        t2 = time()
		@info "I am $(t2-t0) s into the computation.\n Finished $i-th run out of $(N_tot) for (W, m, reptyp) = ($W, $(trm(m)), $reptyp). This run took $(t2-t1) s.\n\n\n\n"
    end
end
function km_walkthrough_repul′(W_range::AbstractVector, m_range::AbstractVector, N::Int64, a::Float64, reptyp_range::IRange{T}, rdisp::ReducedDispersion, saver_o::Saver, opt::Optim.Options, int_rtol::Float64 = 1e-7, limits::Int64 = 200, i₀::Int64 = 1) where {T<:AbstractRepulsionType}
	tups = vec([tup for tup in product(W_range, m_range, reptyp_range)])[1:i₀-1]
    N_tot = length(tups)
    t0 = time()
    for i in eachindex(tups)
        t1 = time()
        W, m, reptyp = tups[i]
		@info "I am currently dealing with $i-th tuple of (W, m, reptyp), which is (W, m, reptyp) = ($W, $(trm(m)), $reptyp).\n\n"
		f_nm, f_seed, bs_seed, τs = get_optimum(W, m, N, a, reptyp, rdisp, opt, int_rtol, limits)
		jldopen("$(saver_o.dirname)/dset_$(saver_o.iter).jld2", "w") do file
			@stash!(file, W, m, N, a, reptyp, rdisp, opt, int_rtol, limits, f_nm, f_seed, bs_seed, τs)
		end
		saver_o.iter += 1
        t2 = time()
		@info "I am $(t2-t0) s into the computation.\n Finished $i-th run out of $(N_tot) for (W, m, reptyp) = ($W, $(trm(m)), $reptyp). This run took $(t2-t1) s.\n\n\n\n"
    end
end

trm(m) = m
trm(::Nothing) = "Inf"

export km_walkthrough_repul, km_walkthrough_repul′

extract_number(str::AbstractString) = parse(Int, match(r"\d+", str).match)
function get_files_in_dir(dir::AbstractString = pwd())
	curdir = pwd()
	files = readdir(dir)
	cd(dir)
	files = files |> files -> filter(str -> !isnothing(match(r"jld2", str)), files) |> files -> sort!(files; lt = (str1, str2)->isless(extract_number(str1), extract_number(str2)))
	cd(curdir)
	return map(file -> "$dir/$file", files)
end

ext: jl


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    PositivityPreservingLimiterZhangShu(; threshold, variables)

The fully-discrete positivity-preserving limiter of
- Zhang, Shu (2011)
  Maximum-principle-satisfying and positivity-preserving high-order schemes
  for conservation laws: survey and new developments
  [doi: 10.1098/rspa.2011.0153](https://doi.org/10.1098/rspa.2011.0153)
The limiter is applied to all scalar `variables` in their given order
using the associated `thresholds` to determine the minimal acceptable values.
The order of the `variables` is important and might have a strong influence
on the robustness.
"""
struct PositivityPreservingLimiterZhangShu{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}}
  thresholds::Thresholds
  variables::Variables
end

function PositivityPreservingLimiterZhangShu(; thresholds, variables)
  PositivityPreservingLimiterZhangShu(thresholds, variables)
end

function (limiter!::PositivityPreservingLimiterZhangShu)(
    u_ode, integrator, semi::AbstractSemidiscretization, t)
  u = wrap_array(u_ode, semi)
  @trixi_timeit timer() "positivity-preserving limiter" limiter_zhang_shu!(
    u, limiter!.thresholds, limiter!.variables, mesh_equations_solver_cache(semi)...)
end

# Iterate over tuples in a type-stable way using "lispy tuple programming",
# similar to https://stackoverflow.com/a/55849398:
# Iterating over tuples of different functions isn't type-stable in general
# but accessing the first element of a tuple is type-stable. Hence, it's good
# to process one element at a time and replace iteration by recursion here.
# Note that you shouldn't use this with too many elements per tuple since the
# compile times can increase otherwise - but a handful of elements per tuple
# is definitely fine.
function limiter_zhang_shu!(u, thresholds::NTuple{N,<:Real}, variables::NTuple{N,Any},
                            mesh, equations, solver, cache) where {N}
  threshold = first(thresholds)
  remaining_thresholds = Base.tail(thresholds)
  variable = first(variables)
  remaining_variables = Base.tail(variables)

  limiter_zhang_shu!(u, threshold, variable, mesh, equations, solver, cache)
  limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations, solver, cache)
  return nothing
end

# terminate the type-stable iteration over tuples
function limiter_zhang_shu!(u, thresholds::Tuple{}, variables::Tuple{},
                            mesh, equations, solver, cache)
  nothing
end

include("positivity_zhang_shu_dg1d.jl")
include("positivity_zhang_shu_dg2d.jl")
include("positivity_zhang_shu_dg3d.jl")

end # @muladd


text:
# Julia interface when not defined in source files

eltype{T}(x::DistMultiVec{T}) = T

function size(A::ElementalMatrix, i::Integer)
    if i < 1
        error("dimension out of range")
    elseif i == 1
        return Int(height(A))
    elseif i == 2
        return Int(width(A))
    else
        return 1
    end
end

size(A::ElementalMatrix) = (size(A, 1), size(A, 2))

(*){T<:ElementalMatrix}(A::T, B::T)      = A_mul_B!(one(eltype(A)), A, B, zero(eltype(A)), similar(A, (size(A, 1), size(B, 2))))
(*){T}(A::DistSparseMatrix{T}, B::DistMultiVec{T}) = A_mul_B!(one(T), A, B, zero(T), similar(B, (size(A, 1), size(B, 2))))
Ac_mul_B{T<:ElementalMatrix}(A::T, B::T) = Ac_mul_B!(one(eltype(A)), A, B, zero(eltype(A)), similar(A, (size(A, 2), size(B, 2))))
Ac_mul_B{T}(A::DistSparseMatrix{T}, B::DistMultiVec{T}) = Ac_mul_B!(one(T), A, B, zero(T), similar(B, (size(A, 2), size(B, 2))))

# Spectral
svd(A::ElementalMatrix) = svd!(copy(A))
svd(A::ElementalMatrix, ctrl::SVDCtrl) = svd!(copy(A), ctrl)
svdvals(A::ElementalMatrix, ctrl::SVDCtrl) = svdvals!(copy(A), ctrl)

# conversions to and from julia arrays

# function copy!{T}(dest::Matrix{T}, src::Base.VecOrMat{T})
#     m, n = size(src, 1), size(src, 2)
#     resize!(dest, m, n)
#     Base.unsafe_copy!(pointer(dest), pointer(src), m*n)
#     return dest
# end
# function copy!{T}(dest::Base.Matrix{T}, src::Matrix{T})
#     m, n = size(dest)
#     if m != size(src, 1) || n != size(src, 2)
#         throw(DimensionMisMatch("source and destination must have same shape"))
#     end
#     Base.unsafe_copy!(pointer(dest), pointer(src), m*n)
#     return dest
# end

function copy!{T}(dest::DistMatrix{T}, src::Base.VecOrMat)
    m, n = size(src, 1), size(src, 2)
    zeros!(dest, m, n)
    if MPI.commRank(comm(dest)) == 0
        for j = 1:n
            for i = 1:m
                queueUpdate(dest, i, j, src[i,j])
            end
        end
    end
    processQueues(dest)
    return dest
end

function convert{T}(::Type{Matrix{T}}, A::Base.VecOrMat{T})
    m, n = size(A, 1), size(A, 2)
    B = Matrix(T)
    resize!(B, m, n)
    Base.unsafe_copy!(pointer(B), pointer(A), m*n)
    return B
end
function convert{T}(::Type{Base.Matrix{T}}, A::Matrix{T})
    m, n = size(A)
    B = Base.Matrix{T}(m, n)
    Base.unsafe_copy!(pointer(B), pointer(A), m*n)
    return B
end

function convert{T}(::Type{DistMatrix{T}}, A::Base.VecOrMat{T})
    m, n = size(A, 1), size(A, 2)
    B = DistMatrix(T)
    zeros!(B, m, n)
    if MPI.commRank(comm(B)) == 0
        for j = 1:n
            for i = 1:m
                queueUpdate(B, i, j, A[i,j])
            end
        end
    end
    processQueues(B)
    return B
end

function convert{T}(::Type{DistMultiVec{T}}, A::Base.VecOrMat{T})
    m, n = size(A, 1), size(A, 2)
    B = DistMultiVec(T)
    zeros!(B, m, n)
    if MPI.commRank(comm(B)) == 0
        for j = 1:n
            for i = 1:m
                queueUpdate(B, i, j, A[i,j])
            end
        end
    end
    processQueues(B)
    return B
end

function convert{T}(::Type{DistMatrix{T}}, A::DistMultiVec{T})
    m, n = size(A)
    B = DistMatrix(T)
    zeros!(B, m, n)
    for j = 1:n
        for i = 1:localHeight(A)
            xij = getLocal(A, i, j)
            queueUpdate(B, globalRow(A, i), j, xij)
        end
    end
    processQueues(B)
    return B
end

function norm(x::ElementalMatrix)
    if size(x, 2) == 1
        return nrm2(x)
    else
        return twoNorm(x)
    end
end

# Multiplication
(*){T}(A::DistMatrix{T}, B::Base.VecOrMat{T}) = A*convert(DistMatrix{T}, B)
(*){T}(A::DistMultiVec{T}, B::Base.VecOrMat{T}) = convert(DistMatrix{T}, A)*convert(DistMatrix{T}, B)
(*){T}(A::DistSparseMatrix{T}, B::Base.VecOrMat{T}) = A*convert(DistMultiVec{T}, B)
Ac_mul_B{T}(A::DistMatrix{T}, B::Base.VecOrMat{T}) = Ac_mul_B(A, convert(DistMatrix{T}, B))
Ac_mul_B{T}(A::DistMultiVec{T}, B::Base.VecOrMat{T}) = Ac_mul_B(convert(DistMatrix{T}, A), convert(DistMatrix{T}, B))
Ac_mul_B{T}(A::DistSparseMatrix{T}, B::Base.VecOrMat{T}) = Ac_mul_B(A, convert(DistMultiVec{T}, B))


text:
export command!,
    component!,
    modal!,
    select!

"""
    command!(
        f::Function
        c::Client
        name::AbstractString
        description::AbstractString;
        kwargs...
    )

Adds a handler for INTERACTION CREATE gateway events where the InteractionData's `name` field matches `name`.
Adds this command to `c.commands` or `c.guild_commands` based on the presence of `guild`.
The `f` parameter signature should be:
```
    (ctx::Context, args...) -> Any
```
Where args is a list of all the Command Options

For example a command that takes a user u and a number n as input should have this signature:
```
(ctx::Context, u::User, n::Int) -> Any
```
and the arguments would automatically get converted.

**Note:** The argument names **must** match the Option names. The arguments can be ordered in any way. If no type is specified, no conversion will be performed, so Discord objects will be `Snowflake`s.
"""
function command!(c::Client, f::Function, name::AbstractString, desc::AbstractString, legacy::Bool; auto_ack=true, kwargs...)
    haskey(kwargs, :options) && check_option.(kwargs[:options])
    namecheck(name, r"^[\w-]{1,32}$", 32, "ApplicationCommand Name")
    namecheck(name, 100, "ApplicationCommand Description")
    if !auto_ack push!(c.no_auto_ack, name) end
    legacy ? add_handler!(c, OnInteractionCreate(f; name=name)) : add_handler!(c, OnInteractionCreate(generate_command_func(f); name=name))
end

function command!(f::Function, c::Client, name::AbstractString, description::AbstractString; legacy::Bool=true, auto_ack=true, kwargs...)
    command!(c, f, name, description, legacy; auto_ack=auto_ack, kwargs...)
    add_command!(c; name=name, description=description, kwargs...)
end
function command!(f::Function, c::Client, g::Number, name::AbstractString, description::AbstractString; legacy::Bool=true, auto_ack=true, kwargs...)
    command!(c, f, name, description, legacy; auto_ack=auto_ack, kwargs...)
    add_command!(c, Snowflake(g); name=name, description=description, kwargs...)
end
command!(f::Function, c::Client, g::String, name::AbstractString, description::AbstractString; kwargs...) = command!(f, c, parse(Int, g), name, description; kwargs...)

function modal!(f::Function, c::Client, custom_id::String, int::Interaction; kwargs...)
    add_handler!(c, OnInteractionCreate(generate_command_func(f), custom_id=custom_id))
    push!(c.no_auto_ack, custom_id)
    respond_to_interaction_with_a_modal(c, int.id, int.token; custom_id=custom_id, kwargs...)
end
modal!(f::Function, c::Client, custom_id::String, ctx::Context; kwargs...) = modal!(f, c, custom_id, ctx.interaction; compkwfix(; kwargs...)...)

function check_option(o::ApplicationCommandOption)
    namecheck(o.name, r"^[\w-]{1,32}$", 32, "ApplicationCommandOption Name")
    namecheck(o.description, 100, "ApplicationCommandOption Description")
end

namecheck(val::String, patt::Regex, len::Int, of::String) = (!occursin(patt, val) || length(val) > len) && throw(NamingError(val, of))
namecheck(val::String, patt::Regex, of::String) = namecheck(val, patt, 32, of)
namecheck(val::String, len::Int, of::String) = namecheck(val, r"", len, of)
namecheck(val::String, of::String) = namecheck(val, 32, of)

function generate_command_func(f::Function)
    args = handlerargs(f)
    @debug "Generating typed args" args=args
    g = (ctx::Context) -> begin
        a = makeargs(ctx, args)
        f(ctx, a...)
    end
    g
end

function generate_select_command_func(f::Function)
    g = (ctx::Context) -> begin
        f(ctx, ctx.interaction.data.values)
    end
    g
end

function makeargs(ctx::Context, args)
    v = []
    args = args[2:end]
    if length(args) == 0
        return v
    end
    o = opt(ctx)
    for x = args
        t = last(x)
        arg = get(o, string(first(x)), nothing)
        if arg === nothing
            push!(v, arg)
        else
            push!(v, conv(t, arg, ctx))
        end
    end
    v
end

function conv(::Type{T}, arg, ctx::Context) where T
    if T in [String, Int, Bool, Any]
        return arg
    end
    id = snowflake(arg)
    if T == User
        return ctx.interaction.data.resolved.users[id]
    elseif T == Member
        u = ctx.interaction.data.resolved.users[id]
        m = ctx.interaction.data.resolved.members[id]
        m.user = u
        return m
    elseif T == Role
        return ctx.interaction.data.resolved.roles[id]
    elseif T == DiscordChannel
        return ctx.interaction.data.resolved.channels[id]
    end
    arg
end

"""
    component!(
        f::Function
        c::Client
        custom_id::AbstractString
        kwargs...
    )

Adds a handler for INTERACTION CREATE gateway events where the InteractionData's `custom_id` field matches `custom_id`.
The `f` parameter signature should be:
```
(ctx::Context) -> Any
```
"""
function component!(f::Function, c::Client, custom_id::AbstractString; auto_ack::Bool=true, auto_update_ack::Bool=true, kwargs...)
    add_handler!(c, OnInteractionCreate(f; custom_id=custom_id))
    if !auto_ack push!(c.no_auto_ack, custom_id) end
    if auto_update_ack push!(c.auto_update_ack, custom_id) end
    return Component(custom_id=custom_id; kwargs...)
end

"""
    select!(
        f::Function
        c::Client
        custom_id::AbstractString
        args::Vector{Tuple}
        kwargs...
    )

Adds a handler for INTERACTION CREATE gateway events where the InteractionData's `custom_id` field matches `custom_id`.
The `f` parameter signature should be:
```
(ctx::Context, choices::Vector{String}) -> Any
```
"""
function select!(f::Function, c::Client, custom_id::AbstractString, args...; auto_ack::Bool=true, auto_update_ack::Bool=true, kwargs...)
    add_handler!(c, OnInteractionCreate(generate_select_command_func(f); custom_id=custom_id))
    if !auto_ack push!(c.no_auto_ack, custom_id) end
    if auto_update_ack push!(c.auto_update_ack, custom_id) end
    if length(args) == 0 throw(FieldRequired("options", "SelectOption")) end
    return Component(custom_id=custom_id; options=[SelectOption(a...) for a in args], type=3, kwargs...)
end

"""
    handle(
        c::Client
        handlers,
        data
    )
Determines whether handlers are appropriate to call and calls them if so.
Creates an `AbstractContext` based on the event using the `data` provided and passes it to the handler.
"""
function handle(c::Client, handlers::Vector{Handler}, data::Dict, t::Symbol)
    ctx = context(t, data::Dict)
    hasproperty(ctx, :interaction) && handle_ack(c, ctx)
    isvalidcommand = (h) -> return (!hasproperty(h, :name) || (!ismissing(ctx.interaction.data.name) && h.name == ctx.interaction.data.name))
    isvalidcomponent = (h) -> return (!hasproperty(h, :custom_id) || (!ismissing(ctx.interaction.data.custom_id) && h.custom_id == ctx.interaction.data.custom_id))
    isvalid = (h) -> return isvalidcommand(h) && isvalidcomponent(h)
    @debug "$(repr.(handlers))" currrent = repr(ctx)
    for hh in handlers
        isvalid(hh) && (runhandler(c, hh, ctx, t))
    end
end

function handle_ack(c::Client, ctx::Context)
    iscomp = !ismissing(ctx.interaction.data.custom_id)
    !ismissing(ctx.interaction.data.type) && ctx.interaction.data.type == 5 && return
    if !iscomp || !(ctx.interaction.data.custom_id in c.no_auto_ack)
        @debug "Acking interaction" name=ctx.interaction.data.name custom_id=ctx.interaction.data.custom_id
        if iscomp && (ctx.interaction.data.custom_id in c.auto_update_ack)
            update_ack_interaction(c, ctx.interaction.id, ctx.interaction.token)
        elseif (!ismissing(ctx.interaction.data.name) && !(ctx.interaction.data.name in c.no_auto_ack)) || iscomp
            ack_interaction(c, ctx.interaction.id, ctx.interaction.token)
        end
    end
end

"""
Runs a handler with given context
"""
function runhandler(c::Client, h::Handler, ctx::Context, t::Symbol)
    @debug "Running handler" h=Handler type=t
    @spawn begin try h.f(ctx) catch e
            kws = logkws(c; handler=t, exception=(e, catch_backtrace()))
            @warn "Handler errored" kws...
    end end
end

"""
Calls handle with a list of all found handlers.
"""
function handle(c::Client, t::Symbol, data::Dict)
    haskey(c.handlers, t) ? handle(c, c.handlers[t], data, t) : return nothing
end


text:
"""
gpreg_all_threads(t, beta, st, gpreg=gpmodels ; nt = 16, bt = 1, n = size(beta, 2))

Wrapper to run a Gaussian Process regression function `gp_reg ∈ {gp_reg_matern52, gp_reg_const, gp_reg_linear}` with multiple threads.

Allows balance of threads to run regression on independent probes in parallel vs threads used by BLAS library for parallelising linear algebra used in regression of a single probe.

Currently, recommend setting BLAS threads `bt = 1` and parallelising over multiple probes.

"""
function gpreg_all_threads(t, beta, st, gpreg=gpmodels ; nt = 16, bt = 1, n = size(beta, 2))
        btc = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
        BLAS.set_num_threads(bt)

    dft = DataFrame(const_mll    = Vector{Float64}(undef, n),
                    const_σn     = Vector{Float64}(undef, n),
                    const_σf     = Vector{Float64}(undef, n),
                    linear_mll   = Vector{Float64}(undef, n),
                    linear_σn    = Vector{Float64}(undef, n),
                    linear_σf    = Vector{Float64}(undef, n),
                    linear_ℓ     = Vector{Float64}(undef, n),
                    mat52_mll    = Vector{Float64}(undef, n),
                    mat52_σn     = Vector{Float64}(undef, n),
                    mat52_σf     = Vector{Float64}(undef, n),
                    mat52_ℓ      = Vector{Float64}(undef, n))
    μ = Matrix{Float64}(undef, length(st), n)
    v = Matrix{Float64}(undef, length(st), n)

    p = Progress(n, barlen=75)
    Threads.@threads for i = 1:n
        params, gpu, gpv = gpreg(t, beta[:, i], st)
        dft[i, :] = params
        μ[:, i] .= gpu
        v[:, i] .= gpv
        next!(p)
    end

    BLAS.set_num_threads(btc)
    dft, μ, v
end

filtmissing(t, y) = t, y
function filtmissing(t, y::Vector{Union{T, Missing}}) where {T}
    ind = .!ismissing.(y)
    t[ind], Vector{T}(y[ind])
end

"""
    gpmodels(t, y, st)

Runs GP regression for constant, linear and Matern52 kernels.
"""
function gpmodels(td, yd, st)

    ### for now hardcode const, linear and mat52 regression
    ### improvement would be to specify vector of kernels and loop over

    ### remove missings if necessary
    t, y = filtmissing(td, yd)

    gpc = gpreg_const(t, y, st)
    gpl = gpreg_linear(t, y, st)
    gpm = gpreg_matern52(t, y, st)

    params = (const_mll    = gpc.gp.mll,
     const_σn    = exp(gpc.gp.logNoise.value),
     const_σf    = sqrt(gpc.gp.kernel.σ2),
     linear_mll  = gpl.gp.mll,
     linear_σn   = exp(gpl.gp.logNoise.value),
     linear_σf   = sqrt(gpl.gp.kernel.kleft.σ2),
     linear_ℓ    = sqrt(gpl.gp.kernel.kright.ℓ2),
     mat52_mll   = gpm.gp.mll,
     mat52_σn    = exp(gpm.gp.logNoise.value),
     mat52_σf    = sqrt(gpm.gp.kernel.σ2),
     mat52_ℓ     = gpm.gp.kernel.ℓ)

     params, gpm.μ,  gpm.v
end

"""
    gpmodels_prior(t, y, st)

Runs GP regression for constant, linear and Matern52 kernels with priors on hyperparameters for constant and linear kernels
"""
function gpmodels_prior(td, yd, st)

    ### for now hardcode const, linear and mat52 regression
    ### improvement would be to specify vector of kernels and loop over

    ### remove missings if necessary
    t, y = filtmissing(td, yd)

    gpc = gpreg_const(t, y, st)
    gpl = gpreg_linear(t, y, st)
    gpm = gpreg_matern52_prior(t, y, st)

    params = (const_mll    = gpc.gp.mll,
     const_σn    = exp(gpc.gp.logNoise.value),
     const_σf    = sqrt(gpc.gp.kernel.σ2),
     linear_mll  = gpl.gp.mll,
     linear_σn   = exp(gpl.gp.logNoise.value),
     linear_σf   = sqrt(gpl.gp.kernel.kleft.σ2),
     linear_ℓ    = sqrt(gpl.gp.kernel.kright.ℓ2),
     mat52_mll   = gpm.gp.mll,
     mat52_σn    = exp(gpm.gp.logNoise.value),
     mat52_σf    = sqrt(gpm.gp.kernel.σ2),
     mat52_ℓ     = gpm.gp.kernel.ℓ)

     params, gpm.μ,  gpm.v
end

"""
    gpreg(t, y, st, kernel, kernparams=(1.0, 1.0), logNoise=log(std(y)))

Generic GP stationary zero mean regression with Gaussian likelihood. Params:
  - `t` vector of time points (in this case ages)
  - `y` vector of signal values (in this case prob methylation)
  - `st` vector of time points to evaluate GP `mean` and `var`
  - `kernel` GP `kernel` with inital params e.g. `Mat52Iso`
  - `logσn` log(σn)

Returns:
Named tuple (gp=gp, μ=μ, v=v)
"""
function gpreg(t, y, st, kernel, logσn=log(std(y)))
    gp = GP(t, y, MeanZero(), kernel, logσn)
    optimize!(gp)
    μ, v = predict_y(gp, st)
    (gp=gp, μ=μ, v=v)
end

"""
    gpreg_matern52_prior(t, y, st, kernparams=(2.0, log(mean(y))), logσn=log(std(y)), kernpriors=[Normal(4.7, 1.2), Normal(-1, 6)])

GP regression with `Mat52Iso(ℓ, σf)` kernel.

`kernparms = (logℓ, log(σf))`

`ℓ`  - initial lengthscale

`σf` - initial function variance
"""
function gpreg_matern52_prior(t, y, st, kernparams=(2.0, log(mean(y))), logσn=log(std(y)), kernpriors=[Normal(4.7, 1.2), Normal(-1, 6)])
    logℓ, logσf = kernparams
    kern = Mat52Iso(logℓ, logσf)
    set_priors!(kern, kernpriors)
    gpreg(t, y, st, kern, logσn)
end

"""
    gpreg_matern52(t, y, st, kernparams=(1.0, 1.0) , log(std(y))))

GP regression with `Mat52Iso(ℓ, σf)` kernel.

`kernparms = (logℓ, log(σf))`

`ℓ`  - initial lengthscale

`σf` - initial function variance
"""
function gpreg_matern52(t, y, st, kernparams=(1.0, log(mean(y))), logσn=log(std(y)))
    logℓ, logσf = kernparams
    gpreg(t, y, st, Mat52Iso(logℓ, logσf), logσn)
end

"""
    gpreg_const(t, y, st, kernparams=(1.0) , log(std(y))))

GP regression with `Const(σf)`.

`kernparms = (log(σf))`

`σf` - initial function variance
"""
function gpreg_const(t, y, st, kernparams=(log(mean(y))), logσn=log(std(y)))
    logσf = kernparams
    gpreg(t, y, st, Const(logσf), logσn)
end

"""
    gpreg_linear(t, y, st, kernparams=(1.0) , log(std(y))))

GP regression with `Const(σf) + LinIso(ℓ)` kernel.

`kernparms = (logℓ, log(σf))`

`ℓ`  - initial lengthscale

`σf` - initial function variance
"""
function gpreg_linear(t, y, st, kernparams=(1.0, log(mean(y))), logσn=log(std(y)))
    logℓ, logσf = kernparams
    gpreg(t, y, st, Const(logσf) + LinIso(logℓ), logσn)
end

"""
    gpreg_linear_prior(t, y, st, kernparams=(1.0, log(mean(y))), logσn=log(std(y)), kernpriors=[Normal(5.7, 2.2)])

GP regression with `Const(σf) + LinIso(ℓ)` kernel.

`kernparms = (logℓ, log(σf))`

`ℓ`  - initial lengthscale

`σf` - initial function variance
"""
function gpreg_linear_prior(t, y, st, kernparams=(1.0, log(mean(y))), logσn=log(std(y)), kernpriors=[Normal(5.7, 2.2)])
    logℓ, logσf = kernparams
    kc = Const(logσf)
    kl = LinIso(logℓ)
    set_priors!(kl, kernpriors)
    k = kc + kl
    gpreg(t, y, st, k, logσn)
end


text:
module SDEModels
using Parameters
using StaticArrays
using LinearAlgebra

abstract type AbstractSDE{D,M} end
abstract type JumpSDE{D,M} <: AbstractSDE{D,M} end
abstract type StateIndependentDiffusion{D,M} <: AbstractSDE{D,M} end

abstract type AbstractScheme end
abstract type ExplicitScheme <: AbstractScheme end
abstract type ImplicitScheme <: AbstractScheme end
abstract type ConditionalScheme{T} <: ExplicitScheme end
abstract type UnconditionalScheme <: ExplicitScheme end

include("schemes/schemes.jl")
include("models/models.jl")
include("simulation.jl")
include("multilevel.jl")

export @sde_model, dim, drift, diffusion, jump, corrected_drift, variables, subdivide
export sample, sample!, simulate, simulate!
export MultilevelScheme, npaths
export transition, logtransition
export StaticArrays

end # module


text:
# tests default choice model

@generator SCMGen begin # prefix SCM (for Sampler Choice Model) to avoid type name clashes
    start() = reps(a,1+1,2+2) # non-literal arguments to allow us to pass range to godelnumber during testing
    a() = choose(Int,-2-2,1+1)
    a() = choose(Float64,-2.0-2.0,1.0+1.0)
    a() = choose(Bool)
end

@generator SCMRepsGen begin
    start() = reps(a,1+1,2+2) # non-literal arguments to allow us to pass range to godelnumber during testing
    a() = 'a'
end

@generator SCMChooseBoolGen begin
    start() = choose(Bool) #
end

@generator SCMChooseIntGen begin
    start() = choose(Int,-2-2,1+1) # non-literal arguments to allow us to pass range to godelnumber during testing
end

@generator SCMChooseFloat64Gen begin
    start() = choose(Float64,-2.0-2.0,1.0+1.0) # non-literal arguments to allow us to pass range to godelnumber during testing
end

@generator SCMRuleGen begin
    start() = x()
    x() = 'a'
    x() = 'b'
    x() = 'c'
    x() = 'd'
end

# internally will use range of choice points, so convenient for testing choice model when more than one choice point
@generator SCMChooseStringGen begin
    start() = choose(String,"a(b|c)d+ef?")
end

@generator SCMEstimateParamsGen begin
	start() = a()
	a() = choose(Bool)
	a() = mult("x")
end

@testset "sampler choice model" begin

	@testset "constructor" begin

	    gn = SCMGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)
	    @test typeof(cm) == DataGenerators.SamplerChoiceModel

	end

	@testset "set/get parameters and ranges" begin

		gn = SCMGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)

		@testset "paramranges" begin
		    ranges = paramranges(cm)
		    @test typeof(ranges) <: Vector{Tuple{Float64,Float64}}
		    @test length(ranges) == 5
		end

		@testset "getparams" begin
		    params = getparams(cm)
		    @test typeof(params) <: Vector{Float64}
		    @test length(params) == 5
		end

		@testset "setparams" begin
		    newparams = [(paramrange[1]+paramrange[2])/2 for paramrange in paramranges(cm)]
		    setparams!(cm, newparams)
		    @test length(getparams(cm)) == 5 # can't check equality of params since some adjustment can be made by the cm
		end

	end

	@testset "rule choice point" begin

		gn = SCMRuleGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		# type ChoiceContext
		# 	derivationstate::DerivationState
		# 	cptype::Symbol
		# 	cpid::Uint64
		# 	datatype::DataType
		# 	lowerbound::Real
		# 	upperbound::Real
		# end
		cc = DataGenerators.ChoiceContext(DummyDerivationState(), :rule, cpids[1], Int, 1, 4)

		@mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		    gnum, trace = DataGenerators.godelnumber(cm, cc)
		    @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		    @mtest_values_are [1,2,3,4] gnum
		end

		@testset "sampler" begin

			@testset "uses a Categorical distribution" begin
				@test length(cm.samplers) == 1
				sampler = first(values(cm.samplers))
				@test typeof(sampler) <: DataGenerators.CategoricalSampler
			end

			@testset "get and set parameters" begin
				sampler = first(values(cm.samplers))
				@test getparams(cm) == [0.25,0.25,0.25,0.25,]
				@test paramranges(cm) == [(0.0,1.0),(0.0,1.0),(0.0,1.0),(0.0,1.0),]
				setparams!(cm,[0.03,0.03,0.01,0.01])
				@test round(getparams(cm),4) == [0.375,0.375,0.125,0.125,] # round to avoid precision errors
			end

		end

	end

	@testset "sequence choice point" begin

		gn = SCMRepsGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		@testset "small finite range" begin

			cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int, 0, 2)

			@mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
			    gnum, trace = DataGenerators.godelnumber(cm, cc)
			    @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
			    @test 0 <= gnum <= 2 # default choice model restricts sequence lengths to a maximum of 3 more than minimum
			    @mtest_values_are [0,1,2] gnum
			end

		end

		@testset "large finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int, 11, 16)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 11 <= gnum <= 16 # default choice model restricts sequence lengths to a maximum of 3 more than minimum
		        @mtest_values_include [11,13] gnum # 16 isn't very likely to occur
		    end

		end

		@testset "semi-finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int, 1, typemax(Int))

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 1 <= gnum <= typemax(Int64)
		        @mtest_values_vary gnum
		    end

		end

		@testset "sampler" begin

			cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int64, 7, typemax(Int64))

			@mtestset "sampler consistent with a offset geometric distribution" reps=Main.REPS alpha=Main.ALPHA begin
				gnum, trace = DataGenerators.godelnumber(cm, cc)
				@mtest_distributed_as Geometric(0.5) (gnum-7)
			end

			@testset "get and set parameters" begin
				sampler = first(values(cm.samplers))
				@test getparams(cm) == [0.5,]
				@test paramranges(cm) == [(0.0,1.0),]
				setparams!(cm,[0.6])
				@test getparams(cm) == [0.6,]
			end

		end

	end

	@testset "Bool value choice point" begin

		gn = SCMChooseBoolGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Bool, false, true)

		@mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		    gnum, trace = DataGenerators.godelnumber(cm, cc)
		    @test convert(Bool,gnum) != nothing  # raises exception if value can't be converted
		    @mtest_values_are [false,true] gnum
		end

		@testset "sampler" begin

			@mtestset "sampler consistent with a Bernoulli distribution" reps=Main.REPS alpha=Main.ALPHA begin
				gnum, trace = DataGenerators.godelnumber(cm, cc)
				@mtest_distributed_as Bernoulli(0.5) gnum
			end

			@testset "get and set parameters" begin
				sampler = first(values(cm.samplers))
				@test getparams(cm) == [0.5,]
				@test paramranges(cm) == [(0.0,1.0),]
				setparams!(cm,[0.6])
				@test getparams(cm) == [0.6,]
			end

		end

	end

	@testset "Int value choice point" begin

		gn = SCMChooseIntGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		@testset "small finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, -1, 2)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test -1 <= gnum <= 2
		        @mtest_values_are [-1,0,1,2] gnum
		    end

		end

		@testset "large finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, 11, 16)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 11 <= gnum <= 16
		        @mtest_values_include [11,13,16] gnum # just a selection of possible values including end points
		    end

		end

		@testset "semi-finite range (upper)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, 128, typemax(Int))

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 128 <= gnum <= typemax(Int)
		        @mtest_values_vary gnum
		    end

		end

		@testset "semi-finite range (lower)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, typemin(Int), 128)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test typemin(Int) <= gnum <= 128
		        @mtest_values_vary gnum
		    end

		end

		@testset "infinite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, typemin(Int), typemax(Int))

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test typemin(Int) <= gnum <= typemax(Int)
		        @mtest_values_vary gnum
		    end

		end

		@testset "sampler" begin

			cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, 29, 301)

			@mtestset "sampler consistent with a discrete uniform distribution" reps=Main.REPS alpha=Main.ALPHA begin
				gnum, trace = DataGenerators.godelnumber(cm, cc)
				@mtest_distributed_as DiscreteUniform(29,301) gnum
			end

			@testset "get and set parameters" begin
				sampler = first(values(cm.samplers))
				@test getparams(cm) == []
				@test paramranges(cm) == []
			end

		end

	end

	@testset "Float64 value choice point" begin

		gn = SCMChooseFloat64Gen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		@testset "finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -42.2, -8.7)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @test -42.2 <= gnum <= -8.7
		        @mtest_that_sometimes round(gnum) != gnum
		        @mtest_values_vary gnum
		    end

		end

		@testset "semi-finite range (upper)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, 450001.6, Inf)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @test 450001.6 <= gnum
		        @mtest_values_vary gnum
		    end

		end

		@testset "semi-finite range (lower)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -Inf, 450001.6)

		    @mtestset "valid Godel numbers returned" begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @test gnum <= 450001.6
		        @mtest_values_vary gnum
		    end

		end

		@testset "infinite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -Inf, Inf)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @mtest_values_vary gnum
		    end

		end

		@testset "sampler" begin

			cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -180.7, 123.728)

			@mtestset "sampler consistent with a uniform distribution" reps=Main.REPS alpha=Main.ALPHA begin
				gnum, trace = DataGenerators.godelnumber(cm, cc)
				@mtest_distributed_as Uniform(-180.7,123.728) gnum
			end

			@testset "get and set parameters" begin
				sampler = first(values(cm.samplers))
				@test getparams(cm) == []
				@test paramranges(cm) == []
			end

		end

	end

	@testset "generate for model with multiple choice points" begin

		gn = SCMChooseStringGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    cm = choicemodel(gn)

		@mtestset "full range of values generated using choice model" reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn, choicemodel=cm)
		    @test ismatch(r"^a(b|c)d+ef?$", td)
		    @mtest_values_include [1,2,3] count(x->x=='d', td)
		    @mtest_values_are ['b','c'] td[2]
		    @mtest_values_are ['e','f'] td[end]
		end

	end

	@testset "non-default mapping" begin

		gn = SCMRuleGen()

		function nondefaultmapping(info::Dict)
			cptype = info[:type]
			if cptype == :rule
				sampler = DataGenerators.DiscreteUniformSampler()
			elseif cptype == :sequence
				sampler = DataGenerators.GeometricSampler()
			elseif cptype == :value
				datatype = info[:datatype]
				if datatype <: Bool
					sampler = DataGenerators.BernoulliSampler()
				elseif datatype <: Integer # order of if clauses matters here since Bool <: Integer
					sampler = DataGenerators.DiscreteUniformSampler()
				else # floating point, but may also be a rational type
					sampler = DataGenerators.UniformSampler()
				end
			else
				error("unrecognised choice point type when creating non-default choice model")
			end
		end

		setchoicemodel!(gn, SamplerChoiceModel(gn, choicepointmapping=nondefaultmapping))
	    cm = choicemodel(gn)

		@test length(cm.samplers) == 1
		sampler = first(values(cm.samplers))
		@test typeof(sampler) <: DataGenerators.DiscreteUniformSampler

	end

	@testset "estimate parameters" begin

		gn = SCMEstimateParamsGen()
		setchoicemodel!(gn, SamplerChoiceModel(gn))
	    scm1 = choicemodel(gn)
		scm2 = deepcopy(scm1)

		# parameters should be: 1:2 categorical for rule choice; 3: bernoulli for choose(Bool); 4: geometric for mult()
		params = [0.3, 0.7, 0.33, 0.45]
		setparams!(scm1, params)

		otherparams = [0.6, 0.4, 0.58, 0.77]
		setparams!(scm2, otherparams)

		cmtraces = map(1:100) do i
			result, state = generate(gn)
			state.cmtrace
		end

		estimateparams!(scm2, cmtraces)

		# to check, we create samplers from the groups of parameters
		# (could also simply 'look inside' the choice model, but this would be less robust to code changes)
		scm2params = getparams(scm2)
		cc = dummyChoiceContext()

		cat2 = DataGenerators.CategoricalSampler(2, scm2params[1:2])
	    @mtestset "consistent with categorical" reps=Main.REPS alpha=Main.ALPHA begin
        	x, trace = DataGenerators.sample(cat2, (0,1), cc)
	        @mtest_distributed_as Categorical(params[1:2]) x
	    end

		bern2 = DataGenerators.BernoulliSampler(scm2params[3:3])
	    @mtestset "consistent with Bernoulli" reps=Main.REPS alpha=Main.ALPHA begin
        	x, trace = DataGenerators.sample(bern2, (0,1), cc)
	        @mtest_distributed_as Bernoulli(params[3]) x
	    end

		geom2 = DataGenerators.GeometricSampler(scm2params[4:4])
	    @mtestset "consistent with geometric" reps=Main.REPS alpha=Main.ALPHA begin
        	x, trace = DataGenerators.sample(geom2, (0,1), cc)
	        @mtest_distributed_as Geometric(params[4]) x
	    end

	end

end


text:
# create a timer object, used throughout this package, see TimerOutputs.jl
global const to = TimerOutput();
global const npml = 20

# struct to store FD params that are used in @initialize
mutable struct FD_params
    ndims::Int64 # 2D or 3D
    order::Int64 # stencil order
    use_gpu::Bool #
    datatype::DataType
end
# default finite-difference stencil order
const _fd = FD_params(2, 2, false, Float32)

"""
Return axis names of 1D, 2D or 3D fields
"""
function dim_names(N, prefix = "", suffix = ""; order = 1)
    if (N == 3)
        names = (order == 2) ? [:zz, :yy, :xx, :xz, :xy, :yz] : [:z, :y, :x]
    elseif (N == 2)
        names = (order == 2) ? [:zz, :xx, :xz] : [:z, :x]
    elseif (N == 1)
        names = [:z]
    else
        error("invalid dim num")
    end
    return [Symbol(prefix, string(m), suffix) for m in names]
end


text:
module ModuleB

using LinearAlgebra

export myPower, myDeflation, myPowerMethod, myTridiag, myTridiagX, myTridiagG, myTridEigQR,
mySymEigQR, stedc!, myJacobi, myLanczos, myBidiagX, myBidiagY, myBidiagG, myJacobiR

function myPower(A::Matrix,x::Vector,tol::Real)
    y=A*x
    ν=x⋅y
    steps=1
    while norm(y-ν*x)>tol
        x=y/norm(y)
        y=A*x
        ν=x⋅y
        steps+=1
    end
    ν, y/norm(y), steps
end

function myDeflation(A::Matrix,x::Vector)
    X,R=qr(x[:,:])
    # To make sure the returned matrix symmetric use
    # full(Symmetric(X[:,2:end]'*A*X[:,2:end]))
    X[:,2:end]'*A*X[:,2:end]
end

function myPowerMethod(A::Matrix{T}, tol::Real) where T
    n=size(A,1)
    S =  T==Int ? Float64 : T
    λ=Vector{S}(undef,n)
    for i=1:n
        λ[i],x,steps=myPower(A,ones(n-i+1),tol)
        A=myDeflation(A,x)
    end
    λ
end

function myTridiag(A::Matrix)
    # Normalized Householder vectors are stored in the lower
    # triangular part of A below the first subdiagonal
    n=size(A,1)
    T=Float64
    A=map(T,A)
    v=Vector{T}(undef,n)
    Trid=SymTridiagonal(zeros(n),zeros(n-1))
    for j = 1 : n-2
        μ = sign(A[j+1,j])*norm(A[j+1:n, j])
        if μ != zero(T)
            β =A[j+1,j]+μ
            v[j+2:n] = A[j+2:n,j] / β
        end
        A[j+1,j]=-μ
        A[j,j+1]=-μ
        v[j+1] = one(T)
        γ = -2 / (v[j+1:n]⋅v[j+1:n])
        w = γ* A[j+1:n, j+1:n]*v[j+1:n]
        q = w + γ * v[j+1:n]*(v[j+1:n]⋅w) / 2
        A[j+1:n, j+1:n] = A[j+1:n,j+1:n] + v[j+1:n]*q' + q*v[j+1:n]'
        A[j+2:n, j] = v[j+2:n]
    end
    SymTridiagonal(diag(A),diag(A,1)), tril(A,-2)
end

# Extract X
function myTridiagX(H::Matrix)
    n=size(H,1)
    T=Float64
    X = Matrix{T}(I,n,n)
    v=Vector{T}(undef,n)
    for j = n-2 : -1 : 1
        v[j+1] = one(T)
        v[j+2:n] = H[j+2:n, j]
        γ = -2 / (v[j+1:n]⋅v[j+1:n])
        w = γ * X[j+1:n, j+1:n]'*v[j+1:n]
        X[j+1:n, j+1:n] = X[j+1:n, j+1:n] + v[j+1:n]*w'
    end
    X
end

# Tridiagonalization using Givens rotations
function myTridiagG(A::Matrix)
    n=size(A,1)
    T=Float64
    X=Matrix{T}(I,n,n)
    for j = 1 : n-2
        for i = j+2 : n
            G,r=givens(A,j+1,i,j)
            A=(G*A)*G'
            X*=G'
        end
    end
    SymTridiagonal(diag(A),diag(A,1)), X
end

function myTridEigQR(A1::SymTridiagonal)
    A=deepcopy(A1)
    n=length(A.dv)
    T=Float64
    λ=Vector{T}(undef,n)
    B=Matrix{T}
    if n==1
        return map(T,A.dv)
    end
    if n==2
        τ=(A.dv[end-1]-A.dv[end])/2
        μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
        # Only rotation
        B=A[1:2,1:2]
        G,r=givens(B-μ*I,1,2,1)
        B=(G*B)*G'
        return diag(B)[1:2]
    end
    steps=1
    k=0
    while k==0 && steps<=10
        # Shift
        τ=(A.dv[end-1]-A.dv[end])/2
        μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
        # First rotation
        B=A[1:3,1:3]
        G,r=givens(B-μ*I,1,2,1)
        B=(G*B)*G'
        A.dv[1:2]=diag(B)[1:2]
        A.ev[1:2]=diag(B,-1)
        bulge=B[3,1]
        # Bulge chasing
        for i = 2 : n-2
            B=A[i-1:i+2,i-1:i+2]
            B[3,1]=bulge
            B[1,3]=bulge
            G,r=givens(B,2,3,1)
            B=(G*B)*G'
            A.dv[i:i+1]=diag(B)[2:3]
            A.ev[i-1:i+1]=diag(B,-1)
            bulge=B[4,2]
        end
        # Last rotation
        B=A[n-2:n,n-2:n]
        B[3,1]=bulge
        B[1,3]=bulge
        G,r=givens(B,2,3,1)
        B=(G*B)*G'
        A.dv[n-1:n]=diag(B)[2:3]
        A.ev[n-2:n-1]=diag(B,-1)
        steps+=1
        # Deflation criterion
        k=findfirst(abs.(A.ev) .< sqrt.(abs.(A.dv[1:n-1].*A.dv[2:n]))*eps(T))
        k=k==nothing ? 0 : k
    end
    λ[1:k]=myTridEigQR(SymTridiagonal(A.dv[1:k],A.ev[1:k-1]))
    λ[k+1:n]=myTridEigQR(SymTridiagonal(A.dv[k+1:n],A.ev[k+1:n-1]))
    λ
end

function myTridEigQR(A1::SymTridiagonal,U::Matrix)
    # U is either the identity matrix or the output from myTridiagX()
    A=deepcopy(A1)
    n=length(A.dv)
    T=Float64
    λ=Vector{T}(undef,n)
    B=Matrix{T}
    if n==1
        return map(T,A.dv), U
    end
    if n==2
        τ=(A.dv[end-1]-A.dv[end])/2
        μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
        # Only rotation
        B=A[1:2,1:2]
        G,r=givens(B-μ*I,1,2,1)
        B=(G*B)*G'
        U*=G'
        return diag(B)[1:2], U
    end
    steps=1
    k=0
    while k==0 && steps<=10
        # Shift
        τ=(A.dv[end-1]-A.dv[end])/2
        μ=A.dv[end]-A.ev[end]^2/(τ+sign(τ)*sqrt(τ^2+A.ev[end]^2))
        # First rotation
        B=A[1:3,1:3]
        G,r=givens(B-μ*I,1,2,1)
        B=(G*B)*G'
        U[:,1:3]*=G'
        A.dv[1:2]=diag(B)[1:2]
        A.ev[1:2]=diag(B,-1)
        bulge=B[3,1]
        # Bulge chasing
        for i = 2 : n-2
            B=A[i-1:i+2,i-1:i+2]
            B[3,1]=bulge
            B[1,3]=bulge
            G,r=givens(B,2,3,1)
            B=(G*B)*G'
            U[:,i-1:i+2]=U[:,i-1:i+2]*G'
            A.dv[i:i+1]=diag(B)[2:3]
            A.ev[i-1:i+1]=diag(B,-1)
            bulge=B[4,2]
        end
        # Last rotation
        B=A[n-2:n,n-2:n]
        B[3,1]=bulge
        B[1,3]=bulge
        G,r=givens(B,2,3,1)
        B=(G*B)*G'
        U[:,n-2:n]*=G'
        A.dv[n-1:n]=diag(B)[2:3]
        A.ev[n-2:n-1]=diag(B,-1)
        steps+=1
        # Deflation criterion
        k=findfirst(abs.(A.ev) .< sqrt.(abs.(A.dv[1:n-1].*A.dv[2:n]))*eps())
        k=k==nothing ? 0 : k
    end
    λ[1:k], U[:,1:k]=myTridEigQR(SymTridiagonal(A.dv[1:k],A.ev[1:k-1]),U[:,1:k])
    λ[k+1:n], U[:,k+1:n]=myTridEigQR(SymTridiagonal(A.dv[k+1:n],A.ev[k+1:n-1]),U[:,k+1:n])
    λ, U
end

function mySymEigQR(A::Matrix)
    T,H=myTridiag(A)
    X=myTridiagX(H)
    myTridEigQR(T,X)
end

### DSTEDC
# Part of the preamble of lapack.jl
const liblapack = Base.liblapack_name
import LinearAlgebra.BLAS.@blasfunc
# import ..LinAlg: BlasFloat, Char, BlasInt, LAPACKException,
    # DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare
import LinearAlgebra.BlasInt
function chklapackerror(ret::BlasInt)
    if ret == 0
        return
    elseif ret < 0
        throw(ArgumentError("invalid argument #$(-ret) to LAPACK call"))
    else # ret > 0
        throw(LAPACKException(ret))
    end
end

for (stedc, elty) in
    ((:dstedc_,:Float64),
    (:sstedc_,:Float32))
    @eval begin
        """
        COMPZ is CHARACTER*1
          = 'N':  Compute eigenvalues only.
          = 'I':  Compute eigenvectors of tridiagonal matrix also.
          = 'V':  Compute eigenvectors of original dense symmetric
                  matrix also.  On entry, Z contains the orthogonal
                  matrix used to reduce the original matrix to
                  tridiagonal form.
        """
        function stedc!(compz::Char, dv::Vector{$elty}, ev::Vector{$elty},
                Z::Array{$elty})
            n = length(dv)
            ldz=n
            if length(ev) != n - 1
                throw(DimensionMismatch("ev has length $(length(ev))
                        but needs one less than dv's length, $n)"))
            end
            w = deepcopy(dv)
            u = deepcopy(ev)
            lwork=5*n^2
            work = Array{$elty}(undef,lwork)
            liwork=6+6*n+5*n*round(Int,ceil(log(n)/log(2)))
            iwork = Array{BlasInt}(undef,liwork)
            # info = Array{BlasInt}(undef,5)
            info = Ref{BlasInt}()
            ccall((@blasfunc($stedc), liblapack), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                    Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                    Ptr{BlasInt}), compz, n, w, u, Z, ldz, work,
                    lwork, iwork, liwork, info)
                chklapackerror(info[])
            return w,Z
        end
    end
end

function myJacobi(A1::Array{T}) where T<:Real
    A=deepcopy(A1)
    n=size(A,1)
    U=Matrix{T}(I,n,n)
    # Tolerance for rotation
    tol=sqrt(map(T,n))*eps(T)
    # Counters
    p=n*(n-1)/2
    sweep=0
    pcurrent=0
    # First criterion is for standard accuracy, second one is for relative accuracy
    # while sweep<30 && vecnorm(A-diagm(diag(A)))>tol
    while sweep<30 && pcurrent<p
        sweep+=1
        # Row-cyclic strategy
        for i = 1 : n-1
            for j = i+1 : n
                # Check for the tolerance - the first criterion is standard,
                # the second one is for relative accuracy for PD matrices
                # if A[i,j]!=zero(T)
                if abs(A[i,j])>tol*sqrt(abs(A[i,i]*A[j,j]))
                    # Compute c and s
                    τ=(A[i,i]-A[j,j])/(2*A[i,j])
                    t=sign(τ)/(abs(τ)+sqrt(1+τ^2))
                    c=1/sqrt(1+t^2)
                    s=c*t
                    G=LinearAlgebra.Givens(i,j,c,s)
                    # A=G*A
                    lmul!(G,A)
                    # A*=G'
                    rmul!(A,adjoint(G))
                    A[i,j]=zero(T)
                    A[j,i]=zero(T)
                    # U*=G'
                    rmul!(U,adjoint(G))
                    pcurrent=0
                else
                    pcurrent+=1
                end
            end
        end
    end
    diag(A), U
end

function myLanczos(A::Array{T}, x::Vector{T}, k::Int) where T
    n=size(A,1)
    X=Array{T}(undef,n,k)
    dv=Array{T}(undef,k)
    ev=Array{T}(undef,k-1)
    X[:,1]=x/norm(x)
    for i=1:k-1
        z=A*X[:,i]
        dv[i]=X[:,i]⋅z
        # Three-term recursion
        if i==1
            z=z-dv[i]*X[:,i]
        else
            # z=z-dv[i]*X[:,i]-ev[i-1]*X[:,i-1]
            # Full reorthogonalization - once or even twice
            z=z-sum([(z⋅X[:,j])*X[:,j] for j=1:i])
            z=z-sum([(z⋅X[:,j])*X[:,j] for j=1:i])
        end
        μ=norm(z)
        if μ==0
            Tr=SymTridiagonal(dv[1:i-1],ev[1:i-2])
            return eigvals(Tr), X[:,1:i-1]*eigvecs(Tr), X[:,1:i-1], μ
        else
            ev[i]=μ
            X[:,i+1]=z/μ
        end
    end
    # Last step
    z=A*X[:,end]
    dv[end]=X[:,end]⋅z
    z=z-dv[end]*X[:,end]-ev[end]*X[:,end-1]
    μ=norm(z)
    Tr=SymTridiagonal(dv,ev)
    eigvals(Tr), X*eigvecs(Tr), X, μ
end

# Extract X
function myBidiagX(H::Matrix)
    m,n=size(H)
    T=typeof(H[1,1])
    X = Matrix{T}(I,m,n)
    v = Array{T}(undef,m)
    for j = n : -1 : 1
        v[j] = one(T)
        v[j+1:m] = H[j+1:m, j]
        γ = -2 / (v[j:m]⋅v[j:m])
        w = γ * X[j:m, j:n]'*v[j:m]
        X[j:m, j:n] = X[j:m, j:n] + v[j:m]*w'
    end
    X
end

# Extract Y
function myBidiagY(H::AbstractMatrix)
    n,m=size(H)
    T=typeof(H[1,1])
    Y = Matrix{T}(I,n,n)
    v = Array{T}(undef,n)
    for j = n-2 : -1 : 1
        v[j+1] = one(T)
        v[j+2:n] = H[j+2:n, j]
        γ = -2 / (v[j+1:n]⋅v[j+1:n])
        w = γ * Y[j+1:n, j+1:n]'*v[j+1:n]
        Y[j+1:n, j+1:n] = Y[j+1:n, j+1:n] + v[j+1:n]*w'
    end
    Y
end

# Bidiagonalization using Givens rotations
function myBidiagG(A1::Matrix)
    A=deepcopy(A1)
    m,n=size(A)
    T=typeof(A[1,1])
    X=Matrix{T}(I,m,m)
    Y=Matrix{T}(I,n,n)
    for j = 1 : n
        for i = j+1 : m
            G,r=givens(A,j,i,j)
            # Use the faster in-place variant
            # A=G*A
            lmul!(G,A)
            # X=G*X
            lmul!(G,X)
        end
        for i=j+2:n
            G,r=givens(A',j+1,i,j)
            # A*=adjoint(G)
            rmul!(A,adjoint(G))
            # Y*=adjoint(G)
            rmul!(Y,adjoint(G))
        end
    end
    X',Bidiagonal(diag(A),diag(A,1),'U'), Y
end

function myJacobiR(A1::AbstractMatrix)
    A=deepcopy(A1)
    m,n=size(A)
    T=typeof(A[1,1])
    V=Matrix{T}(I,n,n)
    # Tolerance for rotation
    tol=sqrt(map(T,n))*eps(T)
    # Counters
    p=n*(n-1)/2
    sweep=0
    pcurrent=0
    # First criterion is for standard accuracy, second one is for relative accuracy
    # while sweep<30 && vecnorm(A-diagm(diag(A)))>tol
    while sweep<30 && pcurrent<p
        sweep+=1
        # Row-cyclic strategy
        for i = 1 : n-1
            for j = i+1 : n
                # Compute the 2 x 2 sumbatrix of A'*A
                F=A[:,[i,j]]'*A[:,[i,j]]
                # Check the tolerance - the first criterion is standard,
                # the second one is for relative accuracy
                # if A[i,j]!=zero(T)
                #
                if abs(F[1,2])>tol*sqrt(F[1,1]*F[2,2])
                    # Compute c and s
                    τ=(F[1,1]-F[2,2])/(2*F[1,2])
                    t=sign(τ)/(abs(τ)+sqrt(1+τ^2))
                    c=1/sqrt(1+t^2)
                    s=c*t
                    G=LinearAlgebra.Givens(i,j,c,s)
                    A*=G'
                    # rmul!(A,adjoint(G))
                    V*=G'
                    # rmul!(A,adjoint(G))
                    pcurrent=0
                else
                    pcurrent+=1
                end
            end
        end
    end
    σ=[norm(A[:,k]) for k=1:n]
    for k=1:n
        A[:,k]./=σ[k]
    end
    A, σ, V
end

end

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/ModuleB.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
module TicketToRide
using Reexport
@reexport using Yao, YaoExtensions

include("utils.jl")
include("noise.jl")
include("fake.jl")
include("prune.jl")

end # module


text:
using CompilerPluginTools
using Test

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

@testset "interpreter" begin
    include("interp.jl")
end

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

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

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


text:
module RINEX

    using FileIO, Unitful, Dates, DataFrames

    export RINEXHeader, rinex_read

    include("NavMessageFile.jl")
    include("RINEXHeader.jl")
    include("RINEXData.jl")
    include("io.jl")
    include("read.jl")

    const EOH = "END OF HEADER       \n"
end


text:
using NeuroEvolution
using Cambrian
using Test

function test_identical(neur1::Neuron, neur2::Neuron)
    @test neur1.position == @test neur2.position
    @test neur1.activ_func == @test neur2.activ_func
    @test neur1.input == @test neur2.input
    @test neur1.output == @test neur2.output
end

function test_identical(g1::Gene, g2::Gene)
    same::Bool=true
    same &= g1.inno_nb == g2.inno_nb
    same &= g1.origin == g2.origin
    same &= g1.destination == g2.destination
    same &= g1.weight == g2.weight
    same &= g1.enabled == g2.enabled
    same
end

function test_identical(ind::NEATIndiv, ind2::NEATIndiv)
    same::Bool=true
    same &= length(ind.genes) == length(ind2.genes)
    same &= length(ind.neuron_pos) == length(ind2.neuron_pos)
    same &= all(ind.fitness .== ind2.fitness)

    if !same
        return false
    end
    for i in eachindex(ind.genes)
        same &= test_identical(ind.genes[i], ind2.genes[i])
    end

    for i in eachindex(ind.neuron_pos)
        same &=  ind.neuron_pos == ind2.neuron_pos
    end

    same &=  all(ind.fitness .== ind2.fitness)
    same
end

include("save_load.jl")
include("test_xor.jl")
include("individual.jl")

include("network.jl")
include("mutation.jl")
include("crossover.jl")

include("species.jl")
include("run.jl")

include("HyperNEAT.jl")


text:
module Model

export PlanktonModel
export TimeStep!

using CUDA
using StructArrays
using LinearAlgebra

using PlanktonIndividuals.Architectures
using PlanktonIndividuals.Grids
using PlanktonIndividuals.Fields
using PlanktonIndividuals.Biogeochemistry
using PlanktonIndividuals.Parameters
using PlanktonIndividuals.Plankton
using PlanktonIndividuals.Diagnostics

using PlanktonIndividuals: AbstractMode, CarbonMode, QuotaMode, MacroMolecularMode

import Base: show

include("timestepper.jl")
include("models.jl")
include("time_step.jl")

end


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

"""
k-forms can be viewed as a linear map from k vectors or a k-vector to a scalar.
Or as the "dual" input to some implicit (k k) tensor.   !me check that this is correct, I kinda just coined this "implicit tensor word" on the spot...
A (k k) tensor is a multi-linear map from k forms + k vectors to a scalar.

"""

export
apply,
𝑑

using Multivectors
using Multivectors: isazero
import Multivectors.kform
import Multivectors: ⋆

"""
    apply(α, v)

Apply k-form α to v.
"""
apply( α::K, v::V ) where { TF<:ZForm,T,N, K<:KVector{TF,N}, V<:KVector{T,N} } =
  mapreduce(*,+,scalar.(α),scalar.(v)) #mapreduce(α⋅reverse(u)

"""
    𝑑(k, jacobian)

Exterior Derivative of KVector of k-forms with current jacobian
"""
𝑑( k::K, ∂ϕ∂x ) where {K<:KVector} = mapreduce(b->𝑑(b, ∂ϕ∂x), +, k)

function 𝑑( b::B, ∂ϕ∂x ) where {T, N, B<:Blade{T,N}}
  bᵢ = factor(b)
  β = mapreduce(one, ∧, bᵢ[3:end])
  k = N-1
  α = bᵢ[1]*bᵢ[2]
  𝑑α = 𝑑(one(bᵢ[2]), ∂ϕ∂x)
  𝑑α∧β + (-1)^k*α∧𝑑(β, ∂ϕ∂x)
end

function 𝑑( b::B, ∂ϕ∂x ) where {T, B<:Blade{T,1}}
  s,eⱼ = factor(b)
  be = basis_1blades(b)
  j = subspace(eⱼ)[1]
  # expand 1-form on scalar out to ds*e₁ + ds*e₂...
  s*sum(Iterators.filter( !isazero, map( i->(∂ϕ∂x[j,i]*be[i])∧one(eⱼ), 1:length(be) )))
end

#== with 0-forms ==#
𝑑( s::R ) where {R<:Real} = zero(s)

"""
    𝑑(k)

Exterior derivative operator on KVector wedged with ZForms.
Differential of k-form.
"""
𝑑( k::K ) where {D, Z<:ZForm{D}, K<:KVector{Z}} = mapreduce(b->𝑑(b), +, k)

"differential of simple k-form"
function 𝑑( b::B ) where {T, N, B<:Blade{T,N}}
  bᵢ = factor(b)
  # 𝑑α∧β + (-1)ᵏα∧𝑑β
  β = mapreduce(e->e(one(T)), ∧, bᵢ[3:end])
  k = N-1
  α = bᵢ[1]
  sum(filter(!isazero, [𝑑(α∧bᵢ[2])∧β, (-1)^k*α∧𝑑(β)]))
end

𝑑( b::B ) where { T<:Real, B<:Blade{T,1} } = zero(T)

"differential of simple 1-form"
function 𝑑( b::B ) where {D, T<:ZForm{D}, B<:Blade{T,1}}
  j = subspace(b)
  be = basis_1blades(b)
  dxⱼ = be[j]
  k = b.x
  sum(Iterators.filter( !isazero, map( i->be[i](ZForm{D}(k.gradient[i], k.hessian[i,:]))∧one(dxⱼ), 1:D )))
end

"""
  𝑑( k, 𝐼 )

differential of 0-form. 𝐼 is the psuedovector of the algebra
"""
function 𝑑( k::Z, 𝐼::B ) where {D, Z<:ZForm{D}, B<:Blade}
  be = basis_1blades(𝐼)
  sum(Iterators.filter(!isazero, map( i->be[i](ZForm{D}(k.gradient[i], k.hessian[i,:])), 1:D)))
end

"""
  𝑑( k )

differential of 0-form
"""
function 𝑑( k::Z ) where {Z<:ZForm}
  @warn "resolving Blade types in top-level module via dual(1)"
  #𝐼 = 1pseudoscalar(Z) #dual(1)
  𝐼 = dual(1)
  𝑑(k, 𝐼)
end

const differential = 𝑑

Base.iszero(b::B) where {D,K,N,B<:KVector{ZForm{D},K,N}} = mapreduce(isazero, max, b)

𝑑(𝑓::F) where F<:Function = α->𝑑(𝑓(α))
⋆(𝑓::F) where F<:Function = α->⋆(𝑓(α))

δ(α) = ⋆(𝑑(⋆α))
Δ(α) = sum(filter(!isazero, [δ(𝑑(α)), 𝑑(δ(α))]))

#Base.:*(α::KF,u::K) where{F<:ZForm,KF<:KVector{F},K<:KVector} = sum(map(kfᵢ->kfᵢ.x.f(scalar.(u)...),α))

# very literal implementation of applying a k-form
"""
    apply(α, u, k)

Apply differential k-form α at point u to k-vector k.
Point u is a 1-vector.
"""
function apply(α::KF, u::P, k::K) where{F<:ZForm, KF<:KVector{F},
                                       U, P<:KVector{U,1},
                                       T, K<:KVector{T}}
  rα = raise.(α)
  αk = [ (eᵢ*reverse(one(eʲ)), j) for (i,eᵢ) in zip(1:length(k),k)
                                           for (j,eʲ) in zip(1:length(rα),untype.(rα)) ]
  mapreduce( ((s,j),)->s*apply(α[j].x, u), +, filter(!isazero∘first, αk))
end

apply(α::B, u, k::KF) where {F<:ZForm, KF<:KVector{F}, B<:Blade} = apply(KVector(α), u, k)

const CN = Union{KVector, Blade}
apply(α::CN, u::CN, k::CN) = apply(KVector(α), KVector(u), KVector(k))

"""
    apply(α, u)

Apply differential 0-form α at point u.
Point u is a 1-vector.
"""
apply(α::F, u::P) where {F<:ZForm, U, P<:KVector{U,1}} = α.f(coords(u)...)

"Constant 0-form.  Same value at any point"
apply(s::T, u::P) where {T<:Real, U, P<:KVector{U,1}} = s

Base.:*(s::T, b::B) where {T<:ZForm, B<:KVector} = B(s*b.k)
Base.:*(b::B, s::T) where {T<:ZForm, B<:KVector} = B(b.k*s)


text:
function value_iteration(
    mountaincar::MountainCar,
    params::Array{Float64};
    threshold = 1e-5,
    gamma = 1.0,
    max_iterations = 1e3,
)
    T = generate_transition_matrix(mountaincar, params)
    R = generate_cost_vector(mountaincar)
    value_iteration(
        mountaincar,
        T,
        R;
        threshold = threshold,
        gamma = gamma,
        max_iterations = max_iterations,
    )
end

function value_iteration(
    mountaincar::MountainCar,
    transition_matrix::Matrix{Int64},
    cost_vector::Vector{Float64};
    threshold = 1e-5,
    gamma = 1.0,
    max_iterations = 1e3,
)
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization + 1
    n_actions = 2
    V_old = copy(cost_vector)
    pi = zeros(Int64, n_states)
    V = 2 * V_old
    Q = zeros((n_states, n_actions))
    error_vec = get_error_vec(V, V_old)
    criterion = maximum(error_vec)
    count = 0
    while criterion >= threshold && count < max_iterations
        count += 1
        V_old = copy(V)
        for a = 1:n_actions
            Q[:, a] = cost_vector .+ gamma .* V_old[transition_matrix[:, a]]
        end
        idxs = argmin(Q, dims = 2)
        for s = 1:n_states
            pi[s] = idxs[s][2]
        end
        V = Q[idxs]
        error_vec = get_error_vec(V, V_old)
        criterion = maximum(error_vec)
    end
    # removing absorbing state
    pi[1:n_states-1], V[1:n_states-1], count != max_iterations
end

function value_iteration(
    mountaincar::MountainCar,
    transition_matrix::Matrix{Int64},
    cost_matrix::Matrix{Float64};
    threshold = 1e-5,
    gamma = 1.0,
    max_iterations = 1e3,
)
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization + 1
    n_actions = 2
    V_old = minimum(cost_matrix, dims=2)
    pi = zeros(Int64, n_states)
    V = 2 * V_old
    Q = zeros((n_states, n_actions))
    error_vec = get_error_vec(V, V_old)
    criterion = maximum(error_vec)
    count = 0
    while criterion >= threshold && count < max_iterations
        count += 1
        V_old = copy(V)
        for a = 1:n_actions
            Q[:, a] = cost_matrix[:, a] .+ gamma .* V_old[transition_matrix[:, a]]
        end
        idxs = argmin(Q, dims = 2)
        for s = 1:n_states
            pi[s] = idxs[s][2]
        end
        V = Q[idxs]
        error_vec = get_error_vec(V, V_old)
        criterion = maximum(error_vec)
    end
    # removing absorbing state
    pi[1:n_states-1], V[1:n_states-1], count != max_iterations
end

function value_iteration(
    mountaincar::MountainCar,
    params::MountainCarParameters,
)
    value_iteration(
        mountaincar,
        vec(params),
    )
end

function get_error_vec(V::Array{Float64}, V_old::Array{Float64})::Array{Float64}
    # abs.(V - V_old) ./ (V_old .+ 1e-50)
    abs.(V - V_old)
end

function generate_transition_matrix(mountaincar::MountainCar, params::Array{Float64})
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization + 1
    n_actions = 2
    T = zeros(Int64, (n_states, n_actions))
    actions = getActions(mountaincar)
    for a = 1:n_actions
        u = actions[a]
        for s = 1:n_states-1
            disc_state = idx_to_disc_state(mountaincar, s)
            if checkGoal(mountaincar, disc_state)
                T[s, a] = n_states  # absorbing state
            else
                s_next = disc_state_to_idx(
                    mountaincar,
                    step(mountaincar, disc_state, u, params)[1],
                )
                T[s, a] = s_next
            end
        end
        T[n_states, a] = n_states  # absorbing state
    end
    T
end

function generate_cost_vector(mountaincar::MountainCar)
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization + 1
    R = zeros(n_states)
    for s = 1:n_states-1
        R[s] = getCost(mountaincar, idx_to_disc_state(mountaincar, s))
    end
    R[n_states] = 0
    R
end

function generate_cost_matrix(mountaincar::MountainCar)
    n_states = mountaincar.position_discretization *
        mountaincar.speed_discretization + 1
    n_actions = 2
    R = zeros(n_states, n_actions)
    for s = 1:n_states-1
        for a = 1:n_actions
            R[s, a] = getCost(mountaincar, idx_to_disc_state(mountaincar, s))
        end
    end
    R[n_states, 1] = 0
    R[n_states, 2] = 0
    R
end

function iterative_policy_evaluation(
    mountaincar::MountainCar,
    params::Array{Float64},
    policy::Array{Int64};
    threshold = 1e-5,
    gamma = 1.0,
    max_iterations = 1e3,
)
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization + 1
    T = generate_transition_matrix(mountaincar, params)
    R = generate_cost_vector(mountaincar)
    V_old = copy(R)
    V = 2 * V_old
    error_vec = get_error_vec(V, V_old)
    criterion = maximum(error_vec)
    count = 0
    while criterion >= threshold && count < max_iterations
        count += 1
        V_old = copy(V)
        next = [T[s, policy[s]] for s in 1:n_states-1]
        push!(next, n_states)  # absorbing state
        V = R .+ gamma .* V_old[next]
        error_vec = get_error_vec(V, V_old)
        criterion = maximum(error_vec)
    end
    V[1:n_states-1]
end

function finite_horizon_value_iteration(
    mountaincar::MountainCar,
    params::Array{Float64},
    horizon::Int64,
)
    n_states = mountaincar.position_discretization * mountaincar.speed_discretization + 1
    n_actions = 2
    T = generate_transition_matrix(mountaincar, params)
    R = generate_cost_vector(mountaincar)
    V = [copy(R) for _ in 1:horizon+1]
    pi = [zeros(Int64, n_states) for _ in 1:horizon]
    for t in horizon:-1:1
        Q = zeros((n_states, n_actions))
        V_next = V[t+1]
        for a = 1:n_actions
            Q[:, a] = R .+ V_next[T[:, a]]
        end
        idxs = argmin(Q, dims = 2)
        for s = 1:n_states
            pi[t][s] = idxs[s][2]
        end
        V[t] = Q[idxs][:, 1]
    end
    pi, V
end


text:
"""
    radians_to_length(radians::Real, units::String="kilometers")::Real

Convert a distance measurement (assuming a spherical Earth) from radians to a more friendly unit.
Valid units: miles, nauticalmiles, inches, yards, meters, metres, kilometers, centimeters, feet
"""
function radians_to_length(radians::Real, units::String="kilometers")::Real
    factor = factors[units]
    if factor == nothing
        throw(error("$(units) is not a valid unit."))
    end
    return radians*factor
end

"""
    length_to_radians(distance::Real, units::String="kilometers")::Real

Convert a distance measurement (assuming a spherical Earth) from a real-world unit to radians.
Valid units: miles, nauticalmiles, inches, yards, meters, metres, kilometers, centimeters, feet
"""
function length_to_radians(distance::Real, units::String="kilometers")::Real
    factor =  factors[units]
    if factor == nothing
        throw(error("$(units) is not a valid unit."))
    end
    return distance / factor
end

"""
    length_to_degrees(distance::Real, units::String="kilometers")::Real

Convert a distance measurement (assuming a spherical Earth) from a real-world unit into degrees.
Valid units: miles, nauticalmiles, inches, yards, meters, metres, centimeters, kilometres, feet
"""
function length_to_degrees(distance::Real, units::String="kilometers")::Real
    return rad2deg(length_to_radians(distance, units))
end

"""
    convert_length(length::Real, originalUnit::String="kilometers", finalUnit::String="kilometers")::Real

Convert a length to the requested unit.
"""
function convert_length(length::Real, originalUnit::String="kilometers", finalUnit::String="kilometers")::Real
    return length >= 0 ? radians_to_length(length_to_radians(length, originalUnit), finalUnit) : error("'length' must be a positive number.")
end

"""
    convert_area(area::Real, originalUnit::String="meters", finalUnit::String="kilometers")::Real

Convert an area to the requested unit.
"""
function convert_area(area::Real, originalUnit::String="meters", finalUnit::String="kilometers")::Real
    if area < 0
        throw(error("'area' must be a positive number."))
    end

    startFactor =  area_factors[originalUnit]
    endFactor =  area_factors[finalUnit]

    return (area / startFactor) * endFactor
end

@inline function allequal(x)
    length(x) < 2 && return true
    e1 = x[1]
    i = 2
    @inbounds for i=2:length(x)
        x[i] == e1 || return false
    end
    return true
end

"""
    to_WGS84(pos::Point)

Convert 900913 x/y values to lon/lat.
"""
function to_WGS84(pos::Point)
    a = 6378137.0

    return [rad2deg(pos.coordinates[1]) / a,
        rad2deg((pi * 0.5) - atan(exp(-pos.coordinates[2] / a)))]
end

"""
    to_mercator(pos::Point)

Convert lon/lat values to 900913 x/y.
"""
function to_mercator(pos::Point)
    a = 6378137.0
    extent = 20037508.342789244

    coords = pos.coordinates

    adjusted = (abs(coords[1]) <= 180) ? coords[1] : (coords[1] - (sign(coords[1] * 360)))
    xy = [deg2rad(a * adjusted), a * log(atan((pi * 0.25) + (deg2rad(0.5 * coords[2]))))]

    xy[1] > extent && (xy[1] = extent)
    xy[1] < -extent && (xy[1] = -extent)
    xy[2] > extent && (xy[2] = extent)
    xy[2] < -extent && (xy[2] = -extent)

    return xy
end

"""
    area(geojson::T) where {T <: Union{AbstractFeatureCollection, AbstractGeometry}}
Take one or more features and return their area in square meters.
"""
function area(geojson::T) where {T <: Union{AbstractFeatureCollection, AbstractGeometry}}
    if geotype(geojson) === :FeatureCollection
        results = []
        for feat in geojson.features
            push!(results, calculateArea(feat.geometry))
        end

        return results
    else
        return calculateArea(geojson)
    end
end

function calculateArea(geom::T) where {T <: AbstractGeometry}
    total = 0
    type = geotype(geom)

    (type === :Point || type === :LineString || type === :MultiLineString) && return 0

    type === :Polygon && return polygonArea(geom.coordinates)

    if type === :MultiPolygon
        for i in 1:length(geom.coordinates)
            total += polygonArea(geom.coordinates[i])
        end

        return total
    end

    return 0
end

function polygonArea(coords)
    total = 0

    if length(coords) > 0
        total += abs(ringArea(coords[1]))

        for i in 2:length(coords)
            total -= abs(ringArea(coords[i]))
        end
    end

    return total
end

function ringArea(coords)
    total = 0
    lowIndex = 1
    midIndex = 1
    upIndex = 1

    if length(coords) > 2
        for i in 1:length(coords) - 1
            if i === length(coords) - 2
                lowIndex = length(coords) - 2
                midIndex = length(coords) - 1
                upIndex = 1
            elseif i === length(coords) - 1
                lowIndex = length(coords) - 1
                midIndex = 1
                upIndex = 2
            else
                lowIndex = i
                midIndex = i + 1
                upIndex = i + 2
            end

            p1 = coords[lowIndex]
            p2 = coords[midIndex]
            p3 = coords[upIndex]

            total += (p3[1] * (pi / 180) - (p1[1] * (pi / 180))) * sin(p2[2] * (pi / 180))
        end

        total = total * 6378137 * 6378137 / 2
    end

    return total
end

"""
    clean(geojson::AbstractGeometry, mutate::Bool=false)

Remove redundant coordinates from any GeoJSON Geometry.

# Examples
```jldoctest
julia> line = LineString([[0, 0], [0, 2], [0, 5], [0, 8], [0, 10], [0, 10]])
LineString(Array{Float64,1}[[0.0, 0.0], [0.0, 2.0], [0.0, 5.0], [0.0, 8.0], [0.0, 10.0], [0.0, 10.0]])

julia> clean(line)
LineString(Array{Float64,1}[[0.0, 0.0], [0.0, 10.0]])
```
"""
function clean(geojson::AbstractGeometry, mutate::Bool=false)

    type = geotype(geojson)
    isequal(type, :Point) && return geojson

    coords = []

    if isequal(type, :LineString) || isequal(type, :MultiPoint)
        coords = cleanline(geojson)
    elseif isequal(type, :Polygon)

        line = polygon_to_line(geojson)
        push!(coords, cleanline(line))

    elseif isequal(type, :MultiPolygon)
        points = []
        for poly in geojson.coordinates
            line = polygon_to_line(Polygon(poly))
            if geotype(line) === :MultiLineString
                for i in eachindex(line.coordinates)
                    push!(points, cleanline(LineString(line.coordinates[i])))
                end
            else
                push!(points, cleanline(line))
            end
        end
        push!(coords, points)
    elseif isequal(type, :MultiLineString)
        for line in geojson.coordinates
            push!(coords, cleanline(LineString(line)))
        end

    else
        throw(error("$(type) is not supported."))

    end

    if mutate
        geojson.coordinates = coords
    else
        geojson = deepcopy(geojson)
        geojson.coordinates = coords
    end

    return geojson
end

"""
    clean!(geojson::AbstractGeometry)

Remove redundant coordinates from any GeoJSON Geometry. It modifies the GeoJSON geometry in place.
"""
clean!(geojson::AbstractGeometry) = clean(geojson, true)

"""
    clean(geojson::AbstractFeature, mutate::Bool=false)

Remove redundant coordinates from any GeoJSON Feature.
"""
clean(geojson::AbstractFeature, mutate::Bool=false) = clean(geojson.geometry, mutate)

"""
    clean!(geojson::AbstractFeature)

Remove redundant coordinates from any GeoJSON Feature. It modifies the GeoJSON geometry in place.
"""
clean!(geojson::AbstractFeature) = clean!(geojson.geometry)

function cleanline(line)
    points = line.coordinates

    (length(points) == 2 && points[1] != points[2]) && return points

    new_points = []
    index = length(points) - 1
    np_length = length(new_points)

    push!(new_points, points[1])

    for i in 2:index
        prev = new_points[length(new_points)]

        if points[i][1] == prev[1] && points[i][2] == prev[2]
            continue
        else
            push!(new_points, points[i])
            np_length = length(new_points)

            if np_length > 2
                if isPointOnSegment(new_points[np_length - 2], new_points[np_length], new_points[np_length- 1])
                    splice!(new_points, length(new_points)-1)
                end
            end
        end
    end

    push!(new_points, points[length(points)])

    np_length = length(new_points)

    (points[1] == points[end] && np_length < 4) && throw(error("invalid polygon."))
    if isPointOnSegment(new_points[np_length - 2], new_points[np_length], new_points[np_length- 1])
        splice!(new_points, length(new_points)-1)
    end

    return new_points
end


text:
########################################################
#
# Robust MDP
#
# Andrew Keith
# 24 Jan 2018
# UNCLASSIFIED, Distribution Unlimited (A)
# MIT license
#
########################################################

############################
# Define constants
###########################

# Epochs
# finite horizon
N = 5

# States
S = [1,2,3]

# Actions
A = [1,2,3]

# Rewards
R = [1.0 10 10; 10 10 10; 2 2 2]

rN = [1; 10.0; 2.0]

# Uncertainty Sets
P = zeros(Float64,3,3,3)
P[:,:,1] = [0.8 0.1 0.1; 0.8 0.1 0.1; 0.8 0.1 0.1]
P[:,:,2] = [0.1 0.8 0.1; 0.1 0.8 0.1; 0.1 0.8 0.1]
P[:,:,3] = [0.1 0.1 0.8; 0.1 0.1 0.8; 0.1 0.1 0.8]

# Bisection Algorithm - Entropy Uncertainty Sets

function sig(lam::Float64,beta::Float64, v::Array{Float64,1}, q::Array{Float64,1})
  lam * log(q' * exp.(v / lam)) + beta * lam
end

function biEntropy(del::Float64, beta::Float64, vimax::Float64, v::Array{Float64,1}, q::Array{Float64,1})
  result = 0.0
  if beta == 0
    result = q' * v
  else
    laml = 0 + del
    lamu = (vimax - q' * v) / beta
    lam = (laml + lamu) / 2

    while lamu - laml > del
      lam = (laml + lamu) / 2
      if sig(lam,beta,v,q) < sig(laml,beta,v,q)
        laml = copy(lam)
      else
        lamu = copy(lam)
      end
    end
    result = sig(lam,beta,v,q)
  end
  result
end

q = P[1,:,1]
v = rN
vimax = maximum(v)

tempsigs = Vector{Float64}(1)

nS = size(S,1)
nA = size(A,1)
eps = 1e-5
del = eps / N
policy = zeros(Int64,nS,N -1)
value = zeros(Float64,nS,N)
value[:,N] = rN
beta = 12.0

for t = (N-1):-1:1
  v = value[:,t + 1]
  vimax = maximum(v)
  for i = 1:nS
    vmin = Inf
    amin = 0
    for j = 1:nA
      sighat = biEntropy(del,beta,vimax,v,P[i,:,j])
      # println("i = $i, j = $j, vmin = $vmin, sighat = $sighat")
      if (R[i,j] + sighat) < vmin
        vmin = R[i,j] + sighat
        amin = j
      end
    end
    value[i,t] = vmin
    policy[i,t] = amin
  end
end

println(policy)
println(value)

max_stars_repo_name: ajkeith/ImpreciseMDP

max_issues_repo_name: ajkeith/ImpreciseMDP

max_forks_repo_name: ajkeith/ImpreciseMDP


text:
module GViz

repmat = repeat

pkgpath(paths...) = joinpath(@__DIR__, "..", paths...)

include("params.jl")

include("RefFrame.jl")

include("render.jl")
include("ticks2.jl")
include("geometry.jl")
include("colors.jl")

include("datatrack.jl")
include("genomicaxis.jl")

include("GenomicFeatureTrack.jl")

include("AbstractGenomicFeature.jl")
include("GTF.jl")

end