text:
# A collection of solvers for submodular maximization problems with partition
# matroid constraints
#
# Some solvers are only applicable to certain problem types. Note type
# constraints.

# Utilities

# Optionally threaded loop
function optionally_threaded(f, collection; threaded = false)
  if threaded
    Threads.@threads for x in collection
      f(x)
    end
  else
    for x in collection
      f(x)
    end
  end
end

#######################################
# Basic solvers (sequential and myopic)
#######################################

# Myopic solver
function solve_myopic(p::PartitionProblem; threaded=false)
  selection = ElementArray(p)(undef, get_num_agents(p))

  optionally_threaded(1:get_num_agents(p), threaded=threaded) do ii
    # Solve given no knowledge of prior decisions
    solution_element = solve_block(p, ii, empty(p))

    selection[ii] = solution_element
  end

  evaluate_solution(p, selection)
end

# sequential solver
function solve_sequential(p::PartitionProblem)
  selection = empty(p)

  for ii = 1:get_num_agents(p)
    solution_element = solve_block(p, ii, selection)

    push!(selection, solution_element)
  end

  evaluate_solution(p, selection)
end

#############################################################
# Solvers that require direct access to the partition matroid
#
# Optimal, worst-case, and randomized solvers
#
# These solvers require explicit representations
#############################################################

# optimal solver
function solve_optimal(p::ExplicitPartitionProblem)
  indices = get_element_indices(p.partition_matroid)

  v0 = evaluate_solution(p, empty())

  # collect because product returns a tuple
  op(s::Solution, b) = partial_max(s, evaluate_solution(p, collect(b)))

  foldl(op, product(indices...); init=v0)
end

# anti-optimal solver
function solve_worst(p::ExplicitPartitionProblem)
  indices = get_element_indices(p.partition_matroid)

  v0 = Solution(Inf, empty())

  # collect because product returns a tuple
  op(s::Solution, b) = partial_min(s, evaluate_solution(p, collect(b)))

  foldl(op, product(indices...), init=v0)
end

# random solver
function solve_random(p::ExplicitPartitionProblem)
  indices = get_element_indices(p.partition_matroid)

  evaluate_solution(p, map(x->rand(x), indices))
end
# In general we require a method to produce random solutions
function sample_block(p::PartitionProblem, ind::Integer)
  error("Please define sample_block method to enable random sampling")
end
function solve_random(p::PartitionProblem)
  block_indices = 1:length(p.partition_matroid)

  evaluate_solution(p, map(x->sample_block(p, x), block_indices))
end

########################################################################
# DAG solver variants
#
# These solvers produce reduced versions of sequential greedy solvers by
# partitioning robots and pruning edges
########################################################################

export DAGSolver, solve_dag, sequence, in_neighbors, deleted_edge_weight,
solve_problem

abstract type DAGSolver end
# in_neighbors(d::DAGSolver, agent_index) = <agent in neighbors>
# sequence(d::DAGSolver) = <sequence of agent ids>

function solve_dag(d::DAGSolver, p::PartitionProblem; threaded=false)
  # The selections is an array of assignments on the matroid
  #
  # Generally, solution elements will line up with the blocks of the matroid,
  # but code should not rely on this fact
  selections = ElementArray(p)(undef, get_num_agents(p))

  for partition in partitions(d)
    function solve(agent_index)
      neighbor_selection::ElementArray(p) =
        map(x->selections[x], in_neighbors(d, agent_index))

      # Index of agent and its solution
      solution_element = solve_block(p, agent_index, neighbor_selection)

      # Store the index of agent's solution
      selections[agent_index] = solution_element
    end

    # Threaded execution is optional
    optionally_threaded(solve, partition; threaded=threaded)
  end

  evaluate_solution(p, selections)
end

solve_problem(d::DAGSolver, p::PartitionProblem; kwargs...) =
  solve_dag(d, p; kwargs...)

function deleted_edge_weight(d::DAGSolver, W::Array{Float64})
  weight = 0.0

  s = sequence(d)
  for ii in 1:length(s)
    agent_index = s[ii]

    nominal_neighbors = s[1:ii-1]
    neighbors = in_neighbors(d, agent_index)

    deleted_edges = setdiff(nominal_neighbors, neighbors)

    weight += reduce(deleted_edges; init=0.0) do w, edge
      w + W[agent_index, edge]
    end
  end

  weight
end

function deleted_edge_weight(d::DAGSolver, p::PartitionProblem)
  weights = compute_weight_matrix(p)

  deleted_edge_weight(d, weights)
end

###########################
# basic partitioned solvers
###########################

export PartitionSolver, generate_by_local_partition_size,
  generate_by_global_partition_size, solve_n_partitions

# helper function to construct the partitions for the solver
function construct_partitions(partition_numbers)
  num_partitions = maximum(partition_numbers)

  partitions = [Int64[] for x in 1:num_partitions]

  for agent_index in 1:length(partition_numbers)
    partition_index = partition_numbers[agent_index]

    push!(partitions[partition_index], agent_index)
  end

  partitions
end

# Generic partition solver
struct PartitionSolver <: DAGSolver
  # Array partitioning agents
  # Inner arrays are blocks and elements are agent ids
  partitions::Array{Array{Int64,1},1}
  # the index of the block in "partitions" containing a given agent_index
  agent_partition_numbers::Array{Int64,1}

  # Construct the dag using just the partition numbers
  PartitionSolver(x) = new(construct_partitions(x), x)
end

sequence(p::PartitionSolver) = vcat(p.partitions...)
partitions(p::PartitionSolver) = p.partitions

solver_rank(p::PartitionSolver) = length(p.agent_partition_numbers)
communication_span(p::PartitionSolver) = length(partitions(p))

function in_neighbors(p::PartitionSolver, agent_index)
  partition_index = p.agent_partition_numbers[agent_index]

  vcat(p.partitions[1:(partition_index-1)]...)
end

# general random partition solver framework from paper
function generate_by_local_partition_size(local_partition_sizes)
  partition_numbers = map(x->rand(1:x), local_partition_sizes)

  PartitionSolver(partition_numbers)
end

function generate_by_global_partition_size(num_agents, partition_size)
  generate_by_local_partition_size(fill(partition_size, num_agents))
end

# fixed number of partitions
function solve_n_partitions(num_partitions, p::PartitionProblem;
                            threaded=false)
  num_agents = length(p.partition_matroid)

  partition_solver = generate_by_global_partition_size(num_agents,
                                                       num_partitions)

  solve_dag(partition_solver, p, threaded=threaded)
end

#######################
# Adaptive partitioning
#######################

export compute_global_num_partitions, compute_local_num_partitions,
  generate_global_adaptive, generate_local_adaptive, solve_global_adaptive,
  solve_local_adaptive

# global adaptive number of partitions
function compute_global_num_partitions(desired_suboptimality,
                                       p::PartitionProblem)
  n = convert(Int64, ceil(total_weight(p)
                          / (get_num_agents(p)*desired_suboptimality)))
  max(1, n)
end

# compute nominal local number of partitions
function compute_local_num_partitions(desired_suboptimality,
                                      p::PartitionProblem)
  W  = compute_weight_matrix(p)

  map(1:length(p.partition_matroid)) do ii
    n = convert(Int64, ceil(sum(W[ii,:]) / (2 * desired_suboptimality)))
    max(1, n)
  end
end

# generate planner using global adaptive number of partitions
function generate_global_adaptive(desired_suboptimality, p::PartitionProblem)
  num_partitions = compute_global_num_partitions(desired_suboptimality, p)

  generate_by_global_partition_size(length(p.partition_matroid), num_partitions)
end

# generate planner using local adaptive number of partitions
function generate_local_adaptive(desired_suboptimality, p::PartitionProblem)
  local_partition_sizes = compute_local_num_partitions(desired_suboptimality, p)

  generate_by_local_partition_size(local_partition_sizes)
end

function solve_global_adaptive(desired_suboptimality, p::PartitionProblem)
  solver = generate_global_adaptive(desired_suboptimality, p)

  solve_dag(solver, p)
end

function solve_local_adaptive(desired_suboptimality, p::PartitionProblem)
  solver = generate_local_adaptive(desired_suboptimality, p)

  solve_dag(solver, p)
end

#######################
# range-limited solvers
#######################

export RangeSolver, solve_communication_range_limit

struct RangeSolver <: DAGSolver
  problem::PartitionProblem
  nominal_solver::DAGSolver
  communication_range::Float64
end

sequence(x::RangeSolver) = sequence(x.nominal_solver)
partitions(x::RangeSolver) = partitions(x.nominal_solver)

function in_neighbors(x::RangeSolver, agent_index)
  agents = x.problem.partition_matroid
  nominal_neighbors = in_neighbors(x.nominal_solver, agent_index)

  neighbors = Int64[]
  for neighbor in nominal_neighbors
    dist = norm(get_center(agents[neighbor]) - get_center(agents[agent_index]))

    if dist < x.communication_range
      push!(neighbors, neighbor)
    end
  end

  neighbors
end

function solve_communication_range_limit(p::PartitionProblem;
                                         num_partitions,
                                         communication_range,
                                         threaded=false)
  num_agents = length(p.partition_matroid)

  partition_solver = generate_by_global_partition_size(num_agents,
                                                       num_partitions)

  range_solver = RangeSolver(p, partition_solver, communication_range)

  solve_dag(range_solver, p, threaded=threaded)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
"""
    MapFun <: AugmentorBase.Operation

Description
--------------

Maps the given function over all individual array elements.

This means that the given function is called with an individual
elements and is expected to return a transformed element that
should take the original's place. This further implies that the
function is expected to be unary. It is encouraged that the
function should be consistent with its return type and
type-stable.

Usage
--------------

    MapFun(fun)

Arguments
--------------

- **`fun`** : The unary function that should be mapped over all
    individual array elements.

See also
--------------

[`AggregateThenMapFun`](@ref), [`ConvertEltype`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor, ColorTypes
input = testpattern()

# subtract the constant RGBA value from each pixel
augment(input, MapFun(px -> px - RGBA(0.5, 0.3, 0.7, 0.0)))

# separate channels to scale each numeric element by a constant value
pl = SplitChannels() |> MapFun(el -> el * 0.5) |> CombineChannels(RGBA)
augment(input, pl)
```
"""
struct MapFun{T} <: Operation
    fun::T
end

@inline supports_lazy(::Type{<:MapFun}) = true

function applyeager(op::MapFun, input::AbstractArray, param)
    maybe_copy(map(op.fun, input))
end

function applylazy(op::MapFun, input::AbstractArray, param)
    mappedarray(op.fun, input)
end

function showconstruction(io::IO, op::MapFun)
    print(io, nameof(typeof(op)), '(', op.fun, ')')
end

function Base.show(io::IO, op::MapFun)
    if get(io, :compact, false)
        print(io, "Map function \"", op.fun, "\" over image")
    else
        print(io, "AugmentorBase.")
        showconstruction(io, op)
    end
end

# --------------------------------------------------------------------

"""
    AggregateThenMapFun <: AugmentorBase.Operation

Description
--------------

Compute some aggregated value of the current image using the
given function `aggfun`, and map that value over the current
image using the given function `mapfun`.

This is particularly useful for achieving effects such as
per-image normalization.

Usage
--------------

    AggregateThenMapFun(aggfun, mapfun)

Arguments
--------------

- **`aggfun`** : A function that takes the whole current image as
    input and which result will also be passed to `mapfun`. It
    should have a signature of `input -> agg`, where `input` will the
    the current image. What type and value `agg` should be is up
    to the user.

- **`mapfun`** : The binary function that should be mapped over
    all individual array elements. It should have a signature of
    `(px, agg) -> new_px` where `px` is a single element of the
    current image, and `agg` is the output of `aggfun`.

See also
--------------

[`MapFun`](@ref), [`ConvertEltype`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor
input = testpattern()

# subtract the average RGB value of the current image
augment(input, AggregateThenMapFun(input -> mean(input), (px, agg) -> px - agg))
```
"""
struct AggregateThenMapFun{A,M} <: Operation
    aggfun::A
    mapfun::M
end

@inline supports_lazy(::Type{<:AggregateThenMapFun}) = true

function applyeager(op::AggregateThenMapFun, input::AbstractArray, param)
    agg = op.aggfun(input)
    maybe_copy(map(x -> op.mapfun(x, agg), input))
end

function applylazy(op::AggregateThenMapFun, input::AbstractArray, param)
    agg = op.aggfun(input)
    mappedarray(x -> op.mapfun(x, agg), input)
end

function showconstruction(io::IO, op::AggregateThenMapFun)
    print(io, nameof(typeof(op)), '(', op.aggfun, ", ", op.mapfun, ')')
end

function Base.show(io::IO, op::AggregateThenMapFun)
    if get(io, :compact, false)
        print(io, "Map result of \"", op.aggfun, "\" using \"", op.mapfun, "\" over image")
    else
        print(io, "AugmentorBase.")
        showconstruction(io, op)
    end
