text:
using Transducers
using DiffEqJump

struct MarginalEnsemble{JP <: DiffEqBase.AbstractJumpProblem,U0}
    jump_problem::JP
    dist::TrajectoryDistribution
    dtimes::Vector{Float64}

    # Handling the initial condition. Can be either a vector that specifies the exact initial
    # condition, or an EmpiricalDistribution specifying a probability distribution over initial
    # conditions.
    u0::U0
end

struct ConditionalEnsemble{JP <: DiffEqBase.AbstractJumpProblem,IX,DX}
    jump_problem::JP
    dist::TrajectoryDistribution
    indep_idxs::IX
    dep_idxs::DX
    dtimes::Vector{Float64}
end

struct SXconfiguration{uTs,uTx,tType}
    s_traj::Trajectory{uTs,tType}
    x_traj::Trajectory{uTx,tType}
end

Base.copy(c::SXconfiguration) = SXconfiguration(copy(c.s_traj), copy(c.x_traj))
struct SRXconfiguration{uTs,Utr,Utx,tType}
    s_traj::Trajectory{uTs,tType}
    r_traj::Trajectory{Utr,tType}
    x_traj::Trajectory{Utx,tType}
end

Base.copy(c::SRXconfiguration) = SRXconfiguration(copy(c.s_traj), copy(c.r_traj), copy(c.x_traj))

ensurevec(a::AbstractVector) = a
ensurevec(a) = SVector(a)

function initial_log_likelihood(ensemble::ConditionalEnsemble, u0, x_traj)
    0.0
end

function initial_log_likelihood(ensemble::MarginalEnsemble, u0::AbstractVector, x_traj::Trajectory)
    s0 = u0[1]
    x0 = x_traj.u[1][1]
    if ensemble.u0 isa EmpiricalDistribution
        # compute log P(x0 | s0) = log P(s0, x0) - log ∑ₓ P(s0, x)
        logpdf(ensemble.u0, [s0, x0]) - log(sum(x -> pdf(ensemble.u0, [s0, x]), ensemble.u0.axes[2]))
    else
        0.0
    end
end

function Base.getindex(conf::SXconfiguration, index)
    merge_trajectories(conf.s_traj, conf.x_traj) |> Map((u, t, i)::Tuple -> (ensurevec(u[index]), t, i)) |> Thin() |> collect_trajectory
end

function Base.getindex(conf::SRXconfiguration, index)
    merge_trajectories(conf.s_traj, conf.r_traj, conf.x_traj) |> Map((u, t, i)::Tuple -> (ensurevec(u[index]), t, i)) |> Thin() |> collect_trajectory
end

sample_initial_condition(u0::AbstractVector) = copy(u0)
sample_initial_condition(u0::EmpiricalDistribution) = rand(u0)

sample_initial_condition(ens::MarginalEnsemble) = ens.u0 isa EmpiricalDistribution ? rand(ens.u0)[[1]] : ens.u0
sample_initial_condition(ens::ConditionalEnsemble) = ens.jump_problem.prob.u0

abstract type JumpNetwork end

"""
    SimpleSystem(input_network, output_network, u0, ps, px, dtimes, dist=nothing)

A `SimpleSystem` consists of two ModelingToolkit `ReactionSystem`s:
One that generates the input trajectories, and one that generates the
output trajectories.
"""
struct SimpleSystem <: JumpNetwork
    sn::ReactionSystem
    xn::ReactionSystem

    u0::Union{AbstractVector,EmpiricalDistribution}

    ps::AbstractVector
    px::AbstractVector

    dtimes

    jump_problem::AbstractJumpProblem
    dist::TrajectoryDistribution
end

function SimpleSystem(sn, xn, u0, ps, px, dtimes, dist=nothing)
    joint = merge(sn, xn)

    tp = (first(dtimes), last(dtimes))
    p = vcat(ps, px)
    dprob = DiscreteProblem(joint, sample_initial_condition(u0), tp, p)
    jprob = JumpProblem(convert(ModelingToolkit.JumpSystem, joint), dprob, Direct(), save_positions=(false, false))

    if dist === nothing
        update_map = build_update_map(joint, xn)
        dist = distribution(joint, p; update_map)
    end

    SimpleSystem(sn, xn, u0, ps, px, dtimes, jprob, dist)
end

struct CompiledSimpleSystem{JP}
    system::SimpleSystem
    marginal_ensemble::MarginalEnsemble{JP}
end

compile(s::SimpleSystem) = CompiledSimpleSystem(s, MarginalEnsemble(s))
marginal_density(csx::CompiledSimpleSystem, algorithm, conf::SXconfiguration) = log_marginal(simulate(algorithm, conf, csx.marginal_ensemble))
conditional_density(csx::CompiledSimpleSystem, algorithm, conf::SXconfiguration) = -energy_difference(conf, csx.marginal_ensemble)

"""
    ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes, dist=nothing; aggregator=Direct())

A `ComplexSystem` consists of three ModelingToolkit `ReactionSystem`s:
One that generates the input trajectories, and one that generates the
output trajectories.
"""
struct ComplexSystem <: JumpNetwork
    sn::ReactionSystem
    rn::ReactionSystem
    xn::ReactionSystem

    u0::Union{AbstractVector,EmpiricalDistribution}

    ps::AbstractVector
    pr::AbstractVector
    px::AbstractVector

    dtimes

    jump_problem
    dist::TrajectoryDistribution
end

function ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes, dist=nothing; aggregator=Direct())
    joint = merge(merge(sn, rn), xn)

    tp = (first(dtimes), last(dtimes))
    p = vcat(ps, pr, px)
    dprob = DiscreteProblem(joint, sample_initial_condition(u0), tp, p)
    jprob = JumpProblem(convert(ModelingToolkit.JumpSystem, joint), dprob, aggregator, save_positions=(false, false))

    if dist === nothing
        update_map = build_update_map(joint, xn)
        dist = distribution(joint, p; update_map)
    end

    ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes, jprob, dist)
end

struct CompiledComplexSystem{JP,JPC,IXC,DXC}
    system::ComplexSystem
    marginal_ensemble::MarginalEnsemble{JP}
    conditional_ensemble::ConditionalEnsemble{JPC,IXC,DXC}
end

compile(s::ComplexSystem) = CompiledComplexSystem(s, MarginalEnsemble(s), ConditionalEnsemble(s))
marginal_density(csrx::CompiledComplexSystem, algorithm, conf::SRXconfiguration) = log_marginal(simulate(algorithm, marginal_configuration(conf), csrx.marginal_ensemble))
conditional_density(csrx::CompiledComplexSystem, algorithm, conf::SRXconfiguration) = log_marginal(simulate(algorithm, conf, csrx.conditional_ensemble))

tspan(sys::JumpNetwork) = (first(sys.dtimes), last(sys.dtimes))

reaction_network(system::SimpleSystem) = merge(system.sn, system.xn)
reaction_network(system::ComplexSystem) = merge(merge(system.sn, system.rn), system.xn)

function _solve(system::SimpleSystem)
    sol = solve(system.jump_problem, SSAStepper())
end

function generate_configuration(system::SimpleSystem)
    joint = reaction_network(system)
    jp = remake(system.jump_problem, u0=sample_initial_condition(system.u0))
    integrator = init(jp, SSAStepper(), tstops=())
    trajectory_iter = SSAIter(integrator)
    trajectory = trajectory_iter |> collect_trajectory

    s_spec = independent_species(system.sn)
    s_idxs = sort(SVector(species_indices(joint, s_spec)...))
    s_traj = sub_trajectory(trajectory, s_idxs)

    x_spec = independent_species(system.xn)
    x_idxs = sort(SVector(species_indices(joint, x_spec)...))
    x_traj = sub_trajectory(trajectory, x_idxs)

    SXconfiguration(s_traj, x_traj)
end

function _solve(system::ComplexSystem)
    sol = solve(system.jump_problem, SSAStepper())
end

function generate_configuration(system::ComplexSystem)
    # we first generate a joint SRX trajectory
    joint = reaction_network(system)
    jp = remake(system.jump_problem, u0=sample_initial_condition(system.u0))
    integrator = init(jp, SSAStepper(), tstops=())
    trajectory = SSAIter(integrator) |> collect_trajectory

    # then we extract the signal
    s_spec = independent_species(system.sn)
    s_idxs = sort(SVector(species_indices(joint, s_spec)...))
    s_traj = sub_trajectory(trajectory, s_idxs)

    # the R trajectory
    r_spec = independent_species(system.rn)
    r_idxs = sort(SVector(species_indices(joint, r_spec)...))
    r_traj = sub_trajectory(trajectory, r_idxs)

    # finally we extract the X part from the SRX trajectory
    x_spec = independent_species(system.xn)
    x_idxs = sort(SVector(species_indices(joint, x_spec)...))
    x_traj = sub_trajectory(trajectory, x_idxs)

    SRXconfiguration(s_traj, r_traj, x_traj)
end

function build_update_map(joint::ReactionSystem, xn::ReactionSystem)
    update_map = Int[]
    spmap = Catalyst.speciesmap(joint)
    mapper = (x, y)::Pair -> spmap[x] => y

    unique_netstoich = unique(map(r -> mapper.(r.netstoich), Catalyst.reactions(xn)))

    for react in Catalyst.reactions(joint)
        new_index = 0
        for (k, un) in enumerate(unique_netstoich)
            if mapper.(react.netstoich) == un
                new_index = k
                break
            end
        end

        push!(update_map, new_index)
    end
    update_map
end

function MarginalEnsemble(system::SimpleSystem)
    joint = reaction_network(system)
    s_idxs = species_indices(joint, Catalyst.species(system.sn))

    dprob = DiscreteProblem(system.sn, sample_initial_condition(system.u0)[s_idxs], tspan(system), system.ps)
    jprob = JumpProblem(convert(ModelingToolkit.JumpSystem, system.sn), dprob, Direct(), save_positions=(false, false))

    u0 = system.u0 isa EmpiricalDistribution ? system.u0 : system.u0[s_idxs]
    MarginalEnsemble(jprob, system.dist, collect(system.dtimes), u0)
end

function MarginalEnsemble(system::ComplexSystem)
    sr_network = merge(system.sn, system.rn)
    joint = merge(sr_network, system.xn)
    sr_idxs = species_indices(joint, Catalyst.species(sr_network))

    dprob = DiscreteProblem(sr_network, sample_initial_condition(system.u0)[sr_idxs], tspan(system), vcat(system.ps, system.pr))
    jprob = JumpProblem(convert(ModelingToolkit.JumpSystem, sr_network), dprob, Direct(), save_positions=(false, false))

    MarginalEnsemble(jprob, system.dist, collect(system.dtimes), system.u0[sr_idxs])
end

function ConditionalEnsemble(system::ComplexSystem)
    joint = merge(merge(system.sn, system.rn), system.xn)
    r_idxs = species_indices(joint, Catalyst.species(system.rn))
    dprob = DiscreteProblem(system.rn, sample_initial_condition(system.u0)[r_idxs], (first(system.dtimes), last(system.dtimes)), system.pr)
    jprob = JumpProblem(convert(ModelingToolkit.JumpSystem, system.rn), dprob, Direct(), save_positions=(false, false))

    indep_species = independent_species(system.rn)
    indep_idxs = species_indices(system.rn, indep_species)

    dep_species = dependent_species(system.xn)
    dep_idxs = indexin(species_indices(system.rn, dep_species), indep_idxs)

    ConditionalEnsemble(jprob, system.dist, indep_idxs, dep_idxs, collect(system.dtimes))
end

# returns a list of species in `a` that also occur in `b`
function intersecting_species(a::ReactionSystem, b::ReactionSystem)
    intersect(Catalyst.species(a), Catalyst.species(b))
end

# returns a list of species in `a` that are not in `b`
function unique_species(a::ReactionSystem, b::ReactionSystem)
    setdiff(Catalyst.species(a), Catalyst.species(b))
end

function species_indices(rs::ReactionSystem, species)
    getindex.(Ref(Catalyst.speciesmap(rs)), species)
end

function independent_species(rs::ReactionSystem)
    i_spec = []
    for r in Catalyst.reactions(rs)
        push!(i_spec, getindex.(r.netstoich, 1)...)
    end
    unique(s for s∈i_spec)
end

function dependent_species(rs::ReactionSystem)
    setdiff(Catalyst.species(rs), independent_species(rs))
end

function sample(configuration::T, system::MarginalEnsemble; θ=0.0)::T where T <: SXconfiguration
    if θ != 0.0
        error("can only use DirectMC with JumpNetwork")
    end
    jprob = remake(system.jump_problem, u0=sample_initial_condition(system))
    integrator = DiffEqBase.init(jprob, SSAStepper(), tstops=(), numsteps_hint=0)
    iter = SSAIter(integrator)
    s_traj = collect_trajectory(iter)
    SXconfiguration(s_traj, configuration.x_traj)
end