end

ext: jl


text:
# This file was generated, do not modify it. # hide
preds = predict(mtm, rows=test)
rmsl(y[test], preds)


text:
module MNIST2

export trainfeatures, testfeatures,
       trainlabel, testlabel,
       traindata, testdata

const IMAGEOFFSET = 16
const LABELOFFSET = 8

const NROWS = 28
const NCOLS = 28

const TRAINIMAGES = joinpath(
    dirname(@__FILE__), "..", "data", "train-images.idx3-ubyte"
)
const TRAINLABELS = joinpath(
    dirname(@__FILE__), "..", "data", "train-labels.idx1-ubyte"
)
const TESTIMAGES = joinpath(
    dirname(@__FILE__), "..", "data", "t10k-images.idx3-ubyte"
)
const TESTLABELS = joinpath(
    dirname(@__FILE__), "..", "data", "t10k-labels.idx1-ubyte"
)

function imageheader(filename::AbstractString)
    io = open(filename, "r")
    magic_number = bswap(read(io, UInt32))
    total_items = bswap(read(io, UInt32))
    nrows = bswap(read(io, UInt32))
    ncols = bswap(read(io, UInt32))
    close(io)
    return (
        magic_number,
        Int(total_items),
        Int(nrows),
        Int(ncols)
    )
end

function labelheader(filename::AbstractString)
    io = open(filename, "r")
    magic_number = bswap(read(io, UInt32))
    total_items = bswap(read(io, UInt32))
    close(io)
    return magic_number, Int(total_items)
end

function getimage(filename::AbstractString, index::Integer)
    io = open(filename, "r")
    seek(io, IMAGEOFFSET + NROWS * NCOLS * (index - 1))
    image_t = read(io, UInt8, (MNIST.NROWS, MNIST.NCOLS))
    close(io)
    return image_t'
end

function getlabel(filename::AbstractString, index::Integer)
    io = open(filename, "r")
    seek(io, LABELOFFSET + (index - 1))
    label = read(io, UInt8)
    close(io)
    return label
end

function trainimage(index::Integer)
    convert(Array{Float64}, getimage(TRAINIMAGES, index))
end

function testimage(index::Integer)
    convert(Array{Float64}, getimage(TESTIMAGES, index))
end

function trainlabel(index::Integer)
    convert(Float64, getlabel(TRAINLABELS, index))
end

function testlabel(index::Integer)
    convert(Float64, getlabel(TESTLABELS, index))
end

trainfeatures(index::Integer) = vec(trainimage(index))

testfeatures(index::Integer) = vec(testimage(index))

function traindata()
    _, nimages, nrows, ncols = imageheader(TRAINIMAGES)
    features = Array{Float64}(nrows * ncols, nimages)
    labels = Array{Float64}(nimages)
    for index in 1:nimages
        features[:, index] = trainfeatures(index)
        labels[index] = trainlabel(index)
    end
    return features, labels
end

function testdata()
    _, nimages, nrows, ncols = imageheader(TESTIMAGES)
    features = Array{Float64}(nrows * ncols, nimages)
    labels = Array{Float64}(nimages)
    for index in 1:nimages
        features[:, index] = testfeatures(index)
        labels[index] = testlabel(index)
    end
    return features, labels
end

end # module


text:
import Base: isless

isless(atomic::Threads.Atomic{Int64}, number::Int64) = isless(atomic.value, number)


text:
# fbp2_example.jl

#using Revise
using Sinograms
using MIRT: image_geom, sino_geom
# todo: ImagePhantoms
using MIRT: ellipse_im_params, ellipse_im, ellipse_sino #, :shepplogan
using MIRTjim: jim
using Plots

# Geometry
down = 2 # down-sampling factor for quicker tests
ig = image_geom(nx=512, ny=504, fov=500)
ig = ig.down(down)

sg = sino_geom(:par, nb=1088, na=984, orbit=360, orbit_start=0, d=541/949, offset=0.25,)
sg = sg.down(down)

# Phantom object
ell = ellipse_im_params(ig, :shepplogan)
xtrue = ellipse_im(ig, ell, oversample = 4, hu_scale=1000)
sino = ellipse_sino(sg, ell, oversample = 4) * 1000 # hu_scale

clim = (1 .+ (-1, 1) .* 0.05) .* 1000
p1 = jim(ig.x, ig.y, xtrue, "Phantom"; clim)
p2 = jim(sg.r, sg.ad, sino, "Sinogram", aspect_ratio=:none)

plan = fbp2(sg, ig)
@info "plan complete"
result,sino_filtered = fbp2(plan, sino)
@info "fbp complete"

p3 = jim(ig.x, ig.y, result, "FBP"; clim = clim)
#p3 = jim(ig.x, ig.y, result)
#p4 = jim(sr.r, sg.ad, sino_filtered, "Filtered Sino")
#p4 = plot(sino_filtered, label="")
p4 = jim(ig.x, ig.y, result-xtrue, "error"; clim=(-1,1).*100)

jim(p1,p2,p3,p4)

#=
#original example:

down = 2
ig = image_geom(nx=512, ny=504, fov=500)
ig = ig.down(down)
sg = sino_geom(:par, nb=888, na=984, orbit=180 , d=541/949, offset=0.25)
sg=sg.down(down)

ell = ellipse_im_params(ig, :shepplogan)
xtrue = ellipse_im(ig, ell, oversample = 4, hu_scale=1000)
sino = ellipse_sino(sg, ell, oversample = 4)

clim = (1 .+ (-1, 1) .* 0.05) .* 1000
p1 = jim(ig.x, ig.y, xtrue; clim)
p2 = jim(sg.r, sg.ad, sino, aspect_ratio=:none)

plan=fbp2(sg, ig)
println("plan complete")
result,sino_filtered=fbp2(plan,sino)
println("fbp complete")

p3 = jim(ig.x, ig.y, result; clim)
p4 = jim(real(sino_filtered))
jim(p1,p2,p3,p4)

=#


text:
"""
    Features can be extracted from QueueRecords.
They always return a FeatureValue, which allows the encoding of discrete / continuous / missing values,
which can also be forced to a Float64.
"""
abstract type AbstractFeature end

module FeatureState
    # good
    const GOOD        = 0 # value is perfectly A-okay

    # caution
    const INSUF_HIST  = 1 # value best-guess was made due to insufficient history (ex, acceleration set to zero due to one timestamp)

    # bad (in these cases fval.v is typically set to NaN as well)
    const MISSING     = 2 # value is missing (no car in front, etc.)
    const CENSORED_HI = 3 # value is past an operating threshold
    const CENSORED_LO = 4 # value is below an operating threshold
end

struct FeatureValue
    v::Float64 # feature value
    i::Int # used to encode

    FeatureValue(v::Float64, i::Int=FeatureState.GOOD) = new(v, i)
end

is_feature_valid(fval::FeatureValue) = fval.i == FeatureState.GOOD || fval.i == FeatureState.INSUF_HIST
Base.convert(::Type{Float64}, fval::FeatureValue) = fval.v

const SYMBOL_TO_FEATURE = Dict{Symbol, AbstractFeature}()

is_symbol_a_feature(sym::Symbol) = haskey(SYMBOL_TO_FEATURE, sym)
allfeatures() = values(SYMBOL_TO_FEATURE)
symbol2feature(sym::Symbol) = SYMBOL_TO_FEATURE[sym]

###

function generate_feature_functions(
    name::String,
    sym::Symbol,
    inherent_type::DataType,
    units::String;
    lowerbound::Float64=-Inf,
    upperbound::Float64=Inf,
    can_be_missing::Bool=false,
    censor_lo::Float64=NaN,
    censor_hi::Float64=NaN,
    history::Int=1,
    )

    for feature in values(SYMBOL_TO_FEATURE)
        @assert(sym != Symbol(feature), "symb: $name -> $feature")
    end
    @assert(lowerbound ≤ upperbound)

    feature_name = Symbol("Feature_" * name)
    const_name   = Symbol(uppercase(name))
    sym_feature  = Meta.quot(sym)

    @eval begin
        export $const_name
        struct $feature_name <: AbstractFeature end
        const       $const_name  = ($feature_name)()
        units(          ::$feature_name)  = $(units)
        inherent_type(  ::$feature_name)  = $(inherent_type)
        lowerbound(     ::$feature_name)  = $(lowerbound)
        upperbound(     ::$feature_name)  = $(upperbound)
        can_be_missing( ::$feature_name)  = $(can_be_missing)
        history(        ::$feature_name)  = $(history)
        Base.Symbol(    ::$feature_name)  = $sym_feature
        SYMBOL_TO_FEATURE[Symbol($const_name)] = $const_name
    end
end

function _get_feature_derivative_backwards(
    f::AbstractFeature,
    rec::EntityQueueRecord{S,D,I},
    roadway::R,
    vehicle_index::Int,
    pastframe::Int=0,
    frames_back::Int=1,
    ) where {S,D,I,R}

    id = rec[pastframe][vehicle_index].id

    retval = FeatureValue(0.0, FeatureState.INSUF_HIST)
    pastframe2 = pastframe - frames_back

    if pastframe_inbounds(rec, pastframe) && pastframe_inbounds(rec, pastframe2)

        veh_index_curr = vehicle_index
        veh_index_prev = findfirst(id, rec[pastframe2])

        if veh_index_prev != nothing
            curr = convert(Float64, get(f, rec, roadway, veh_index_curr, pastframe))
            past = convert(Float64, get(f, rec, roadway, veh_index_prev, pastframe2))
            Δt = get_elapsed_time(rec, pastframe2, pastframe)
            retval = FeatureValue((curr - past) / Δt)
        end
    end

    return retval
end

###

generate_feature_functions("PosFs", :posFs, Float64, "m")
generate_feature_functions("PosFt", :posFt, Float64, "m")
generate_feature_functions("PosFyaw", :posFyaw, Float64, "rad")
generate_feature_functions("Speed", :speed, Float64, "m/s")
generate_feature_functions("VelFs", :velFs, Float64, "m/s")
generate_feature_functions("VelFt", :velFt, Float64, "m/s")

generate_feature_functions("Acc", :acc, Float64, "m/s^2")
function Base.get(::Feature_Acc, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0) where {S,D,I,R}
    _get_feature_derivative_backwards(SPEED, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("AccFs", :accFs, Float64, "m/s²")
function Base.get(::Feature_AccFs, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0) where {S,D,I,R}
    _get_feature_derivative_backwards(VELFS, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("AccFt", :accFt, Float64, "m/s²")
function Base.get(::Feature_AccFt, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0) where {S,D,I,R}
    _get_feature_derivative_backwards(VELFT, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("Jerk", :jerk, Float64, "m/s³")
function Base.get(::Feature_Jerk, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0) where {S,D,I,R}
    _get_feature_derivative_backwards(ACC, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("JerkFs", :jerkFs, Float64, "m/s³")
function Base.get(::Feature_JerkFs, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0) where {S,D,I,R}
    _get_feature_derivative_backwards(ACCFS, rec, roadway, vehicle_index, pastframe)
end
generate_feature_functions("JerkFt", :jerkFt, Float64, "m/s³")
function Base.get(::Feature_JerkFt, rec::EntityQueueRecord{S,D,I}, roadway::R, vehicle_index::Int, pastframe::Int=0) where {S,D,I,R}
    _get_feature_derivative_backwards(ACCFT, rec, roadway, vehicle_index, pastframe)
end


text:
using CatViews
using Base.Test

@testset "CatView tests" begin
    x = CatView([1,2,3,4],[5,6,7,8])

    @test length(x) == 8
    @test ndims(x) == 1
    [ @test x[i] == i for i in 1:length(x) ]
    @test x[end] == 8
    @test_throws BoundsError x[0]
    @test_throws BoundsError x[9]

    a = randn(10,10)
    b = randn(9,9)

    x = CatView(view(a,:,2),view(b,:,3))
    @test length(x) == 19
    [ @test x[i] == a[i,2] for i in 1:size(a,1) ]
    [ @test x[size(a,1)+i] == b[i,3] for i in 1:size(b,1) ]

    ## Test mutates original object
    x[1] = 55
    @test a[1,2] == 55.0
    x[1:3] = [55, 56, 57]
    @test a[1,2] == 55.0
    @test a[2,2] == 56.0
    @test a[3,2] == 57.0

    randn!(x)
    @test isapprox(a[:,2],x[1:10])
    @test isapprox(b[:,3],x[11:end])

    ## Matrix multiply
    c = vcat(a[:,2],b[:,3])
    A = randn(19,19)
    @test isapprox(A*c,A*x)

    ## Test other constructors
    x = CatView([1,2,3,4],[5,6,7,8])
    x = CatView(x,[9,10,11,12])
    @test length(x) == 12
    [ @test x[i] == i for i in 1:length(x) ]

    y = CatView([1,2,3,4],[5,6,7,8],[9,10,11,12])
    @test length(y) == 12
    [ @test y[i] == i for i in 1:length(y) ]

    ## Iteration
    x = CatView([1,2],[3,4])
    @test all(x .== [1,2,3,4])
end

@testset "splitview tests" begin
    @testset "without preallocation" begin

        x,(A,B),s,e = splitview(Int64,(2,3),(4,5))

        # test sizes
        @test length(x) == 26
        @test size(A) == (2,3)
        @test size(B) == (4,5)
        @test eltype(A) == Int64
        @test eltype(B) == Int64

        # test partititions between vectors
        @test s == (1,7)
        @test e == (6,26)

        # set A to 1:6 by mutating x
        [ x[i] = i for i = s[1]:e[1] ]
        @test all(A == [1 3 5; 2 4 6])

        # set B to 1:20 by mutating it directly
        [ B[i] = i for i in 1:length(B) ]
        [ @test x[i+6] == i for i in 1:20 ]

        # test on a higher-order array
        sz1 = (2,2,2)
        sz2 = (2,3,2)
        sz3 = (3,2,2)
        x,(A,B,C),s,e = splitview(Int64,sz1,sz2,sz3)

        @test size(A) == sz1
        @test size(B) == sz2
        @test size(C) == sz3
        @test s == (1, length(A)+1, length(A)+length(B)+1)
        @test e == (length(A), length(A)+length(B), length(A)+length(B)+length(C))
    end

    @testset "with preallocation" begin

        # preallocate x
        x = collect(1:26)
        (A,B),s,e = splitview(x,(2,3),(4,5))

        #test sizes
        @test size(A) == (2,3)
        @test size(B) == (4,5)
        @test s == (1,7)
        @test e == (6,26)

        # test that A and B faithfully view into x
        for i = 1:6
            @test A[i] == i
            @test B[i] == 6+i
        end

        # test on higher-order array
        sz1 = (2,2,2)
        sz2 = (2,3,2)
        sz3 = (3,2,2)
        N = prod(sz1)+prod(sz2)+prod(sz3)
        x = collect(1:N)
        (A,B,C),s,e = splitview(x,sz1,sz2,sz3)

        @test eltype(A) == Int
        @test size(A) == sz1
        @test size(B) == sz2
        @test size(C) == sz3
        @test s == (1, length(A)+1, length(A)+length(B)+1)
        @test e == (length(A), length(A)+length(B), length(A)+length(B)+length(C))

        # test that arrays correctly hold 1:N
        for (i,j,X) in zip(s,e,(A,B,C))
            for (xi,Xi) in zip(i:j,eachindex(X))
                @test x[xi] == X[Xi] == xi
            end
        end

    end
end


text:
## load test dependencies and set paths to testing data
using Dates
using LightGraphs
using NCDatasets
using StaticArrays
using Test
using UnPack
using Wflow
using CSV
using BenchmarkTools
using Base.MathConstants: eulergamma
import Polynomials

# ensure test data is present
testdir = @__DIR__
datadir = joinpath(testdir, "data")
isdir(datadir) || mkdir(datadir)

"Download a test data file if it does not already exist"
function testdata(version, source_filename, target_filename)
    target_path = joinpath(datadir, target_filename)
    base_url = "https://github.com/visr/wflow-artifacts/releases/download"
    url = string(base_url, '/', string('v', version), '/', source_filename)
    isfile(target_path) || download(url, target_path)
    return target_path
end

staticmaps_rhine_path = testdata(v"0.1", "staticmaps.nc", "staticmaps-rhine.nc")
staticmaps_moselle_path = testdata(v"0.2.1", "staticmaps.nc", "staticmaps-moselle.nc")
staticmaps_lahn_path = testdata(v"0.2.1", "staticmaps-lahn.nc", "staticmaps-lahn.nc")
forcing_moselle_path = testdata(v"0.2", "forcing-2000.nc", "forcing-moselle.nc")
forcing_lahn_path = testdata(v"0.2", "forcing-lahn.nc", "forcing-lahn.nc")
instates_moselle_path = testdata(v"0.2.1", "instates-moselle.nc", "instates-moselle.nc")
forcing_sbm_gw_path = testdata(v"0.2.1", "forcing-sbm-groundwater.nc", "forcing-sbm-groundwater.nc")
staticmaps_sbm_gw_path = testdata(v"0.2.1", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")

include("testing_utils.jl")

## run all tests
@testset "Wflow.jl" begin
    include("horizontal_process.jl")
    include("io.jl")
    include("vertical_process.jl")
    include("reservoir_lake.jl")
    include("run_sbm.jl")
    include("run_hbv.jl")
    include("run_sbm_gwf.jl")
    include("run.jl")
    include("groundwater.jl")
    include("utils.jl")
end


text:
using Random, LinearAlgebra, BlockArrays, FillArrays
using Stheno: BS, unbox, are_conformal, chol, ABM, LazyPDMat, Xt_invA_X, Xt_invA_Y

@testset "block_arrays" begin

    # Test construction of `BlockVector` from a vector of vectors. Also test copying.
    let
        rng, N, N′ = MersenneTwister(123456), 5, 6
        x, x′ = randn(rng, N), randn(rng, N′)
        x̂ = BlockVector([x, x′])
        @test length(x̂) == N + N′
        @test x̂ == vcat(x, x′)
        @test getblock(x̂, 1) == x && getblock(x̂, 2) == x′
        @test blocksizes(x̂, 1) == [N, N′]
        @test blocksizes(x̂) == (blocksizes(x̂, 1),)
        @test blocklengths(x̂) == blocksizes(x̂, 1)
    end

    # Test construction of `BlockMatrix` from matrix of matrices. Also test copying.
    let
        rng, P1, P2, P3, Q1, Q2 = MersenneTwister(123456), 2, 3, 6, 4, 5
        X11, X12 = randn(rng, P1, Q1), randn(rng, P1, Q2)
        X21, X22 = randn(rng, P2, Q1), randn(rng, P2, Q2)
        X31, X32 = randn(rng, P3, Q1), randn(rng, P3, Q2)
        X = BlockMatrix(reshape([X11, X21, X31, X12, X22, X32], 3, 2))
        @test size(X, 1) == P1 + P2 + P3
        @test size(X, 2) == Q1 + Q2
        @test getblock(X, 1, 1) == X11
        @test getblock(X, 1, 2) == X12
        @test getblock(X, 2, 1) == X21
        @test getblock(X, 2, 2) == X22

        @test blocksizes(X, 1) == [P1, P2, P3]
        @test blocksizes(X, 2) == [Q1, Q2]
        @test blocksizes(X) == (blocksizes(X, 1), blocksizes(X, 2))

        @test BlockMatrix([X11, X21, X31, X12, X22, X32], 3, 2) == X

        Y = BlockMatrix([X11, X21, X31])
        @test size(Y, 1) == P1 + P2 + P3
        @test size(Y, 2) == Q1
        @test getblock(Y, 1, 1) == X11
        @test getblock(Y, 2, 1) == X21
        @test blocksizes(Y, 1) == [P1, P2, P3]
        @test blocksizes(Y, 2) == [Q1]
    end

    # Test Symmetric block matrix construction and util.
    let
        rng, P1, P2 = MersenneTwister(123456), 3, 4
        A11, A12, A22 = randn(rng, P1, P1), randn(rng, P1, P2), randn(rng, P2, P2)
        A_ = BlockMatrix([A11, zeros(Float64, P2, P1), A12, A22], 2, 2)
        A = Symmetric(A_)
        @test unbox(A) === A_
        @test Symmetric(A) === A
        @test nblocks(A) == nblocks(A_)
        @test nblocks(A, 1) == nblocks(A_, 1)
        @test nblocks(A, 2) == nblocks(A_, 2)
        @test blocksize(A, 1, 1) == blocksize(A_, 1, 1) &&
            blocksize(A, 1, 2) == blocksize(A_, 1, 2) &&
            blocksize(A, 2, 1) == blocksize(A_, 2, 1) &&
            blocksize(A, 2, 2) == blocksize(A_, 2, 2)
        @test size(A) == size(A_)
        @test getindex(A, 1, 2) == getindex(A_, 1, 2)
        @test eltype(A) == eltype(A_)

        @test getblock(A, 1, 1) === Symmetric(getblock(A_, 1, 1))
        @test getblock(A, 2, 1) == transpose(getblock(A_, 1, 2))
        @test getblock(A, 1, 2) === getblock(A_, 1, 2)
        @test getblock(A, 2, 2) == Symmetric(getblock(A_, 2, 2))

    end

    # Test triangular basics.
    let
        rng, P1, P2 = MersenneTwister(123456), 3, 4
        A11, A12, A22 = randn(rng, P1, P1), randn(rng, P1, P2), randn(rng, P2, P2)
        A = BlockMatrix([A11, zeros(Float64, P2, P1), A12, A22], 2, 2)
        U, L = UpperTriangular(A), LowerTriangular(A)

        # UpperTriangular functionality.
        @test unbox(U) === A
        @test blocksize(U, 1, 1) == blocksize(A, 1, 1)
        @test blocksize(U, 2, 1) == blocksize(A, 2, 1)
        @test blocksizes(U, 1) == blocksizes(A, 1)
        @test blocksizes(U, 2) == blocksizes(A, 2)
        @test blocksizes(U) == blocksizes(A)

        @test getblock(U, 1, 1) == UpperTriangular(getblock(A, 1, 1))
        @test getblock(U, 2, 1) === Zeros{Float64}(P2, P1)
        @test getblock(U, 2, 2) == UpperTriangular(getblock(A, 2, 2))

        @test BlockMatrix(U) isa AbstractBlockMatrix
        @test Matrix(BlockMatrix(U)) == Matrix(U)

        # LowerTriangular functionality.
        @test unbox(L) === A
        @test blocksize(L, 1, 1) == blocksize(A, 1, 1)
        @test blocksize(L, 2, 1) == blocksize(A, 2, 1)
        @test blocksizes(L, 1) == blocksizes(A, 1)
        @test blocksizes(L, 2) == blocksizes(A, 2)
        @test blocksizes(L) == blocksizes(A)

        @test getblock(L, 1, 1) == LowerTriangular(getblock(A, 1, 1))
        @test getblock(L, 1, 2) === Zeros{Float64}(P1, P2)
        @test getblock(L, 2, 2) == LowerTriangular(getblock(A, 2, 2))

        @test BlockMatrix(L) isa AbstractBlockMatrix
        @test Matrix(BlockMatrix(L)) == Matrix(L)

        # Note the subtle difference here. This is correct.
        @test getblock(U, 1, 2) === getblock(A, 1, 2)
        @test getblock(L, 2, 1) == getblock(A, 2, 1)
    end

    # Test copying.
    let
        # BlockVector
        rng, N, N′ = MersenneTwister(123456), 5, 6
        x, x′ = randn(rng, N), randn(rng, N′)
        x̂ = BlockVector([x, x′])

        @test copy(x̂) isa BlockVector
        @test copy(x̂) == x̂

        # BlockMatrix
        rng, P1, P2 = MersenneTwister(123456), 2, 3
        X11, X12 = randn(rng, P1, P1), randn(rng, P1, P2)
        X21, X22 = randn(rng, P2, P1), randn(rng, P2, P2)
        X = BlockMatrix(reshape([X11, X21, X12, X22], 2, 2))
        @test copy(X) isa BlockMatrix
        @test copy(X) == X

        # block-symmetric
        B = Symmetric(X)
        @test copy(B) isa Symmetric{T, <:AbstractBlockMatrix{T}} where T
        @test copy(B) == B

        # LowerTriangular
        L = LowerTriangular(X)
        @test copy(L) isa LowerTriangular{T, <:AbstractBlockMatrix{T}} where T
        @test copy(L) == L

        # UpperTriangular
        U = UpperTriangular(X)
        @test copy(U) isa UpperTriangular{T, <:AbstractBlockMatrix{T}} where T
        @test copy(U) == U

    end

    # Transposition
    let
        # Regular block matrices
        rng, P1, P2, P3, Q1, Q2 = MersenneTwister(123456), 2, 3, 6, 4, 5
        X11, X12 = randn(rng, P1, Q1), randn(rng, P1, Q2)
        X21, X22 = randn(rng, P2, Q1), randn(rng, P2, Q2)
        X31, X32 = randn(rng, P3, Q1), randn(rng, P3, Q2)
        X = BlockMatrix(reshape([X11, X21, X31, X12, X22, X32], 3, 2))

        for foo in [adjoint, transpose]

            @test Matrix(foo(X)) == foo(Matrix(X))

            @test nblocks(foo(X), 1) == nblocks(X, 2)
            @test nblocks(foo(X), 2) == nblocks(X, 1)
            @test nblocks(foo(X)) == reverse(nblocks(X))

            @test blocksizes(foo(X), 1) == blocksizes(X, 2)
            @test blocksizes(foo(X), 2) == blocksizes(X, 1)

            @test size(foo(X), 1) == size(X, 2)
            @test size(foo(X), 2) == size(X, 1)
            @test size(foo(X)) == reverse(size(X))

            @test foo(foo(X)) == X
        end

        # Symmetric matrices should be invariant under transposition
        rng, P1, P2 = MersenneTwister(123456), 2, 3
        X11, X12 = randn(rng, P1, P1), randn(rng, P1, P2)
        X21, X22 = randn(rng, P2, P1), randn(rng, P2, P2)
        X = BlockMatrix(reshape([X11, X21, X12, X22], 2, 2))
        B = Symmetric(X)

        for foo in [adjoint, transpose]
            @test foo(B) isa BS
            @test foo(B) === B
        end

        # Triangular block matrices
        for foo in [adjoint, transpose]
            U = UpperTriangular(X)
            # @test foo(U) isa LowerTriangular{T, <:AbstractBlockMatrix{T}} where T
            @test Matrix(foo(U)) == foo(Matrix(U))

            L = LowerTriangular(X)
            # @test foo(L) isa UpperTriangular{T, <:AbstractBlockMatrix{T}} where T
            @test Matrix(foo(L)) == foo(Matrix(L))
        end
    end

    # # # Test transposition of block vectors.
    # # let
    # #     # @test transpose(getblock(x̂, 1)) == getblock(transpose(x̂), 1)
    # #     # @test adjoint(getblock(x̂, 1)) == getblock(adjoint(x̂), 1)
    # #     # @test blocklengths(x̂) == blocklengths(transpose(x̂))
    # #     # @test blocklengths(x̂) == blocklengths(adjoint(x̂))
    # # end

    # Test multiplication of a `BlockMatrix` by a `BlockVector`.
    let
        rng, P1, P2, Q1, Q2 = MersenneTwister(123456), 2, 3, 4, 5
        X11, X12 = randn(rng, P1, Q1), randn(rng, P1, Q2)
        X21, X22 = randn(rng, P2, Q1), randn(rng, P2, Q2)
        x1, x2 = randn(rng, Q1), randn(rng, Q2)
        A, x = BlockMatrix([X11, X21, X12, X22], 2, 2), BlockVector([x1, x2])

        @test length(A * x) == P1 + P2
        @test are_conformal(A, x)
        @test A * x isa AbstractBlockVector
        @test A * x ≈ Matrix(A) * Vector(x)

        x̂1, x̂2 = randn(rng, P1), randn(rng, P2)
        x̂ = BlockVector([x̂1, x̂2])

        @test are_conformal(A', x̂)
        @test A' * x̂ isa AbstractBlockVector
        @test Vector(A' * x̂) ≈ Matrix(A)' * Vector(x̂)

        @test are_conformal(transpose(A), x̂)
        @test transpose(A) * x̂ isa AbstractBlockVector
        @test Vector(transpose(A) * x̂) ≈ transpose(Matrix(A)) * Vector(x̂)

        @test_throws AssertionError A * Vector(x)
        @test_throws AssertionError transpose(A) * Vector(x̂)
        @test_throws AssertionError A' * Vector(x̂)
    end

    # Test multiplication of a `BlockMatrix` by a `BlockMatrix`.
    let
        rng, P1, P2, Q1, Q2, R1, R2, R3 = MersenneTwister(123456), 2, 3, 4, 5, 6, 7, 8
        A11, A12 = randn(rng, P1, Q1), randn(rng, P1, Q2)
        A21, A22 = randn(rng, P2, Q1), randn(rng, P2, Q2)
        B11, B12, B13 = randn(rng, Q1, R1), randn(rng, Q1, R2), randn(rng, Q1, R3)
        B21, B22, B23 = randn(rng, Q2, R1), randn(rng, Q2, R2), randn(rng, Q2, R3)
        A = BlockMatrix([A11, A21, A12, A22], 2, 2)
        B = BlockMatrix([B11, B21, B12, B22, B13, B23], 2, 3)
        @test size(A * B, 1) == P1 + P2
        @test size(A * B, 2) == R1 + R2 + R3

        @test are_conformal(A, B)
        @test A * B isa AbstractBlockMatrix
        @test A * B ≈ Matrix(A) * Matrix(B)

        for foo in [adjoint, transpose]
            @test foo(B) * foo(A) isa AbstractBlockMatrix
            @test foo(B) * foo(A) ≈ foo(Matrix(B)) * foo(Matrix(A))

            @test foo(B) * B isa AbstractBlockMatrix
            @test foo(B) * B ≈ foo(Matrix(B)) * Matrix(B)

            @test B * foo(B) isa AbstractBlockMatrix
            @test B * foo(B) ≈ Matrix(B) * foo(Matrix(B))
        end
    end

    # Test `chol`, logdet, and backsolving.
    let
        rng, P1, P2, P3 = MersenneTwister(123456), 3, 4, 5
        tmp = randn(rng, P1 + P2 + P3, P1 + P2 + P3)
        A_ = tmp * tmp' + 1e-9I
        A = BlockArray(A_, [P1, P2, P3], [P1, P2, P3])
        @assert A_ == A

        # Compute chols and compare.
        U_, U = cholesky(Symmetric(A_)).U, cholesky(Symmetric(A)).U
        @test U isa UpperTriangular{<:Real, <:ABM}
        @test U_ ≈ U
        @test U_ ≈ Matrix(U)

        # Test `logdet`.
        @test logdet(U) ≈ logdet(U_)

        # Test backsolving for block vector.
        x1, x2, x3 = randn(rng, P1), randn(rng, P2), randn(rng, P3)
        x = BlockVector([x1, x2, x3])

        @test U \ x isa AbstractBlockVector
        @test size(U \ x) == size(U_ \ Vector(x))
        @test U \ x ≈ U_ \ Vector(x)

        @test U' \ x isa AbstractBlockVector
        @test typeof(U') <: Adjoint{<:Real, <:UpperTriangular{<:Real, <:ABM}}
        @test size(U' \ x) == size(U_' \ Vector(x))
        @test U' \ x ≈ U_' \ Vector(x)

        @test transpose(U) \ x isa AbstractBlockVector
        @test typeof(transpose(U)) <: Transpose{<:Real, <:UpperTriangular{<:Real, <:ABM}}
        @test size(transpose(U) \ x) == size(U_' \ Vector(x))
        @test transpose(U) \ x ≈ U_' \ Vector(x)

        # Test backsolving for block matrix
        Q1, Q2 = 7, 6
        X11, X12 = randn(rng, P1, Q1), randn(rng, P1, Q2)
        X21, X22 = randn(rng, P2, Q1), randn(rng, P2, Q2)
        X31, X32 = randn(rng, P3, Q1), randn(rng, P3, Q2)
        X = BlockMatrix([X11, X21, X31, X12, X22, X32], 3, 2)

        @test U \ X isa AbstractBlockMatrix
        @test size(U \ X) == size(U_ \ Matrix(X))
        @test U \ X ≈ U_ \ Matrix(X)

        @test U' \ X isa AbstractBlockMatrix
        @test size(U' \ X) == size(U_' \ Matrix(X))
        @test U' \ X ≈ U_' \ Matrix(X)
    end

    # Test UniformScaling interaction.
    let
        rng, N = MersenneTwister(123456), 5
        X = BlockArray(randn(rng, N, N), [2, 3], [2, 3])
        @test X + I isa AbstractBlockMatrix
        @test I + X isa AbstractBlockMatrix
        @test Matrix(X + I) == Matrix(X) + I
        @test Matrix(I + X) == I + Matrix(X)

        Xs = Symmetric(X)
        @test Xs + I isa Symmetric{T, <:AbstractBlockMatrix{T}} where T
        @test I + Xs isa Symmetric{T, <:AbstractBlockMatrix{T}} where T
        @test Matrix(Xs + I) == Matrix(Xs) + I
        @test Matrix(I + Xs) == I + Matrix(Xs)
    end

    # Run covariance matrix tests with a BlockMatrix.
    let
        rng, P1, P2, Q1, Q2, R1, R2 = MersenneTwister(123456), 2, 3, 4, 5, 6, 7
        P, Q, R = P1 + P2, Q1 + Q2, R1 + R2
        B = randn(rng, P, P)
        A_ = BlockArray(B' * B + UniformScaling(1e-6), [P1, P2], [P1, P2])
        A = LazyPDMat(Symmetric(A_), 1e-12)
        x = BlockArray(randn(rng, P), [P1, P2])
        X = BlockArray(randn(rng, P, Q), [P1, P2], [Q1, Q2])
        X′ = BlockArray(randn(rng, P, R), [P1, P2], [R1, R2])

        # Check utility functionality.
        @test size(A) == size(A_)
        @test size(A, 1) == size(A_, 1)
        @test size(A, 2) == size(A_, 2)

        @test Matrix(A) == A_
        @test A == A

        # Test unary operations.
        @test logdet(A) ≈ logdet(Matrix(A_))
        @test chol(A) isa UpperTriangular{T, <:AbstractBlockMatrix{T}} where T
        @test chol(A) ≈ chol(Matrix(A_) + A.ϵ * I)

        # Test binary operations.
        @test x + x isa AbstractBlockVector
        @test Vector(x + x) == Vector(x) + Vector(x)
        @test x - x isa AbstractBlockVector
        @test Vector(x - x) == Vector(x) - Vector(x)

        @test A_ + A_ isa AbstractBlockMatrix
        @test Matrix(A + A) == A_ + A_
        @test A_ - A_ isa AbstractBlockMatrix
        @test Matrix(A - A) == A_ - A_

        @test chol(A).data isa AbstractBlockMatrix
        @test x' * (Matrix(A_) \ x) ≈ Xt_invA_X(A, x)
        @test Xt_invA_X(A, X) isa LazyPDMat
        @test chol(Xt_invA_X(A, X)).data isa AbstractBlockMatrix
        @test maximum(abs.(Matrix(X') * (Matrix(A_) \ Matrix(X)) - Matrix(Xt_invA_X(A, X)))) < 1e-6
        @test Xt_invA_Y(X′, A, X) isa AbstractBlockMatrix
        @test maximum(abs.(Matrix(X′)' * (Matrix(A_) \ Matrix(X)) - Xt_invA_Y(X′, A, X))) < 1e-6
        @test A \ X isa AbstractBlockMatrix
        @test maximum(abs.(Matrix(A_) \ X - A \ X)) < 1e-6
    end
end


text:
filter_ext(path, ext::AbstractString) = filter(fn -> splitext(fn)[2] == ext, readdir(path))
filter_ext(path, exts::Union{Tuple,AbstractArray}) =
    filter(fn -> splitext(fn)[2] in exts, readdir(path))
filter_ext(path, ext::Nothing) = readdir(path)

maybewindow2indices(A, window::Tuple) = maybewindow2indices(A, dims(A), window::Tuple)
maybewindow2indices(A, dims::Tuple, window::Tuple) =
    window == () ? () : to_indices(A, DD.dims2indices(dims, window))

# Read from the paraent dataset, using the window indices if they exist
# We try to load as little data from disk as possible.
readwindowed(A::AbstractGeoArray, window::Tuple{}) = GeoArray(A)
readwindowed(A, window::Tuple{}) = Array(A)
readwindowed(A, window::Tuple{}, I...) = A[I...]
readwindowed(A, window::Tuple, I...) = A[Base.reindex(window, I)...]
readwindowed(A, window::Tuple) = readwindowed(A, window...)
readwindowed(A, i, I...) = A[i, I...]
readwindowed(A) = Array(A)

# Get a metadata field
getmeta(A::AbstractGeoArray, key, fallback) = getmeta(metadata(A), key, fallback)
getmeta(m::Metadata, key, fallback) = get(val(m), key, fallback)
getmeta(m::NoMetadata, key, fallback) = fallback

# Check that array order matches expectation
checkarrayorder(A, order::Order) = map(d -> checkarrayorder(d, order), dims(A))
checkarrayorder(A, order::Tuple) = map(checkarrayorder, dims(A), order)
checkarrayorder(dim::Dimension, order::Order) =
    arrayorder(dim) == order || @warn "Array order for `$(DD.basetypeof(order))` is `$(arrayorder(dim))`, usualy `$order`"

checkindexorder(A, order::Order) = map(d -> checkindexorder(d, order), dims(A))
checkindexorder(A, order::Tuple) = map(checkindexorder, dims(A), order)
checkindexorder(dim::Dimension, order::Order) =
    indexorder(dim) == order || @warn "Array order for `$(DD.basetypeof(order))` is `$(indexorder(dim))`, usualy `$order`"

cleankeys(keys) = Tuple(Symbol.(keys))

function _to_sampled(A)
    if any(m -> !(m isa AbstractSampled), mode(A))
        rebuild(A; data=data(A), dims=map(_to_sampled, dims(A)))
    end
end
_to_sampled(mode, dim) = dim
_to_sampled(mode::AbstractSampled, dim::SpatialDim) = dim
_to_sampled(mode, dim::SpatialDim) = rebuild(dim; mode=Sampled(mode(dim)))


text:
include("../covariatefun.jl")

using Ipopt, JuMP
import JuMP.MathOptInterface
const MOI = MathOptInterface

function opt_UCB(cdf, lx, ux)
    cdf_fixed(v) = cdf(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
    d = length(lx)
    fun(x) = x[1] #x\in R^d+1
    #cdf_wrapper(x) = cdf_fixed(collect(x))
    cdf_wrapper(x) = cdf_fixed(x)
    model = Model(Ipopt.Optimizer)
    register(model, :constraint, 3, cdf_wrapper; autodiff = true)
    register(model, :objective, 3, fun; autodiff = true)
    @variable(model, lx[i] <=au[i=1:3] <= ux[i])
    @NLobjective(model, Min, objective(au))
    @NLconstraint(model, constraint(au) == 0.25)
    JuMP.optimize!(model)
    return value.(model)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module FastCGI

using Sockets

import Base: show, close, write, flush, ==
export FCGIServer, FCGIClient, FCGIRequest, process, isrunning, stop, set_server_param, set_server_runner

include("bufferedoutput.jl")
include("types.jl")
include("runners.jl")
include("requesthandler.jl")
include("server.jl")
include("client.jl")

end # module


text:
using Test
using ViennaRNA
using ViennaRNA: LibRNA
using Unitful

@testset "FoldCompound" begin
    for s in ["GGGAAACCC",
              ["GGG-A-CC-", "G-GGC-G-G"],
              ["GGG-A-CC-", "G-GGC-G-G", "CGC----CG"]]
        fc = FoldCompound(s)
        @test length(fc) == 9
        @test sum(size(fc)) == length(fc)

        # test FoldCompound properties
        @test fc.circular         == LibRNA.VRNA_MODEL_DEFAULT_CIRC
        @test fc.dangles          == LibRNA.VRNA_MODEL_DEFAULT_DANGLES
        @test fc.gquadruplex      == LibRNA.VRNA_MODEL_DEFAULT_GQUAD
        @test fc.has_exp_matrices == false
        @test fc.log_ML           == LibRNA.VRNA_MODEL_DEFAULT_LOG_ML
        @test fc.min_loop_size    == LibRNA.TURN
        @test fc.no_GU_basepairs  == LibRNA.VRNA_MODEL_DEFAULT_NO_GU
        @test fc.no_GU_closure    == LibRNA.VRNA_MODEL_DEFAULT_NO_GU_CLOSURE
        @test fc.no_lonely_pairs  == LibRNA.VRNA_MODEL_DEFAULT_NO_LP
        @test fc.nstrands         == 1
        @test fc.params_name      == "RNA - Turner 2004"
        @test fc.special_hairpins == LibRNA.VRNA_MODEL_DEFAULT_SPECIAL_HP
        @test fc.temperature      == 37u"°C"
        @test fc.uniq_ML          == LibRNA.VRNA_MODEL_DEFAULT_UNIQ_ML

        # options
        fc = FoldCompound(s; options=ViennaRNA.LibRNA.VRNA_OPTION_MFE)
        @test length(fc) == 9
        # use different energy parameter sets
        fc = FoldCompound(s; params=:RNA_Turner1999)
        @test length(fc) == 9
        @test_throws ArgumentError FoldCompound(s; params=:UNKNOWN_ENERGY_PARAMS)
        # use different temperatures to evaluate energies
        fc = FoldCompound(s; temperature=55u"°C")
        @test length(fc) == 9
        @test fc.temperature == 55u"°C"
        fc = FoldCompound(s; temperature=310u"K")
        @test length(fc) == 9
        @test fc.temperature == 310u"K"
        @test_throws Unitful.DimensionError FoldCompound(s; temperature=100u"m")

        # circular
        fc = FoldCompound(s; circular=true)
        @test length(fc) == 9
        @test fc.circular == true
        # dangles
        fc = FoldCompound(s; dangles=1)
        @test length(fc) == 9
        @test fc.dangles == 1
        @test_throws ArgumentError FoldCompound(s; dangles=-1)
        @test_throws ArgumentError FoldCompound(s; dangles=4)
        # gquadruplex
        fc = FoldCompound(s; gquadruplex=true)
        @test length(fc) == 9
        @test fc.gquadruplex == true
        # log_ML
        fc = FoldCompound(s; log_ML=true)
        @test length(fc) == 9
        @test fc.log_ML == true
        # min_loop_size
        fc = FoldCompound(s; min_loop_size=2)
        @test length(fc) == 9
        @test fc.dangles == 2
        @test_throws ArgumentError FoldCompound(s; min_loop_size=-1)
        # no_GU_basepairs
        fc = FoldCompound(s; no_GU_basepairs=true)
        @test length(fc) == 9
        @test fc.no_GU_basepairs == true
        # no_GU_closure
        fc = FoldCompound(s; no_GU_closure=true)
        @test length(fc) == 9
        @test fc.no_GU_closure == true
        # no_lonely_pairs
        fc = FoldCompound(s; no_lonely_pairs=true)
        @test length(fc) == 9
        @test fc.no_lonely_pairs == true
        # special_hairpins
        fc = FoldCompound(s; special_hairpins=false)
        @test length(fc) == 9
        @test fc.special_hairpins == false
        # uniq_ML
        fc = FoldCompound(s; uniq_ML=true)
        @test length(fc) == 9
        @test fc.uniq_ML == true
        # show
        buf = IOBuffer()
        show(buf, MIME("text/plain"), fc)
        @test length(String(take!(buf))) > 0
    end

    # multiple strands
    fc = FoldCompound("GGG&AAA&CCC")
    @test length(fc) == 9
    @test size(fc) == (3, 3, 3)
    @test fc.nstrands == 3

    # multistrand alifold (not supported yet)
    @test_throws ArgumentError fc = FoldCompound(["GGG&AAA&CCC", "G-G&AA-&CCA"])
    # @test length(fc) == 9
    # @test size(fc) == (3, 3, 3)
    # @test fc.nstrands == 3

    # multiple strands don't work in comparative (alifold) mode
    @test_throws ArgumentError FoldCompound(["GGGA", "GCGG", "G"])
    @test_throws ArgumentError FoldCompound(["GG&C", "AA&U"])

    # test setting of min_loop_size with single and multiple strands
    # - this didn't work in ViennaRNA-2.5.0 with multiple strands
    #   (e.g. "GGG&CCC") and h=3 (min_loop_size was always set to 0 by
    #   ViennaRNA). Works as of ViennaRNA-2.5.1.
    for h = 0:5
        fc = FoldCompound("GGG"; min_loop_size=h)
        @test fc.min_loop_size == h
        fc = FoldCompound("GGG&CCC"; min_loop_size=h)
        @test fc.min_loop_size == h
        fc = FoldCompound("GGG&CCC&AAAA"; min_loop_size=h)
        @test fc.min_loop_size == h
        fc = FoldCompound(["GGG", "G-C"]; min_loop_size=h)
        @test fc.min_loop_size == h
        fc = FoldCompound(["GGG", "G-C", "-AG"]; min_loop_size=h)
        @test fc.min_loop_size == h
    end
end

@testset "Pairtable" begin
    pt = Pairtable("(((...)))")
    @test length(pt) == 9
    # getindex
    @test pt[1] == 9
    @test pt[2] == 8
    @test pt[3] == 7
    @test pt[4] == 0
    @test pt[5] == 0
    @test pt[6] == 0
    @test pt[7] == 3
    @test pt[8] == 2
    @test pt[9] == 1
    @test_throws BoundsError pt[0]
    @test_throws BoundsError pt[length(pt) + 1]
    # setindex!
    pt[3] = 4
    @test pt[3] == 4
    @test pt[4] == 3
    @test pt[7] == 0
    @test_throws ArgumentError pt[1] = 1
    @test_throws ArgumentError pt[0] = 0
    @test_throws ArgumentError pt[length(pt) + 1] = 0
    @test_throws ArgumentError pt[1] = -1
    @test_throws ArgumentError pt[1] = length(pt) + 1
    # show
    buf = IOBuffer()
    show(buf, MIME("text/plain"), pt)
    @test length(String(take!(buf))) > 0
    # basepairs
    @test Set(basepairs(Pairtable("(((...)))"))) == Set([(1,9), (2,8), (3,7)])
end

@testset "bp_distance" begin
    @test bp_distance("((..))", "(....)") == 1
    @test_throws ArgumentError bp_distance("()", ".")
end

@testset "tree_edit_dist" begin
    @test tree_edit_dist("...", "()") == 3.0f0
end

@testset "mean_bp_distance" begin
    seq = "GGGAAACCC"
    fc = FoldCompound(seq)
    @test_throws ArgumentError mean_bp_distance(fc)
    partfn(fc)
    @test mean_bp_distance(fc) isa AbstractFloat
    @test mean_bp_distance(seq) isa AbstractFloat
end

@testset "ensemble_defect" begin
    seq = "GGGAAACCCC"
    str = "(((....)))"
    fc = FoldCompound(seq)
    partfn(fc)
    @test ensemble_defect(fc, str) isa AbstractFloat
    @test ensemble_defect(fc, Pairtable(str)) isa AbstractFloat
    @test_throws ArgumentError ensemble_defect(fc, ".")
    @test_throws ArgumentError ensemble_defect(fc, Pairtable("."))
    @test ensemble_defect(seq, Pairtable(str)) isa AbstractFloat
    @test ensemble_defect(seq, str) isa AbstractFloat
end

@testset "prob_of_structure" begin
    seq = "GGGAAACCCC"
    str = "(((....)))"
    fc = FoldCompound(seq)
    partfn(fc)
    @test prob_of_structure(fc, str) isa AbstractFloat
    @test_throws ArgumentError prob_of_structure(fc, ".")
    @test prob_of_structure(seq, str) isa AbstractFloat
end

@testset "energy" begin
    seq = "GGGAAAACCCC"
    str = "(((....)))."
    fc = FoldCompound(seq)
    @test energy(fc, str) isa Unitful.Quantity
    @test energy(seq, str) isa Unitful.Quantity
    redirect_stdout(devnull) do
        # suppress stdout output from ViennaRNA for verbose option
        @test energy(fc, str; verbose=true) isa Unitful.Quantity
        @test energy(seq, str; verbose=true) isa Unitful.Quantity
    end
    @test_throws ArgumentError energy(fc, ".")
    @test_throws ArgumentError energy("(...)", ".")
    # test different energy parameter sets
    @test energy(FoldCompound(seq; params=:RNA_Turner1999), str) ≈ -2.90u"kcal/mol" atol=1e-3u"kcal/mol"
    @test energy(FoldCompound(seq; params=:RNA_Andronescu2007), str) ≈ -2.15u"kcal/mol" atol=1e-3u"kcal/mol"
    # test different temperatures
    @test energy(FoldCompound(seq; temperature=30u"°C"), str) ≈ -3.52u"kcal/mol" atol=1e-3u"kcal/mol"
    @test energy(FoldCompound(seq; temperature=37u"°C"), str) ≈ -2.90u"kcal/mol" atol=1e-3u"kcal/mol"
    @test energy(FoldCompound(seq; temperature=55u"°C"), str) ≈ -1.24u"kcal/mol" atol=1e-3u"kcal/mol"
    @test energy(FoldCompound(seq; temperature=323.15u"K"), str) ≈ -1.69u"kcal/mol" atol=1e-3u"kcal/mol" # 323.15u"K" is 50u"°C"
end

@testset "mfe" begin
    seq = "GGGAAACCCC"
    fc = FoldCompound(seq)
    @test mfe(fc) isa Tuple{String,Unitful.Quantity}
    @test mfe(seq) isa Tuple{String,Unitful.Quantity}
end

@testset "partfn" begin
    seq = "GGGAAACCCC"
    fc = FoldCompound(seq)
    @test partfn(fc) isa Tuple{String,Unitful.Quantity}
    @test partfn(seq) isa Tuple{String,Unitful.Quantity}
end

@testset "bpp" begin
    seq = "GGGAAACCCC"
    n = length(seq)

    fc = FoldCompound(seq)
    @test_throws ArgumentError bpp(fc)
    partfn(fc)
    p = bpp(fc)
    @test eltype(p) <: AbstractFloat
    @test size(p) == (n,n)

    p = bpp(seq)
    @test eltype(p) <: AbstractFloat
    @test size(p) == (n,n)
end

@testset "pbacktrack" begin
    seq = "GGGGGAAAAACCCCCCCCAUUCA"
    fc = FoldCompound(seq; uniq_ML=true)
    partfn(fc)
    @test length(pbacktrack(fc)) == 1
    @test length(pbacktrack(seq)) == 1
    s = pbacktrack(fc; num_samples=10)
    @test length(s) == 10
    s = pbacktrack(fc; num_samples=5,
                   options=ViennaRNA.LibRNA.VRNA_PBACKTRACK_NON_REDUNDANT)
    @test length(s) == 5
end

@testset "mea" begin
    seq = "GGGGGAAAAACCCCCCCCAUUCA"
    fc = FoldCompound(seq)
    partfn(fc)
    @test mea(fc) isa Tuple{String,AbstractFloat}
    @test mea(fc; gamma=1.0) isa Tuple{String,AbstractFloat}
    @test mea(fc; gamma=0.5) isa Tuple{String,AbstractFloat}
    @test mea(seq) isa Tuple{String,AbstractFloat}
    @test mea(seq; gamma=1.0) isa Tuple{String,AbstractFloat}
    @test mea(seq; gamma=0.5) isa Tuple{String,AbstractFloat}
    fc = FoldCompound(seq)
    @test_throws ArgumentError mea(fc)
end

@testset "centroid" begin
    seq = "GGGGGAAAAACCCCCCCCAUUCA"
    fc = FoldCompound(seq)
    partfn(fc)
    @test centroid(fc) isa Tuple{String,AbstractFloat}
    @test centroid(seq) isa Tuple{String,AbstractFloat}
    fc = FoldCompound(seq)
    @test_throws ArgumentError centroid(fc)
end

@testset "subopt" begin
    seq = "GGGGGAAAAACCCCCCCCAUUCA"
    fc = FoldCompound(seq; uniq_ML=true)
    s = subopt(fc; delta=5u"kcal/mol")
    @test s isa Vector{Tuple{String,Unitful.Quantity}}
    s = subopt(fc; delta=5u"kcal/mol", sorted=true)
    @test s isa Vector{Tuple{String,Unitful.Quantity}}
    s = subopt(seq; delta=5u"kcal/mol")
    @test s isa Vector{Tuple{String,Unitful.Quantity}}
    s = subopt(seq; delta=5u"kcal/mol", sorted=true)
    @test s isa Vector{Tuple{String,Unitful.Quantity}}
end

@testset "subopt_zuker" begin
    seq = "GGGGAAAACCCC"
    # TODO: the following sequences generate warnings
    # seq = "GGGUAAACCCAUUCAC"
    # seq = "GGGGGAAAAACCCCCCCCAUUCA"
    fc = FoldCompound(seq)
    s = subopt_zuker(fc)
    @test s isa Vector{Tuple{String,Unitful.Quantity}}
    s = subopt_zuker(seq)
    @test s isa Vector{Tuple{String,Unitful.Quantity}}
end

@testset "inverse_fold" begin
    s = "((((....))))"
    @test inverse_fold("A"^length(s), s) isa Tuple{String,AbstractFloat}
    @test_throws ArgumentError inverse_fold("A", "()")
    @test inverse_pf_fold("A"^length(s), s) isa Tuple{String,Unitful.Quantity}
    @test_throws ArgumentError inverse_pf_fold("A", "()")
end

@testset "neighbors" begin
    fc = FoldCompound("GGGAAACCC")
    pt = Pairtable(".((...)).")
    @test neighbors(fc, pt) == [[(-2, -8)], [(-3, -7)], [(1, 9)]]
end

@testset "plot_coords" begin
    function test_plot_xy(s, x, y)
        @test x isa Vector{Float32}
        @test length(x) == length(s)
        @test y isa Vector{Float32}
        @test length(y) == length(s)
    end
    s = "(((...)))"
    pt = Pairtable(s)

    # test: plot_coords(::String), plot_coords(::Pairtable)
    x, y = plot_coords(s)
    test_plot_xy(s, x, y)
    x, y = plot_coords(pt)
    test_plot_xy(pt, x, y)
    for plot_type in (:simple, :naview, :circular, :turtle, :puzzler)
        x, y = plot_coords(s; plot_type)
        test_plot_xy(s, x, y)
        x, y = plot_coords(pt; plot_type)
        test_plot_xy(pt, x, y)
    end
    @test_throws ArgumentError x, y = plot_coords(s; plot_type = :unknown)
    @test_throws ArgumentError x, y = plot_coords(pt; plot_type = :unknown)
    # zero-sized inputs
    x, y = plot_coords("")
    test_plot_xy("", x, y)
    x, y = plot_coords(Pairtable(""))
    test_plot_xy(Pairtable(""), x, y)
end

@testset "heat_capacity" begin
    s = "GGGAAACCC"
    fc = FoldCompound(s)
    Tmin = 10u"°C"
    Tmax = 60u"°C"
    Tincrement = 1u"°C"
    mpoints = 2
    # TODO: this range call doesn't work without ustrip
    n = length(range(start=ustrip(Tmin), stop=ustrip(Tmax), step=ustrip(Tincrement)))

    hcs = heat_capacity(fc, Tmin, Tmax, Tincrement; mpoints)
    @test hcs isa Vector{Tuple{typeof(1.0f0u"°C"),typeof(1.0f0u"kcal/mol/K")}}
    @test length(hcs) == n
    hcs = heat_capacity(fc, Tmin, Tmax; mpoints)
    @test hcs isa Vector{Tuple{typeof(1.0f0u"°C"),typeof(1.0f0u"kcal/mol/K")}}
    @test length(hcs) == n
    hcs = heat_capacity(s, Tmin, Tmax, Tincrement; mpoints)
    @test hcs isa Vector{Tuple{typeof(1.0f0u"°C"),typeof(1.0f0u"kcal/mol/K")}}
    @test length(hcs) == n
    hcs = heat_capacity(s, Tmin, Tmax; mpoints)
    @test hcs isa Vector{Tuple{typeof(1.0f0u"°C"),typeof(1.0f0u"kcal/mol/K")}}
    @test length(hcs) == n
end

include("utils.jl")
include("plot_structure.jl")


text:
using BinaryBuilder

name = "JSBSim"
version = v"1.1.5"
julia_version = v"1.5.3"

# Collection of sources required to build JSBSim
sources = [
    GitSource("https://github.com/JSBSim-Team/jsbsim.git",
              "9bc866713431d3296441697a91517b32c9c017d3"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/jsbsim
sed -ri 's/PROJECT_VERSION_PATCH\s\"6\"/PROJECT_VERSION_PATCH \"5\"/g' CMakeLists.txt
mkdir build && cd build
FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    FLAGS+=(-DCMAKE_CXX_FLAGS_RELEASE="-D_POSIX_C_SOURCE")
fi
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_DOCS=OFF \
    -DBUILD_PYTHON_MODULE=OFF \
    -DBUILD_JULIA_PACKAGE=ON \
    -DJulia_PREFIX=${prefix} \
    "${FLAGS[@]}" \
    ..
cmake --build . --target JSBSimJL -- -j${nproc}
install_license $WORKSPACE/srcdir/jsbsim/COPYING
cp julia/*JSBSimJL*.$dlext $libdir/.
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)

filter!(p -> libc(p) != "musl" && !Sys.isfreebsd(p), platforms) # muslc is not supported

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    Dependency("libcxxwrap_julia_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8",
    julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Distributions
import Statistics, Plots #Warning: unrecommened - to the personal standard- Reason: has a huge set of other packages- once initial deguggin finished -Remove
import AbstractPlotting.MakieLayout
#--- Bernoulli

"""
Returns a Bernoulli Probability Density Function (PDF) that is;

1. Discrete process
2. I.I.D

"""
function pdf(ntimes)

  return Bernoulli(ntimes)
end

function cdf(pdf, x)
  return sum(pdf(x)) #discrete sum
end
#function #later

function sign(value = 0.5, comparisonFunction = rand(), trueValue = -1, falseValue = 1)

  multiplier += ifelse(comparisonFunction < value, trueValue, falseValue)
  return multiplier
end

"""
A function useful for benchmarking julia's statistical functions

  tested
"""
function _globalRandomWalker()
  @time begin
    pos = 0
    count = 0
    _mean = 0
    Mean = 0
    data = []
    numsteps = 10^4
    numwalkers = 10^4
    final_square_positions = Int[]

    for i in 1:numwalkers
      for j in 1:numsteps
        pos += ifelse(rand() < 0.5, -1, +1)
        count += 1
      end
      push!(final_square_positions, pos^2)
      push!(data, pos)
    end
    println("Total Displacement " * string(pos))
    _mean = pos / count #adding  a naive arithmetic mean
    Mean = Statistics.mean(final_square_positions)
    println("Mean Displacement " * string(_mean))
    println("Mean Squared Displacement " * string(Mean))

    return pos, _mean, Mean, data
  end
end

pos, _mean, Mean, data = _globalRandomWalker()
#root-squared mean
trueMean = Statistics.mean(data)
n = size(data)[1]
rsMean = sqrt(Mean)
#=Q. is mean ≈ rsMean? if not , is it dispered by how much
ended up with 3 means:
1.trueMean (based on Statistics.jl)
2.rsMean (square root of squared values) [hugely large , Disqualified] - say, what if this is the only one that truely expresses the mean?
3._mean (non-weighted, naive mean )
=#
error2 = max(_mean, rsMean) - min(_mean, rsMean) # rsMean did not work  here very smoothly
error = max(_mean, trueMean) - min(_mean, trueMean)
#=by only comparing the 2 results above, only a correct Mean would be close to the true one i.e. calculated naive  _mean =#
#then we can visualize error  with y

x = 1:n; #preferable to use linspace
y = data; # These are the plotting data
Plots.plot(x, y)

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using Clang

const SFML_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include")
const CONFIG_H = joinpath(SFML_INCLUDE, "SFML", "Config.h")

# create a work context
ctx = DefaultContext()

# parse headers
parse_headers!(ctx, [CONFIG_H], includes=[SFML_INCLUDE, CLANG_INCLUDE])

# settings
ctx.options["is_struct_mutable"] = false  # for nested struct

# write output
for trans_unit in ctx.trans_units
    root_cursor = getcursor(trans_unit)
    push!(ctx.cursor_stack, root_cursor)
    header = spelling(root_cursor)
    @info "wrapping header: $header ..."
    # loop over all of the child cursors and wrap them, if appropriate.
    ctx.children = children(root_cursor)
    for (i, child) in enumerate(ctx.children)
        child_name = name(child)
        child_header = filename(child)
        ctx.children_index = i
        # choose which cursor to wrap
        startswith(child_name, "__") && continue  # skip compiler definitions
        child_name in keys(ctx.common_buffer) && continue  # already wrapped
        child_header != header && continue

        wrap!(ctx, child)
    end
end

# write "common" definitions: types, typealiases, etc.
common_file = joinpath(@__DIR__, "..", "src", "config.jl")
open(common_file, "w") do f
    println(f, "# Automatically generated using Clang.jl\n")
    print_buffer(f, dump_to_buffer(ctx.common_buffer))
end


text:
module TestHetero
using GaussianProcesses, ScikitLearnBase
using GaussianProcesses: set_params!, get_params, num_params,
                         update_cK!, update_mll!, update_target_and_dtarget!
using Test, Random

Random.seed!(1)
@testset "GPE heteroscedastic" begin
    d, n = 3, 10

    X = 2π * rand(d, n)
    y = [sum(sin, view(X, :, i)) / d for i in 1:n]
    mZero = MeanZero()
    kern = SE(0.0, 0.0)
    logNoise = range(0.0, stop=1.0, length=n)

    ntest = 5
    Xtest = randn(d, ntest)

    @testset "VectorParam" begin
        vp = GaussianProcesses.wrap_param(logNoise)
        @test all(get_params(vp) .== logNoise)
        @test num_params(vp) == n
        @test length(vp) == n
    end

    @testset "GPE constructors" begin
        gp1 = GPE(X, y, mZero, kern, logNoise)
        gp2 = GPE(X, y, mZero, kern, GaussianProcesses.wrap_param(logNoise))
    end

    gp = GP(X, y, mZero, kern, logNoise)

    @testset "update_cK!" begin
        update_cK!(gp)
    end
    @testset "update_mll" begin
        update_mll!(gp)
        @test isfinite(gp.mll)
    end

    @testset "update_dmll!" begin
        update_target_and_dtarget!(gp; noise=false)
        @test isfinite(gp.mll)
        @test isfinite(gp.target)
        @test all(isfinite.(gp.dmll))
        @test all(isfinite.(gp.dtarget))
    end

    @testset "num_params" begin
        @test num_params(gp) == 0+2+n
    end

    @testset "show" begin
        show(devnull, gp) # doesn't crash
    end

    @testset "Predictive mean at test locations" begin
        f_pred, sig = predict_f(gp, Xtest)
        f_pred, sig = predict_f(gp, Xtest; full_cov=true)
    end

    # ScikitLearn interface test
    @testset "ScikitLearn interface" begin
        gp_sk = ScikitLearnBase.fit!(GPE(), X', y)
        f_pred = ScikitLearnBase.predict(gp_sk, X')
    end

    # Modify kernel and update
    @testset "Update" begin
        gp.kernel.ℓ2 = 4.0
        X_pred = 2π * rand(d, n)

        GaussianProcesses.update_target!(gp)
        f_pred, sig = predict_f(gp, X_pred)
    end

    #Check that the rand function works
    @testset "Random GP sampling" begin
        X_test = 2π * rand(d, n)
        samples = rand(gp, X_test)
    end

    @testset "params round trip" begin
        params_1 = deepcopy(get_params(gp))
        set_params!(gp, params_1)
        params_2 = get_params(gp)
        @test params_1 ≈ params_2
    end

end
end


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

name = "OpenCV"
version = v"0.5.0"

# Collection of sources required to complete build
sources = [
    GitSource(
      "https://github.com/opencv/opencv.git",
      "d5fd2f0155ffad366f9ac912dfd6d189a7a6a98e",
    ),
]
# Bash recipe for building across all platforms
script = raw"""
Julia_PREFIX=$prefix

mkdir build
cd build
cmake -DCMAKE_FIND_ROOT_PATH=$prefix \
      -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      ../opencv/
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
cd ..
install_license opencv/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
 Linux(:i686, libc=:glibc),
 Linux(:x86_64, libc=:glibc),
 Linux(:aarch64, libc=:glibc),
 Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
 Linux(:powerpc64le, libc=:glibc),
 Linux(:i686, libc=:musl),
 Linux(:x86_64, libc=:musl),
 Linux(:aarch64, libc=:musl),
 Linux(:armv7l, libc=:musl, call_abi=:eabihf),
 MacOS(:x86_64),
 #FreeBSD(:x86_64), <- fails,
 Windows(:i686),
 Windows(:x86_64),
] |> expand_cxxstring_abis

# The products that we will ensure are always built
products = [
    LibraryProduct(["libopencv_calib3d", "libopencv_calib3d450"], :libopencv_calib3d),
    LibraryProduct(["libopencv_core", "libopencv_core450"], :libopencv_core),
    LibraryProduct(["libopencv_dnn", "libopencv_dnn450"], :libopencv_dnn),
    LibraryProduct(["libopencv_features2d", "libopencv_features2d450"], :libopencv_features2d),
    LibraryProduct(["libopencv_flann", "libopencv_flann450"], :libopencv_flann),
    LibraryProduct(["libopencv_gapi", "libopencv_gapi450"], :libopencv_gapi),
    LibraryProduct(["libopencv_highgui", "libopencv_highgui450"], :libopencv_highgui),
    LibraryProduct(["libopencv_imgcodecs", "libopencv_imgcodecs450"], :libopencv_imgcodecs),
    LibraryProduct(["libopencv_imgproc", "libopencv_imgproc450"], :libopencv_imgproc),
    LibraryProduct(["libopencv_ml", "libopencv_ml450"], :libopencv_ml),
    LibraryProduct(["libopencv_objdetect", "libopencv_objdetect450"], :libopencv_objdetect),
    LibraryProduct(["libopencv_photo", "libopencv_photo450"], :libopencv_photo),
    LibraryProduct(["libopencv_stitching", "libopencv_stitching450"], :libopencv_stitching),
    LibraryProduct(["libopencv_video", "libopencv_video450"], :libopencv_video),
    LibraryProduct(["libopencv_videoio", "libopencv_videoio450"], :libopencv_videoio),
]

# Dependencies that must be installed before this package can be built
dependencies = []

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


text:
module JLLisp
export T, Atom, Number_, List, Null, Nil

abstract type T end
abstract type Atom <: T end
abstract type Number_ <: Atom end
abstract type List <: T end
struct Null <: List end
const Nil = Null()
Base.string(x::Null) = "NIL"

module Cons_
    import ..JLLisp
    mutable struct Cons <: JLLisp.List
        car::JLLisp.T
        cdr::JLLisp.T
    end

    function Cons()
        Cons(JLLisp.Null(), JLLisp.Null())
    end

    function Base.string(cons::Cons)
        str = ""
        list = cons
        str *= "("
        while true
            str *= string(list.car)
            if list.cdr == JLLisp.Nil
                str *= ")"
                break
            elseif isa(list.cdr, JLLisp.Atom)
                str *= " . " * string(list.cdr) * ")"
                break
            else
                str *= " "
                list = list.cdr
            end
        end
        return str
    end
end # Cons_

module Eval
    import ..JLLisp
    const maxstacksize = 65536
    const stack = Vector{JLLisp.T}(undef, maxstacksize)
    stackP = 1

    function eval_(form::JLLisp.T)
        if isa(form, JLLisp.Symbols.Symbol_)
            try
                global symbolvalue = form.value
            catch
                error("Unbound Variable Error: $(form)")
            end
            return symbolvalue
        end

        if isa(form, JLLisp.Null) return form end
        if isa(form, JLLisp.Atom) return form end
        car = form.car
        isa(car, JLLisp.Symbols.Symbol_) || error("Not a Symbol: $(car)")
        try
            global fun = car.fn
        catch
            error("Undefined Function Error: $(car)")
        end

        # system functions
        if isa(fun, JLLisp.Function_.Func)
            argumentlist = form.cdr
            return JLLisp.Function_.funcall(fun, argumentlist)
        end

        # evaluate S expression
        if isa(fun, JLLisp.Cons_.Cons)
            cdr = fun.cdr
            lambdalist = cdr.car
            body = cdr.cdr
            lambdalist == JLLisp.Null() && return evalbody(body)

            return bindevalbody(lambdalist, body, form.cdr)
        end
        error("Not a Function: $(fun)")
    end

    # bindevalbody を実装
    function bindevalbody(lambda::JLLisp.Cons_.Cons, body::JLLisp.Cons_.Cons, form::JLLisp.Cons_.Cons)
        oldstackP = stackP
        while true
            ret = eval_(form.car)
            global stack[stackP] = ret
            global stackP += 1
            form.cdr == JLLisp.Nil && break
            form = form.cdr
        end

        arglist = lambda
        sp = oldstackP
        while true
            sym = arglist.car
            sym.value, stack[sp] = stack[sp], sym.value
            sp += 1
            arglist.cdr == JLLisp.Nil && break
            arglist = arglist.cdr
        end

        ret = evalbody(body)

        arglist = lambda
        global stackP = oldstackP
        while true
            sym= arglist.car
            sym.value = stack[oldstackP]
            oldstackP += 1
            arglist.cdr == JLLisp.Nil && break
            arglist = arglist.cdr
        end

        return ret
    end

    # bindbody
    function evalbody(body::JLLisp.Cons_.Cons)
        while true
            body.cdr == JLLisp.Nil && return JLLisp.Eval.eval_(body.car)
            body = body.cdr
        end
    end
end # Eval

module Symbols
    import ..JLLisp

    mutable struct Symbol_ <: JLLisp.Atom
        name::String
        value::JLLisp.T
        fn::JLLisp.T

        function Symbol_(name::String)
            x = new()
            x.name = name
            x.value = JLLisp.Nil
            x.fn = JLLisp.Nil
            x
        end
    end

    Base.string(x::Symbol_) = x.name
    const symboltable = Dict{String, Symbol_}()
    function symbol_(name::String)
        get!(symboltable, name, Symbol_(name))
    end

    const symbolT = symbol_("T")
    symbolT.value = symbolT
    const symbolQuit = symbol_("QUIT")

end # Symbols_

module Integer__
    import ..JLLisp
    struct Integer_ <: JLLisp.Number_
        value::Integer
    end

    lisptf(bool) = bool ? JLLisp.Symbols.symbolT : JLLisp.Null()
    add(a::Integer_, b::Integer_) = Integer_(+(a.value, b.value))
    sub(a::Integer_, b::Integer_) = Integer_(-(a.value, b.value))
    mul(a::Integer_, b::Integer_) = Integer_(*(a.value, b.value))
    Base.div(a::Integer_, b::Integer_) = Integer_(div(a.value, b.value))
    ge(a::Integer_, b::Integer_) = lisptf(>=(a.value, b.value))
    le(a::Integer_, b::Integer_) = lisptf(<=(a.value, b.value))
    gt(a::Integer_, b::Integer_) = lisptf(>(a.value, b.value))
    lt(a::Integer_, b::Integer_) = lisptf(<(a.value, b.value))
    numberequal(a::Integer_, b::Integer_) = lisptf(a.value == b.value)
    Base.string(x::Integer_) = string(x.value)
end # Integer__

module Function_
    import ..JLLisp
    abstract type Func <: JLLisp.Atom end

    function funcall(fn::Func, arguments::JLLisp.List)
        return JLLisp.Null()
    end

    function regist(name::String, fn::Func)
        sym = JLLisp.Symbols.symbol_(name)
        sym.fn = fn
    end

    struct Car <: Func end
    struct Cdr <: Func end
    struct FunCons <: Func end
    struct Eq <: Func end
    struct Add <: Func end
    struct Sub <: Func end
    struct Mul <: Func end
    struct Div <: Func end
    struct Ge <: Func end
    struct Le <: Func end
    struct Gt <: Func end
    struct Lt <: Func end
    struct NumberEqual <: Func end
    struct Quote <: Func end
    struct Setq <: Func end
    struct Defun <: Func end
    struct If <: Func end
    struct TypeOf <: Func end
    struct SymbolFunction <: Func end

    function funcall(fn::Car, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        return arg1 == JLLisp.Null() ? JLLisp.Null() : arg1.car
    end

    function funcall(fn::Cdr, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        return arg1 == JLLisp.Null() ? JLLisp.Null() : arg1.cdr
    end

    function funcall(fn::FunCons, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Cons_.Cons(arg1, arg2)
    end

    function funcall(fn::Eq, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        if arg1 == arg2
            return JLLisp.Symbols.symbolT
        else
            return JLLisp.Nil
        end
    end

    function funcall(fn::Add, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.add(arg1, arg2)
    end

    function funcall(fn::Sub, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.sub(arg1, arg2)
    end

    function funcall(fn::Mul, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.mul(arg1, arg2)
    end

    function funcall(fn::Div, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.div(arg1, arg2)
    end

    function funcall(fn::Ge, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.ge(arg1, arg2)
    end

    function funcall(fn::Le, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.le(arg1, arg2)
    end

    function funcall(fn::Gt, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.gt(arg1, arg2)
    end

    function funcall(fn::Lt, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.lt(arg1, arg2)
    end

    function funcall(fn::NumberEqual, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        return JLLisp.Integer__.numberequal(arg1, arg2)
    end

    function funcall(fn::Quote, arguments::JLLisp.List)
        return arguments.car
    end

    function funcall(fn::Setq, arguments::JLLisp.List)
        arg1 = arguments.car
        arg2 = JLLisp.Eval.eval_(arguments.cdr.car)
        sym = arg1
        value = JLLisp.Eval.eval_(arg2)
        sym.value = value
        return value
    end

    function funcall(fn::Defun, arguments::JLLisp.List)
        arg1 = arguments.car
        args = arguments.cdr
        fun = arg1
        lambda = JLLisp.Cons_.Cons()
        lambda.car = JLLisp.Symbols.symbol_("LAMBDA")
        lambda.cdr = args
        fun.fn = lambda
        return fun
    end

    function funcall(fn::If, arguments::JLLisp.List)
        arg1 = arguments.car
        args = arguments.cdr
        arg2 = args.car
        arg3 = args.cdr == JLLisp.Nil ? JLLisp.Nil : args.cdr.car
        if JLLisp.Eval.eval_(arg1) != JLLisp.Nil
            return JLLisp.Eval.eval_(arg2)
        else
            return JLLisp.Eval.eval_(arg3)
        end
    end

    function funcall(fn::TypeOf, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        type = typeof(arg1).name |> string |> uppercase
        return JLLisp.Symbols.symbol_(type)
    end

    function funcall(fn::SymbolFunction, arguments::JLLisp.List)
        arg1 = JLLisp.Eval.eval_(arguments.car)
        return arg1.fn
    end

    function registSystemFunctions()
        regist("CAR", Car())
        regist("CDR", Cdr())
        regist("CONS", FunCons())
        regist("EQ", Eq())
        regist("+", Add())
        regist("-", Sub())
        regist("*", Mul())
        regist("/", Div())
        regist(">=", Ge())
        regist("<=", Le())
        regist(">", Gt())
        regist("<", Lt())
        regist("=", NumberEqual())
        regist("QUOTE", Quote())
        regist("SETQ", Setq())
        regist("DEFUN", Defun())
        regist("IF", If())
        regist("TYPE-OF", TypeOf())
        regist("SYMBOL-FUNCTION", SymbolFunction())
    end
    registSystemFunctions()
end # Function_

module Reader
    using ..JLLisp
    const charbuffersize = 256
    charbuff = Vector{Char}()
    ch = ' '
    line = ""
    indexofline = 1
    linelength = 0

    function read(;ex = "")
        if isempty(ex)
            global line = readline()
        else # For test
            global line = ex
        end
        prepare()
        return getSexp()
    end

    function prepare()
        global indexofline = 1
        global linelength = length(line)
        global charbuff = Vector{Char}(line)
        push!(charbuff, '\0')
        getchar()
    end

    function getchar()
        global ch = charbuff[indexofline]
        global indexofline += 1
    end

    function skipspace()
        while isspace(ch)
            getchar()
        end
    end

    function getSexp()
        while true
            skipspace()
            ch == '('   &&  return makelist()
            ch == '\''  &&  return makequote()
            ch == '-'   &&  return makeminusnumber()
            isdigit(ch) &&  return makenumber()
            return makesymbol()
        end
    end

    function makelist()
        getchar()
        skipspace()
        if ch == ')'
            getchar()
            return JLLisp.Null()
        end
        top = JLLisp.Cons_.Cons()
        list = top
        while true
            list.car = getSexp()
            skipspace()
            indexofline > linelength && JLLisp.Null() # 不等号要注意
            ch == ')' && break
            if ch == '.'
                getchar()
                list.cdr = getSexp()
                skipspace()
                getchar()
                return top
            end
            list.cdr = JLLisp.Cons_.Cons()
            list = list.cdr
        end
        getchar()
        return top
    end

    function makequote()
        top = JLLisp.Cons_.Cons()
        list = top
        list.car = JLLisp.Symbols.symbol_("QUOTE")
        list.cdr = JLLisp.Cons_.Cons()
        list = list.cdr
        getchar()
        list.car = getSexp()
        return top
    end

    function makeminusnumber()
        nch = charbuff[indexofline]
        isdigit(nch) || return makesymbolinternal(Char[ch])
        return makenumber()
    end

    function makenumber()
        str = Vector{Char}()
        if ch == '-'
            push!(str, ch)
            getchar()
        end
        while indexofline <= linelength + 1
            ch == '(' || ch == ')' && break
            isspace(ch) && break
            if !isdigit(ch)
                global indexofline -= 1
                return makesymbolinternal(str)
            end
            push!(str, ch)
            getchar()
        end
        value = parse(Int, String(str))
        return JLLisp.Integer__.Integer_(value)
    end

    function makesymbol()
        global ch = uppercase(ch)
        str = Char[ch]
        return makesymbolinternal(str)
    end

    function makesymbolinternal(str::Vector{Char})
        while indexofline < linelength + 1
            getchar()
            ch in ('(', ')') && break
            isspace(ch) && break
            global ch = uppercase(ch)
            push!(str, ch)
        end
        symstr = String(str)

        symstr == "NIL" && return JLLisp.Null()
        return JLLisp.Symbols.symbol_(symstr)
    end
end # Reader

module TopLevel
    export repl
    using ..JLLisp
    function repl()
        println("Welcome to JLLisp! (2019-8-11)")
        println("> Copyright (C) goropikari 2019.")
        println("> Type quit and hit Enter for leaving JLLisp.")
        JLLisp.Function_.registSystemFunctions()
        while true
            try
                print("> ")
                sexp = JLLisp.Reader.read()
                sexp == JLLisp.Symbols.symbolQuit && break
                @async if eof(stdin) # exit when Ctrl-d is typed
                    println("\nBye!")
                    exit(0)
                end
                ret = JLLisp.Eval.eval_(sexp)
                println(string(ret))
            catch e
                println(typeof(e))
            end
        end
        println("Bye!")
    end
end

using .TopLevel
export repl

end # JLLisp module


text:
@testset "QUBO creation" begin
   # TODO fill
end

@testset "QUBO operations" begin
   # TODO fill
end


text:
function event_occurred(obj::FieldAbstractDetector, arg0::FieldSpacecraftState, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (FieldSpacecraftState, jboolean), arg0, arg1)
end

function g(obj::FieldAbstractDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_handler(obj::FieldAbstractDetector)
    return jcall(obj, "getHandler", FieldEventHandler, ())
end

function get_max_check_interval(obj::FieldAbstractDetector)
    return jcall(obj, "getMaxCheckInterval", RealFieldElement, ())
end

function get_max_iteration_count(obj::FieldAbstractDetector)
    return jcall(obj, "getMaxIterationCount", jint, ())
end

function get_threshold(obj::FieldAbstractDetector)
    return jcall(obj, "getThreshold", RealFieldElement, ())
end

function init(obj::FieldAbstractDetector, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate), arg0, arg1)
end

function is_forward(obj::FieldAbstractDetector)
    return jcall(obj, "isForward", jboolean, ())
end

function reset_state(obj::FieldAbstractDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "resetState", FieldSpacecraftState, (FieldSpacecraftState,), arg0)
end

function with_handler(obj::FieldAbstractDetector, arg0::FieldEventHandler)
    return jcall(obj, "withHandler", FieldEventDetector, (FieldEventHandler,), arg0)
end

function with_max_check(obj::FieldAbstractDetector, arg0::RealFieldElement)
    return jcall(obj, "withMaxCheck", FieldEventDetector, (RealFieldElement,), arg0)
end

function with_max_iter(obj::FieldAbstractDetector, arg0::jint)
    return jcall(obj, "withMaxIter", FieldEventDetector, (jint,), arg0)
end

function with_threshold(obj::FieldAbstractDetector, arg0::RealFieldElement)
    return jcall(obj, "withThreshold", FieldEventDetector, (RealFieldElement,), arg0)
end


text:
@testset "GHC" begin
  Z = [ones(10,10) 2ones(10,10); 3ones(10,10) 4ones(10,10)]
  𝒮 = georef((Z=Z,))
  p = partition(𝒮, GHC(4, 1.0))
  @test length(p) == 4
  @test all(nelements.(p) .== 100)
  p1, p2, p3, p4 = p
  @test mean(coordinates(centroid(p1, ind)) for ind in 1:nelements(p1)) == [5.0,5.0]
  @test mean(coordinates(centroid(p2, ind)) for ind in 1:nelements(p2)) == [15.0,5.0]
  @test mean(coordinates(centroid(p3, ind)) for ind in 1:nelements(p3)) == [5.0,15.0]
  @test mean(coordinates(centroid(p4, ind)) for ind in 1:nelements(p4)) == [15.0,15.0]

  C = cluster(𝒮, GHC(4, 1.0))
  @test C.cluster == vec(Z')

  𝒮′ = georef(values(𝒮), centroid.(domain(𝒮)))
  C′ = cluster(𝒮′, GHC(4, 1.0))
  @test C.cluster == C′.cluster

  𝒮 = georef((z=[√(i^2+j^2) for i in 1:50, j in 1:50],))
  p = partition(𝒮, GHC(50, 1.0))
  @test length(p) == 50

  if visualtests
    @test_reference "data/ghc.png" plot(p)
  end
end


text:
import JSON

struct Emitter
  io::IO
  cache::Cache
end

function make_emitter(io, verbose::Bool)
  let cache = verbose ? NoopCache() : RollingCache()
    Emitter(io, cache)
  end
end

function emit_raw(e::Emitter, s::AbstractString)
  print(e.io, s)
end

function emit_tag(e::Emitter, x::AbstractString)
  emit(e, "~$x", true)
end

function emit(e::Emitter, x::AbstractString, cacheable::Bool)
  let value = iscacheable(x, cacheable) ? write!(e.cache, x) : x
      print(e.io, JSON.json(value))
  end
end

function emit(e::Emitter, x::Integer)
  print(e.io, JSON.json(x))
end

function emit(e::Emitter, x::Bool)
  print(e.io, JSON.json(x))
end

function emit_null(e::Emitter, askey::Bool)
  askey ? emit_tag(e, "_") : emit_raw(e, "null")
end

function emit_array_start(e::Emitter)
  print(e.io, "[")
end

function emit_array_end(e::Emitter)
  print(e.io, "]")
end

function emit_map_start(e::Emitter)
  print(e.io, "{")
end

function emit_map_sep(e::Emitter, i=2)
  emit_array_sep(e, i)
end

function emit_key_sep(e::Emitter)
  print(e.io, ":")
end

function emit_map_end(e::Emitter)
  print(e.io, "} ")
end

function emit_array_sep(e::Emitter, i=2)
  if i != 1
    print(e.io, ",")
  end
end


text:
function FGradDiv2VecDSS!(grad1CG,grad2CG,v1CG,v2CG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;

vCon = Param.CacheC1
DvCon = Param.CacheC2
vC1 = Param.CacheC3
D1vC1 = vCon
D2vC1 = DvCon
grad1TempCG = Param.CacheC3
grad2TempCG = Param.CacheC4
grad1 = Param.Cache1
grad2 = Param.Cache2
JC = Param.cache.JC

vCon .= v1CG.*Param.dXdxIC11 .+ v2CG.*Param.dXdxIC12;
mul!(reshape(vC1,OP,OP*NF*nz),CG.DS,reshape(vCon,OP,OP*nz*NF))
vCon .= v1CG.*Param.dXdxIC21 .+ v2CG.*Param.dXdxIC22;
mul!(reshape(PermutedDimsArray(DvCon,(2,1,3,4)),OP,OP*NF*nz),CG.DS,reshape(PermutedDimsArray(vCon,(2,1,3,4)),OP,OP*nz*NF))
vC1 .= vC1 .+ DvCon

mul!(reshape(D1vC1,OP,OP*NF*nz),CG.DW,reshape(vC1,OP,OP*nz*NF))
mul!(reshape(PermutedDimsArray(D2vC1,(2,1,3,4)),OP,OP*NF*nz),CG.DW,reshape(PermutedDimsArray(vC1,(2,1,3,4)),OP,OP*nz*NF))

grad1TempCG .= (Param.dXdxIC11.*D1vC1 .+ Param.dXdxIC21.*D2vC1)./ JC;
grad2TempCG .= (Param.dXdxIC12.*D1vC1 .+ Param.dXdxIC22.*D2vC1)./ JC;

grad1 .= 0
grad2 .= 0
for iM=1:size(CG.FaceGlob,1)
  i = CG.FaceGlob[iM].Ind
  arr = reshape(CG.Glob[:,i,:],OP*OP*size(i,1))
  grad1[arr,:,:] .= grad1[arr,:,:] .+
    reshape(grad1TempCG[:,:,i,:]
    ,OP*OP*size(i,1),nz);
  grad2[arr,:,:] .= grad2[arr,:,:] .+
    reshape(grad2TempCG[:,:,i,:]
    ,OP*OP*size(i,1),nz);
end

grad1.=grad1./CG.M;
grad2.=grad2./CG.M;
grad1CG .= reshape(grad1[reshape(CG.Glob,OP*OP*NF,1),:]
  ,OP,OP,NF,nz);
grad2CG .= reshape(grad2[reshape(CG.Glob,OP*OP*NF,1),:]
  ,OP,OP,NF,nz);

end

function FGradDiv2VecDSS(v1CG,v2CG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
dXdxIC = Param.cache.dXdxIC
JC = Param.cache.JC
vC1=reshape(
  CG.DS*reshape(v1CG.*dXdxIC[:,:,:,:,1,1] +
  v2CG.*dXdxIC[:,:,:,:,1,2]
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz) +
  permute(
  reshape(
  CG.DS*reshape(
  permute(v1CG.*dXdxIC[:,:,:,:,2,1] +
  v2CG.*dXdxIC[:,:,:,:,2,2]
  ,[2 1 3 4])
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz)
  ,[2 1 3 4]);

D1cCG=reshape(
  CG.DW*reshape(vC1,OP,OP*NF*nz)
  ,OP,OP,NF,nz);
D2cCG=permute(reshape(
  CG.DW*reshape(
  permute(vC1
  ,[2 1 3 4])
  ,OP,OP*NF*nz)
  ,OP,OP,NF,nz)
  ,[2 1 3 4]);

gradCG=zeros(OP,OP,NF,nz,2);
gradCG[:,:,:,:,1]=(dXdxIC[:,:,:,:,1,1].*D1cCG +
  dXdxIC[:,:,:,:,2,1].*D2cCG)./JC;
gradCG[:,:,:,:,2]=(dXdxIC[:,:,:,:,1,2].*D1cCG +
  dXdxIC[:,:,:,:,2,2].*D2cCG)./JC;

grad=zeros(CG.NumG,nz,2);
for iM=1:size(CG.FaceGlob,1)
  i = CG.FaceGlob[iM].Ind
  arr = reshape(CG.Glob[:,i,:],OP*OP*size(i,1))
  grad[arr,:,:] = grad[arr,:,:] .+
    reshape(gradCG[:,:,i,:,:],OP*OP*size(i,1),nz,2);
end
grad[:,:,1]=grad[:,:,1]./CG.M;
grad[:,:,2]=grad[:,:,2]./CG.M;
grad1CG=reshape(grad[reshape(CG.Glob,OP*OP*NF,1),:,1]
  ,OP,OP,NF,nz);
grad2CG=reshape(grad[reshape(CG.Glob,OP*OP*NF,1),:,2]
  ,OP,OP,NF,nz);
return (grad1CG,grad2CG)
end