function collect_samples(initial::SXconfiguration, ensemble::MarginalEnsemble, num_samples::Int)
    result = Array{Float64,2}(undef, length(ensemble.dtimes), num_samples)
    for result_col ∈ eachcol(result)
        u0 = sample_initial_condition(ensemble)
        jprob = remake(ensemble.jump_problem, u0=u0)
        integrator = DiffEqBase.init(jprob, SSAStepper(), tstops=(), numsteps_hint=0)
        iter = SSAIter(integrator) |> Map((u, t, i)::Tuple -> (u, t, 0))
        cumulative_logpdf!(result_col, ensemble.dist, merge_trajectories(iter, initial.x_traj), ensemble.dtimes)
        result_col .+= initial_log_likelihood(ensemble, u0, initial.x_traj)
    end

    result
end

function propagate(conf::SXconfiguration, ensemble::MarginalEnsemble, u0, tspan::Tuple)
    jprob = remake(ensemble.jump_problem, u0=u0, tspan=tspan)
    integrator = DiffEqBase.init(jprob, SSAStepper(), tstops=(), numsteps_hint=0)
    ix1 = max(searchsortedfirst(conf.x_traj.t, tspan[1]) - 1, 1)
    iter = SSAIter(integrator) |> Map((u, t, i)::Tuple -> (u, t, 0))

    log_weight = trajectory_energy(ensemble.dist, iter |> MergeWith(conf.x_traj, ix1), tspan=tspan)

    copy(integrator.u), log_weight
end

function energy_difference(configuration::SXconfiguration, ensemble::MarginalEnsemble)
    log_p0 = initial_log_likelihood(ensemble, configuration.s_traj.u[1], configuration.x_traj)
    - cumulative_logpdf(ensemble.dist, configuration.s_traj |> MergeWith(configuration.x_traj), ensemble.dtimes) .- log_p0
end

function sample(configuration::SRXconfiguration, system::ConditionalEnsemble; θ=0.0)
    if θ != 0.0
        error("can only use DirectMC with JumpNetwork")
    end
    driven_jp = DrivenJumpProblem(system.jump_problem, configuration.s_traj)
    integrator = init(driven_jp)
    iter = SSAIter(integrator)
    rtraj = sub_trajectory(iter, system.indep_idxs)
    SRXconfiguration(configuration.s_traj, rtraj, configuration.x_traj)
end

function collect_samples(initial::SRXconfiguration, ensemble::ConditionalEnsemble, num_samples::Int)
    driven_jp = DrivenJumpProblem(ensemble.jump_problem, initial.s_traj)

    result = Array{Float64,2}(undef, length(ensemble.dtimes), num_samples)
    for result_col ∈ eachcol(result)
        integrator = init(driven_jp)
        iter = SSAIter(integrator) |> Map((u, t, i)::Tuple -> (u, t, 0))
        cumulative_logpdf!(result_col, ensemble.dist, merge_trajectories(iter, initial.x_traj), ensemble.dtimes)
    end

    result
end

function simulate(algorithm::DirectMCEstimate, initial::Union{SXconfiguration,SRXconfiguration}, system)
    samples = collect_samples(initial, system, algorithm.num_samples)
    DirectMCResult(samples)
end

function create_integrator(conf::SRXconfiguration, ensemble::ConditionalEnsemble, u0, tspan::Tuple)
    jprob = remake(ensemble.jump_problem, u0=u0, tspan=tspan)
    driven_jp = DrivenJumpProblem(jprob, conf.s_traj)
    init(driven_jp)
end

function propagate(conf::SRXconfiguration, ensemble::ConditionalEnsemble, u0, tspan::Tuple)
    integrator = create_integrator(conf, ensemble, u0, tspan)
    iter = SSAIter(integrator) |> Map((u, t, i)::Tuple -> (u, t, 0))
    ix1 = max(searchsortedfirst(conf.x_traj.t, tspan[1]), 1)

    log_weight = trajectory_energy(ensemble.dist, iter |> MergeWith(conf.x_traj, ix1), tspan=tspan)

    copy(integrator.u), log_weight
end

function energy_difference(configuration::SRXconfiguration, ensemble::ConditionalEnsemble)
    traj = merge_trajectories(configuration.s_traj, configuration.r_traj, configuration.x_traj)
    -cumulative_logpdf(ensemble.dist, traj, ensemble.dtimes)
end

function marginal_configuration(conf::SRXconfiguration)
    new_s = conf.s_traj |> MergeWith(conf.r_traj) |> Map((u, t, i)::Tuple -> (SVector(u...), t, i)) |> collect_trajectory
    SXconfiguration(new_s, conf.x_traj)
end

# MCMC Moves in trajectory space

mutable struct TrajectoryChain{Ensemble} <: MarkovChain
    ensemble::Ensemble
    # interaction parameter
    θ::Float64

    # to save statistics
    last_regrowth::Float64
    accepted_list::Vector{Float64}
    rejected_list::Vector{Float64}
end

chain(ensemble; θ::Real=1.0) = TrajectoryChain(ensemble, θ, 0.0, Float64[], Float64[])

# reset statistics
function reset(pot::TrajectoryChain)
    resize!(pot.accepted_list, 0)
    resize!(pot.rejected_list, 0)
end

function accept(pot::TrajectoryChain)
    push!(pot.accepted_list, pot.last_regrowth)
end

function reject(pot::TrajectoryChain)
    push!(pot.rejected_list, pot.last_regrowth)
end

function energy(conf::SXconfiguration, chain::TrajectoryChain; θ=chain.θ)
    if θ > zero(θ)
        -θ * trajectory_energy(chain.ensemble.dist, conf.s_traj |> MergeWith(conf.x_traj))
    else
        0.0
    end
end

function energy(conf::SRXconfiguration, chain::TrajectoryChain; θ=chain.θ)
    if θ > zero(θ)
        traj = merge_trajectories(conf.s_traj, conf.r_traj, conf.x_traj)
        θ * logpdf(chain.ensemble.dist, traj)
    else
        0.0
    end
end

function propose!(new_conf, old_conf, chain::TrajectoryChain)
    chain.last_regrowth = propose!(new_conf, old_conf, chain.ensemble)
    new_conf
end
propose!(new_conf::SXconfiguration, old_conf::SXconfiguration, ensemble::MarginalEnsemble) = propose!(new_conf.s_traj, old_conf.s_traj, ensemble)

function propose!(new_traj::Trajectory, old_traj::Trajectory, ensemble::MarginalEnsemble)
    jump_problem = ensemble.jump_problem

    regrow_duration = rand() * duration(old_traj)

    if rand(Bool)
        shoot_forward!(new_traj, old_traj, jump_problem, old_traj.t[end] - regrow_duration)
    else
        shoot_backward!(new_traj, old_traj, jump_problem, old_traj.t[begin] + regrow_duration)
    end

    regrow_duration
end

function shoot_forward!(new_traj::Trajectory, old_traj::Trajectory, jump_problem::DiffEqBase.AbstractJumpProblem, branch_time::Real)
    branch_value = old_traj(branch_time)
    branch_point = searchsortedfirst(old_traj.t, branch_time)
    tspan = (branch_time, old_traj.t[end])

    empty!(new_traj.u)
    empty!(new_traj.t)
    empty!(new_traj.i)
    append!(new_traj.u, @view old_traj.u[begin:branch_point - 1])
    append!(new_traj.t, @view old_traj.t[begin:branch_point - 1])
    append!(new_traj.i, @view old_traj.i[begin:branch_point - 1])

    jump_problem = DiffEqBase.remake(jump_problem; u0=branch_value, tspan=tspan)
    integrator = DiffEqBase.init(jump_problem, SSAStepper())
    iter = SSAIter(integrator)
    new_branch = collect_trajectory(iter)

    append!(new_traj.t, new_branch.t)
    append!(new_traj.u, new_branch.u)
    append!(new_traj.i, new_branch.i)
    nothing
end

function shoot_backward!(new_traj::Trajectory, old_traj::Trajectory, jump_problem::DiffEqBase.AbstractJumpProblem, branch_time::Real)
    branch_value = old_traj(branch_time)
    branch_point = searchsortedfirst(old_traj.t, branch_time)
    tspan = (old_traj.t[begin], branch_time)

    jump_problem = DiffEqBase.remake(jump_problem; u0=branch_value, tspan=tspan)
    integrator = DiffEqBase.init(jump_problem, SSAStepper())
    iter = SSAIter(integrator)
    new_branch = collect_trajectory(iter)

    empty!(new_traj.u)
    empty!(new_traj.t)
    empty!(new_traj.i)

    append!(new_traj.u, @view new_branch.u[end - 1:-1:begin])
    append!(new_traj.u, @view old_traj.u[branch_point:end])
    append!(new_traj.i, @view new_branch.i[end - 1:-1:begin])
    append!(new_traj.i, @view old_traj.i[branch_point:end])

    for rtime in @view new_branch.t[end:-1:begin + 1]
        push!(new_traj.t, branch_time - rtime)
    end
    append!(new_traj.t, @view old_traj.t[branch_point:end])
    nothing
end


text:
# Test get/set objective coefficients in LP

using Gurobi, Test

@testset "LP 03" begin
	env = Gurobi.Env()
	setparam!(env, "OutputFlag", 0)

	# original model
	#
	#   maximize  2x + 2y
	#
	#	s.t. 0.2 <= x, y <= 1
	#

	model = gurobi_model(env;
		name="lp_03",
		sense=:maximize,
		f=[2.0, 2.0],
		lb=[0.2, 0.2],
		ub=[1.0, 1.0])

	lb_ = lowerbounds(model)
	ub_ = upperbounds(model)
	c_ = objcoeffs(model)

	@test lb_ == [0.2, 0.2]
	@test ub_ == [1.0, 1.0]
	@test c_ == [2.0, 2.0]

	optimize(model)

	@test get_solution(model) == [1.0, 1.0]
	@test get_objval(model) == 4.0

	# change objective (warm start)
	#
	#	maximize x - y
	#
	#	s.t. 0.2 <= x, y <= 1
	#

	set_objcoeffs!(model, [1, -1])
	update_model!(model)

	c_ = objcoeffs(model)
	@test c_ == [1.0, -1.0]

	optimize(model)

	@test get_solution(model) == [1.0, 0.2]
	@test get_objval(model) == 0.8

	GC.gc()
end

ext: jl


text:
using Printf
using CUDA
using Oceananigans
using Oceananigans.TurbulenceClosures: Explicit, VerticallyImplicit
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary
using LinearAlgebra

function boundary_clustered(N, L, ini)
    Δz(k)   = k < N / 2 + 1 ? 2 / (N - 1) * (k - 1) + 1 : - 2 / (N - 1) * (k - N) + 1
    z_faces = zeros(N+1)
    for k = 2:N+1
        z_faces[k] = z_faces[k-1] + Δz(k-1)
    end
    z_faces = z_faces ./ z_faces[end] .* L .+ ini
    return z_faces
end

function center_clustered(N, L, ini)
    Δz(k)   = k < N / 2 + 1 ? 2 / (N - 1) * (k - 1) + 1 : - 2 / (N - 1) * (k - N) + 1
    z_faces = zeros(N+1)
    for k = 2:N+1
        z_faces[k] = z_faces[k-1] + 3 - Δz(k-1)
    end
    z_faces = z_faces ./ z_faces[end] .* L .+ ini
    return z_faces
end

grid = RectilinearGrid(GPU(), size=(512, 256),
                       x = (-10, 10),
                       z = (0, 5),
                topology = (Periodic, Flat, Bounded))

# Gaussian bump of width "1"
bump(x, y, z) = z < exp(-x^2)

@inline show_name(t) = t isa Explicit ? "explicit" : "implicit"

grid_with_bump = ImmersedBoundaryGrid(grid, GridFittedBoundary(bump))

# Tidal forcing
tidal_forcing(x, y, z, t) = 1e-4 * cos(t)

for time_stepper in (ExplicitTimeDiscretization, VerticallyImplicitTimeDiscretization)

    model = HydrostaticFreeSurfaceModel(grid = grid_with_bump,
                                        momentum_advection = CenteredSecondOrder(),
                                        free_surface = ExplicitFreeSurface(gravitational_acceleration=10),
                                        closure = ScalarDiffusivity(time_stepper, ν=1e-2, κ=1e-2),
                                        tracers = :b,
                                        buoyancy = BuoyancyTracer(),
                                        coriolis = FPlane(f=sqrt(0.5)),
                                        forcing = (u = tidal_forcing,))

    # Linear stratification
    set!(model, b = (x, y, z) -> 4 * z)

    progress_message(s) = @info @sprintf("[%.2f%%], iteration: %d, time: %.3f, max|w|: %.2e",
                                100 * s.model.clock.time / s.stop_time, s.model.clock.iteration,
                                s.model.clock.time, maximum(abs, model.velocities.w))

    gravity_wave_speed = sqrt(model.free_surface.gravitational_acceleration * grid.Lz)

    Δt = CUDA.@allowscalar 0.1 * minimum(grid.Δxᶜᵃᵃ) / gravity_wave_speed

    simulation = Simulation(model, Δt = Δt, stop_time = 50000Δt)

    serialize_grid(file, model) = file["serialized/grid"] = model.grid.grid

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                        schedule = TimeInterval(0.1),
                                                        prefix = "internal_tide_$(show_name(time_stepper))",
                                                        init = serialize_grid,
                                                        force = true)

    simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(10))

    run!(simulation)

    @info """
        Simulation complete.
        Output: $(abspath(simulation.output_writers[:fields].filepath))
    """
end

using JLD2
using Plots
ENV["GKSwstype"] = "100"

function nice_divergent_levels(c, clim; nlevels=20)
    levels = range(-clim, stop=clim, length=nlevels)
    cmax = maximum(abs, c)
    clim < cmax && (levels = vcat([-cmax], levels, [cmax]))
    return (-clim, clim), levels
end

function nan_solid(x, z, u, bump)
    Nx, Nz = size(u)
    x2 = reshape(x, Nx, 1)
    z2 = reshape(z, 1, Nz)
    u[bump.(x2, 0, z2)] .= NaN
    return nothing
end

function visualize_internal_tide_simulation(prefix)

    filename = prefix * ".jld2"
    file = jldopen(filename)

    grid = adapt(CPU(), file["serialized/grid"])

    bump(x, y, z) = z < exp(-x^2)

    xu, yu, zu = nodes((Face, Center, Center), grid)
    xw, yw, zw = nodes((Center, Center, Face), grid)
    xb, yb, zb = nodes((Center, Center, Center), grid)

    b₀ = file["timeseries/b/0"][:, 1, :]

    iterations = parse.(Int, keys(file["timeseries/t"]))

    anim = @animate for (i, iter) in enumerate(iterations)

        @info "Plotting iteration $iter of $(iterations[end])..."

        u = file["timeseries/u/$iter"][:, 1, :]
        w = file["timeseries/w/$iter"][:, 1, :]
        b = file["timeseries/b/$iter"][:, 1, :]
        t = file["timeseries/t/$iter"]

        b′ = b .- b₀

        wlims, wlevels = nice_divergent_levels(w, 1e-4)
        ulims, ulevels = nice_divergent_levels(u, 1e-3)
        blims, blevels = nice_divergent_levels(b′, 1e-4)

        nan_solid(xu, zu, u, bump)
        nan_solid(xw, zw, w, bump)
        nan_solid(xb, zb, b, bump)

        u_title = @sprintf("x velocity, t = %.2f", t)

        u_plot = contourf(xu, zu, u'; title = u_title,                  color = :balance, aspectratio = :equal, linewidth = 0, levels = ulevels, clims = ulims)
        w_plot = contourf(xw, zw, w'; title = "z velocity",             color = :balance, aspectratio = :equal, linewidth = 0, levels = wlevels, clims = wlims)
        b_plot = contourf(xb, zb, b′'; title = "buoyancy perturbation", color = :balance, aspectratio = :equal, linewidth = 0, levels = blevels, clims = blims)

        plot(u_plot, w_plot, b_plot, layout = (3, 1), size = (1200, 1200))
    end

    mp4(anim, prefix * ".mp4", fps = 16)

    close(file)
end

function plot_implicit_explicit_difference(filename)

    file_explicit = jldopen(filename * "_explicit.jld2")
    file_implicit = jldopen(filename * "_implicit.jld2")

    iterations = parse.(Int, keys(file_explicit["timeseries/t"]))

    comparison_u = zeros(length(iterations))
    comparison_w = zeros(length(iterations))
    comparison_b = zeros(length(iterations))

    for (i, iter) in enumerate(iterations)

        u_explicit = file_explicit["timeseries/u/$iter"][:, 1, :]
        w_explicit = file_explicit["timeseries/w/$iter"][:, 1, :]
        b_explicit = file_explicit["timeseries/b/$iter"][:, 1, :]

        u_implicit = file_implicit["timeseries/u/$iter"][:, 1, :]
        w_implicit = file_implicit["timeseries/w/$iter"][:, 1, :]
        b_implicit = file_implicit["timeseries/b/$iter"][:, 1, :]

        comparison_u[i] = norm(u_explicit .- u_implicit)
        comparison_w[i] = norm(w_explicit .- w_implicit)
        comparison_b[i] = norm(b_explicit .- b_implicit)
    end

    kwargs = (linewidth = 2, foreground_color_legend = nothing, legendfontsize = 12, legend = :right, grid = false,
              xtickfontsize = 12, ytickfontsize=12, xlabel = "time", ylabel = "norm of difference")

     plot(iterations, comparison_u, label = "u"; kwargs...)
    plot!(iterations, comparison_w, label = "w"; kwargs...)
    plot!(iterations, comparison_b, label = "b"; kwargs...)

    savefig(filename * "_comparison_implicit_explicit.png")

    close(file_explicit)
    close(file_implicit)
end

visualize_internal_tide_simulation("internal_tide_explicit")
visualize_internal_tide_simulation("internal_tide_implicit")
plot_implicit_explicit_difference("internal_tide")


text:
"""
    function share_sgf!(d::Dict)
`sgf`: State Government Finance data.

Note: D.C. is not included in the original data set, so assume its SGFs equal Maryland's.
"""
function share_sgf!(d::Dict)
    println("  Calculating regional shares of State Government Finance data")
    d[:sgf] = d[:sgf] / transform_over(d[:sgf], :r)
    verify_over(d[:sgf],:r) !== true && @error("SGF shares don't sum to 1.")
    return dropnan!(d[:sgf])
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using BrainFlow

# enable all possible logs from all three libs
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)
BrainFlow.enable_dev_logger(BrainFlow.DATA_HANDLER)
BrainFlow.enable_dev_logger(BrainFlow.ML_MODULE)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)
sampling_rate = BrainFlow.get_sampling_rate(BrainFlow.SYNTHETIC_BOARD)
nfft = BrainFlow.get_nearest_power_of_two(sampling_rate)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_board_data(board_shim)
BrainFlow.release_session(board_shim)

eeg_channels = BrainFlow.get_eeg_channels(BrainFlow.SYNTHETIC_BOARD)

bands = BrainFlow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
feature_vector = vcat(bands[1], bands[2])

# calc concentration
model_params = BrainFlowModelParams(metric = "concentration", classifier = "KNN")
BrainFlow.prepare(model_params)
print(BrainFlow.predict(feature_vector, model_params))
BrainFlow.release(model_params)

# calc relaxation
model_params = BrainFlowModelParams(metric = "relaxation", classifier = "regression")
BrainFlow.prepare(model_params)
print(BrainFlow.predict(feature_vector, model_params))
BrainFlow.release(model_params)


text:
using Test, SpatioTemporalSystems, DrWatson
# Test:

p = (dt = 0.05,)
sts = STS{:bocf}(N = 100, T = 10, Δt = 1, S = 500, ic = "spiral", p=p)
d = makesim(sts)

U = d["U"]

using PyPlot
close("all")
fig, axs = subplots(1, 3, figsize = (10, 4))

for i in 1:3
    axs[i].imshow(U[i])
end


text:
function DLDA(ALPHA,BETA)
      A=[-0.0410   -0.0410   -0.0420   -0.0400   -0.0430   -0.0440   -0.0430
         -0.0520   -0.0530   -0.0530   -0.0520   -0.0490   -0.0480   -0.0490
         -0.0530   -0.0530   -0.0520   -0.0510   -0.0480   -0.0480   -0.0470
         -0.0560   -0.0530   -0.0510   -0.0520   -0.0490   -0.0470   -0.0450
         -0.0500   -0.0500   -0.0490   -0.0480   -0.0430   -0.0420   -0.0420
         -0.0560   -0.0510   -0.0490   -0.0480   -0.0420   -0.0410   -0.0370
         -0.0820   -0.0660   -0.0430   -0.0420   -0.0420   -0.0200   -0.0030
         -0.0590   -0.0430   -0.0350   -0.0370   -0.0360   -0.0280   -0.0130
         -0.0420   -0.0380   -0.0260   -0.0310   -0.0250   -0.0130   -0.0100
         -0.0380   -0.0270   -0.0160   -0.0260   -0.0210   -0.0140   -0.0030
         -0.0270   -0.0230   -0.0180   -0.0170   -0.0160   -0.0110   -0.0070
         -0.0170   -0.0160   -0.0140   -0.0120   -0.0110   -0.0100   -0.0080]
      #
      S= 0.2 * ALPHA
      K= floor(S)+3
      if(K <= 1)
         K= 2
      end
      if(K >=  12)
         K=  11
      end
      DA= S - K +3
      L = K + fix( 1.1*sign(DA) )
      #
      S= 0.1 * BETA
      M= floor(S)+4
      if(M <= 1)
         M= 2
      end
      if(M >=  7)
         M=  6
      end
      DB= S - M + 4
      N= M + fix( 1.1*sign(DB) )
      #
      K=convert(Int64, K)
      L=convert(Int64, L)
      M=convert(Int64, M)
      N=convert(Int64, N)
      T= A[K,M]
      U= A[K,N]
      V= T + abs(DA) * (A[L,M] - T)
      W= U + abs(DA) * (A[L,N] - U)
      return dldav = V + (W-V)  * abs(DB)
end


text:
using FixedPointNumbers, Colors, ColorVectorSpace, SimpleTraits, ImageAxes, ImageMetadata, OffsetArrays
using Test
import Dates: now
using Unitful: m
import AxisArrays
using AxisArrays: AxisArray, axisnames, (..)

@testset "indexing" begin
    # 1d images
    for A = (rand(3),
             view(rand(6), 1:2:5),
             view(rand(6), [1,2,4]),
             rand(RGB{N0f8}, 3),
             rand(Gray{N0f8}, 3))
        img = ImageMeta(A, prop1=1, prop2=[1,2,3])
        @test eltype(img) == eltype(A)
        @test IndexStyle(img) == IndexStyle(A)
        @test ndims(img) == 1
        @test size(img) == (3,)
        @test data(img) === A
        for i = 1:3
            @test @inferred(img[i]) === A[i]
        end
        for I in eachindex(img)
            @test @inferred(img[I]) === A[I]
        end
        k = 0
        for a in img
            @test a == A[k+=1]
        end
        img[2] = zero(eltype(img))
        @test A[2] == zero(eltype(A))
        img[3] = oneunit(eltype(img))
        @test A[3] == oneunit(eltype(A))
        @test_throws BoundsError img[0]
        @test_throws BoundsError img[4]
        @test img["prop1"] == 1
        @test img["prop2"] == [1,2,3]
        img["prop1"] = -1
        @test img["prop1"] == -1
    end
    for A in (rand(3,5),
              view(rand(4,6), 1:3, 1:5),
              view(rand(4,5), [1,2,4], :),
              reshape(1:15, 3, 5),
              rand(RGB{N0f8}, 3, 5),
              rand(Gray{Float32}, 3, 5))
        img = ImageMeta(A, prop1=1, prop2=[1,2,3]) # TODO: add @inferred (see julia #17719)
        @test eltype(img) == eltype(A)
        @test IndexStyle(img) == IndexStyle(A)
        @test ndims(img) == 2
        @test size(img) == (3,5)
        @test data(img) === A
        for j = 1:5, i = 1:3
            @test @inferred(img[i,j]) === A[i,j]
        end
        for k = 1:15
            @test @inferred(img[k]) === A[k]
        end
        for I in eachindex(img)
            @test @inferred(img[I]) === A[I]
        end
        k = 0
        for a in img
            @test a == A[k+=1]
        end
        if !isa(A, typeof(reshape(1:15, 3, 5)))
            img[2,3] = zero(eltype(img))
            @test A[2,3] == zero(eltype(A))
            img[4] = oneunit(eltype(img))
            @test A[4] == oneunit(eltype(A))
        end
        @test_throws BoundsError img[0,0]
        @test_throws BoundsError img[4,1]
        @test_throws BoundsError img[1,6]
        @test img["prop1"] == 1
        @test img["prop2"] == [1,2,3]
        img["prop1"] = -1
        @test img["prop1"] == -1
        # vector-indexing
        @test isa(@inferred(img[:,:]), ImageMeta) && img[:,:] == img
        @test isa(@inferred(img[:]), ImageMeta) && img[:] == A[:]
        @test isa(@inferred(img[1:2,1:2]), ImageMeta) && img[1:2,1:2] == A[1:2,1:2]
        @test isa(@inferred(img[1:2,:]), ImageMeta) && img[1:2,:] == A[1:2,:]
        @test isa(@inferred(img[:,1:2]), ImageMeta) && img[:,1:2] == A[:,1:2]
        @test isa(@inferred(view(img,:,:)), ImageMeta) && view(img,:,:) == img
        @test isa(@inferred(view(img,:)), ImageMeta) && view(img,:) == A[:]
        @test isa(@inferred(view(img,1:2,1:2)), ImageMeta) && view(img,1:2,1:2) == A[1:2,1:2]
        @test isa(@inferred(view(img,1:2,:)), ImageMeta) && view(img,1:2,:) == A[1:2,:]
        @test isa(@inferred(view(img,:,1:2)), ImageMeta) && view(img,:,1:2) == A[:,1:2]
    end
    # Test bounds-checking removal by @inbounds
    if Base.JLOptions().check_bounds != 1 && Base.JLOptions().can_inline == 1
        set5!(x) = @inbounds x[5] = 1.234
        get5(x) = @inbounds x[5]
        aa = zeros(3)
        sizehint!(aa, 10)  # make sure we don't cause a segfault
        a = ImageMeta(aa)
        @test_throws BoundsError a[5]
        set5!(a)
        val = get5(a)
        @test val == 1.234
        a = zeros(3,5)
    end

    # Issue #10
    A = AxisArray(rand(3,5), :y, :x)
    B = ImageMeta(A, info="blah")
    Broi = B[2:3, 2:3]
    @test isa(Broi, ImageMeta)
    @test axisnames(Broi) == (:y, :x)
    A1, B1 = A[2:7], B[2:7]
    @test isa(B1, ImageMeta) && A1 == B1
    Broi = view(B, 2:3, 2:3)
    @test isa(Broi, ImageMeta)
    @test axisnames(Broi) == (:y, :x)

    # Non-1 indexing
    Ao = OffsetArray(rand(2,3), 0:1, -1:1)
    A = ImageMeta(Ao, info="blah")
    @test axes(A) === axes(Ao)
end

@testset "convert" begin
    A = rand(3,5)
    M = convert(ImageMeta, A)
    @test isa(M, ImageMeta)
    @test convert(ImageMeta, M) === M
    @test convert(ImageMeta{Float64}, M) === M
    @test eltype(convert(ImageMeta{Gray{Float64}}, M)) == Gray{Float64}
    @test eltype(convert(ImageMeta{Gray}, M)) == Gray{Float64}
    @test convert(ImageMeta{Gray}, M) == M
    @test convert(ImageMeta{Gray}, A) == A
end

@testset "reinterpret" begin
    # It's possible that reinterpretc shouldn't be defined for ImageMeta, but...
    A = rand(Float32, 4, 5)
    M = ImageMeta(A, meta=true)
    Mr = reinterpretc(Gray, M)
    @test eltype(Mr) == Gray{Float32}
    @test Mr["meta"] = true
    # Ensure that it gets defined for the formerly un-reinterpretable
    A = zeros(Float32, 4, 5)
    M = ImageMeta(view(A, 1:2:3, 1:4), meta=true)
    Mr = reinterpretc(Gray, M)
    Mr = reinterpretc(Gray{Float32}, M)
    @test eltype(Mr) == Gray{Float32}
    @test Mr["meta"] = true
    Mr[:] .= Gray{Float32}(0.5)
    @test all(A[1:2:3,1:4] .== 0.5)
end

@testset "copy/similar" begin
    img = ImageMeta(rand(3,5); prop1 = 1, prop2 = [1,2,3])
    img2 = copy(img)
    @test img2.data == img.data
    img2[2,2] = -1
    @test img2[2,2] < 0
    @test img[2,2] >= 0
    img2["prop2"][2] = -2
    @test img2["prop2"] == [1,-2,3]
    @test img["prop2"] == [1,2,3]
    img2 = similar(img)
    @test img2["prop1"] == 1
    @test img2["prop2"] == [1,2,3]
    @test img2.data != img.data
    img2["prop3"] = 7
    @test !haskey(img, "prop3")
    img2 = similar(img, RGB{Float16}, (Base.OneTo(5),))
    @test img2["prop1"] == 1
    @test img2["prop2"] == [1,2,3]
    @test eltype(img2) == RGB{Float16}
    @test size(img2) == (5,)
    @test img2.data != img.data
    img2["prop3"] = 7
    @test !haskey(img, "prop3")

    A = AxisArray(rand(3,5), :y, :x)
    B = ImageMeta(A, info="blah")
    C = similar(B)
    @test isa(C, ImageMeta)
    @test isa(C.data, AxisArray)
    @test eltype(C) == Float64
    C = similar(B, RGB{Float16})
    @test isa(C, ImageMeta)
    @test isa(C.data, AxisArray)
    @test eltype(C) == RGB{Float16}
end

@testset "copy/shareproperties/viewim" begin
    img = ImageMeta(rand(3,5); prop1 = 1, prop2 = [1,2,3])
    @test !isempty(properties(img))
    v = view(img, 1:2, 1:2)
    c = img[1:2, 1:2]
    @test v["prop1"] == 1
    @test c["prop1"] == 1
    img2 = copyproperties(img, reshape(1:15, 5, 3))
    @test size(img2) == (5,3)
    img2["prop1"] = -1
    @test img["prop1"] == 1
    img2 = shareproperties(img, reshape(1:15, 5, 3))
    @test size(img2) == (5,3)
    img2["prop1"] = -1
    @test img["prop1"] == -1
    @test v["prop1"] == -1
    @test c["prop1"] == 1
    imgb = ImageMeta(rand(RGB{N0f8}, 2, 2), propa = "hello", propb = [1,2])
    copy!(img, imgb, "propa", "propb")
    @test img["propa"] == "hello"
    @test img["propb"] == [1,2]
    img["propb"][2] = 10
    @test img["propb"] == [1,10]
    @test imgb["propb"] == [1,2]
    delete!(img, "propb")
    @test  haskey(img, "propa")
    @test !haskey(img, "propb")

    img = ImageMeta(AxisArray(rand(3,5,8),
                              Axis{:x}(1:3),
                              Axis{:y}(1:5),
                              Axis{:time}(0.1:0.1:0.8));
                    prop1 = 1, prop2 = [1,2,3])
    v = view(img, Axis{:time}(0.25..0.5))
    c = img[Axis{:time}(0.25..0.5)]
    @test v["prop1"] == 1
    @test c["prop1"] == 1
end

@testset "views" begin
    for A1 in (rand(Gray{N0f8}, 4, 5), rand(RGB{Float32}, 4, 5))
        t1 = now()
        M1 = ImageMeta(A1, date=t1)
        vM1 = channelview(M1)
        @test isa(vM1, ImageMeta)
        @test vM1["date"] == t1
        @test data(vM1) == channelview(A1)
        @test isa(rawview(vM1), ImageMeta)
        if ndims(rawview(vM1)) == 3
            @test rawview(vM1)[1,2,1] === rawview(channelview(A1))[1,2,1]
        else
            @test rawview(vM1)[1,2] === rawview(channelview(A1))[1,2]
        end
    end
    for (A1,C) in ((rand(UInt8, 4, 5), Gray), (rand(UInt8, 3, 4, 5), RGB))
        M1 = ImageMeta(A1)
        vM1 = normedview(M1)
        @test isa(vM1, ImageMeta)
        @test data(vM1) == normedview(A1)
        cvM = colorview(C, vM1)
        @test isa(cvM, ImageMeta)
        @test cvM[1,2] === colorview(C, normedview(A1))[1,2]
    end
    A = AxisArray(rand(RGB{N0f8}, 3, 5), :y, :x)
    t = now()
    M = ImageMeta(A, date=t)
    vM = channelview(M)
    @test isa(vM, ImageMeta)
    @test colordim(vM) == 1
    pvM = permutedims(vM, (2,3,1))
    @test colordim(pvM) == 3
    @test pvM["date"] == t
    sleep(0.1)
    pvM["date"] = now()
    @test M["date"] == t && pvM["date"] != t
    vM = permuteddimsview(M, (2,1))
    @test vM["date"] == t
    @test axisnames(vM) == (:x, :y)
end

@testset "meta-axes" begin
    A = AxisArray(rand(3,5), :y, :x)
    M = ImageMeta(A)
    @test AxisArrays.HasAxes(M) == AxisArrays.HasAxes{true}()
    @test AxisArrays.axes(M) == (Axis{:y}(1:3), Axis{:x}(1:5))
    @test AxisArrays.axes(M, 2) == Axis{:x}(1:5)
    @test axisdim(M, Axis{:y}) == 1
    @test size(M, Axis{:y}) == 3
    @test size(M, Axis{:x}()) == 5
    @test axes(M, Axis{:y}) === Base.OneTo(3)
    @test axes(M, Axis{:x}()) == Base.OneTo(5)
    @test axisdim(M, Axis{:x}) == 2
    @test_throws ErrorException axisdim(M, Axis{:z})
    @test axisnames(M) == (:y, :x)
    @test axisvalues(M) == (1:3, 1:5)
    @test M[Axis{:y}(2:3)] == A[Axis{:y}(2:3)]
    @test view(M, Axis{:y}(2:3)) == A[Axis{:y}(2:3)]
    M[Axis{:y}(2:3), Axis{:x}(1)] .= -5
    @test all(A[2:3,1] .== -5)
end

@testset "show" begin
    for supp in (Set(["prop3"]), "prop3")
        img = ImageMeta(rand(3,5); prop1 = 1, prop2 = [1,2,3], suppress = supp, prop3 = "hide")
        str = string(img)
        @test occursin("ImageMeta with", str)
        @test occursin("prop1", str)
        @test occursin("prop2", str)
        @test occursin("$([1,2,3])", str)
        @test occursin("prop3", str)
        @test !occursin("hide", str)
        @test occursin("<suppressed>", str)
        @test !occursin("suppress:", str)
        io = IOBuffer()
        show(io, MIME("text/plain"), img)
        str2 = String(take!(io))
        @test str == str2
    end
end

@testset "traits" begin
    img = ImageMeta(AxisArray(rand(3,5,8),
                              Axis{:x}(1:3),
                              Axis{:y}(1:5),
                              Axis{:time}(0.1:0.1:0.8)))
    @test @inferred(timedim(img)) == 3
    @test @inferred(pixelspacing(img)) === (1,1)
    @test spacedirections(img) === ((1,0),(0,1))
    img["spacedirections"] = "broken"
    @test spacedirections(img) == "broken"
    @test sdims(img) == 2
    @test @inferred(coords_spatial(img)) == (1,2)
    @test spatialorder(img) == (:x, :y)
    @test nimages(img) == 8
    @test @inferred(size_spatial(img)) == (3,5)
    @test @inferred(indices_spatial(img)) == (Base.OneTo(3), Base.OneTo(5))
    assert_timedim_last(img)
    @test timeaxis(img) == Axis{:time}(0.1:0.1:0.8)
end

@testset "spatialprops" begin
    img = ImageMeta(rand(3,5),
                    spatialproperties=Set(["vector","matrix","tuple"]),
                    vector=[1,2],
                    matrix=[1 3; 2 4],
                    tuple=(1,2))
    for imgp in (img', permutedims(img, (2,1)), permuteddimsview(img, (2,1)))
        @test imgp.data == img.data'
        @test imgp["vector"] == [2,1]
        @test imgp["matrix"] == [4 2; 3 1]
        @test imgp["tuple"]  == (2,1)
        @test img["vector"] == [1,2]
        @test img["matrix"] == [1 3; 2 4]
        @test img["tuple"]  == (1,2)
    end
end

@testset "dimchange" begin
    M = ImageMeta([1,2,3,4])
    Mp = M'
    @test ndims(Mp) == 2 && isa(Mp, ImageMeta)

    M = ImageMeta([1,2,3,4],
                  spatialproperties=["vector"],
                  vector=[1])
    @test_throws ErrorException M'
end

@testset "Inference" begin
    # Heavily-nested types tax Julia's inference capabilities
    a = rand(N2f14, 3, 5, 2)
    a1 = a[1,1,1]
    aa = AxisArray(a, Axis{:y}(0:2), Axis{:x}(0:4), Axis{:z}(0:1))
    cv = colorview(RGB, aa, zeroarray, aa)
    @test @inferred(cv[1,1,1]) == RGB(a1, zero(a1), a1)
    au = AxisArray(a, Axis{:y}(0m:1m:2m), Axis{:x}(0m:1m:4m), Axis{:z}(0m:1m:1m))
    cv = colorview(RGB, au, zeroarray, au)
    @test @inferred(cv[1,1,1]) == RGB(a1, zero(a1), a1)
    am = ImageMeta(au)
    cv = colorview(RGB, am, zeroarray, am)
    @inferred cv[1,1,1]
end

@testset "AxisArray_CartesianIndex" begin
    M = reshape([1,2,3,4], 2,2)
    M = AxisArray(M)
    img = ImageMeta(M)
    @test 1 == img[1, CartesianIndex(1)] #Int and cartesianindex
    @test maximum(img,dims=2) == reshape([3,4],2,1)
    @test mean(img,dims=1) == reshape([1.5,3.5], 1,2)
end

@testset "AbstractDicts" begin
    img = ImageMeta(AxisArray(rand(3,5,8),
                              Axis{:x}(1:3),
                              Axis{:y}(1:5),
                              Axis{:time}(0.1:0.1:0.8)),
                    IdDict{String,Any}("pixelspacing" => (1,1)))
    @test @inferred(pixelspacing(img)) === (1,1)

    rand_img = rand(10,10)
    dict = Dict("attr1" => 3, "attr2" => 4)
    @test ImageMeta(rand_img, dict)[1, 1] == rand_img[1, 1]
end

nothing


text:
handle("workspace") do mod
  workspace(mod)
end

function workspace(mod)
  mod = getmodule(mod)
  ns = Symbol.(CodeTools.filtervalid(names(mod; all = true)))
  filter!(ns) do n
    !Base.isdeprecated(mod, n) &&
    string(n) != last(split(string(mod), '.'))
  end
  contexts = [d(:context => string(mod), :items => map(n -> wsitem(mod, n), ns))]
  isdebugging() && prepend!(contexts, JunoDebugger.contexts())
  contexts
end

wsitem(mod, name) = begin
  val = getfield′(mod, name)
  wsitem(mod, name, val)
end
wsitem(mod, name, @nospecialize(val)) = begin
  d(:name       => name,
    :value      => render′(Inline(), val),
    :nativetype => nativetype(mod, name, val),
    :type       => wstype(mod, name, val),
    :icon       => wsicon(mod, name, val))
end

nativetype(mod, name, @nospecialize(val)) = DocSeeker.determinetype(mod, name)
nativetype(mod, name, ::Undefined) = "Undefined"

#=
@NOTE: `wstype` and `wsicon` are also used for completions / docs
=#

wstype(mod, name, @nospecialize(val)) = isconst(mod, name) ? "constant" : "variable"
wstype(mod, name, f::Function) = ismacro(f) ? "snippet" : "function"
wstype(mod, name, ::Type) = "type"
wstype(mod, name, ::Module) = "module"
wstype(mod, name, ::Expr) = "mixin"
wstype(mod, name, ::Symbol) = "tag"
wstype(mod, name, ::Exception) = "mixin"
wstype(mod, name, ::Undefined) = "ignored"

wsicon(mod, name, @nospecialize(val)) = isconst(mod, name) ? "c" : "v"
wsicon(mod, name, f::Function) = ismacro(f) ? "icon-mention" : "λ"
wsicon(mod, name, ::Type) = "T"
wsicon(mod, name, ::Module) = "icon-package"
wsicon(mod, name, ::Number) = "n"
wsicon(mod, name, ::AbstractVector) = "icon-list-ordered"
wsicon(mod, name, ::AbstractArray) = "icon-file-binary"
wsicon(mod, name, ::AbstractDict) = "icon-list-unordered"
wsicon(mod, name, ::AbstractString) = "icon-quote"
wsicon(mod, name, ::Regex) = "icon-quote"
wsicon(mod, name, ::Expr) = "icon-code"
wsicon(mod, name, ::Symbol) = "icon-code"
wsicon(mod, name, ::Exception) = "icon-bug"
wsicon(mod, name, ::Undefined) = "icon-circle-slash"


text:
module GradientBoost

# Load sources files
include("util.jl")
include("loss.jl")
include("gb.jl")
include("gb_bl.jl")
include("gb_dt.jl")
include("ml.jl")

end # module


text:
#
# from: Equivalence scales: rationales, uses and assumptions
# Jenny Chanfreau and Tania Burchardt Table 1 p.5
# FIXME we need complete table 1 omplementation
# missing mclem-ahc, whatevere pse is
#

@enum Equivalence_Scale_Type oxford modified_oecd eq_square_root mcclements eq_per_capita
@enum Equiv_Person_Type eq_head eq_spouse_of_head eq_other_adult eq_dependent_child eq_hbai_ahc

struct EQ_Person
      age::Integer
      ptype::Equiv_Person_Type
end

"""
 from: Equivalence scales: rationales, uses and assumptions
 Jenny Chanfreau and Tania Burchardt Table 1 p.5
 https://www2.gov.scot/Topics/Statistics/Browse/Social-Welfare/IncomePoverty/equivalence-scales-paper

 FIXME we need complete table 1 implementation
 missing mclem-ahc, whatevere pse is
"""
function get_equivalence_scales( people::Vector{EQ_Person} )::Dict{Equivalence_Scale_Type,Real}

      eq = 0.0
      scales = Dict{Equivalence_Scale_Type,Real}()
      np = size(people)[1]
      # choose someone to be the head of the unit
      pos_of_head = 1
      positions_of_head = findall((p -> (p.ptype == eq_head) & (p.age >= 16)), people)
      num_heads = size(positions_of_head)[1]
      if num_heads > 0
            pos_of_head = positions_of_head[1]
      end
      scales[eq_square_root] = sqrt(np)
      scales[eq_per_capita] = np
      ox = 0.0
      oecd = 0.0
      mclem = 0.0
      n_extra_adults = 0
      hbai = 0.0
      for pno in 1:np
            if pno == pos_of_head
                  ox += 1
                  oecd += 1
                  mclem += 1
                  hbai += 1
            else
                  if people[pno].age <= 14
                        ox += 0.5
                        oecd += 0.3
                  else
                        ox += 0.7
                        oecd += 0.5
                  end
                  if people[pno].ptype == eq_dependent_child
                        if people[pno].age in 0:1
                              mclem += 0.148
                        elseif people[pno].age in 2:4
                              mclem += 0.295
                        elseif people[pno].age in 5:7
                              mclem += 0.344
                        elseif people[pno].age in 8:10
                              mclem += 0.377
                        elseif people[pno].age in 11:12
                              mclem += 0.410
                        elseif people[pno].age in 13:15
                              mclem += 0.443
                        else
                              mclem += 0.590
                        end
                        if people[pno].age <= 14
                              hbai += 0.34
                        else
                              hbai += 0.72
                        end
                  elseif people[pno].ptype == eq_spouse_of_head
                        mclem += 0.640
                        hbai += 0.72
                  else # other adult
                        n_extra_adults += 1
                        if n_extra_adults == 1
                              mclem += 0.75
                              hbai += 0.72
                        else
                              mclem += 0.59
                              hbai += 0.72
                        end
                  end
            end # person not the head
      end # person loop
      @assert mclem > 0
      @assert oecd > 0
      @assert ox > 0
      scales[oxford] = ox
      scales[modified_oecd] = oecd
      scales[mcclements] = mclem
      scales[eq_hbai_ahc] = hbai
      scales
end # get_equivalence_scales


text:
@testset "Bitmap Complex" begin
    # empty complex
    cplx = BitmapComplex()
    @test size(cplx) == (0,)
    @test size(cplx,0) == 0
    @test length(cplx) == 0
    @test dim(cplx) == 0
    @test eltype(cplx) == Float64
    cplx = BitmapComplex(Int)
    @test eltype(cplx) == Int

    # 2D complex
    img = [
        7 6  5 ;
        8 20 4;
        1 4  6 ;
    ]
    cplx = BitmapComplex(img)
    @test size(cplx) == (16, 24, 9)
    @test size(cplx,0) == 16
    @test length(cplx) == 49
    @test dim(cplx) == 2
    @test eltype(cplx) == Int

    @test ComputationalHomology.topcubes(cplx) == [9, 11, 13, 23, 25, 27, 37, 39, 41]

    @test ComputationalHomology.cellindex(cplx, 1, 1) == 9

    cdc  = [ComputationalHomology.celldim(cplx, i) for i in 1:length(cplx)]
    cdt = vec([0 1 0 1 0 1 0 1 2 1 2 1 2 1 0 1 0 1 0 1 0 1 2 1 2 1 2 1 0 1 0 1 0 1 0 1 2 1 2 1 2 1 0 1 0 1 0 1 0])
    @test cdc == cdt

    @test faces(cplx, 9) == [2, 16, 10, 8]
    @test cofaces(cplx, 3) == [10, 2, 4]

    flt = filtration(cplx)
    @test eltype(flt) == BitmapComplex{Int}
    @test valtype(flt) == Float32
    @test minimum(flt) == 1
    @test maximum(flt) == 20

    ∂ = boundary(flt)
    @test length(∂) == length(complex(flt))
    @test count(!iszero, sparse(∂)) == 84

    @testset "PH Intervals " for (d, itrs) in diagram(flt)
        titrs = d == 0 ? diagram(1.0=>Inf) : diagram(8.0=>20.0)
        for itr in itrs
            @test itr ∈ titrs
        end
    end

    # Perseus format
    str = "2\n3\n3\n1\n4\n6\n8\n20\n4\n7\n6\n5"
    io = IOBuffer(str)
    cplx = read(io, BitmapComplex, Int)
    @test size(cplx) == (16, 24, 9)
    @test dim(cplx) == 2
    @test eltype(cplx) == Int

    # img = reshape([0,0,2,2,2,3,0,0,0,2,0,0,1,3,0,0,0,0,2,1,2,0,0,0,0,0,0,4,1,0,0,0,0,0,1,1,1,1,0,0,0,0,2,0,0,2,1,0,0,0,2,0,0,0,2,0,0,0,1,2,2,3,0,0],8,8)
    # flt2 = filtration(BitmapComplex(maximum(img) .- img))
    # diagram(flt2)
    # z = maximum(img) .- img
    # join(mapslices(v->join(v,"\t"), map(v->v == 4 ? "" : string(v), z), dims=2), "\n") |> println

end


text:
#=
P.4 Complex step Method

for analytic t i = sqrt(-i)

    f'(x) ≈ Im f(x + i*h)/h

Error O(h^2)
h not restricted by rounding errors
can take h = 10^100

- Approximate automatic Differntiation
- Alg for f must not emply complex arithmetic

Paper:
-Lyness & Moler (1967)
- Squire & Trapp (1998)

f(x) = atan(x) / (1+ exp(-x^2)) at x=2:

Prof Higham: "derivative is correct to 15 figures"

Accurate derivative
0.274623728154858 gets all digits correct

Complex step with h = le-100. #gets all digits correct
0.274623728154858

Forward difference with h = 8*u.
0.250000000000000

Forward difference with h = sqrt(8*u).
0.274623729288578

#takeaway: with optimal choice, you get only Half of the digits correct

"The shortest and best way between two
truths of the real domain often passes
through the imaginary one."
Jacques Hadamard (1945)
log(z) is the Alogorithm for which Im log z E (-1r, 1r]

log(z1 + log(z2)    (9.1)

log z1 * z2 = log(abs(z1*z2)) + i*arg(z1*z2)
            = Log(abs(z1) + Log(z2) + i*arg(z1)+ i*arg(z2)
            = log(z1) + log(z2)

Difference between Log vs log:

book: an introduction to complex analysis R.P Agarwal 2011
turning the page, 9.1 holds if you
take the appropriate
"Branch of the log (in each of those occurences of the log )
(not an actual identity to the principal Log)
So, the formula is not , in general, true
Q.1how can we understand this?
Q2. how can we compute this in practice ?
Goes on say it must take appropriate Brach for
    each occurence of "log"

Book: NIST Handbook (Olver et al. 2010)
a very widely-used mathematical functions
arcsin, Arccos are "general values" of inverse sine, inverse cosine
9
4.24(iii) Addition Formulas

Arcsin(u)  ± Arcsin(v)
= - Arcsin(u*(l - v^2)^1/2 ± v(l - 1t2 )^1/2) 4.24.13
Arccos(u) = Arcos(v)
= Arccos(u*wi ∓ ((1-u^2))*(1 - u^2)*(1-v^2))^1/2) 4.24.14

Arcsin, Arccos are "general values" of inverse sine, inverse
cosine

Each occuruce of the "Multi-valued function" could be from Different Branch
Q. if I have a particular u & v
for my principal branch, I want to know : what can i say
 the Author's Understanding (topic)v.s. Reader Understanding(topic )

 #---P.12. unwinding Number

 U(z) = z- log(exp(x)^z)/2*pi*i

 Note: log(exp(x)^z)+ 2*pi*i*U(z)
 papers:
 -Corless,Hare & Jeffery (1996)
 -Apostol (1974):special case

 #---P.13.
 unwinding number is an integer (an alternative formula)

 U(z) = z-log(exp(x)^z)/2*pi*i = ceil(Im(z)-pi/2*pi)

 when is log exp^z = z  #NOT A USUAL ^
 U(z) = 0 iff Im(z) ∈ (i*pi,pi} #semi-open interval
me: reimann surface of logarithm forms a helix variable color(whose value vary)

comment: turns out
1.unwinding number is an integer
2.reimann Surface: counts which branch of Reimann Surface (you are on)

 when log(e^z)=z  [me: no errors!] #NOT USUAL ^
 it will be as long as imaginary part (Im(z))
 is in a VErtical strip (horizontal strip) - i digress
 between [-pi,pi]

 #--- Roundtrip Relations
another identity:
cos(acos(z)) = z [minimal error ]

WARNING: AVOID AMBIGUITY WITH ^
SEE:https://github.com/JuliaLang/julia/blob/master/base/complex.jl
there are 2 imaginary number specifications :
im , imag
 =#

z1 = imag #generic function with 19 methods
z2 = im::Complex
z1 != z2
z = z2
cos(acos(z)) == z #false
cos(acos(z)) ≈ z  #true! [already an identity]
# acos(cos(z)) = z # Syntax Error cos(z) is not a valid function

function log10(z::Complex)
    a = log(z)
    a / log(oftype(real(a), 10))
end
function log2(z::Complex)
    a = log(z)
    a / log(oftype(real(a), 2))
end

function exp(z::Complex)
    zr, zi = reim(z)
    if isnan(zr)
        Complex(zr, zi == 0 ? zi : zr)
    elseif !isfinite(zi)
        if zr == Inf
            Complex(-zr, oftype(zr, NaN))
        elseif zr == -Inf
            Complex(-zero(zr), copysign(zero(zi), zi))
        else
            Complex(oftype(zr, NaN), oftype(zi, NaN))
        end
    else
        er = exp(zr)
        if iszero(zi)
            Complex(er, zi)
        else
            s, c = sincos(zi)
            Complex(er * c, er * s)
        end
    end
end

function expm1(z::Complex{T}) where {T<:Real}
    Tf = float(T)
    zr, zi = reim(z)
    if isnan(zr)
        Complex(zr, zi == 0 ? zi : zr)
    elseif !isfinite(zi)
        if zr == Inf
            Complex(-zr, oftype(zr, NaN))
        elseif zr == -Inf
            Complex(-one(zr), copysign(zero(zi), zi))
        else
            Complex(oftype(zr, NaN), oftype(zi, NaN))
        end
    else
        erm1 = expm1(zr)
        if zi == 0
            Complex(erm1, zi)
        else
            er = erm1 + one(erm1)
            if isfinite(er)
                wr = erm1 - 2 * er * (sin(convert(Tf, 0.5) * zi))^2
                return Complex(wr, er * sin(zi))
            else
                s, c = sincos(zi)
                return Complex(er * c, er * s)
            end
        end
    end
end

function log1p(z::Complex{T}) where {T}
    zr, zi = reim(z)
    if isfinite(zr)
        isinf(zi) && return log(z)
        # This is based on a well-known trick for log1p of real z,
        # allegedly due to Kahan, only modified to handle real(u) <= 0
        # differently to avoid inaccuracy near z==-2 and for correct branch cut
        u = one(float(T)) + z
        u == 1 ? convert(typeof(u), z) : real(u) <= 0 ? log(u) : log(u) * z / (u - 1)
    elseif isnan(zr)
        Complex(zr, zr)
    elseif isfinite(zi)
        Complex(T(Inf), copysign(zr > 0 ? zero(T) : convert(T, pi), zi))
    else
        Complex(T(Inf), T(NaN))
    end
end

function exp2(z::Complex{T}) where {T}
    z = float(z)
    er = exp2(real(z))
    theta = imag(z) * log(convert(float(T), 2))
    s, c = sincos(theta)
    Complex(er * c, er * s)
end

function exp10(z::Complex{T}) where {T}
    z = float(z)
    er = exp10(real(z))
    theta = imag(z) * log(convert(float(T), 10))
    s, c = sincos(theta)
    Complex(er * c, er * s)
end

# _cpow helper function to avoid method ambiguity with ^(::Complex,::Real)
function _cpow(z::Union{T,Complex{T}}, p::Union{T,Complex{T}}) where {T}
    z = float(z)
    p = float(p)
    Tf = float(T)
    if isreal(p)
        pᵣ = real(p)
        if isinteger(pᵣ) && abs(pᵣ) < typemax(Int32)
            # |p| < typemax(Int32) serves two purposes: it prevents overflow
            # when converting p to Int, and it also turns out to be roughly
            # the crossover point for exp(p*log(z)) or similar to be faster.
            if iszero(pᵣ) # fix signs of imaginary part for z^0
                zer = flipsign(copysign(zero(Tf), pᵣ), imag(z))
                return Complex(one(Tf), zer)
            end
            ip = convert(Int, pᵣ)
            if isreal(z)
                zᵣ = real(z)
                if ip < 0
                    iszero(z) && return Complex(Tf(NaN), Tf(NaN))
                    re = power_by_squaring(inv(zᵣ), -ip)
                    im = -imag(z)
                else
                    re = power_by_squaring(zᵣ, ip)
                    im = imag(z)
                end
                # slightly tricky to get the correct sign of zero imag. part
                return Complex(re, ifelse(iseven(ip) & signbit(zᵣ), -im, im))
            else
                return ip < 0 ? power_by_squaring(inv(z), -ip) : power_by_squaring(z, ip)
            end
        elseif isreal(z)
            # (note: if both z and p are complex with ±0.0 imaginary parts,
            #  the sign of the ±0.0 imaginary part of the result is ambiguous)
            if iszero(real(z))
                return pᵣ > 0 ? complex(z) : Complex(Tf(NaN), Tf(NaN)) # 0 or NaN+NaN*im
            elseif real(z) > 0
                return Complex(real(z)^pᵣ, z isa Real ? ifelse(real(z) < 1, -imag(p), imag(p)) : flipsign(imag(z), pᵣ))
            else
                zᵣ = real(z)
                rᵖ = (-zᵣ)^pᵣ
                if isfinite(pᵣ)
                    # figuring out the sign of 0.0 when p is a complex number
                    # with zero imaginary part and integer/2 real part could be
                    # improved here, but it's not clear if it's worth it…
                    return rᵖ * complex(cospi(pᵣ), flipsign(sinpi(pᵣ), imag(z)))
                else
                    iszero(rᵖ) && return zero(Complex{Tf}) # no way to get correct signs of 0.0
                    return Complex(Tf(NaN), Tf(NaN)) # non-finite phase angle or NaN input
                end
            end
        else
            rᵖ = abs(z)^pᵣ
            ϕ = pᵣ * angle(z)
        end
    elseif isreal(z)
        iszero(z) && return real(p) > 0 ? complex(z) : Complex(Tf(NaN), Tf(NaN)) # 0 or NaN+NaN*im
        zᵣ = real(z)
        pᵣ, pᵢ = reim(p)
        if zᵣ > 0
            rᵖ = zᵣ^pᵣ
            ϕ = pᵢ * log(zᵣ)
        else
            r = -zᵣ
            θ = copysign(Tf(π), imag(z))
            rᵖ = r^pᵣ * exp(-pᵢ * θ)
            ϕ = pᵣ * θ + pᵢ * log(r)
        end
    else
        pᵣ, pᵢ = reim(p)
        r = abs(z)
        θ = angle(z)
        rᵖ = r^pᵣ * exp(-pᵢ * θ)
        ϕ = pᵣ * θ + pᵢ * log(r)
    end

    if isfinite(ϕ)
        return rᵖ * cis(ϕ)
    else
        iszero(rᵖ) && return zero(Complex{Tf}) # no way to get correct signs of 0.0
        return Complex(Tf(NaN), Tf(NaN)) # non-finite phase angle or NaN input
    end
end

#=Roundtrip Relations
Theorem

if Re z != k*pi ; k \in Z then

    acos(cos(z)) #not existant per se
    = z - 2*pi*U(i*z)*sign(z-2*pi*U(i*z)))
Corollary
acos(cos(z)) = z iff Re z z∈(0,pi)

=#
#=Accurate Difference

how to compute Accurate Value of
log \lambda , λ1,λ2 ∈ C
=#
# λ1, λ2 = #rand(im,1)# no rand for complex
#unwinding number
U(z) = (z - log(exp(z))) / 2 * pi * im # takeaway:the use of generic functon promotes lazy processing - easy on memory, fast in response

λ1 = 0.75;
λ2 = 0.25;
lhs = log(λ2) - log(λ1)

z = (λ2 - λ1) / (λ2 + λ1)

#TODO: check / Recheck is required
#z = log(exp(z)) + 2 * pi * im * U(z) #no method matching exp(::Float64)
#=
log(λ2) - log(λ1)
someway, i can rewrite the expression , for more accurate Evaluation
by introducing the z:
#TODO: check error1 = abs(lhs - rhs1)
=#

rhs1 = log(λ2 / λ1) + 2 * pi * im * U(log(λ2) - log(λ1))# my bad: U(z) * (log(λ2) - log(λ1) # now, it's corrrect #checked

rhs1_5 = log((1 + z) / (1 - z)) + 2 * pi * im * U(log(λ2) - log(λ1)) #this value is also sound
#result clustering rhs1, rhs1_5
error_initial = rhs1 - rhs1_5 # error is the same  (same values ) # logical error: comparing 2 RHSs # TODO: compare each RHS with its respective LHS

error1 = abs(rhs1) - abs(lhs) # 2.220446049250313e-16 < \epsilon≈ 0

#= "The Difference of Logs as a log Ratio
problem: need a correction term
solution: that's where the unwinding number steps in  "
=log((1+z)/(1-z) + 2*π*i*u(log(λ2)-log(λ1))

=assuming we have a good atanh function
=#
rhs2 = atanh(z) + 2 * pi * im * U(log(λ2) - log(λ1)) #used in state of art matrix log codes

#=you can be sure that won't give you an accurate estimation
(even if lambdas are close )
used in state of the art matrix log codes

emphasize: z still has λ2 - λ1 (in it)
(but they're subtracting original numbers, that was given)

=#

error = abs(z - rhs2)
#--- the Significance (& the Plateau)
#= one value isn't & won't be enough for correctly identifying the Significance of values
    Thus, it is in my pure interest to develop a 'Gray Area' around 0.05
    My prine concern is to build a more valuable plateau around the Significant value (rather than relying on an arbitrary line in sand)
    This would render results as one of three: Significant, insignificant, indeterminant
    The hardest thing is failing to admit our indeterminance
    Our determenance is defined by the thickness of the indeterminance, the plateau must always be smaller than the Significance itself
    Let's say we have a line of Significance at 0.05 , then there must be a free-zone revolving around it, that is equal to each other- as above, so below
    0.05 ± 0.01 : 0.01 thus, absolute thickness would make up 0.02

    =#
isSignificatant = error < 0.05 == true || error >= 0.05 == false

max_stars_repo_path: Tutorial/Nick Higham/log.jl

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_path: Tutorial/Nick Higham/log.jl

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_path: Tutorial/Nick Higham/log.jl

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using HTTP, Random, JSON

const EXPIRY = 30
const CHALLENGE_EXPIRY = 10
const PORT = 8888

const ID_TOKEN = Random.randstring(100)
const TOKEN = Ref(Dict())

challenge_response_map = Dict()
challenge_timeout = Dict()
response_challenge_map = Dict()
function challenge_handler(req)
    println("challenge_handler")
    challenge = String(req.body)
    response = Random.randstring(64)
    challenge_timeout[challenge] = time() + CHALLENGE_EXPIRY

    response_challenge_map[response] = challenge

    return HTTP.Response(200, response)
end

refresh_token_map = Dict()
function response_handler(req)
    println("response_handler")

    response = HTTP.URIs.URI(req.target).query
    @show response
    if haskey(response_challenge_map, response)
        challenge = response_challenge_map[response]
        @show challenge
        challenge_response_map[challenge] = response
    else
        return HTTP.Response(400)
    end

    refresh_token = Random.randstring(10)
    TOKEN[] = Dict(
        "user_name" => "firstname lastname",
        "user_email" => "user@email.com",
        "id_token" => ID_TOKEN,
        "refresh_token" => refresh_token,
        "refresh_url" => "http://localhost:$(PORT)/auth/renew/token.toml/v2/",
        "expires_in" => EXPIRY,
        "expires_at" => round(Int, time()) + EXPIRY
    )

    @show response

    return HTTP.Response(200)
end

function claimtoken_handler(req)
    println("claimtoken_handler")
    payload = JSON.parse(String(req.body))
    @show payload
    @show challenge_response_map
    if haskey(challenge_response_map, payload["challenge"]) &&
       challenge_response_map[payload["challenge"]] == payload["response"]

        delete!(challenge_response_map, payload["challenge"])
        delete!(response_challenge_map, payload["response"])
        @show JSON.json(TOKEN[])
        return HTTP.Response(200, JSON.json(Dict(
            "token" => TOKEN[]
        )))
    else
        expires_in = round(Int, challenge_timeout[payload["challenge"]] - time())
        @show expires_in
        return HTTP.Response(200, JSON.json(Dict(
            "expiry" => expires_in
        )))
    end

end

function renew_handler(req)
    println("renew_handler")
    auth = HTTP.header(req, "Authorization")
    @show auth
    refresh_token = match(r"Bearer (.+)", auth)[1]

    @assert refresh_token == TOKEN[]["refresh_token"]

    TOKEN[]["refresh_token"] = Random.randstring(10)
    TOKEN[]["expires_at"] = ceil(Int, time() + EXPIRY)

    return HTTP.Response(200, JSON.json(Dict(
        "token" => TOKEN[]
    )))
end

function check_validity(req)
    payload = JSON.parse(String(req.body))

    return HTTP.Response(200, payload == TOKEN[])
end

router = HTTP.Router()
HTTP.@register(router, "POST", "/auth/challenge", challenge_handler)
HTTP.@register(router, "GET", "/auth/response", response_handler)
HTTP.@register(router, "POST", "/auth/claimtoken", claimtoken_handler)
HTTP.@register(router, "GET", "/auth/renew/token.toml/v2", renew_handler)
HTTP.@register(router, "POST", "/auth/isvalid", check_validity)

function run()
    println("starting server")
    HTTP.serve(router, "127.0.0.1", PORT)
    readline()
end

run()


text:
module MIDIInterface # PushInterface

# https://github.com/Ableton/push-interface/blob/master/doc/AbletonPush2MIDIDisplayInterface.asc#midi-interface

end # module PushInterface.MIDIInterface


text:
#
# Superblock
#

const SUPERBLOCK_SIGNATURE = htol(0x0a1a0a0d46444889) # UInt8[0o211, 'H', 'D', 'F', '\r', '\n', 0o032, '\n']

# https://www.hdfgroup.org/HDF5/doc/H5.format.html#FileMetaData
# Superblock (Version 2)
# struct Superblock
#     file_consistency_flags::UInt8
#     base_address::Int64
#     superblock_extension_address::RelOffset
#     end_of_file_address::Int64
#     root_group_object_header_address::RelOffset
# end

struct Superblock
    version::UInt8
    base_address::UInt64
    size_of_offsets::UInt8
    size_of_lengths::UInt8
    end_of_file_address::UInt64
    data
end

jlsizeof(::Union{Type{Superblock},Superblock}) =
    12+jlsizeof(RelOffset)*4+4

function jlread(io::IO, ::Type{Superblock})
    cio = begin_checksum_read(io)

    # Signature
    signature = jlread(cio, UInt64)
    signature == SUPERBLOCK_SIGNATURE || throw(InvalidDataException())

    # Version
    version = jlread(cio, UInt8)
    if version == 0
        version_free_space_storage = jlread(cio, UInt8) # has to be zero
        version_root_group_symbol_table_enty = jlread(cio, UInt8) # has to be zero
        jlread(cio, UInt8)
        version_share_header_msg_format = jlread(cio, UInt8) # has to be zero
        size_of_offsets = jlread(cio, UInt8)
        OffsetType = uintofsize(size_of_offsets)
        size_of_lengths = jlread(cio, UInt8)
        jlread(cio, UInt8)
        group_leaf_node_k = jlread(cio, UInt16) # must be greater than zero
        group_internal_node_k = jlread(cio, UInt16) # must be greater than zero
        # Unused File consistency flags
        jlread(cio, UInt32)
        #indexed_storage_internal_node_k = jlread(cio, UInt16) # must be greater than zero
        #jlread(cio, UInt16)
        base_address = jlread(cio, OffsetType) # base adress for offsets within file (also absolute address of superblock)
        adress_free_space_info = jlread(cio, OffsetType)  # Undefined Adress
        end_of_file_address = jlread(cio, OffsetType) # absolute adress of first byte past end of data
        driver_info_block_adress = jlread(cio, OffsetType) # undefined of relative adress of driver info block
        #root_group_symbol_table_entry = jlread(cio, UInt32) # symbol table entry of root group

        link_name_offset = jlread(cio, OffsetType)
        obj_header_adress = jlread(cio, OffsetType)
        cachetype = jlread(cio, UInt32)
        reserved = jlread(cio, UInt32)
        scratchspace = jlread(cio, UInt128)

        # Discard Checksum
        end_checksum(cio)

        Superblock(0, base_address, size_of_offsets, size_of_lengths, end_of_file_address,
                    (;#version_free_space_storage,
                    #version_root_group_symbol_table_enty,
                    #version_share_header_msg_format,
                    group_leaf_node_k,
                    group_internal_node_k,
                    #indexed_storage_internal_node_k,
                    #adress_free_space_info,
                    driver_info_block_adress,
                    #root_group_symbol_table_entry
                    link_name_offset,
                    obj_header_adress,
                    cachetype,
                       ))
    elseif version == 2 || version == 3

        # Size of offsets and size of lengths
        size_of_offsets = jlread(cio, UInt8)
        size_of_lengths = jlread(cio, UInt8)
        (size_of_offsets == 8 && size_of_lengths == 8) || throw(UnsupportedFeatureException())

        # File consistency flags
        file_consistency_flags = jlread(cio, UInt8)

        # Addresses
        base_address = jlread(cio, Int64)
        superblock_extension_address = jlread(cio, RelOffset)
        end_of_file_address = jlread(cio, Int64)
        root_group_object_header_address = jlread(cio, RelOffset)

        # Checksum
        cs = end_checksum(cio)
        jlread(io, UInt32) == cs || throw(InvalidDataException())

        Superblock(version, base_address, size_of_offsets, size_of_lengths, end_of_file_address,
                    (;
                    superblock_extension_address,
                    root_group_object_header_address))
    else
        throw(UnsupportedVersionException("superblock version $version is not supported."))
    end
end

function jlwrite(io::IO, s::Superblock)
    cio = begin_checksum_write(io, 8+4+4*jlsizeof(RelOffset))
    jlwrite(cio, SUPERBLOCK_SIGNATURE::UInt64)    # Signature
    jlwrite(cio, UInt8(2))                        # Version
    jlwrite(cio, UInt8(8))                        # Size of offsets
    jlwrite(cio, UInt8(8))                        # Size of lengths
    jlwrite(cio, s.file_consistency_flags::UInt8)
    jlwrite(cio, s.base_address::Int64)
    jlwrite(cio, s.superblock_extension_address::RelOffset)
    jlwrite(cio, s.end_of_file_address::Int64)
    jlwrite(cio, s.root_group_object_header_address::RelOffset)
    jlwrite(io, end_checksum(cio))
end

function find_superblock(f)#::JLDFile)
    # Search at 0, 512, 1024, 2048 ...
    for offset in (0, 512, 1024, 2048, 4096)
        seek(f.io, offset)
        # Signature
        signature = jlread(f.io, UInt64)
        if signature == SUPERBLOCK_SIGNATURE
            @info "Found Superblock at $(offset)"
            seek(f.io, offset)
            return jlread(f.io, Superblock)
        end
    end
    throw(InvalidDataException("Did not find a Superblock"))
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libiio_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libiio")
JLLWrappers.@generate_main_file("libiio", UUID("8da1a5bb-246c-59c5-9ba1-a6804f8e92f0"))
end  # module libiio_jll


text:
@with_kw struct rtgM3{A,B,C} <: RTGmodel
    model::A = rtgM3_model()
    protein_lookup = get_protein_lookup(model)
    p::C = init_p(model; idx_hill_coefs=[])
    u::B = init_u(model, protein_lookup)
end

"""
This model is modified version of rtgM2.
Major changes:
1. Change input to be michaelis menten
2. Change the translocation to be Michaelis menten: (which is not ture). See [supplemental material](https://www.pnas.org/content/pnas/suppl/2006/06/30/0604085103.DC1/04085SuppText.pdf)
3. Rtg1p and Rtg3p enter nucleus separately

References
-----------

1. Hao, N., Budnik, B. A., Gunawardena, J., & O'Shea, E. K. (2013). Tunable signal processing through modular control of transcription factor translocation. Science, 339(6118), 460-464. [[Article](https://science.sciencemag.org/content/339/6118/460), [SOM](https://www.pnas.org/content/pnas/suppl/2006/06/30/0604085103.DC1/04085SuppText.pdf)]

2. Mitochondrial Signaling: The Retrograde Response. 2004. (Reaction of Bmh-Mks on Rtg13 dimer) [[Article](https://www.cell.com/molecular-cell/pdf/S1097-2765(04)00179-0.pdf)]
"""
function rtgM3_model()
    rtgM3= @reaction_network  begin
        #=
        Input
        =#
        (1), s→s

        #=
        Modulation Layer
        =#
        (mm(s, ksV, ksD), k2I), Rtg2_ina_c ↔ Rtg2_act_c #[1]
        (k2M, kn2M), Rtg2_act_c + Mks ↔ Rtg2Mks_c
        (kBM, knBM), Bmh + Mks ↔ BmhMks

        #=
        Activation Model:
        =#
        # NLSact: Nuclear Localization Signal (NLS) activation
        (k13I + mm(BmhMks, k13IV, k13ID)), Rtg13_a_c → Rtg13_i_c
        (k3A_c, k3I_c), Rtg3_i_c ↔ Rtg3_a_c
        (k3I_n), Rtg3_a_n → Rtg3_i_n

        # Rtg1Binding: Rtg13 formation
        (k13_c, kn13_c), Rtg1_c + Rtg3_a_c ↔ Rtg13_a_c
        (k13_c, kn13_c), Rtg1_c + Rtg3_i_c ↔ Rtg13_i_c

        (k13_n, kn13_n), Rtg1_n + Rtg3_a_n ↔ Rtg13_a_n #[3]
        (k13_n, kn13_n), Rtg1_n + Rtg3_i_n ↔ Rtg13_i_n #[3]

        # Translocation
        (k1in , k1out), Rtg1_c ↔ Rtg1_n #[2]
        (k3inA, k3outA), Rtg3_a_c ↔ Rtg3_a_n #[2]
        (k3inI, k3outI), Rtg3_i_c ↔ Rtg3_i_n

        #(k4, kn4), Rtg13_a_c ↔ Rtg13_a_n # Deleted [2]
    end ksV ksD k2I k2M kn2M kBM knBM k13I k13IV k13ID k3A_c k3I_c k3I_n k13_c kn13_c k13_n kn13_n k1in k1out k3inA k3outA k3inI k3outI

    ##
    return rtgM3
end


text:
"""
    pcaeigen(X, weights = ones(size(X, 1)); nlv)
PCA by Eigen factorization.
* `X` : X-data.
* `weights` : Weights of the observations.
* `nlv` : Nb. principal components (PCs).

Noting D a (n, n) diagonal matrix of weights for the observations (rows of `X`),
the function makes an Eigen factorization of X' * D * X.

`X` is internally centered.
"""
function pcaeigen(X, weights = ones(size(X, 1)); nlv)
    pcaeigen!(copy(X), weights; nlv = nlv)
end

function pcaeigen!(X, weights = ones(size(X, 1)); nlv)
    X = ensure_mat(X)
    n, p = size(X)
    nlv = min(nlv, n, p)
    weights = mweights(weights)
    sqrtw = sqrt.(weights)
    xmeans = colmeans(X, weights)
    center!(X, xmeans)
    zX = Diagonal(sqrtw) * X
    res = eigen!(Symmetric(zX' * zX); sortby = x -> -abs(x))
    P = res.vectors[:, 1:nlv]
    eig = res.values[1:min(n, p)]
    eig[eig .< 0] .= 0
    sv = sqrt.(eig)
    T = X * P
    Pca(T, P, sv, xmeans, weights, nothing, nothing)
end

"""
    pcaeigenk(X, weights = ones(size(X, 1)); nlv)
PCA by Eigen factorization of the kernel form (XX') for wide matrices.
* `X` : X-data.
* `weights` : Weights of the observations.
* `nlv` : Nb. principal components (PCs).

Noting D a (n, n) diagonal matrix of weights for the observations (rows of X),
the function makes an Eigen factorization of D^(1/2) * X * X' D^(1/2).
This is the "kernel cross-product trick" version of the PCA algorithm (Wu et al. 1997).
For wide matrices (n << p, where p is the nb. columns) and n not too large,
this algorithm can be much faster than the others.

`X` is internally centered.

## References

Wu, W., Massart, D.L., de Jong, S., 1997. The kernel PCA algorithms for wide data. Part I: Theory and algorithms.
Chemometrics and Intelligent Laboratory Systems 36, 165-172. https://doi.org/10.1016/S0169-7439(97)00010-5
"""
function pcaeigenk(X, weights = ones(size(X, 1)); nlv)
    pcaeigenk!(copy(X), weights; nlv = nlv)
end

function pcaeigenk!(X, weights = ones(size(X, 1)); nlv)
    X = ensure_mat(X)
    n, p = size(X)
    nlv = min(nlv, n, p)
    weights = mweights(weights)
    sqrtw = sqrt.(weights)
    xmeans = colmeans(X, weights)
    center!(X, xmeans)
    zX = Diagonal(sqrtw) * X
    res = eigen!(Symmetric(zX * zX'); sortby = x -> -abs(x))
    eig = res.values[1:min(n, p)]
    eig[eig .< 0] .= 0
    sv = sqrt.(eig)
    zV = scale(res.vectors[:, 1:nlv], sv[1:nlv])
    P = zX' * zV
    T = X * P
    Pca(T, P, sv, xmeans, weights, nothing, nothing)
end


text:
# Transformations between types
using SparseArrays

export taylor_expansion_pep
export shift_and_scale
export mobius_transform

"""
    struct ShiftScaledNEP <: NEP
    ShiftScaleNEP(orgnep::NEP[,shift=0][,scale=1])

Transforms a nep (orgnep) M(λ)v to a new nep T(λ)=M(scale*λ+shift). This can be used if the method does not have an easy implementation of shift and scaling. Usage of this transformation can slow down the algorithm.
"""
struct ShiftScaledNEP <: NEP
    shift::Number
    scale::Number
    orgnep::NEP
    function ShiftScaledNEP(orgnep;shift=0,scale=1)
        return new(shift,scale,orgnep)
    end
end
"""
    shift_and_scale(orgnep::NEP;shift=0,scale=1)

Transforms the orgnep by defining a new NEP from the relation
T(λ)=M(scale * λ+shift) where M is the orgnep. This function tries
 to preserve the NEP type, e.g., a shift_and_scale operation on
an SPMF-object, return an SPMF object. If it cannot preserve
the type, it will return a nep of the struct `ShiftScaledNEP`.

# Example
```julia-repl
julia> nep0=nep_gallery("pep0")
julia> σ=3; α=10;
julia> nep1=shift_and_scale(nep0,shift=σ,scale=α)
julia> opnorm(compute_Mder(nep0,α*(4+4im)+σ)-compute_Mder(nep1,4+4im))
8.875435870738592e-12
```
"""
function shift_and_scale(orgnep::NEP;shift=0,scale=1)
    return ShiftScaledNEP(orgnep,shift=shift,scale=scale);
end

# Native implementation for PEPs by transformation of the coefficients. This is not necessarily fast.
function shift_and_scale(orgnep::PEP;shift=0,scale=1)
    At=copy(shift*scale*orgnep.A); # Allocate correct type (not so fast)
    m=size(At,1)-1
    for j=0:m
        AA = zero(At[j+1])
        for i=j:m
            factor=((scale^j)*(shift^(i-j))*factorial(i)/factorial(i-j))/factorial(j)
            AA+=orgnep.A[i+1]*factor;

        end
        At[j+1]=AA;
    end
    return PEP(At)
end

# A shift and rescale of a DEP is DEP where the matrix coefficients and the delays are rescaled and another matrix coefficient with delay zero is added.
function shift_and_scale(orgnep::DEP;shift=0,scale=1)
    J = Matrix{eltype(orgnep.A[1])}(I, size(orgnep.A[1]))
    return DEP([broadcast(*,orgnep.A,exp.(-orgnep.tauv*shift)/scale); [-shift/scale * J] ],[orgnep.tauv*scale; 0])
end

# Native implementation for SPMF. Only for generic SPMF (not AbstractSPMF)
function shift_and_scale(orgnep::SPMF_NEP;shift=0,scale=1)
    orgfv=get_fv(orgnep);
    m=size(orgfv,1);
    fv = Vector{Function}(undef, m)
    # create anonymous functions corresponding to the
    # shift and scaled problem
    for i=1:m
        fv[i] = S -> orgfv[i](scale*S + shift*one(S))::((S isa Number) ? Number : Matrix)
    end
    # create a new SPMF with transformed functions
    return SPMF_NEP(get_Av(orgnep),fv;Schur_fact=orgnep.Schur_factorize_before);
end

# Size and issparse implemented by delegation
function size(nep::ShiftScaledNEP,dim=-1)
    return size(nep.orgnep,dim)
end
function issparse(nep::ShiftScaledNEP)
    return issparse(nep.orgnep)
end

function compute_Mder(nep::ShiftScaledNEP,λ::Number,i::Integer=0)
    # Just scale the orgnep call. Chain rule for differentiation.
    return (nep.scale^i)*compute_Mder(nep.orgnep,nep.scale*λ+nep.shift,i);
end

function compute_MM(nep::ShiftScaledNEP,S,V)
    # Just call orgnep with a different S
    return compute_MM(nep.orgnep,S*nep.scale+eye(S)*nep.shift,V)
end

function compute_Mlincomb(nep::ShiftScaledNEP,λ::Number,V::AbstractVecOrMat,
                          a::Vector=ones(size(V,2)))
    # Multiply the coefficient matrix V with a diagonal matrix
    # corresponding to the chain rule
    p=size(V,2);
    z=nep.scale.^Array{eltype(V),1}(0:p-1)
    W=V * diagm(0 => z) # not very fast
    return compute_Mlincomb(nep.orgnep,nep.scale*λ+nep.shift,W,a);
end

"""
    mobius_transform(orgnep::NEP,[,a=1][,b=0][,c=0][,d=1])

Transforms a nep (orgnep) ``M(λ)v`` to a new nep ``T(λ)=M((aλ+b)/(cλ+d))``.
This function tries to preserve the type such that `T`
and `M` are of the same NEP-type (see `shift_and_scale()`).
If it cannot be preserved it will return a `MobiusTransformedNEP`.
It is in general advised to try to preserve the type, and
the use of `MobiusTransformedNEP` can
considerably slow down NEP-access.

# Example

```julia-repl
julia> nep0=nep_gallery("pep0")
julia> a=1; b=3; c=4; d=5;
julia> nep1=mobius_transform(nep0,a=a,b=b,c=c,d=d);
julia> s=3;
julia> opnorm(compute_Mder(nep0,(a*s+b)/(c*s+d))-compute_Mder(nep1,s))
0.0
```
"""
function mobius_transform(orgnep::NEP;a=1,b=0,c=0,d=1)
    return MobiusTransformedNEP(orgnep;a=a,b=b,c=c,d=d)
end

function mobius_transform(orgnep::SPMF_NEP;a=1,b=0,c=0,d=1)
    orgfv=get_fv(orgnep);
    m=size(orgfv,1);
    fv = Vector{Function}(undef, m)
    # create anonymous functions corresponding to the
    # möbius transformed problem
    for i=1:m
        fv[i] = S -> orgfv[i]((a*S + b*one(S)) / (c*S + d*one(S)))::((S isa Number) ? Number : Matrix)
    end
    # create a new SPMF with transformed functions
    return SPMF_NEP(get_Av(orgnep),fv;Schur_fact=orgnep.Schur_factorize_before);
end

"""
Hidden type representing a transformed NEP
"""
struct MobiusTransformedNEP <: NEP
    a::Number
    b::Number
    c::Number
    d::Number
    orgnep::NEP
    function MobiusTransformedNEP(orgnep;a=1,b=0,c=0,d=1)
        return new(a,b,c,d,orgnep)
    end
end

# Size and issparse implemented by delegation
function size(nep::MobiusTransformedNEP,dim=-1)
    return size(nep.orgnep,dim)
end
function issparse(nep::MobiusTransformedNEP)
    return issparse(nep.orgnep)
end

function compute_Mder(nep::MobiusTransformedNEP,λ::Number,i::Integer=0)
    # I did not found a better way
    return compute_Mder_from_MM(nep,λ,i);
end

function compute_MM(nep::MobiusTransformedNEP,S,V)
    # Just call orgnep with a different S
    return compute_MM(nep.orgnep, (nep.a*S + nep.b*I) / (nep.c*S + nep.d*I), V)
end

function compute_Mlincomb(nep::MobiusTransformedNEP,λ::Number,V::AbstractVecOrMat
                          ,a::Vector=ones(size(V,2)))
    # I did not found a better way
    return compute_Mlincomb_from_MM(nep,λ,V,a)
end

# consider renaming this function
"""
    taylor_expansion_pep(orgnep::NEP[,d=2])

Compute the truncated (with `d` terms) Taylor series of the NEP.
The output is a [`PEP`](@ref).
"""
function taylor_expansion_pep(nep::NEP,d::Integer=2)
    A=Vector{Matrix}(undef,d+1)
    for i=0:d
        A[i+1]=compute_Mder(nep,0,i)/factorial(i);
    end
    return PEP(A)
end


text:
using NormalSplines
using Random
using Test

@testset "NormalSplines.jl" begin

include("RK_Wx.jl")

include("RK_Hx.jl")

end


text:
using Robotlib
using Robotlib: T2t

T_RB_T0, T_TF_S0, POSES, MEASURED, LEDs = Robotlib.simulateCalibration_AXYB(100, σt = 0.1)

T_RB_T, T_TF_S, normDist = calibAXYB(POSES, MEASURED, refine=false)
@test T_RB_T ≈ T_RB_T0 rtol = 1e-6
@test T_TF_S ≈ T_TF_S0 rtol = 1e-6
@test length(LEDs) == 3
@test size(LEDs[1]) == (3, 100)
L = leds_to_tf(LEDs)
@test LEDs[1][:,1] ≈ T2t(MEASURED[1])
@test T2t(L[1]) ≈ T2t(MEASURED[1])
@test maximum(norm.(L .- MEASURED)) < 1e-6

## With noise

T_RB_T0, T_TF_S0, POSES, MEASURED, LEDs = Robotlib.simulateCalibration_AXYB(100, σt = 0.1, σy = 0.05)

L = leds_to_tf(LEDs)
T_RB_T, T_TF_S, normDist = calibAXYB(POSES, L, refine=false)
@test T2t(T_RB_T) ≈ T2t(T_RB_T0) rtol = 3e-1
@test T2t(T_TF_S) ≈ T2t(T_TF_S0) rtol = 3e-1

@test Rangle(T_RB_T, T_RB_T0) < deg2rad(5)
@test Rangle(T_TF_S, T_TF_S0) < deg2rad(25)

T_RB_T2, T_TF_S2, normDist2 = calibAXYB(POSES, L, refine=true)
@test T2t(T_RB_T2) ≈ T2t(T_RB_T0) rtol = 3e-1
@test T2t(T_TF_S2) ≈ T2t(T_TF_S0) rtol = 3e-1

@test Rangle(T_RB_T2, T_RB_T0) < deg2rad(5)
@test Rangle(T_TF_S2, T_TF_S0) < deg2rad(25)

# @test norm(T2t(T_RB_T2) - T2t(T_RB_T0)) < norm(T2t(T_RB_T) - T2t(T_RB_T0)) # these are not stable
# @test norm(T2t(T_TF_S2) - T2t(T_TF_S0)) < norm(T2t(T_TF_S) - T2t(T_TF_S0)) # these are not stable

# @test Rangle(T_RB_T2, T_RB_T0) < Rangle(T_RB_T, T_RB_T0) # these are not stable
# @test Rangle(T_TF_S2, T_TF_S0) < Rangle(T_TF_S, T_TF_S0) # these are not stable

@test sum(normDist2) < sum(normDist)

max_forks_repo_name: baggepinnen/Robotlib


text:
# Returns the sparsity patern matrix for matrix A,
# e.g., SparseMatrixCSC{Bool, Int}
function GetMatrixSparsity(A::AbstractVecOrMat)
    if !(A isa SparseMatrixCSC{Bool, Int})
        if A isa AbstractVector
            n = length(A)
            m = 1
        else
            (n,m) = size(A)
        end
        if n == 0 || m == 0
            return sparse(Matrix{Bool}(undef, 0, 0))
        else
            rows  = zeros(length(A))
            cols  = zeros(length(A))
            count = 0
            for r in 1:n
                for c in 1:m
                    if A[r,c] != 0.0
                        count += 1
                        rows[count] = r
                        cols[count] = c
                    end
                end
            end
            return sparse(rows[1:count],cols[1:count],[true for _ in 1:count])
        end
    else # A is already a sparsity pattern matrix
        return A
    end
end


text:
function __init__()
    global _registry_extensibleunion_to_genericfunctions
    global _registry_extensibleunion_to_members
    global _registry_genericfunctions_to_extensibleunions
    empty!(_registry_extensibleunion_to_genericfunctions)
    empty!(_registry_extensibleunion_to_members)
    empty!(_registry_genericfunctions_to_extensibleunions)
    return nothing
end


text:
### Constants

const BLANK = 0
const ONES = 0x3fe #0x1111111110

max_stars_repo_name: Teo-ShaoWei/SudokuSolverSG

max_issues_repo_name: Teo-ShaoWei/SudokuSolverSG

max_forks_repo_name: Teo-ShaoWei/SudokuSolverSG


text:
if VERSION < v"0.4-"
  import AudioIO
  const pa = AudioIO
else
  import PortAudio, LibSndFile
end
import GR

function play()
  if VERSION < v"0.4-"
    f = pa.open("Monty_Python.wav")
    data = read(f)
    close(f)

    pa.Pa_Initialize()
    stream = pa.Pa_OpenDefaultStream(0, 1, pa.paInt16, 44100.0, 1024)
    pa.Pa_StartStream(stream)
  else
    data = LibSndFile.load("Monty_Python.wav")
    stream = PortAudio.PortAudioStream(0, 1, blocksize=1024)
  end
  spectrum = zeros(Int32, 250, 250)

  start = 1
  while start + 1024 < length(data)
    amplitudes = data[start:start+1024]
    start += 1024

    if VERSION < v"0.4-"
      pa.Pa_WriteStream(stream, amplitudes)
      amplitudes /= 22050.0
    else
      PortAudio.write(stream, amplitudes)
    end

    power = log(abs(fft(amplitudes)) + 1) * 50
    spectrum[:, 1] = round(Int, power[1:250]) + 1000
    spectrum = circshift(spectrum, [0, -1])

    produce(spectrum)
  end
end

t = Task(play)
start = time_ns()

while !istaskdone(t)
  spectrum = consume(t)

  if time_ns() - start > 20 * 1000000   # 20ms
    GR.clearws()
    GR.setcolormap(-113)
    GR.setviewport(0, 1, 0, 1)
    GR.cellarray(0, 1, 0, 1, 250, 250, reshape(rotr90(spectrum), 250 * 250))
    GR.updatews()

    start = time_ns()
  end
end