text:
#=
    ---------------------------------------------------------------------------
    Implements a few collections that the main mcmc sampler from the file
    `mcmc.jl` uses to keep track of the current state and history of where it
    was. The following structures are defined:
    - AccptTracker   : for tracking historical acceptance rate
    - ParamHistory   : for storing the parameter chain
    - ActionTracker  : keeps track of what to do on a given iteration
    - Workspace      : main workspace for `mcmc` function from `mcmc.jl`
    - ParamUpdtDefn  : defines a single parameter update step
    - GibbsDefn      : defines an entire gibbs sweep of parameter updates
    ---------------------------------------------------------------------------
=#

import Base: last, getindex, length, display, eltype

"""
    AccptTracker

Tracks historical acceptance rate of the chain
"""
mutable struct AccptTracker
    accpt_imp::Int64            # a number of accepted imputation samples
    prop_imp::Int64             # a number of proposed imputation samples
    accpt_updt::Vector{Int64}   # a number of accepted parameter updates (one per parameter)
    prop_updt::Vector{Int64}    # a number of proposed parameter updates (one per parameter)
    updt_len::Int64             # total number of parameter update steps

    """
        AccptTracker(setup::MCMCSetup)

    Initialise tracker of acceptance rate according to a `setup` of the mcmc
    sampler
    """
    function AccptTracker(setup::MCMCSetup)
        updt_len = length(setup.updt_coord)
        accpt_imp = 0
        prop_imp = 0
        accpt_updt = [0 for i in 1:updt_len]
        prop_updt = [0 for i in 1:updt_len]
        new(accpt_imp, prop_imp, accpt_updt, prop_updt, updt_len)
    end
end

"""
    update!(at::AccptTracker, ::ParamUpdate, i, accepted::Bool)

Update acceptance tracker: increment parameter update info on the `i`th
coordinate with 1 proposed and---if `accepted` is true---then also 1 accepted
sample
"""
function update!(at::AccptTracker, ::ParamUpdate, i, accepted::Bool)
    at.prop_updt[i] += 1
    at.accpt_updt[i] += 1*accepted
end

"""
    update!(at::AccptTracker, ::Imputation, accepted::Bool)

Update acceptance tracker: increment imputation info with 1 proposed and---if
`accepted` is true---then also 1 accepted sample
"""
function update!(at::AccptTracker, ::Imputation, accepted::Bool)
    at.prop_imp += 1
    at.accpt_imp += 1*accepted
end

"""
    accpt_rate(at::AccptTracker, ::ParamUpdate)

Return current acceptance rate for parameter updates
"""
accpt_rate(at::AccptTracker, ::ParamUpdate) = at.accpt_updt./at.prop_updt

"""
    accpt_rate(at::AccptTracker, ::Imputation)

Return current acceptance rate for imputation
"""
accpt_rate(at::AccptTracker, ::Imputation) = at.accpt_imp/at.prop_imp

"""
    display(at::AccptTracker)

Show the acceptance rates
"""
function display(at::AccptTracker)
    print("Imputation acceptance rate: ", accpt_rate(at, Imputation()),
          ".\nParameter update acceptance rate: ",
          accpt_rate(at, ParamUpdate()), ".\n")
end

"""
    ParamHistory{T}

Stores information about parameter history
"""
mutable struct ParamHistory{T}
    θ_chain::Vector{T}  # parameter history
    counter::Int64      # index of the most recently accepted parameter vector

    """
        ParamHistory(setup::MCMCSetup)

    Initialise tracker of parameter history according to a `setup` of the mcmc
    sampler
    """
    function ParamHistory(setup::MCMCSetup)
        N, n = setup.num_mcmc_steps, setup.warm_up
        updt_len = length(setup.updt_coord)

        θ = params(setup.P˟)
        T = typeof(θ)
        θ_chain = Vector{T}(undef, (N-n)*updt_len+1)
        θ_chain[1] = copy(θ)
        new{T}(θ_chain, 1)
    end
end

"""
    update!(ph::ParamHistory, θ)

Update parameter history with a new accepted sample θ
"""
function update!(ph::ParamHistory, θ)
    ph.counter += 1
    ph.θ_chain[ph.counter] = copy(θ)
end

"""
    last(ph::ParamHistory)

Return a copy of the most recently accepted paramter vector
"""
last(ph::ParamHistory) = copy(ph.θ_chain[ph.counter])

"""
    ActionTracker

Keeps track of what to do on a given iteration
"""
struct ActionTracker{T,S,R,U}
    save_iter::T    # Save the path every ... iteration
    verb_iter::S   # Print progress message to console every ... iteration
    warm_up::R      # Number of steps of the chain in which no param update is made
    readjust::U
    param_updt::Bool    # Flag for whether to update parameters at all

    """
        ActionTracker(setup::MCMCSetup)

    Initialise tracker of what to do on a given iteration according to a `setup`
    of the mcmc sampler
    """
    function ActionTracker(setup::MCMCSetup)
        i1, i2, i3 = setup.save_iter, setup.verb_iter, setup.warm_up
        i4 = setup.pCN_readjust_param.step
        s1, s2, s3, s4 = typeof(i1), typeof(i2), typeof(i3), typeof(i4)
        @assert (s1 <: Number) && (s2 <: Number) && (s3 <: Number) && (s4 <: Number)
        new{s1,s2,s3,s4}(i1, i2, i3, i4, setup.param_updt)
    end
end

"""
    act(::SavePath, at::ActionTracker, i)

Determine whether to save path on a given iteration, indexed `i`
"""
function act(::SavePath, at::ActionTracker, i)
    (i > at.warm_up) && (i % at.save_iter == 0)
end

"""
    act(::Verbose, at::ActionTracker, i)

Determine whether to print out information to a console on a given iteration,
indexed `i`
"""
act(::Verbose, at::ActionTracker, i) = (i % at.verb_iter == 0)

"""
    act(::ParamUpdate, at::ActionTracker, i)

Determine whether to update parameters on a given iteration, indexed `i`
"""
act(::ParamUpdate, at::ActionTracker, i) = at.param_updt && (i > at.warm_up)

act(::Readjust, at::ActionTracker, i) = (i % at.readjust == 0) && (i > at.warm_up)

mutable struct SingleElem{T} val::T end

set!(x::SingleElem{T}, y::T) where T = (x.val = y)

const RhoInfoType = NamedTuple{(:step, :scale, :minδ, :maxρ, :trgt, :offset),
                               Tuple{Int64, Float64, Float64, Float64, Float64, Int64}}

"""
    Workspace{ObsScheme,S,TX,TW,R,ST}

The main container of the `mcmc` function from `mcmc.jl` in which most data
pertinent to sampling is stored
"""
struct Workspace{ObsScheme,B,ST,S,TX,TW,R,TP,TZ}# ,Q, where Q = eltype(result)
    Wnr::Wiener{S}         # Wiener, driving law
    XXᵒ::Vector{TX}        # Diffusion proposal paths
    XX::Vector{TX}         # Accepted diffusion paths
    WWᵒ::Vector{TW}        # Driving noise of proposal
    WW::Vector{TW}         # Driving noise of the accepted paths
    Pᵒ::Vector{R}          # Guided proposals parameterised by proposal param
    P::Vector{R}           # Guided proposals parameterised by accepted param
    fpt::Vector            # Additional information about first passage times
    ρ::Vector{Vector{Float64}}      # Memory parameter of the precond Crank-Nicolson scheme
    recompute_ODEs::Vector{Bool}    # Info on whether to recompute H,Hν,c after resp. param updt
    accpt_tracker::AccptTracker     # Object for tracking acceptance rate
    θ_chain::ParamHistory           # Object for tracking parameter history
    action_tracker::ActionTracker   # Object for tracking steps to perform on a given iteration
    skip_for_save::Int64            # Thining parameter for saving path
    paths::Vector                   # Storage with historical, accepted paths
    time::Vector{Float64}           # Storage with time axis
    blocking::B
    blidx::Int64
    x0_prior::TP
    z::SingleElem{TZ}
    pCN_readjust_param::RhoInfoType
    #result::Vector{Q} #TODO come back to later
    #resultᵒ::Vector{Q} #TODO come back to later

    """
        Workspace(setup::MCMCSetup{ObsScheme})

    Initialise workspace of the mcmc sampler according to a `setup` variable
    """
    function Workspace(setup::MCMCSetup{ObsScheme}) where ObsScheme
        # just to make sure that nothing gets messed up if the user decides
        # to later modify `setup` use deepcopies
        x0_prior, Wnr = deepcopy(setup.x0_prior), deepcopy(setup.Wnr)
        XX, WW = deepcopy(setup.XX), deepcopy(setup.WW)
        P, fpt = deepcopy(setup.P), deepcopy(setup.fpt)
        updt_coord = deepcopy(setup.updt_coord)
        pCN_readjust = deepcopy(setup.pCN_readjust_param)

        # forcedSolve defines type by the starting point, make sure it matches
        x0_guess = eltype(eltype(XX))(setup.x0_guess)
        TW, TX, S, R = eltype(WW), eltype(XX), valtype(Wnr), eltype(P)
        ST, TP = typeof(setup.solver), typeof(x0_prior)

        m = length(P)

        y = copy(x0_guess)
        for i in 1:m
            WW[i] = Bridge.samplepath(P[i].tt, zero(S))
            sample!(WW[i], Wnr)
            WW[i], XX[i] = forcedSolve(Euler(), y, WW[i], P[i])    # this will enforce adherence to domain
            while !checkFpt(ObsScheme(), XX[i], fpt[i])
                sample!(WW[i], Wnr)
                forcedSolve!(Euler(), XX[i], y, WW[i], P[i])    # this will enforce adherence to domain
            end
            y = XX[i].yy[end]
        end
        y = x0_guess
        ll = logpdf(x0_prior, y)
        ll += path_log_likhd(ObsScheme(), XX, P, 1:m, fpt, skipFPT=true)
        ll += lobslikelihood(P[1], y)

        XXᵒ = deepcopy(XX)
        WWᵒ = deepcopy(WW)
        Pᵒ = deepcopy(P)

        # compute the white noise that generates x0_guess under the initial posterior
        z = inv_start_pt(y, x0_prior, P[1])
        TZ = typeof(z)
        z = SingleElem{TZ}(z)

        #TODO come back to gradient initialisation
        skip = setup.skip_for_save
        _time = collect(Iterators.flatten(p.tt[1:skip:end-1] for p in P))

        θ_history = ParamHistory(setup)

        blocking = set_blocking(setup.blocking, setup.blocking_params, P)
        ρ = prepare_mem_param(setup.ρ, blocking)
        display(blocking)
        B = typeof(blocking)

        (workspace = new{ObsScheme,B,ST,S,TX,TW,R,TP,TZ}(Wnr, XXᵒ, XX, WWᵒ, WW,
                                                         Pᵒ, P, fpt, ρ,
                                                         check_if_recompute_ODEs(setup),
                                                         AccptTracker(setup),
                                                         θ_history,
                                                         ActionTracker(setup),
                                                         skip, [], _time,
                                                         blocking, 1, x0_prior,
                                                         z, pCN_readjust),
         ll = ll, θ = last(θ_history))
    end

    # NOTE this constructor is no longer in use, can be removed
    """
        Workspace(ws::Workspace{ObsScheme,S,TX,TW,R}, new_ρ::Float64)

    Copy constructor of `workspace`. Keeps everything the same as passed `ws`
    with the exception of new memory parameter for the preconditioned
    Crank-Nicolson scheme, which is changed to `new_ρ`.
    """
    function Workspace(ws::Workspace{ObsScheme,B,ST,S,TX,TW,R,TP,TZ}, new_ρ::Vector{Vector{Float64}}
                       ) where {ObsScheme,B,ST,S,TX,TW,R,TP,TZ}
        new{ObsScheme,B,ST,S,TX,TW,R,TP,TZ}(ws.Wnr, ws.XXᵒ, ws.XX, ws.WWᵒ,
                                            ws.WW, ws.Pᵒ, ws.P, ws.fpt, new_ρ,
                                            ws.recompute_ODEs, ws.accpt_tracker,
                                            ws.θ_chain, ws.action_tracker,
                                            ws.skip_for_save, ws.paths, ws.time,
                                            ws.blocking, ws.blidx, ws.x0_prior,
                                            ws.z, ws.pCN_readjust_param)
    end

    function Workspace(ws::Workspace{ObsScheme,B,ST,S,TX,TW,R̃,TP,TZ},
                       P::Vector{R}, Pᵒ::Vector{R}, idx
                       ) where {ObsScheme,B,ST,S,TX,TW,R̃,R,TP,TZ}
        new{ObsScheme,B,ST,S,TX,TW,R,TP,TZ}(ws.Wnr, ws.XXᵒ, ws.XX, ws.WWᵒ,
                                            ws.WW, Pᵒ, P, ws.fpt, ws.ρ,
                                            ws.recompute_ODEs, ws.accpt_tracker,
                                            ws.θ_chain, ws.action_tracker,
                                            ws.skip_for_save, ws.paths, ws.time,
                                            ws.blocking, idx, ws.x0_prior, ws.z,
                                            ws.pCN_readjust_param)
    end
end

eltype(::SamplePath{T}) where T = T
eltype(::Type{SamplePath{T}}) where T = T
solver_type(::Workspace{O,B,ST}) where {O,B,ST} = ST

next_set_of_blocks(ws::Workspace{O,NoBlocking}) where O = ws

"""
    next(𝔅::ChequeredBlocking, XX, θ)

Switch the set of blocks that are being updated. `XX` is the most recently
sampled (accepted) path. `θ` can be used to change parametrisation.
"""
function next_set_of_blocks(ws::Workspace{O,<:ChequeredBlocking}) where O
    XX, P, Pᵒ, 𝔅 = ws.XX, ws.P, ws.Pᵒ, ws.blocking
    idx = (ws.blidx % 2) + 1
    θ = params(P[1].Target)

    vs = find_end_pts(𝔅, XX, idx)
    Ls = 𝔅.Ls[idx]
    Σs = 𝔅.Σs[idx]
    ch_pts = 𝔅.change_pts[idx]
    aux_flags = 𝔅.aux_flags[idx]

    P_new = [GuidPropBridge(P[i], Ls[i], vs[i], Σs[i], ch_pts[i], θ, aux_flags[i])
             for (i,_) in enumerate(P)]
    Pᵒ_new = [GuidPropBridge(Pᵒ[i], Ls[i], vs[i], Σs[i], ch_pts[i], θ, aux_flags[i])
              for (i,_) in enumerate(Pᵒ)]
    Workspace(ws, P_new, Pᵒ_new, idx)
end

prepare_mem_param(ρ::Number, ::NoBlocking) = [[ρ]]

function prepare_mem_param(ρ::Number, blocking::ChequeredBlocking)
    [[ρ for _ in block_seq] for block_seq in blocking.accpt_tracker.accpt]
end

"""
    act(action, ws::Workspace, i)

Determine whether to perform `action` on a given iteration, indexed `i`
"""
act(action, ws::Workspace, i) = act(action, ws.action_tracker, i)

#NOTE temporary
act(action::Readjust, ws::Workspace, i) = (typeof(ws.blocking) <:ChequeredBlocking) && act(action, ws.action_tracker, i)

"""
    savePath!(ws, wsXX, bXX)

Save the entire path spanning all segments in `XX`. Only 1 in every `ws.skip`
points is saved to reduce storage space. To-be-saved `XX` is set to `wsXX` or
`bXX` depending on whether blocking is used.
"""
function save_path!(ws)
    skip = ws.skip_for_save
    push!(ws.paths, collect(Iterators.flatten(ws.XX[i].yy[1:skip:end-1]
                                               for i in 1:length(ws.XX))))
end

sigmoid(x, a=1.0) = 1.0 / (1.0 + exp(-a*x))
logit(x, a=1.0) = (log(x) - log(1-x))/a

function readjust_pCN!(ws, mcmc_iter)
    at = ws.blocking.short_term_accpt_tracker
    p = ws.pCN_readjust_param
    δ = max(p.minδ, p.scale/sqrt(max(1.0, mcmc_iter/p.step-p.offset)))
    accpt_rates = acceptance(at)
    for (i, a_i) in enumerate(accpt_rates)
        for (j, a_ij) in enumerate(a_i)
            ws.ρ[i][j] = min(sigmoid(logit(ws.ρ[i][j]) - (2*(a_ij > p.trgt)-1)*δ), p.maxρ)
        end
    end
    display_acceptance_rate(ws.blocking, true)
    print_pCN(ws)
    reset!(at)
end

#NOTE _print_info defined in `blocking_schedule`
function print_pCN(ws)
    print("\nρ parameter:\n----------------------\n")
    _print_info(ws.ρ[1])
    _print_info(ws.ρ[2])
end

"""
    ParamUpdtDefn

For a given, single parameter update step defines transition kernels, priors,
which coordinates are updated etc.
"""
struct ParamUpdtDefn{R,S,T,U}
    updt_type::R         # The type of update (Metropolis-Hastings/conjugate etc)
    updt_coord::S        # Which coordinates to update
    t_kernel::T          # Transition kernel for a given parameter update
    priors::U            # Prior over updated parameters
    recompute_ODEs::Bool # Whether given param updt calls for recomputing H,Hν,c

    """
        ParamUpdtDefn(updt_type::R, updt_coord::S, t_kernel::T, priors::U,
                      recompute_ODEs::Bool)

    Initialisation of the complete definition of the parameter update step
    """
    function ParamUpdtDefn(updt_type::R, updt_coord::S, t_kernel::T, priors::U,
                           recompute_ODEs::Bool
                           ) where {R<:ParamUpdateType,S,T,U}
        new{R,S,T,U}(updt_type, updt_coord, t_kernel, priors, recompute_ODEs)
    end
end

"""
    GibbsDefn

Definition of the entire Gibbs sweep.
"""
struct GibbsDefn{N}
    updates::NTuple{N,ParamUpdtDefn}

    """
        GibbsDefn(setup)

    Initialises Gibbs sweep according to the `setup`
    """
    function GibbsDefn(setup)
        recompute_ODEs = check_if_recompute_ODEs(setup)

        updates = [ParamUpdtDefn(ut, uc, tk, pr, ro) for (ut, uc, tk, pr, ro)
                   in zip(setup.updt_type, setup.updt_coord, setup.t_kernel,
                          setup.priors, recompute_ODEs)]
        new{length(updates)}(Tuple(updates))
    end
end

"""
    getindex(g::GibbsDefn, i::Int)

Return `i`th definition of parameter update
"""
getindex(g::GibbsDefn, i::Int) = g.updates[i]

"""
    length(g::GibbsDefn{N})

Return the total number of parameter updates in a single Gibbs sweep
"""
length(g::GibbsDefn{N}) where N = N

"""
    init_adaptation!(adpt::Adaptation{Val{false}}, ws::Workspace)

Nothing to do when no adaptation needs to be done
"""
init_adaptation!(adpt::Adaptation{Val{false}}, ws::Workspace) = nothing

"""
    init_adaptation!(adpt::Adaptation{Val{true}}, ws::Workspace)

Resize internal container with paths in `adpt` to match the length of imputed
paths
"""
function init_adaptation!(adpt::Adaptation{Val{true}}, ws::Workspace)
    m = length(ws.XX)
    resize!(adpt, m, [length(ws.XX[i]) for i in 1:m])
end

"""
    update!(adpt::Adaptation{Val{false}}, ws::Workspace{ObsScheme}, yPr, i, ll,
            solver::ODESolverType)

Nothing to be done for no adaptation
"""
function update!(adpt::Adaptation{Val{false}}, ws::Workspace{ObsScheme}, i,
                 ll) where ObsScheme
    adpt, ll
end

"""
    update!(adpt::Adaptation{Val{false}}, ws::Workspace{ObsScheme}, yPr, i, ll,
            solver::ODESolverType)

Update the proposal law according to the adaptive scheme and the recently saved
history of the imputed paths
"""
function update!(adpt::Adaptation{Val{true}}, ws::Workspace{ObsScheme,B,ST},
                 i, ll) where {ObsScheme,B,ST}
    if i % adpt.skip == 0
        if adpt.N[2] == adpt.sizes[adpt.N[1]]
            X_bar = mean_trajectory(adpt)
            m = length(ws.P)
            for j in 1:m
                Pt = recentre(ws.P[j].Pt, ws.XX[j].tt, X_bar[j])
                update_λ!(Pt, adpt.λs[adpt.N[1]])
                ws.P[j] = GuidPropBridge(ws.P[j], Pt)

                Ptᵒ = recentre(ws.Pᵒ[j].Pt, ws.XX[j].tt, X_bar[j])
                update_λ!(Ptᵒ, adpt.λs[adpt.N[1]])
                ws.Pᵒ[j] = GuidPropBridge(ws.Pᵒ[j], Ptᵒ)
            end

            solve_back_rec!(NoBlocking(), ws, ws.P)
            #solveBackRec!(NoBlocking(), ws.Pᵒ, ST())
            y = ws.XX[1].yy[1]
            z = inv_start_pt(y, ws.x0_prior, ws.P[1])
            set!(ws.z, z)

            for j in 1:m
                inv_solve!(Euler(), ws.XX[j], ws.WW[j], ws.P[j])
            end
            ll = logpdf(ws.x0_prior, y)
            ll += path_log_likhd(ObsScheme(), ws.XX, ws.P, 1:m, ws.fpt)
            ll += lobslikelihood(ws.P[1], y)
            adpt.N[2] = 1
            adpt.N[1] += 1
        else
            adpt.N[2] += 1
        end
    end
    adpt, ll
end


text:
using IncrementalInference
# using KernelDensityEstimate, Gadfly # for vstack

fg = emptyFactorGraph()

N=100

doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]

v1 = addNode!(fg,:x1,ContinuousScalar,N=N)
f1  = addFactor!(fg,[v1], Obsv2(doors, cov', [1.0]))

# tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

v3=addNode!(fg,:x3,ContinuousScalar, N=N)
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))

v4=addNode!(fg,:x4,ContinuousScalar, N=N)
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

v5=addNode!(fg,:x5,ContinuousScalar, N=N)
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

l1=addNode!(fg, :l1, ContinuousScalar, N=N)
addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))

l2=addNode!(fg, :l2, ContinuousScalar, N=N)
addFactor!(fg, [v4,l2], Ranged([16.0],[0.5],[1.0]))

l3=addNode!(fg, :l3, ContinuousScalar, N=N)
addFactor!(fg, [v5,l3], Ranged([16.0],[0.5],[1.0]))

# writeGraphPdf(fg);

tree = prepBatchTree!(fg,drawpdf=true);
# run(`evince bt.pdf`)

# do belief propagation inference over tree once
# inferOverTreeR!(fg, tree, N=100)

ext: jl


text:
# rotation matrix (2D)
function rotation_matrix(x)
	[cos(x) -sin(x); sin(x) cos(x)]
end

# signed distance for a box
function sd_box(p, dim)
	q = abs.(p) - dim
	norm(max.(q, 0.0)) + min(maximum(q), 0.0)
end

function sd_2d_box(p, pose, dim, rnd)
	x, y, θ = pose
	R = rotation_matrix(-θ)
	p_rot = R * (p - pose[1:2])

	return sd_box(p_rot, dim) - rnd
end

# box parameters
dim = [0.1, 0.1]
rnd = 0.0
dim_rnd = dim .- rnd

# setup
# p = [0.5 * sqrt(2.0), 0.5 * sqrt(2.0)]
p = [0.0999, 0.0]
px = 0.0
py = 0.0
θ = 0.0 * π
pose = [px, py, θ]

# problem
@show sdf = sd_2d_box(p, pose, dim_rnd, rnd)

sd_pose(x) = sd_2d_box(p, x, dim_rnd, rnd)
@show Npose = ForwardDiff.gradient(sd_pose, pose)
ForwardDiff.gradient(x->norm(x), zeros(3))

sd_p(x) = sd_2d_box(x, pose, dim_rnd, rnd)
@show Np = ForwardDiff.gradient(sd_p, p)

n_dir = Np[1:2] ./ norm(Np[1:2])
t_dir = rotation_matrix(0.5 * π) * n_dir

r = p - pose[1:2]
m = cross([r; 0.0], [t_dir; 0.0])[3]

@show Ppose = [t_dir[1], t_dir[2], cross([p; 0.0] - [pose[1:2]; 0.0], [t_dir; 0.0])[3]]

@show Pp = -t_dir

max_issues_repo_name: thowell/DirectMotionPlanning


text:
module ClimateTools

# External modules
using Reexport
@reexport using ClimateBase
using NetCDF
@reexport using NCDatasets
using Shapefile
using AxisArrays
using NaNMath
const axes = Base.axes
using ArgCheck
using DataFrames
using Interpolations
using ProgressMeter
using Polynomials
using IterTools
using Statistics
using Random
using Dates
using GeoStats
using InverseDistanceWeighting
using Extremes
using Distances
import Base.vcat
import Base.getindex
import Base.show
import Base.size
import Base.setindex!
import Base.similar
import Base.write
import Statistics.minimum
import Statistics.maximum
import Statistics.std
import Statistics.var
import Statistics.mean
import Base: +
import Base: -
import Base: *
import Base: /
import Base.findmax

# Included files
include("functions.jl")
include("indices.jl")
include("indicators.jl")
include("extract.jl")
include("interface.jl")

include("biascorrect.jl")
include("export.jl")
include("time.jl")
include("spatial.jl")
include("analysis.jl")

export inpoly, inpolygrid, meshgrid, inpolyvec, ndgrid
export findmax, findmin
export frostdays, summerdays, icingdays, tropicalnights
export daysabove10 #, daysbelow0, degdaysabove, degdaysbelow
export customthresover, customthresunder, annualmax, annualmin
export annualmean, annualsum, prcp1
export drought_dc
export ensemble_mean, ensemble_std, ensemble_max, ensemble_min
export load, load2D
export regrid, applymask
export shapefile_coords, shapefile_coords_poly
export resample, spatialsubset
export qqmap, qqmaptf
export biascorrect_extremes
export permute_west_east
export getdim_lat, getdim_lon, isdefined, extractpoly
export polyfit, polyval
export @isdefined
export merge, vaporpressure, approx_surfacepressure
export wbgt, diurnaltemperature, meantemperature
export minimum, maximum, std, var, mean
export daymean, daysum
export monthmean, monthsum, temporalmean
export yearmonthdayhour
export write, findmindist

end #module


text:
using CUDA
using Oceananigans.Architectures: AbstractCPUArchitecture, AbstractGPUArchitecture

import Base: zeros

zeros(FT, ::AbstractCPUArchitecture, N...) = zeros(FT, N...)
zeros(FT, ::AbstractGPUArchitecture, N...) = CUDA.zeros(FT, N...)

zeros(arch::AbstractArchitecture, grid, N...) = zeros(eltype(grid), arch, N...)


text:
module ArrayToolsTests

using Test, Random
using ArrayTools, ArrayTools.AnnotatedArrays, ArrayTools.PseudoArrays
using .AnnotatedArrays:
    propertyname,
    propertynames,
    DynamicallyAnnotatedArray,
    StaticallyAnnotatedArray
using Base: @propagate_inbounds
import Base: getindex, setindex!, checkbounds

function samevalues(A::AbstractArray, B::AbstractArray)
    @assert has_standard_indexing(A, B)
    @assert size(A) == size(B)
    for i ∈ eachindex(A, B)
        A[i] != B[i] && return false
    end
    return true
end

function maxabsdif(A::AbstractArray, B::AbstractArray)
    @assert has_standard_indexing(A, B)
    @assert size(A) == size(B)
    T = promote_type(eltype(A), eltype(B))
    res = zero(T)
    for i ∈ eachindex(A, B)
        res = max(res, abs(A[i] - B[i]))
    end
    return res
end

slice(A::AbstractArray{T,N}, i::Integer) where {T,N} =
    A[colons(Val(N-1))..., i]

generate(::Type{T}, dims::Integer...) where {T} = generate(T, dims)
generate(::Type{T}, dims::NTuple{N,Integer}) where {T,N} =
    generate!(Array{T,N}(undef, dims))
generate!(A::AbstractArray{T,N}) where {T,N} = begin
    k = 0
    @inbounds for i in eachindex(A)
        k += 1
        A[i] = k
    end
    return A
end

const T = Float32
const dims = (3, 4, 5, 6)
const N = length(dims)
A = generate(T, dims)
Vf = view(A, :, :, :, 3:4)   # a flat sub-array
Va = view(A, 2:3, :, :, 3:4) # a non-flat sub-array
V = view(A, :, 2:3, :, :)
S = 1:2:70  # StepRange
U = 3:50    # UnitRange
atol = 1e-6

@testset "Miscellaneous " begin
    # Promotion of element types.
    for (T1,T2,T3) in ((Float32,Float64,Int), (Int16,Int32,Int64))
        @test promote_eltype(zeros(T1,1)) == promote_type(T1)
        @test promote_eltype(AbstractVector{T2}) == promote_type(T2)
        @test promote_eltype(zeros(T2,4), zeros(T3,2,3)) == promote_type(T2,T3)
        @test promote_eltype(zeros(T2,2), Array{T3}) == promote_type(T2,T3)
        @test promote_eltype(zeros(T1,3), zeros(T2,5), AbstractVector{T3}) ==
            promote_type(T1,T2,T3)
        @test promote_eltype(zeros(T1,3), zeros(T2,5), ones(T3,2)) ==
            promote_type(T1,T2,T3)
        @test promote_eltype(DenseMatrix{T1}, zeros(T2,5), AbstractVector{T3}) ==
            promote_type(T1,T2,T3)
    end
    # Conversions.
    @test to_type(Int32, 3) === Int32(3)
    @test to_type(Float32, 3) === Float32(3)
    @test to_type(Vector{Float32}, [1,2,3]) == Float32[1,2,3]
    @test isa(to_type(Vector{Float32}, [1,2,3]), Vector{Float32})
    # Dimensions.
    @test to_size(()) === ()
    @test to_size(5) === (5,)
    @test to_size(Int16(5)) === (5,)
    @test isa(to_size(Int16(5)), Tuple{Int})
    @test standard_size(A) === size(A)
    @test to_size(dims) === dims
    @test to_size(dims...) === dims
    @test to_size(UInt16(dims[1]), dims[2:end]...) === dims
    @test to_size((UInt16(dims[1]), dims[2:end]...)) === dims
    @test check_size(()) == 1
    @test check_size((1,0,2)) == 0
    @test check_size((1,2,3)) == 6
    @test check_size(0xff) == 255
    @test_throws ArgumentError check_size((1,0,-1))
    @test_throws ArgumentError check_size(-1)
    # Axis and axes.
    @test to_axis(1:3) === 1:3
    @test to_axis(UInt16(1):UInt16(3)) === 1:3
    @test to_axis(5) === Base.OneTo(5)
    @test to_axis(UInt16(5)) === Base.OneTo(5)
    @test to_axes(()) === ()
    @test to_axes((Base.OneTo(3), 0:4)) === (Base.OneTo(3), 0:4)
    @test to_axes(5) === (to_axis(5),)
    @test to_axes(2,UInt16(3),0:4) === (to_axis(2),to_axis(3),to_axis(0:4),)
    @test to_axes((2,UInt16(3),0:4)) === (to_axis(2),to_axis(3),to_axis(0:4),)
    #
    # Tests for `allof`, `anyof` and `noneof`.
    #
    W = [true, true, true]
    X = [true, false]
    Y = [false, false]
    for i in randperm(length(X)) # prevent compilation-time optimization
        @test anyof(X[i]) == X[i]
    end
    @test allof(true) == true
    @test anyof(true) == true
    @test noneof(true) == false
    @test allof(false) == false
    @test anyof(false) == false
    @test noneof(false) == true
    @test allof(W) == true
    @test anyof(W) == true
    @test noneof(W) == false
    @test allof(Tuple(W)) == allof(W)
    @test anyof(Tuple(W)) == anyof(W)
    @test noneof(Tuple(W)) == noneof(W)
    @test allof(X) == false
    @test anyof(X) == true
    @test noneof(X) == false
    @test allof(Tuple(X)) == allof(X)
    @test anyof(Tuple(X)) == anyof(X)
    @test noneof(Tuple(X)) == noneof(X)
    @test allof(Y) == false
    @test anyof(Y) == false
    @test noneof(Y) == true
    @test allof(Tuple(Y)) == allof(Y)
    @test anyof(Tuple(Y)) == anyof(Y)
    @test noneof(Tuple(Y)) == noneof(Y)
    for (x, y) in ((W,W),(X,X),(Y,Y),(W,X),(X,Y),(Y,W))
        @test anyof(x, y) == (anyof(x) || anyof(y))
        @test allof(x, y) == (allof(x) && allof(y))
        @test noneof(x, y) == (noneof(x) && noneof(y))
    end
    for x in (W, X, Y)
        @test allof(x) == (allof(minimum, x, x) && allof(maximum, x))
        @test noneof(x) == (noneof(minimum, x, x) && noneof(maximum, x))
        @test anyof(x) == (anyof(minimum, x, x) || anyof(maximum, x))
    end
    #
    # Tests for `cartesian_indices`.
    #
    @test cartesian_indices(A) === CartesianIndices(axes(A))
    @test cartesian_indices(A) === cartesian_indices(size(A))
    @test cartesian_indices(A) === cartesian_indices(axes(A))
    @test cartesian_indices(A) === cartesian_indices(cartesian_indices(A))
    I1 = CartesianIndex(1,2,3)
    I2 = CartesianIndex(5,7,9)
    @test cartesian_indices(I1,I2) ===
        CartesianIndices(([I1[k]:I2[k] for k in 1:length(I1)]...,))
    @test cartesian_indices(3:4) === CartesianIndices((3:4,))
    @test cartesian_indices(3:4,2:5) === CartesianIndices((3:4,2:5))
    @test cartesian_indices((3:4,2:5)) === CartesianIndices((3:4,2:5))
    #
    # Tests for `axis_limits`.
    #
    @test axis_limits(Base.OneTo(7)) === (1,7)
    @test axis_limits(7:16) === (7,16)
    @test axis_limits(Int16(7):Int16(1):Int16(16)) === (7,16)
    @test axis_limits(16:-1:7) === (7,16)
    @test axis_limits(7:-1:16) === (8,7)
    @test_throws ArgumentError axis_limits(7:3:16)
    #
    # Tests for `all_match`.
    #
    @test all_match(nothing, identity)
    @test all_match(1, length, π)
    @test all_match(1, length, π, 2, "e")
    #
    # Tests for `same_axes`.
    #
    let A = Array{Char}(undef, (2,3)),
        B = Array{Int}(undef, 2, 3),
        C = Array{Float32}(undef, (2,3,4))
        @test same_axes(A) === axes(A)
        @test same_axes(B) === axes(B)
        @test same_axes(A,B) == axes(A)
        @test same_axes(A,B) == axes(B)
        @test_throws DimensionMismatch same_axes(A,B,C)
    end
    #
    # Tests for `all_indices`.
    #
    B = rand(T, dims[1:2])
    C = rand(T, dims[1:end-2]..., dims[end], dims[end-1])
    X = rand(T, 6)
    @test all_indices(X) === eachindex(X)
    @test all_indices(A) === eachindex(A)
    @test all_indices(Va) === eachindex(Va)
    @test_throws DimensionMismatch all_indices(A,Va)
    @test eachindex(A,C) === eachindex(A)
    @test_throws DimensionMismatch all_indices(A,C)
    @test all_indices(A, rand(T, dims), rand(T, dims)) === eachindex(A)
    Y = rand(T, dims[1], dims[2]+2, dims[3:end]...)
    Z = view(Y, :, 2:dims[2]+1, colons(length(dims)-2)...)
    @test IndexStyle(Z) === IndexCartesian()
    @test all_indices(A, rand(T, dims), Z) === eachindex(IndexCartesian(), A)
    #
    # Other stuff.
    #
    @test reversemap(x -> x^2, dims) === reverse(map(x -> x^2, dims))
    @test_throws ArgumentError ArrayTools.throw_non_standard_indexing()
end

#
# Tests for rubber indices.
#
@testset "Rubber Indices" begin
    I1 = CartesianIndex(2)
    I2 = CartesianIndex(2,3)
    I3 = CartesianIndex(2,3,4)
    I4 = CartesianIndex(1,2,3,2)
    I5 = CartesianIndex(1,2,3,2,4)
    @test_throws ArgumentError colons(-1)
    for d ∈ 0:12
        tup = ntuple(x -> Colon(), d)
        @test colons(d) === tup
        @test colons(Val(d)) === tup
    end
    for k ∈ 1:dims[end]
        @test samevalues(slice(A, k), A[:,:,:,k])
    end
    @test A[..]           === A
    @test A[..]           == A[:,:,:,:]
    @test A[..,3]         == A[:,:,:,3]
    @test A[2,..]         == A[2,:,:,:]
    @test A[..,2:4,5]     == A[:,:,2:4,5]
    @test A[2:3,..,1,2:4] == A[2:3,:,1,2:4]
    @test A[:,2,..,2:4]   == A[:,2,:,2:4]
    @test A[:,2:3,..]     == A[:,2:3,:,:]
    @test A[..,2:3,:]     == A[:,:,2:3,:]
    @test A[I1,..]        == A[I1,:,:,:]
    @test A[I2,..]        == A[I2,:,:]
    @test A[..,I1]        == A[:,:,:,I1]
    @test A[..,I2]        == A[:,:,I2]
    @test A[I1,..,I2]     == A[I1,:,I2]
    @test A[I2,..,I1]     == A[I2,:,I1]
    @test A[I1,..,I3]     == A[I1,I3]
    @test A[1,2:end,..]   == A[1,2:end,:,:]
    @test A[1,2:end-1,..] == A[1,2:end-1,:,:]
    # The `end` keyword can only appear **before** the rubber index.
    @test_broken A[2,..,3:end] == A[2,:,:,3:end]
    A[1,..] .= 3
    @test all(isequal(3), A[1,:,:,:])
    A[1,..] = A[2,..]
    @test A[1,:,:,:] == A[2,:,:,:]
    A[1,..,2:3] .= A[2,..,3:2:5]
    @test A[1,:,:,2:3] == A[2,:,:,3:2:5]
    # Tests with a 9 dimensional array.
    siz = (3, 4, 2, 3, 4, 2, 3, 4, 2)
    n = length(siz)
    B = generate(Int, siz)
    C = copy(B)
    inds = (2:3, 4, Colon(), CartesianIndex(2),
            CartesianIndex(3), 1:2, 3, 2:3, CartesianIndex(1))
    for k in 0:n, l in 0:n-k
        J1 = inds[1:k]         # leading indices
        J2 = inds[n-l+1:n]     # trailing indices
        J = (J1..., .., J2...,)# with a rubber index
        K = (J1..., colons(n - length(J1) - length(J2))..., J2...,)
        R = B[K...]            # extract region before change
        @test R == B[J...]     # extract values with `getindex`
        B[J...] .*= -1         # change values with `dotview`
        @test R == -B[J...]    # extract values with `getindex`
        B[J...] = R            # restore values with `setindex!`
        @test B == C
    end
end

@testset "Storage       " begin
    B = to_flat_array(Float32, A)
    C = to_flat_array(Float32, Va)
    @test StorageType() === AnyStorage()
    @test StorageType("a") === AnyStorage()
    @test StorageType(A) === FlatStorage()
    @test StorageType(Vf) === FlatStorage()
    @test StorageType(Va) === AnyStorage()
    @test is_flat_array() == false
    @test is_flat_array("a") == false
    @test is_flat_array(S) == false
    @test is_flat_array(A) == true
    @test is_flat_array(Vf) == true
    @test is_flat_array(Va) == false
    @test is_flat_array(A,B,C) == true
    @test is_flat_array(A,B,C,Vf) == true
    @test is_flat_array(A,B,C,Va) == false
    @test to_flat_array(A) === A
    @test to_flat_array(eltype(A), A) === A
    @test pointer(A) != pointer(to_flat_array(Va))
    @test samevalues(S, to_flat_array(S))
    @test eltype(S) === eltype(to_flat_array(S))
    @test maxabsdif(A, B) ≤ atol
    @test maxabsdif(Va, C) ≤ atol
    Da = to_flat_array(Va)
    @test is_flat_array(Da) == true
    @test Da == Va
    @test Da !== Va
    Df = to_flat_array(Vf)
    @test is_flat_array(Df) == true
    @test Df == Vf
    @test Df === Vf
    for n in 1:5
        K = rand(Float64, ntuple(x -> 3, n))
        L = view(K, 2:3, colons(n-1)...)
        M = view(K, colons(n-1)..., 2:3)
        @test StorageType(K) == FlatStorage()
        @test StorageType(L) == (n == 1 ? FlatStorage() : AnyStorage())
        @test StorageType(M) == FlatStorage()
    end
    L = ((nothing,               false),
         ("a",                   false),
         ((),                    false),
         (1,                     false),
         (A,                     true ),
         (view(A, :, :, :, 2:3), true ),
         (view(A, :, :, 2:2, :), false),
         (view(A, :, :, 2, :),   false),
         (view(A, :, :, 2:2, 3), true ),
         (view(A, :, :, 2, 3),   true ),
         (view(A, :, :, 2, 3:3), false))
    for i in randperm(length(L)) # prevent compilation-time optimization
        x, b = L[i]
        @test is_flat_array(x) == b
    end
    @test is_flat_array(A, view(A, :, :, 2, 3), view(A, :, :, 2:2, 3)) == true
    @test is_flat_array(A, view(A, :, :, 2:2, :), view(A, :, :, 2:2, 3)) == false
end

@testset "Indexing      " begin
    # Prior to version 1.1, Base.has_offset_axes does not work for tuples.
    for Q in (A,Va,S,101,Colon(),(),(1,),("e",2,))
        if isa(Q, Tuple) && VERSION < v"1.1"
            @test has_standard_indexing(Q) == true
        else
            @test has_standard_indexing(Q) == !Base.has_offset_axes(Q)
        end
    end
    B = to_fast_array(Float32, A)
    C = to_fast_array(Float32, Va)
    @test has_standard_indexing(A,Va) == (has_standard_indexing(A) &&
                                          has_standard_indexing(Va))
    @test IndexingType(A) === FastIndexing()
    @test IndexingType(Va) === AnyIndexing()
    @test IndexingType("a") === AnyIndexing()
    @test is_fast_array() == false
    @test is_fast_array(S) == true
    @test samevalues(Va, to_fast_array(Va))
    @test samevalues(A, to_fast_array(A))
    @test is_fast_array(A) == true
    @test is_fast_array(A,B,C) == true
    @test is_fast_array(Va) == false
    @test is_fast_array(A,B,C,Va) == false
    @test is_fast_array(to_fast_array(Va)) == true
    @test pointer(A) == pointer(to_fast_array(A))
    @test pointer(A) != pointer(to_fast_array(Va))
    @test samevalues(Va, to_fast_array(Va))
    @test samevalues(A, to_fast_array(A))
    @test maxabsdif(A, B) ≤ atol
    @test maxabsdif(Va, C) ≤ atol
end

@testset "Broadcasting  " begin
    # Check broadcasting of dimensions.
    @test_throws DimensionMismatch bcastsize(3, 4)
    for (a,b) in (((), ()),
                  ((Int16(5), Int32(6), Int64(7)), (5,6,7)),)
        @test bcastsize(a) === b
        @test bcastsize((), a) === b
        @test bcastsize(a, ()) === b
        @test bcastsize(a, (a..., 3,)) === (b..., 3,)
        @test bcastsize(a, (a..., 3,), (a..., 1, 2)) === (b..., 3, 2)
        @test bcastsize((a..., 3,), a) === (b..., 3,)
        @test bcastsize((a..., 1,), (a..., 3,)) === (b..., 3,)
        @test bcastsize((a..., 3,), (a..., 1,)) === (b..., 3,)
    end
    # Check `bcastcopy`.
    T = Int32
    dims1 = (1,2,3)
    dims2 = (3,2,1,4)
    dims3 = bcastsize(dims1, dims2)
    A1 = generate(T, dims1)
    A2 = generate(T, dims2)
    @test A1 .+ zeros(T, dims2) == bcastcopy(A1, dims3...)
    @test A2 .+ zeros(T, dims1) == bcastcopy(A2, map(Int16, dims3))
    C1 = bcastcopy(A, eltype(A), size(A))
    @test C1 !== A && C1 == A
    C2 = bcastcopy(A, eltype(A), map(Int16, size(A))...)
    @test C2 !== A && C2 == A
    C4 = bcastcopy(A, size(A))
    @test C4 !== A && C4 == A
    C5 = bcastcopy(A, map(Int16, size(A))...)
    @test C5 !== A && C5 == A
    # Check that `bcastlazy attempt` to yield the same object.
    @test bcastlazy(A, eltype(A), size(A)) === A
    @test bcastlazy(A, eltype(A), map(Int16, size(A))...) === A
    @test bcastlazy(A, size(A)) === A
    @test bcastlazy(A, map(Int16, size(A))...) === A
    B1 = bcastlazy(A, Int32, size(A))
    @test B1 !== A && B1 == A
    B2 = bcastlazy(Float32(1), dims)
    @test eltype(B2) === Float32
    @test size(B2) === dims
    @test B2 == ones(Float32, dims)
    B3 = bcastlazy(1, Float16, dims...)
    @test eltype(B3) === Float16
    @test size(B3) === dims
    @test B3 == ones(Float16, dims)
    # Check slices of arrays given by `bcastcopy`/`bcastlazy`
    X = generate(Float64, (dims[1], 1, dims[3:end]...))
    X1 = bcastlazy(X, dims)
    X2 = bcastcopy(X, dims)
    for i in 1:dims[2]
        @test X[:,1,:,:] == X1[:,i,:,:]
        @test X[:,1,:,:] == X2[:,i,:,:]
    end
end

# UnfinishedArray does not extend Base.parent()
struct UnfinishedArray{T,N,A<:AbstractArray{T,N},S} <: PseudoArray{T,N,S}
    arr::A
    cnt::Int
end
UnfinishedArray(arr::A, cnt::Integer=0) where {T,N,A<:AbstractArray{T,N}} = begin
    S = typeof(IndexStyle(arr))
    return UnfinishedArray{T,N,A,S}(arr,cnt)
end

# DummyArray does extend Base.parent()
struct DummyArray{T,N,A<:AbstractArray{T,N},S} <: PseudoArray{T,N,S}
    arr::A
    cnt::Int
end
DummyArray(arr::A, cnt::Integer=0) where {T,N,A<:AbstractArray{T,N}} = begin
    S = typeof(IndexStyle(arr))
    return DummyArray{T,N,A,S}(arr,cnt)
end
Base.parent(A::DummyArray) = A.arr

@testset "Custom arrays " begin
    inds = map(n -> Base.OneTo(n), dims)
    T = Float32
    N = length(dims)
    D1 = Dict("units" => "photons", "Δx" => 0.20, "Δy" => 0.15)
    @test isa(D1, Dict{String,Any})
    K2 = (x=true, y=1.8, units="µm")
    D2 = Dict(pairs(K2)...)
    @test isa(D2, Dict{Symbol,Any})
    D3 = Dict(:x => 1, :y => 2, :z => 3)
    @test isa(D3, Dict{Symbol,Int})
    D4 = Dict(:x => 1, Float64 => 3.1, "string" => "hello")
    @test isa(D4, Dict{Any,Any})
    D5 = Dict(1.0 => 1, 2.0 => :2, 3.0 => "3")
    @test isa(D5, Dict{Float64,Any})
    G = AnnotatedArray(zeros(T, dims), pairs(D1)...)
    F = AnnotatedArray{T}(parent(G), Dict{String,Any}())
    H = AnnotatedArray{T,N}(undef, dims, Dict{Symbol,Float32}())
    @test (AnnotatedArrays.keytype(D1) ===
           AnnotatedArrays.keytype(typeof(D1)) ===
           keytype(D1) === keytype(typeof(D1)) === String)
    @test (AnnotatedArrays.valtype(D1) ===
           AnnotatedArrays.valtype(typeof(D1)) ===
           valtype(D1) === valtype(typeof(D1)) === Any)
    @test (AnnotatedArrays.keytype(D2) ===
           AnnotatedArrays.keytype(typeof(D2)) ===
           keytype(D2) === keytype(typeof(D2)) === Symbol)
    @test (AnnotatedArrays.valtype(D2) ===
           AnnotatedArrays.valtype(typeof(D2)) ===
           valtype(D2) === valtype(typeof(D2)) === Any)
    @test (AnnotatedArrays.keytype(K2) ===
           AnnotatedArrays.keytype(typeof(K2)) === Symbol)
    @test (AnnotatedArrays.valtype(K2) ===
           AnnotatedArrays.valtype(typeof(K2)) === Any)
    @test (AnnotatedArrays.keytype(D3) ===
           AnnotatedArrays.keytype(typeof(D3)) ===
           keytype(D3) === keytype(typeof(D3)) === Symbol)
    @test (AnnotatedArrays.valtype(D3) ===
           AnnotatedArrays.valtype(typeof(D3)) ===
           valtype(D3) === valtype(typeof(D3)) === Int)
    @test (AnnotatedArrays.keytype(D4) ===
           AnnotatedArrays.keytype(typeof(D4)) ===
           keytype(D4) === keytype(typeof(D4)) === Any)
    @test (AnnotatedArrays.valtype(D4) ===
           AnnotatedArrays.valtype(typeof(D4)) ===
           valtype(D4) === valtype(typeof(D4)) === Any)
    @test (AnnotatedArrays.keytype(D5) ===
           AnnotatedArrays.keytype(typeof(D5)) ===
           keytype(D5) === keytype(typeof(D5)) === Float64)
    @test (AnnotatedArrays.valtype(D5) ===
           AnnotatedArrays.valtype(typeof(D5)) ===
           valtype(D5) === valtype(typeof(D5)) === Any)

    # Forbidden key types.
    @test_throws ErrorException AnnotatedArray{T}(undef, dims, Dict{Any,Any}())
    @test_throws ErrorException AnnotatedArray{T}(undef, dims, Dict{Int32,Any}())
    @test_throws ErrorException AnnotatedArray{T}(undef, dims, Dict{CartesianIndex,Any}())
    @test_throws ErrorException AnnotatedArray{T}(undef, dims, Dict{UnitRange{Int},Any}())
    @test_throws ErrorException AnnotatedArray{T}(undef, dims, Dict{Colon,Any}())

    # Common errors.
    @test_throws ErrorException AnnotatedArray(undef, dims)
    @test_throws ErrorException AnnotatedArray(undef, dims...)

    # Annotated arrays with immutable properties.
    let X = AnnotatedArray(A, (a = 1,)),
        Y = AnnotatedArray{T}(A, (a = "1", b = "2"))
        @test isa(X, StaticallyAnnotatedArray)
        @test (AnnotatedArrays.keytype(X) === keytype(X) ===
               AnnotatedArrays.keytype(typeof(X)) === keytype(typeof(X)) === Symbol)
        @test (AnnotatedArrays.valtype(X) === valtype(X) ===
               AnnotatedArrays.valtype(typeof(X)) === valtype(typeof(X)) === Int)
        @test nkeys(X) == 1
        @test propertyname(typeof(X), :foo) === :foo
        @test isa(Y, StaticallyAnnotatedArray)
        @test (AnnotatedArrays.keytype(Y) === keytype(Y) ===
               AnnotatedArrays.keytype(typeof(Y)) === keytype(typeof(Y)) === Symbol)
        @test (AnnotatedArrays.valtype(Y) === valtype(Y) ===
               AnnotatedArrays.valtype(typeof(Y)) === valtype(typeof(Y)) === String)
        @test nkeys(Y) == 2
        @test_throws ErrorException X.field = 40
        @test_throws ErrorException X[:field] = 40
        @test :a ∈ propertynames(X)
        @test :b ∈ propertynames(Y)
        @test haskey(X, :a) == true
        @test haskey(X, :b) == false
        @test haskey(X, "b") == false
        @test getkey(X, :a, :x) == :a
        @test getkey(X, :b, :x) == :x
        @test getkey(X, "b", :x) == :x
        @test get(X, :a, :x) == X.a
        @test get(X, :b, :x) == :x
        @test get(X, "b", :x) == :x
        @test X[:a] == X.a
        @test_throws ErrorException X.unknown
        @test_throws ErrorException delete!(X, :a)
        @test_throws ErrorException get!(X, :a, 1)
        @test_throws ErrorException pop!(X, :a)
        @test_throws ErrorException pop!(X, :a, 1)
        @test_throws ErrorException merge!(X, Y)
        @test_throws ErrorException merge!(+, X, Y)
    end
    let X = AnnotatedArray{T}(undef, dims, (a = 1.0, b = 2.0, c = 3.0))
        @test isa(X, StaticallyAnnotatedArray)
        @test eltype(X) === T
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Float64
        @test nkeys(X) == 3
    end
    let X = AnnotatedArray{T,N}(A, (a = 1, b = 2, c = 3, d = 4))
        @test isa(X, StaticallyAnnotatedArray)
        @test eltype(X) === T
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Int
        @test nkeys(X) == 4
    end
    let X = AnnotatedArray{T,N}(undef, dims,
                                (a = 1, b = 2, c = 3, d = 4, e = 5))
        @test isa(X, StaticallyAnnotatedArray)
        @test eltype(X) === T
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Int
        @test nkeys(X) == 5
    end

    # Annotated arrays with mutable properties.
    let X = AnnotatedArray{T,N}(A)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T,N}(A, D3)
        @test isa(X, DynamicallyAnnotatedArray)
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Int
    end
    let X = AnnotatedArray{T,N}(A, pairs(D3)...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    @test_throws ErrorException AnnotatedArray{T,N}(A, D4)
    @test_throws ErrorException AnnotatedArray{T,N}(A, pairs(D4))
    @test_throws ErrorException AnnotatedArray{T,N}(A, pairs(D4)...)
    let X = AnnotatedArray{T,N}(undef, dims)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
        X.foo = 42
        X.bar = π
        @test X.foo === X[:foo] === 42
        @test X.bar === X[:bar] === π
    end
    let X = AnnotatedArray{T,N}(undef, dims, D2)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T,N}(undef, dims, pairs(D2)...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T,N}(A; K2...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T,N}(undef, dims; K2...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(A)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(A, D1)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === String
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
        @test propertyname(typeof(X), :somekey) == "somekey"
    end
    let X = AnnotatedArray{T}(A, pairs(D1)...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === String
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(undef, dims)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(undef, dims, D2)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(undef, dims, pairs(D2)...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(undef, dims; K2...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray(A)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray(A, D1)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === String
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray(A, pairs(D2)...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray(A; K2...)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test ndims(X) == N
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray(A, π => "π", sqrt(2) => "√2", 1.3 => :bar)
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === eltype(A)
        @test axes(X) === axes(A)
        @test nkeys(X) == 3
        @test AnnotatedArrays.keytype(X) === keytype(X) === Real
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
        @test_throws ErrorException propertyname(typeof(X), :foo)
    end
    let X = AnnotatedArray{T}(undef, dims, π=>"π", sqrt(2)=>"√2")
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test size(X) === dims
        @test nkeys(X) == 2
        @test AnnotatedArrays.keytype(X) === keytype(X) === Real
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
    end
    let X = AnnotatedArray{T}(undef, dims...; a="π", b="√2")
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
        @test propertyname(typeof(X), :foo) === :foo
    end
    let X = AnnotatedArray{T,N}(undef, dims...; a="π", b="√2")
        @test isa(X, DynamicallyAnnotatedArray)
        @test eltype(X) === T
        @test size(X) === dims
        @test AnnotatedArrays.keytype(X) === keytype(X) === Symbol
        @test AnnotatedArrays.valtype(X) === valtype(X) === Any
        @test X.a == "π"
        @test X.b == "√2"
        @test :a ∈ propertynames(X)
    end

    Q = UnfinishedArray(Va)
    @test_throws ErrorException Q[1]
    R = DummyArray(Va)

    copyto!(G, rand(dims...)) # exercise setindex! for linear indices
    copyto!(R, rand(size(Va)...)) # exercise setindex! for Cartesian indices
    sum1 = 0.0
    for val in G
        sum1 += val
    end
    @test samevalues(F, G)
    @test sum(F) == sum(G) ≈ sum1
    sum2 = 0.0
    for val in R
        sum2 += val
    end
    @test samevalues(R, Va)
    @test sum(R) == sum(Va) ≈ sum2
    @test pointer(parent(F)) == pointer(parent(G))
    @test length(F) == length(G) == prod(dims)
    @test eltype(F) == eltype(G) == T
    @test ndims(F) == ndims(G) == N
    @test size(F) == size(G) == dims
    @test ntuple(d -> size(F,d), N) == ntuple(d -> size(G,d), N) == dims
    @test axes(F) == axes(G) == inds
    @test ntuple(d -> axes(F,d), N) == ntuple(d -> axes(G,d), N) == inds
    @test Base.axes1(F) == Base.axes1(G) == inds[1]
    @test Base.elsize(F) == Base.elsize(G) == Base.elsize(parent(F))
    @test Base.elsize(R) == Base.elsize(typeof(R)) == Base.elsize(Va)
    @test sizeof(F) == sizeof(G) == sizeof(parent(F))
    @test IndexStyle(F) == IndexStyle(G) == IndexStyle(parent(F))
    @test_throws ErrorException parent(Q)
    @test IndexStyle(R) == IndexStyle(parent(R)) == IndexCartesian()
    @test pairs(IndexStyle(F), F) == pairs(IndexStyle(G), G) == pairs(IndexStyle(parent(F)), parent(F))
    @test pairs(IndexStyle(R), R) == pairs(IndexStyle(Va), Va)
    @test_throws ErrorException size(F,0)
    @test_throws BoundsError F[0]
    @test keytype(F) === keytype(G) === String
    @test keytype(H) === Symbol
    @test valtype(F) === valtype(G) === Any
    @test valtype(H) === Float32

    @test nkeys(F) == nkeys(properties(F)) == 0 && nkeys(G) == nkeys(properties(G)) == 3
    @test keys(F) == keys(properties(F))
    @test keys(G) == keys(properties(G))
    @test values(G) == values(properties(G))
    @test pairs(G) == pairs(properties(G))
    @test_throws ArgumentError (H["ga"] = π)
    @test haskey(H, "ga") == false
    @test (H[:ga] = π) ≈ π
    @test haskey(H, :ga) == true
    @test keys(H) == keys(properties(H))
    @test getkey(H, :ga, nothing) === :ga
    @test getkey(H, :bu, nothing) === nothing
    @test getkey(H, "ga", nothing) === nothing
    @test get(H, :ga, nothing) ≈ π
    @test get(H, :bu, nothing) === nothing
    @test get(H, "ga", nothing) === nothing
    @test get!(H, :ga, nothing) ≈ π
    @test get!(H, :bu, 55) == 55 && haskey(H, :bu) == true
    @test_throws Exception get!(H, "ga", nothing)
    merge!(F, G)
    @test nkeys(F) == nkeys(G)
    @test keys(F) == keys(G)
    merge!(F, G, Dict("ga" => π))
    @test nkeys(F) == nkeys(G) + 1
    @test pop!(F, "ga") == π
    @test_throws KeyError F["ga"]
    @test_throws KeyError pop!(F, "ga")
    @test pop!(F, "ga", 42) == 42
    @test nkeys(delete!(F, "ga")) == nkeys(G)
    @test nkeys(delete!(F, "units")) == nkeys(G) - 1
    @test merge(G, F) == merge(properties(G), properties(F))
    @test merge(D2, F, D1) == merge(D2, properties(F), D1)
    Dat1 = Dict{String,Any}("a" => 1, "b" => sqrt(2), "c" => 4)
    Dat2 = Dict{String,Any}("a" => -1.0, "b" => π, "c" => sqrt(3), "d" => 11)
    for k in keys(F)
        delete!(F, k)
    end
    for k in keys(G)
        delete!(G, k)
    end
    merge!(F, Dat1)
    merge!(G, Dat2)
    @test merge(+, F, Dat2) == merge(+, Dat1, Dat2)
    @test merge(-, Dat1, G) == merge(-, Dat1, Dat2)
    @test merge(*, F, G) == merge(*, Dat1, Dat2)
    @test properties(merge!(+, F, Dat2)) == merge(+, Dat1, Dat2)
    @test merge!(+, copy(Dat1), G) == properties(F)
end

end # module


text:
using Test
using MDCDL
using FFTW
using Random
using LinearAlgebra

@testset "RNSOLT" begin
    include("testsetGenerator.jl")

    rcsd1D = rnsoltValidConfigSet1D()
    rcsd2D = rnsoltValidConfigSet2D()
    rcsd3D = rnsoltValidConfigSet3D()

    rcsd = [ rcsd1D, rcsd2D, rcsd3D ]

    Random.seed!(9387509284)

    @testset "Constructor" begin

        maxDims = 3
        defaultType = Float64

        for d in 1:maxDims
            allcfgset = [ (crdf.I, crord.I .- 1, crnch.I) for crdf in CartesianIndices(tuple(fill(4,d)...)), crord in CartesianIndices(tuple(fill(6+1,d)...)), crnch in CartesianIndices(tuple(fill(10,2)...)) ]
            cfgsetTypeI  = filter(c -> c[3][1] == c[3][2], allcfgset)
            cfgsetTypeII = filter(c -> c[3][1] != c[3][2], allcfgset)
            cfgset = vcat(randsubseq(cfgsetTypeI, 50 / length(cfgsetTypeI))..., randsubseq(cfgsetTypeII, 50 / length(cfgsetTypeII))...)

            for (df, ord, nch) in cfgset
                if prod(df) > sum(nch)
                    @test_throws AssertionError Rnsolt(df, ord, nch)
                    continue
                end

                if !(cld(prod(df),2) <= nch[1] <= sum(nch) - fld(prod(df),2)) || !(fld(prod(df),2) <= nch[2] <= sum(nch) - cld(prod(df),2))
                    @test_throws AssertionError Rnsolt(df, ord, nch)
                    continue
                end

                if nch[1] != nch[2] && any(isodd.(ord))
                    @test_throws AssertionError Rnsolt(df, ord, nch)
                    continue
                end

                nsolt = Rnsolt(df, ord, nch)

                if nch[1] == nch[2]
                    @test nsolt isa RnsoltTypeI
                    @test istype1(nsolt) == true && istype2(nsolt) == false
                else
                    @test nsolt isa RnsoltTypeII
                    @test istype1(nsolt) == false && istype2(nsolt) == true
                end
            end
        end
    end

    @testset "DefaultValues" begin
        @testset "maximally-sampled" begin
            for d in 1:3
                for n = 1:20
                    df = (rand(1:4, d)...,)
                    ord = (fill(0, d)...,)
                    nch = prod(df)

                    x = rand(df...)
                    A = MDCDL.permdctmtx(df...)
                    expctd = A*vec(x)

                    nsolt = Rnsolt(df, ord, nch)
                    pvx = reshape(x, prod(df), :)
                    y = analyze(nsolt, pvx, (fill(1,d)...,))
                    actual = reshape(y, prod(df))

                    @test actual ≈ expctd
                end
            end
        end
    end

    @testset "FilterSymmetry" begin
        for d in 1:length(rcsd), (df, ord, nch) in rcsd[d]
            nsolt = Rnsolt(df, ord, nch)
            rand!(nsolt)

            afb = MDCDL.analysisbank(nsolt)
            Γ = Diagonal(vcat(fill(1,nch[1]), fill(-1,nch[2])))

            @test afb ≈ Γ * reverse(afb; dims=2)
        end
    end

    @testset "AnalysisSynthesis" begin
        # output mode options for analyzer
        oms = [ Shapes.Separated, Shapes.Combined, Shapes.Vec ]
        @testset "NormalOrder" begin
            for d in 1:length(rcsd), (df, ord, nch) in rcsd[d]
                szx = df .* (ord .+ 1)
                nsolt = Rnsolt(df, ord, nch)
                rand!(nsolt)

                x = rand(Float64, szx...)

                foreach(oms) do om
                    nsop = createTransform(nsolt, om(size(x)))

                    y = analyze(nsop, x)
                    rx = synthesize(nsop, y)

                    @test size(x) == size(rx)
                    @test rx ≈ x
                end
            end
        end
        @testset "DimensionPermutation" begin
            for d in 1:length(rcsd), (df, ord, nch) in rcsd[d]
                szx = df .* (ord .+ 1)
                nsolt = Rnsolt(df, ord, nch, perm=(randperm(d)...,))
                rand!(nsolt)

                x = rand(Float64, szx...)

                foreach(oms) do om
                    nsop = createTransform(nsolt, om(size(x)))

                    y = analyze(nsop, x)
                    rx = synthesize(nsop, y)

                    @test size(x) == size(rx)
                    @test rx ≈ x
                end
            end
        end
    end

    @testset "AnalyzerKernel" begin
        for d in 1:length(rcsd), (df, ord, nch) in rcsd[d]
            szx = df .* (ord .+ 1)
            x = rand(Float64, szx)

            nsolt = Rnsolt(df, ord, nch)
            rand!(nsolt)

            ana = createTransform(nsolt, Shapes.Separated())
            ya = analyze(ana, x)

            afs = analysiskernels(nsolt)
            myfilter(A, h) = begin
                ha = zero(A)
                ha[UnitRange.(1, size(h))...] = h
                real(ifft(fft(A).*fft(ha)))
            end
            offset = df .- 1

            sys = map(afs) do af
                fx = myfilter(x, af)
                dwfx = downsample(fx, df, offset)
                circshift(dwfx, (-1 .* fld.(ord,2)))
            end
            @test all(ya .≈ sys)
        end
    end

    @testset "SynthesizerKernel" begin
        for d in 1:length(rcsd), (df, ord, nch) in rcsd[d]

            nsolt = Rnsolt(df, ord, nch)
            rand!(nsolt)

            y = [ rand(Float64,((ord.+1) .* df)...) for p in 1:sum(nch) ]

            syn = createTransform(nsolt, Shapes.Separated())
            x = synthesize(syn, y)

            sfs = synthesiskernels(nsolt)
            myfilter(A, h) = begin
                ha = zero(A)
                ha[UnitRange.(1, size(h))...] = h
                real(ifft(fft(A).*fft(ha)))
            end
            offset = df .- 1

            subrxs = sum(map((yp,fp)->myfilter( MDCDL.upsample(yp, df), fp ), y, sfs))
            subrx = circshift(subrxs, -1 .* df .* cld.(ord,2))
            @test size(x) == size(subrx)
            @test subrx ≈ x
        end
    end

    @testset "Factorization" begin
        for d in 1:length(rcsd), (df, ord, nch) in rcsd[d]
            src = Rnsolt(df, ord, nch)
            dst = Rnsolt(df, ord, nch)
            rand!(src)

            (angs, mus) = getrotations(src)
            setrotations!(dst, angs, mus)

            @test src.W0 ≈ dst.W0 && src.U0 ≈ dst.U0

            foreach(src.Udks, dst.Udks) do srcUs, dstUs
                @test all(srcUs .≈ dstUs)
            end

            if istype2(src)
                foreach(src.Wdks, dst.Wdks) do srcWs, dstWs
                    @test all(srcWs .≈ dstWs)
                end
            end
        end
    end

end

nothing

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
# Implements a 1-Neighborhod Stochastic Local Search for MAP Inference

"""
    localsearch!(evidence,spn,query)

Performs a local search for the maximum value of a configuration of `query` variables for given some `evidence`.
Returns the value of the last configuration visited by the search (a local optimum or maximum iterations limit).

# Parameters
- `evidence`: vector of values for each variable with initial solution; NaN values denote sum-out variables; the algorithm overwrites the vector with the best solution found.
- `spn`: sum-product network.
- `query`: set of variables to be maximized.
- `maxiterations`: maximum number of iterations [default: 100].
"""
function localsearch!(evidence::AbstractVector,spn::SumProductNetwork,query::AbstractSet,maxiterations=100)
    values = Vector{Float64}(undef,length(spn))
    best = logpdf!(values,spn,evidence)
    vdims = vardims(spn)
    it = 1
    before = time_ns()
    while it < maxiterations
        now = time_ns()
        if (now-before) > 60e9
            # report incumbent result every minute
            printstyled("[$it/$maxiterations]"; color = :lightcyan)
            printstyled(" $best "; color = :normal)
            printstyled("  [$((now-before)*1e-9)s]\n"; color = :light_black)
            before = now
        end
        cur = best
        for var in query
            v = evidence[var]
            for i = 1:vdims[var]
                evidence[var] = i
                prob = logpdf!(values,spn,evidence)
                if prob > best
                    best = prob
                    v = i
                    break
                end
            end
            evidence[var] = v
        end
        if cur == best # local optimum reached
            break
        end
        it += 1
    end
    best
end


text:
# model for FashionMNIST
#
# reference: https://github.com/FluxML/model-zoo/blob/master/mnist/conv.jl

using Flux
using Flux: onehotbatch, argmax, crossentropy, throttle, @epochs
using BSON: @save, @load
using Base.Iterators: partition

using MLDatasets # FashionMNIST
using ColorTypes: N0f8, Gray

const Img = Matrix{Gray{N0f8}}

function prepare_train()
    # load full training set
    train_x, train_y = FashionMNIST.traindata() # 6000
    imgs = Img.([train_x[:,:,i] for i in 1:6000])
    labels = onehotbatch(train_y, 0:9)
    imgs, labels
end

function prepare_test()
    # load full test set
    test_x,  test_y  = FashionMNIST.testdata() # 10_000
    test_imgs = Img.([test_x[:,:,i] for i in 1:10_000])

    # Prepare test set
    tX = cat(4, float.(test_imgs[1:10_000])...) |> gpu
    tY = onehotbatch(test_y[1:10_000], 0:9) |> gpu
    tX, tY
end

imgs, labels = prepare_train()
tX, tY = prepare_test()

# Partition into batches of size 1,000
train = [(cat(4, float.(imgs[i])...), labels[:,i])
         for i in partition(1:6000, 1000)]

train = gpu.(train)

m = Chain(
  Conv((2,2), 1=>16, relu),
  x -> maxpool(x, (2,2)),
  Conv((2,2), 16=>8, relu),
  x -> maxpool(x, (2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(288, 10), softmax) |> gpu

# m(train[1][1])

loss(x, y) = crossentropy(m(x), y)

accuracy(x, y) = mean(argmax(m(x)) .== argmax(y))

evalcb = throttle(() -> @show(accuracy(tX, tY)), 1)
opt = ADAM(params(m))

@epochs 5 Flux.train!(loss, train, opt, cb = evalcb)

#=
INFO: Epoch 1
accuracy(tX, tY) = 0.1
accuracy(tX, tY) = 0.1
INFO: Epoch 2
accuracy(tX, tY) = 0.1
accuracy(tX, tY) = 0.1
INFO: Epoch 3
accuracy(tX, tY) = 0.1
accuracy(tX, tY) = 0.1001
INFO: Epoch 4
accuracy(tX, tY) = 0.1163
accuracy(tX, tY) = 0.1947
INFO: Epoch 5
accuracy(tX, tY) = 0.2864
accuracy(tX, tY) = 0.3105
=#

# @save "fashion_mnist_model.bson" m
# weights = Tracker.data.(params(m)) ;
# @save "fashion_mnist_weights.bson" weights


text:
@testset "NLPAtX" begin
  #Test unconstrained NLPAtX
  uncons_nlp_at_x = NLPAtX(zeros(10))

  @test uncons_nlp_at_x.x == zeros(10)
  @test isnan(uncons_nlp_at_x.fx)
  @test uncons_nlp_at_x.gx == zeros(0)
  @test uncons_nlp_at_x.Hx == zeros(0, 0)
  @test uncons_nlp_at_x.mu == zeros(0)
  @test uncons_nlp_at_x.cx == zeros(0)
  @test uncons_nlp_at_x.Jx == zeros(0, 0)

  @test uncons_nlp_at_x.lambda == zeros(0)
  @test isnan(uncons_nlp_at_x.current_time)
  @test isnan(uncons_nlp_at_x.current_score)

  #check constrained NLPAtX
  cons_nlp_at_x = NLPAtX(zeros(10), zeros(10))

  @test cons_nlp_at_x.x == zeros(10)
  @test isnan(cons_nlp_at_x.fx)
  @test cons_nlp_at_x.gx == zeros(0)
  @test cons_nlp_at_x.Hx == zeros(0, 0)
  @test cons_nlp_at_x.mu == zeros(0)
  @test cons_nlp_at_x.cx == zeros(0)
  @test cons_nlp_at_x.Jx == zeros(0, 0)
  @test (false in (cons_nlp_at_x.lambda .== 0.0)) == false
  @test isnan(cons_nlp_at_x.current_time)
  @test isnan(cons_nlp_at_x.current_score)

  update!(cons_nlp_at_x, Hx = ones(20, 20), gx = ones(2), lambda = zeros(2))
  compress_state!(cons_nlp_at_x, max_vector_size = 5, lambda = zeros(0), gx = true)
  @test cons_nlp_at_x.Hx == zeros(0, 0)
  @test cons_nlp_at_x.x == [0.0]
  @test cons_nlp_at_x.lambda == zeros(0)
  @test cons_nlp_at_x.gx == zeros(0)

  # On vérifie que la fonction update! fonctionne
  update!(uncons_nlp_at_x, x = ones(10), fx = 1.0, gx = ones(10))
  update!(uncons_nlp_at_x, lambda = ones(10), current_time = 1.0)
  update!(uncons_nlp_at_x, Hx = ones(10, 10), mu = ones(10), cx = ones(10), Jx = ones(10, 10))

  @test (false in (uncons_nlp_at_x.x .== 1.0)) == false #assez bizarre comme test...
  @test uncons_nlp_at_x.fx == 1.0
  @test (false in (uncons_nlp_at_x.gx .== 1.0)) == false
  @test (false in (uncons_nlp_at_x.Hx .== 1.0)) == false
  @test uncons_nlp_at_x.mu == ones(10)
  @test uncons_nlp_at_x.cx == ones(10)
  @test (false in (uncons_nlp_at_x.Jx .== 1.0)) == false
  @test (false in (uncons_nlp_at_x.lambda .== 1.0)) == false
  @test uncons_nlp_at_x.current_time == 1.0
  @test isnan(uncons_nlp_at_x.current_score)

  reinit!(uncons_nlp_at_x)
  @test uncons_nlp_at_x.x == ones(10)
  @test isnan(uncons_nlp_at_x.fx)
  reinit!(uncons_nlp_at_x, x = zeros(10))
  @test uncons_nlp_at_x.x == zeros(10)
  @test isnan(uncons_nlp_at_x.fx)
  reinit!(uncons_nlp_at_x, zeros(10))
  @test uncons_nlp_at_x.x == zeros(10)
  @test isnan(uncons_nlp_at_x.fx)
  reinit!(uncons_nlp_at_x, zeros(10), l = zeros(0))
  @test uncons_nlp_at_x.x == zeros(10)
  @test isnan(uncons_nlp_at_x.fx)

  c_uncons_nlp_at_x = copy_compress_state(uncons_nlp_at_x, max_vector_size = 5)

  @test c_uncons_nlp_at_x != uncons_nlp_at_x
  @test c_uncons_nlp_at_x.x == [0.0]
  @test c_uncons_nlp_at_x.lambda == [1.0]

  uncons_nlp_at_x.Hx = zeros(10, 10)
  zip_uncons_nlp_at_x =
    compress_state!(uncons_nlp_at_x, keep = true, save_matrix = true, max_vector_size = 5, Hx = 1)
  zip_uncons_nlp_at_x.Hx == 0.0

  nlp_64 = NLPAtX(ones(10))
  nlp_64.x = ones(10)
  nlp_64.fx = 1.0
  nlp_64.gx = ones(10)

  # nlp_32 = convert_nlp(Float32, nlp_64)
  # @test typeof(nlp_32.x[1]) == Float32
  # @test typeof(nlp_32.fx[1]) == Float32
  # @test typeof(nlp_32.gx[1]) == Float32
  # @test isnan(nlp_32.mu[1])
  # @test isnan(nlp_32.current_time)
  #
  # @test typeof(nlp_64.x[1]) == Float64
  # @test typeof(nlp_64.fx[1]) == Float64
  # @test typeof(nlp_64.gx[1]) == Float64
  # @test isnan(nlp_64.mu[1])
  # @test isnan(nlp_64.current_time)

  #Test the _size_check:
  try
    NLPAtX(ones(5), gx = zeros(4))
    @test false
  catch
    @test true
  end
  try
    NLPAtX(ones(5), mu = zeros(4))
    @test false
  catch
    @test true
  end
  try
    NLPAtX(ones(5), Hx = zeros(4, 4))
    @test false
  catch
    @test true
  end
  try
    NLPAtX(ones(5), zeros(1), cx = zeros(2))
    @test false
  catch
    @test true
  end

  # Test matrix types
  state = NLPAtX(ones(5), ones(2), Jx = spzeros(2, 5), Hx = spzeros(5, 5))
  @test typeof(spzeros(2, 5)) == typeof(state.Jx)
  @test typeof(spzeros(5, 5)) == typeof(state.Hx)
  state = NLPAtX(ones(5), ones(2), Jx = spzeros(2, 5))
  @test typeof(spzeros(2, 5)) == typeof(state.Jx)
  @test typeof(zeros(5, 5)) == typeof(state.Hx)
end


text:
## C types
const Ctime_t = Int

## HDF5 types and constants
const haddr_t  = UInt64
const hbool_t  = Cuint
const herr_t   = Cint
const hid_t    = Int64
const hsize_t  = UInt64
const hssize_t = Int64
const htri_t   = Cint   # pseudo-boolean (negative if error)

const H5Z_filter_t = Cint

# For VLEN
struct hvl_t
    len::Csize_t
    p::Ptr{Cvoid}
end
const HVL_SIZE = sizeof(hvl_t)

# Reference types
struct hobj_ref_t
    buf::UInt64 # H5R_OBJ_REF_BUF_SIZE bytes
end
struct hdset_reg_ref_t
    buf::NTuple{12,UInt8} # H5R_DSET_REG_REF_BUF_SIZE bytes
end
const HOBJ_REF_T_NULL = hobj_ref_t(0x0)
const HDSET_REG_REF_T_NULL = hdset_reg_ref_t(ntuple(_ -> 0x0, Val(12)))

#= TODO: when upgraded to using newer HDF5 v1.12 reference API, can replace both with:
struct H5R_ref_t
    # Element type of UInt64 and not UInt8 to get 8-byte alignment
    buf::NTuple{8,UInt64} # H5R_REF_BUF_SIZE bytes
end
const H5R_REF_T_NULL = H5R_ref_t(ntuple(_ -> 0x0, Val(64)))
=#

# For group information
struct H5G_info_t
    storage_type::Cint # enum H5G_storage_type_t
    nlinks::hsize_t
    max_corder::Int64
    mounted::hbool_t
end

# For objects
struct H5_ih_info_t
    index_size::hsize_t
    heap_size::hsize_t
end
struct H5O_info_t # version 1 type H5O_info1_t
    fileno::Cuint
    addr::haddr_t
    otype::Cint # enum H5O_type_t
    rc::Cuint
    atime::Ctime_t
    mtime::Ctime_t
    ctime::Ctime_t
    btime::Ctime_t
    num_attrs::hsize_t
    #{ inlined struct H5O_hdr_info_t named type
    version::Cuint
    nmesgs::Cuint
    nchunks::Cuint
    flags::Cuint
    total::hsize_t
    meta::hsize_t
    mesg::hsize_t
    free::hsize_t
    present::UInt64
    shared::UInt64
    #}
    #{ inlined anonymous struct named meta_size
    meta_obj::H5_ih_info_t
    meta_attr::H5_ih_info_t
    #}
end

# For links
struct H5L_info_t
    linktype::Cint
    corder_valid::hbool_t
    corder::Int64
    cset::Cint # enum H5T_cset_t
    u::haddr_t # union { haddr_t address, size_t val_size }
end

# For registering filters
struct H5Z_class_t # version 2 type H5Z_class2_t
    version::Cint # = H5Z_CLASS_T_VERS
    id::H5Z_filter_t # Filter ID number
    encoder_present::Cuint # Does this filter have an encoder?
    decoder_present::Cuint # Does this filter have a decoder?
    name::Ptr{UInt8} # Comment for debugging
    can_apply::Ptr{Cvoid} # The "can apply" callback
    set_local::Ptr{Cvoid} # The "set local" callback
    filter::Ptr{Cvoid} # The filter callback
end

# MPI communicators required by H5P
abstract  type Hmpih end
primitive type Hmpih32 <: Hmpih 32 end # MPICH C/Fortran, OpenMPI Fortran: 32 bit handles
primitive type Hmpih64 <: Hmpih 64 end # OpenMPI C: pointers (mostly 64 bit)

# Private function to extract exported global library constants.
# Need to call H5open to ensure library is initalized before reading these constants.
# Although these are runtime initalized constants, in practice their values are stable, so
# we can precompile for improved latency.
let libhdf5handle = Ref(Libdl.dlopen(libhdf5))
    ccall(Libdl.dlsym(libhdf5handle[], :H5open), herr_t, ())
    global _read_const(sym::Symbol) = unsafe_load(cglobal(Libdl.dlsym(libhdf5handle[], sym), hid_t))
end

# iteration order constants
const H5_ITER_UNKNOWN = -1
const H5_ITER_INC     = 0
const H5_ITER_DEC     = 1
const H5_ITER_NATIVE  = 2
const H5_ITER_N       = 3

# indexing type constants
const H5_INDEX_UNKNOWN   = -1
const H5_INDEX_NAME      = 0
const H5_INDEX_CRT_ORDER = 1
const H5_INDEX_N = 2

# dataset constants
const H5D_COMPACT      = 0
const H5D_CONTIGUOUS   = 1
const H5D_CHUNKED      = 2

# allocation times (C enum H5D_alloc_time_t)
const H5D_ALLOC_TIME_ERROR = -1
const H5D_ALLOC_TIME_DEFAULT = 0
const H5D_ALLOC_TIME_EARLY = 1
const H5D_ALLOC_TIME_LATE = 2
const H5D_ALLOC_TIME_INCR = 3

# used to "unset" chunk cache configuration parameters
const H5D_CHUNK_CACHE_NSLOTS_DEFAULT = -1 % Csize_t
const H5D_CHUNK_CACHE_NBYTES_DEFAULT = -1 % Csize_t
const H5D_CHUNK_CACHE_W0_DEFAULT = Cdouble(-1)

# error-related constants
const H5E_DEFAULT      = 0

# file access modes
const H5F_ACC_RDONLY     = 0x0000
const H5F_ACC_RDWR       = 0x0001
const H5F_ACC_TRUNC      = 0x0002
const H5F_ACC_EXCL       = 0x0004
const H5F_ACC_DEBUG      = 0x0008
const H5F_ACC_CREAT      = 0x0010
const H5F_ACC_SWMR_WRITE = 0x0020
const H5F_ACC_SWMR_READ  = 0x0040

# Library versions
const H5F_LIBVER_EARLIEST = 0
const H5F_LIBVER_V18      = 1
const H5F_LIBVER_V110     = 2
const H5F_LIBVER_LATEST   = H5F_LIBVER_V110

# object types
const H5F_OBJ_FILE     = 0x0001
const H5F_OBJ_DATASET  = 0x0002
const H5F_OBJ_GROUP    = 0x0004
const H5F_OBJ_DATATYPE = 0x0008
const H5F_OBJ_ATTR     = 0x0010
const H5F_OBJ_ALL      = (H5F_OBJ_FILE|H5F_OBJ_DATASET|H5F_OBJ_GROUP|H5F_OBJ_DATATYPE|H5F_OBJ_ATTR)
const H5F_OBJ_LOCAL    = 0x0020

# other file constants
const H5F_SCOPE_LOCAL   = 0
const H5F_SCOPE_GLOBAL  = 1
const H5F_CLOSE_DEFAULT = 0
const H5F_CLOSE_WEAK    = 1
const H5F_CLOSE_SEMI    = 2
const H5F_CLOSE_STRONG  = 3

# file driver constants
const H5FD_MPIO_INDEPENDENT    = 0
const H5FD_MPIO_COLLECTIVE     = 1
const H5FD_MPIO_CHUNK_DEFAULT  = 0
const H5FD_MPIO_CHUNK_ONE_IO   = 1
const H5FD_MPIO_CHUNK_MULTI_IO = 2
const H5FD_MPIO_COLLECTIVE_IO  = 0
const H5FD_MPIO_INDIVIDUAL_IO  = 1

# object types (C enum H5Itype_t)
const H5I_FILE         = 1
const H5I_GROUP        = 2
const H5I_DATATYPE     = 3
const H5I_DATASPACE    = 4
const H5I_DATASET      = 5
const H5I_ATTR         = 6
const H5I_REFERENCE    = 7
const H5I_VFL          = 8

# Link constants
const H5L_TYPE_HARD    = 0
const H5L_TYPE_SOFT    = 1
const H5L_TYPE_EXTERNAL= 2

# Object constants
const H5O_TYPE_GROUP   = 0
const H5O_TYPE_DATASET = 1
const H5O_TYPE_NAMED_DATATYPE = 2

# Property constants
const H5P_DEFAULT          = hid_t(0)
const H5P_OBJECT_CREATE    = _read_const(:H5P_CLS_OBJECT_CREATE_ID_g)
const H5P_FILE_CREATE      = _read_const(:H5P_CLS_FILE_CREATE_ID_g)
const H5P_FILE_ACCESS      = _read_const(:H5P_CLS_FILE_ACCESS_ID_g)
const H5P_DATASET_CREATE   = _read_const(:H5P_CLS_DATASET_CREATE_ID_g)
const H5P_DATASET_ACCESS   = _read_const(:H5P_CLS_DATASET_ACCESS_ID_g)
const H5P_DATASET_XFER     = _read_const(:H5P_CLS_DATASET_XFER_ID_g)
const H5P_FILE_MOUNT       = _read_const(:H5P_CLS_FILE_MOUNT_ID_g)
const H5P_GROUP_CREATE     = _read_const(:H5P_CLS_GROUP_CREATE_ID_g)
const H5P_GROUP_ACCESS     = _read_const(:H5P_CLS_GROUP_ACCESS_ID_g)
const H5P_DATATYPE_CREATE  = _read_const(:H5P_CLS_DATATYPE_CREATE_ID_g)
const H5P_DATATYPE_ACCESS  = _read_const(:H5P_CLS_DATATYPE_ACCESS_ID_g)
const H5P_STRING_CREATE    = _read_const(:H5P_CLS_STRING_CREATE_ID_g)
const H5P_ATTRIBUTE_CREATE = _read_const(:H5P_CLS_ATTRIBUTE_CREATE_ID_g)
const H5P_OBJECT_COPY      = _read_const(:H5P_CLS_OBJECT_COPY_ID_g)
const H5P_LINK_CREATE      = _read_const(:H5P_CLS_LINK_CREATE_ID_g)
const H5P_LINK_ACCESS      = _read_const(:H5P_CLS_LINK_ACCESS_ID_g)

# Reference constants
const H5R_OBJECT         = 0
const H5R_DATASET_REGION = 1
const H5R_OBJ_REF_BUF_SIZE      = 8  # == sizeof(hobj_ref_t)
const H5R_DSET_REG_REF_BUF_SIZE = 12 # == sizeof(hdset_reg_ref_t)

# Dataspace constants
const H5S_ALL          = hid_t(0)
const H5S_SCALAR       = hid_t(0)
const H5S_SIMPLE       = hid_t(1)
const H5S_NULL         = hid_t(2)
const H5S_UNLIMITED    = typemax(hsize_t)

# Dataspace selection constants
const H5S_SELECT_SET     = 0
const H5S_SELECT_OR      = 1
const H5S_SELECT_AND     = 2
const H5S_SELECT_XOR     = 3
const H5S_SELECT_NOTB    = 4
const H5S_SELECT_NOTA    = 5
const H5S_SELECT_APPEND  = 6
const H5S_SELECT_PREPEND = 7

# type classes (C enum H5T_class_t)
const H5T_INTEGER      = hid_t(0)
const H5T_FLOAT        = hid_t(1)
const H5T_TIME         = hid_t(2)  # not supported by HDF5 library
const H5T_STRING       = hid_t(3)
const H5T_BITFIELD     = hid_t(4)
const H5T_OPAQUE       = hid_t(5)
const H5T_COMPOUND     = hid_t(6)
const H5T_REFERENCE    = hid_t(7)
const H5T_ENUM         = hid_t(8)
const H5T_VLEN         = hid_t(9)
const H5T_ARRAY        = hid_t(10)

# Character types
const H5T_CSET_ASCII   = 0
const H5T_CSET_UTF8    = 1

# Sign types (C enum H5T_sign_t)
const H5T_SGN_NONE     = Cint(0)  # unsigned
const H5T_SGN_2        = Cint(1)  # 2's complement

# Search directions
const H5T_DIR_ASCEND   = 1
const H5T_DIR_DESCEND  = 2

# String padding modes
const H5T_STR_NULLTERM = 0
const H5T_STR_NULLPAD  = 1
const H5T_STR_SPACEPAD = 2

# Type_id constants (LE = little endian, I16 = Int16, etc)
const H5T_STD_I8LE        = _read_const(:H5T_STD_I8LE_g)
const H5T_STD_I8BE        = _read_const(:H5T_STD_I8BE_g)
const H5T_STD_U8LE        = _read_const(:H5T_STD_U8LE_g)
const H5T_STD_U8BE        = _read_const(:H5T_STD_U8BE_g)
const H5T_STD_I16LE       = _read_const(:H5T_STD_I16LE_g)
const H5T_STD_I16BE       = _read_const(:H5T_STD_I16BE_g)
const H5T_STD_U16LE       = _read_const(:H5T_STD_U16LE_g)
const H5T_STD_U16BE       = _read_const(:H5T_STD_U16BE_g)
const H5T_STD_I32LE       = _read_const(:H5T_STD_I32LE_g)
const H5T_STD_I32BE       = _read_const(:H5T_STD_I32BE_g)
const H5T_STD_U32LE       = _read_const(:H5T_STD_U32LE_g)
const H5T_STD_U32BE       = _read_const(:H5T_STD_U32BE_g)
const H5T_STD_I64LE       = _read_const(:H5T_STD_I64LE_g)
const H5T_STD_I64BE       = _read_const(:H5T_STD_I64BE_g)
const H5T_STD_U64LE       = _read_const(:H5T_STD_U64LE_g)
const H5T_STD_U64BE       = _read_const(:H5T_STD_U64BE_g)
const H5T_IEEE_F32LE      = _read_const(:H5T_IEEE_F32LE_g)
const H5T_IEEE_F32BE      = _read_const(:H5T_IEEE_F32BE_g)
const H5T_IEEE_F64LE      = _read_const(:H5T_IEEE_F64LE_g)
const H5T_IEEE_F64BE      = _read_const(:H5T_IEEE_F64BE_g)
const H5T_C_S1            = _read_const(:H5T_C_S1_g)
const H5T_STD_REF_OBJ     = _read_const(:H5T_STD_REF_OBJ_g)
const H5T_STD_REF_DSETREG = _read_const(:H5T_STD_REF_DSETREG_g)
# Native types
const H5T_NATIVE_B8       = _read_const(:H5T_NATIVE_B8_g)
const H5T_NATIVE_INT8     = _read_const(:H5T_NATIVE_INT8_g)
const H5T_NATIVE_UINT8    = _read_const(:H5T_NATIVE_UINT8_g)
const H5T_NATIVE_INT16    = _read_const(:H5T_NATIVE_INT16_g)
const H5T_NATIVE_UINT16   = _read_const(:H5T_NATIVE_UINT16_g)
const H5T_NATIVE_INT32    = _read_const(:H5T_NATIVE_INT32_g)
const H5T_NATIVE_UINT32   = _read_const(:H5T_NATIVE_UINT32_g)
const H5T_NATIVE_INT64    = _read_const(:H5T_NATIVE_INT64_g)
const H5T_NATIVE_UINT64   = _read_const(:H5T_NATIVE_UINT64_g)
const H5T_NATIVE_FLOAT    = _read_const(:H5T_NATIVE_FLOAT_g)
const H5T_NATIVE_DOUBLE   = _read_const(:H5T_NATIVE_DOUBLE_g)
# Other type constants
const H5T_VARIABLE = reinterpret(UInt, -1)

# Filter constants
const H5Z_FLAG_OPTIONAL = 0x0001
const H5Z_FLAG_REVERSE = 0x0100
const H5Z_CLASS_T_VERS = 1


text:
@testset "NLopt" begin
######################### integration with NLopt.jl #########################

# ================ standard NLopt ================
function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 0
        grad[2] = 2x[2]
    end
    return x[2]^2
end

function myconstraint(x::Vector, grad::Vector, a, b)
    if length(grad) > 0
        grad[1] = 2*(x[1]-1)
        grad[2] = -1
    end
    (x[1]-1)^2 - x[2]
end

opt = Opt(:LD_MMA, 2)
opt.lower_bounds = [-Inf, 0.0]
opt.xtol_rel = 1e-4

opt.min_objective = myfunc
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)

(minf,minx,ret) = NLopt.optimize(opt, [1.234, 5.678])
# ================ END NLopt ================
p = ones(10)
Con = CustomOptimizer() do f, df, c, dc, x0, args...
    opt = Opt(:LD_MMA, length(x0))
    bd = zeros(length(x0)); bd[end-1:end] = [-Inf, 0.0]
    opt.lower_bounds = bd
    opt.xtol_rel = 1e-4
    opt.min_objective = (x,g)->(g[:]= df(x); return f(x)[1])
    inequality_constraint!(opt, (x,g)->( g[:]= dc(x);c(x)[1]), 1e-8)
    (minf,minx,ret) = NLopt.optimize(opt, x0)
    minx
end

# unfortunately, we must ensure only minimization variables are in the current tensorflow graph
# for NLopt to work properly
reset_default_graph()
x = Variable([1.234; 5.678])
y = Variable([1.0;2.0])
loss = x[2]^2 + sum(y^2)
c1 = (x[1]-1)^2 - x[2]
opt = Con(loss, inequalities=[c1],var_to_bounds=Dict(x=>[-1.,1.]))
sess = Session(); init(sess)
opt.minimize(sess)
xmin = run(sess, x)

@test true
end

@testset "Optim" begin
    ######################### integration with Optim.jl #########################

    NonCon = CustomOptimizer() do f, df, c, dc, x0, args...
        @show f, df, c, dc, x0
        res = Optim.optimize(f, df, x0; inplace = false)
        res.minimizer
    end
    reset_default_graph()
    x = Variable(rand(2))
    f = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    f1(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
    res = Optim.optimize(f1, rand(2))
    opt = NonCon(f)
    sess = Session(); init(sess)
    opt.minimize(sess)
    xmin = run(sess, x)
    @test norm(xmin-res.minimizer)<1e-2

end

@testset "newton raphson" begin
@test_skip begin

    A = rand(10,10)
    A = A'*A + I
    rhs = rand(10)
    function newton_raphson_f(θ, u)
        r = A*u - rhs
        r, constant(A)
    end
    ADCME.options.newton_raphson.verbose = true
    nr = newton_raphson(newton_raphson_f, zeros(10), missing)
    nr = run(sess, nr)
    @test nr.x≈A\rhs

    u0 = rand(10)
    function newton_raphson_f2(θ, u)
        r = u^3+u - u0
        r, spdiag(3u^2+1.0)
    end

    ADCME.options.newton_raphson.verbose = true
    ADCME.options.newton_raphson.tol = 1e-5
    nr = newton_raphson(newton_raphson_f2, rand(10), missing)
    nr = run(sess, nr)
    uval = nr.x
    @test norm(uval.^3+uval-u0)<1e-3

    # least square
    u0 = rand(10)
    rs = rand(10)
    function newton_raphson_f3(θ, u)
        r = [u^2;u] - [rs.^2;rs]
        r, [spdiag(2*u); spdiag(10)]
    end
    nr = newton_raphson(newton_raphson_f3, rand(10), missing)
    nr = run(sess, nr)
    uval = nr.x
    @test norm(uval-rs)<1e-3
end
end

@testset "NonlinearConstrainedProblem" begin
    @test_skip begin
    θ = Variable(1.8*ones(3))
    u0 = ones(3)
    function f1(θ, u)
        r = u+u^3-θ
        A = spdiag(3u^2+1.0)
        r, A
    end
    function L1(u)
        return sum((u-u0)^2)
    end
    l, u, g = NonlinearConstrainedProblem(f1,L1,θ,zeros(3))
    # the negative gradients and vars are (-g, θ)
    opt_ = AdamOptimizer(1e-1).apply_gradients([(g, θ)])
    init(sess)
    for i = 1:500
        run(sess, opt_)
    end
    @test norm(run(sess, θ)-2.0*ones(3))<1e-8

    θ = constant(ones(1))
    u0 = ones(3)
    Random.seed!(233)
    A = round.(rand(3, 3), digits=1)
    t = zeros(3); t[1:2].=1.0
    function f1(θ, u)
        # r = u-sum(θ)^2-t*sum(θ)
        # A = spdiag(length(u))
        # r, A
        # -t*sum(θ)
        A*u-[sum(θ);sum(θ);constant(0.0)], constant(A)
    end
    function L1(u)
        return sum(u)
    end
    close("all")

    function value_and_gradients_function(θ)
        l, u, g = NonlinearConstrainedProblem(f1,L1,θ,zeros(3))
        l, g
    end
    results = BFGS!(value_and_gradients_function, zeros(3))
    u = run(sess, results)
    @test u≈[2.;2.;2.]
    end
end

@testset "Custom BFGS!" begin
    reset_default_graph()
    x = Variable(2ones(10))
    z = Variable(1.0)
    w = Variable(ones(20,30,10))
    loss = sum((x-1.0)^2+z^2+w^2)
    grads = [gradients(loss, x), gradients(loss, z), gradients(loss, w)]
    vars = [x,z, w]
    global sess = Session(); init(sess)
    l = BFGS!(sess, loss, grads, vars)
    @test l[end]<1e-5
end

@testset "var_to_bounds" begin
    x = Variable(2.0)
    loss = x^2
    init(sess)
    BFGS!(sess, loss, bounds=Dict(x=>[1.0,3.0]))
    @test run(sess, x)≈1.0
end

#=
@testset "Ipopt" begin
######################### integration with Ipopt.jl #########################
# https://github.com/jainachin/tf-ipopt/blob/master/examples/hs071.py
IPOPT = CustomOptimizer() do f, df, c, dc, x0, x_L, x_U
    n = length(x0)
    nz = length(dc(x0))
    g_L, g_U = [-2e19;0.0], [0.0;0.0]
    function eval_jac_g(x, mode, rows, cols, values)
        if mode == :Structure
            rows[1] = 1; cols[1] = 1
            rows[2] = 1; cols[2] = 1
            rows[3] = 2; cols[3] = 1
            rows[4] = 2; cols[4] = 1
        else
            values[:]=dc(x)
        end
    end
    prob = Ipopt.createProblem(n, x_L, x_U, 2, g_L, g_U, nz, 0,
            f, (x,g)->(g[:]=c(x)), (x,g)->(g[:]=df(x)), eval_jac_g, nothing)
    addOption(prob, "hessian_approximation", "limited-memory")
    addOption(prob, "max_iter", 100)

    prob.x = x0
    status = Ipopt.solveProblem(prob)
    println(Ipopt.ApplicationReturnStatus[status])
    println(prob.x)
    prob.x
end

reset_default_graph()
x = Variable([1.0;1.0])
x1 = x[1]; x2 = x[2];
loss = x2
g = x1
h = x1*x1 + x2*x2 - 1
opt = IPOPT(loss, inequalities=[g], equalities=[h], var_to_bounds=Dict(x=>(-1.0,1.0)))
sess = Session(); init(sess)
minimize(opt, sess)
=#

@testset "newton_raphson_with_grad" begin
    function f(θ, x)
        x^3 - θ, 3spdiag(x^2)
    end

    θ = [2. .^3;3. ^3; 4. ^3]
    x = newton_raphson_with_grad(f, constant(ones(3)), θ)
    @test run(sess, x)≈[2.;3.;4.]

    θ = constant([2. .^3;3. ^3; 4. ^3])
    x = newton_raphson_with_grad(f, constant(ones(3)), θ)
    @test run(sess, x)≈[2.;3.;4.]

    @test run(sess, gradients(sum(x), θ))≈1/3*[2. .^3;3. ^3; 4. ^3] .^(-2/3)
end

# @testset "Constrained Optimizer" begin
#     reset_default_graph() # this is very important. UnconstrainedOptimizer only works with a fresh session
#     x = Variable(2*ones(10))
#     y = constant(ones(10))
#     loss = sum((y-x)^4)
#     init(sess)
#     uo = UnconstrainedOptimizer(sess, loss)

#     @test getInit(uo) ≈ 2ones(10)
#     @test getLoss(uo, 3*ones(10))≈ 160.0
#     @test getLossAndGrad(uo, 3*ones(10))[2] ≈  [32.0, 32.0, 32.0, 32.0, 32.0, 32.0, 32.0, 32.0, 32.0, 32.0]
#     @test getLossAndGrad(uo, 3*ones(10))[1] ≈  160.0
# end

# @testset "optimzers" begin
#     ad = AndersonAcceleration()
#     @test_nowarn begin
#         apply!(ad, rand(10), rand(10))
#     end
# end


text:
println("Imports")
using AdventOfCode2018
using AdventOfCode2018.Day4

println("getting input")
input = get_aoc_input(4)

println("start stuff")
guard_info = parse_input(input)
sleep_map = build_sleep_map(guard_info)

println("Solution 1: ", part_1(sleep_map))
println("Solution 2: ", part_2(sleep_map))


text:
using Base.Test
using DataFrames
import WCS

@testset "images" begin

    # A lot of tests are in a single function to avoid having to reload
    # the full image multiple times.
    images = SampleData.get_sdss_images(3900, 6, 269)
    cat_entries = SampleData.get_sdss_catalog(3900, 6, 269)

    ea = make_elbo_args(images, cat_entries, patch_radius_pix=1e-6)

    # Just check some basic facts about the catalog.
    @test length(cat_entries) > 500
    @test 0 < sum([ce.is_star for ce in cat_entries]) < 805

    # Find an object near the middle of the image.
    ctr = WCS.pix_to_world(images[3].wcs, [images[3].H / 2, images[3].W / 2])
    dist = [sqrt((ce.pos[1] - ctr[1])^2 + (ce.pos[2] - ctr[2])^2)
                    for ce in cat_entries]
    obj_index = findmin(dist)[2]    # index of closest object
    obj_loc = cat_entries[obj_index].pos    # location of closest object

    # Test get_source_psf at point while we have the images loaded.
    test_b = 3
    img = ea.images[test_b]
    ea_obj = make_elbo_args(images, cat_entries[obj_index:obj_index])
    vs_obj = DeterministicVI.catalog_init_source(cat_entries[obj_index])
    pixel_loc = WCS.world_to_pix(img.wcs, obj_loc)
    original_psf_val = img.psfmap(pixel_loc[1], pixel_loc[2])

    original_psf_celeste =
        PSF.fit_raw_psf_for_celeste(original_psf_val, ea.psf_K)[1]
    fit_original_psf_val = PSF.get_psf_at_point(original_psf_celeste)

    obj_psf = PSF.get_source_psf(vs_obj[ids.pos], img, ea.psf_K)[1]
    obj_psf_val = PSF.get_psf_at_point(obj_psf)

    # The fits should match exactly.
    @test isapprox(obj_psf_val, fit_original_psf_val, atol=1e-6)

    # The raw psf will not be as good.
    @test isapprox(obj_psf_val, original_psf_val, atol=1e-2)

    cat_several = [cat_entries[1]; cat_entries[obj_index]]
    ea_several = make_elbo_args(ea_obj.images, cat_several)

    # The second set of vp is the object of interest
    point_patch_psf = PSF.get_psf_at_point(ea_several.patches[2, test_b].psf)
    # The threshold for the test below was formerly 1e-6.
    # Is it a problem I needed to increase it?
    @test isapprox(obj_psf_val, point_patch_psf, atol=5e-4)
end


text:
"""
    MulAddStyle

indicates that an `Applied` object should be materialised
via `ArrayLayouts.MulAdd`.
"""
struct MulAddStyle <: AbstractArrayApplyStyle end

ApplyStyle(::typeof(*), ::Type{α}, ::Type{A}, ::Type{B}) where {α<:Number, A<:AbstractMatrix,B<:AbstractVector} =
    MulAddStyle()
ApplyStyle(::typeof(*), ::Type{α}, ::Type{A}, ::Type{B}) where {α<:Number, A<:AbstractMatrix,B<:AbstractMatrix} =
    MulAddStyle()
ApplyStyle(::typeof(*), ::Type{α}, ::Type{A}, ::Type{B}) where {α<:Number, A<:AbstractVector,B<:AbstractMatrix} =
    MulAddStyle()
ApplyStyle(::typeof(+), ::Type{<:Applied{MulAddStyle,typeof(*)}}, ::Type{<:Applied{<:Any,typeof(*)}}) = MulAddStyle() # TODO: simpler second arg
ApplyStyle(::typeof(+), ::Type{<:Applied{MulAddStyle,typeof(*)}}, ::Type{<:AbstractArray}) = MulAddStyle()
ApplyStyle(::typeof(+), ::Type{<:Applied{MulStyle,typeof(*)}}, ::Type{<:Applied{<:Any,typeof(*)}}) = MulAddStyle() # TODO: simpler second arg
ApplyStyle(::typeof(+), ::Type{<:Applied{MulStyle,typeof(*)}}, ::Type{<:AbstractArray}) = MulAddStyle()

_αAB(M::Applied{<:Any,typeof(*),<:Tuple{<:AbstractArray,<:AbstractArray}}, ::Type{T}) where T = tuple(scalarone(T), M.args...)
_αAB(M::Applied{<:Any,typeof(*),<:Tuple{<:Number,<:AbstractArray,<:AbstractArray}}, ::Type{T}) where T = M.args
_αABβC(M::Applied{<:Any,typeof(*)}, ::Type{T}) where T = tuple(_αAB(M, T)..., scalarzero(T), fillzeros(T,axes(M)))

_βC(M::Applied{<:Any,typeof(*)}, ::Type{T}) where T = M.args
_βC(M::AbstractArray, ::Type{T}) where T = (scalarone(T), M)

_αABβC(M::Applied{MulAddStyle,typeof(+)}, ::Type{T}) where T =
    tuple(_αAB(M.args[1], T)..., _βC(M.args[2], T)...)

MulAdd(M::Applied) = MulAdd(_αABβC(M, eltype(M))...)

similar(M::Applied{MulAddStyle}, ::Type{T}) where T = similar(MulAdd(M), T)
copy(M::Applied{MulAddStyle}) = copy(MulAdd(M))
@inline copyto!(dest::AbstractArray, M::Applied{MulAddStyle}) = copyto!(dest, MulAdd(M))


text:
const FATORES_RISCO = [
    :JJ1
    :JM1
    :JM1
    :JM2
    :JM3
    :JM4
    :JM5
    :JM7
    :JM9
    :JT1
    :JT2
    :JT3
    :JT9
    :JI1
    :JI2
    :JI3
    :JI9
    :JP1
    :JP2
    :ME1
    :ME2
    :ME3
    :ME4
    :ME5
    :ME6
    :ME7
    :ME9
    :AA1
    :AA2
    :AA3
    :AA4
    :AA5
    :AA6
    :AA7
    :AA9
    :MC1
    :FF1
    Symbol("998")
    Symbol("999")
]

const LOCAIS_REGISTRO = [:onshore_clearing, :onshore_sem_clearing, :offshore]
const DECODE_LOCAIS_REGISTRO = Dict{Symbol, String}([
    :onshore_clearing     => "01"
    :onshore_sem_clearing => "02"
    :offshore             => "03"
])

const CARTEIRAS = [:trading, :banking]
const DECODE_CARTEIRAS = Dict{Symbol, String}([
    :trading => "01"
    :banking => "02"
])

const VERTICES = [
    Symbol("01")
    Symbol("02")
    Symbol("03")
    Symbol("04")
    Symbol("05")
    Symbol("06")
    Symbol("07")
    Symbol("08")
    Symbol("09")
    Symbol("10")
    Symbol("11")
    Symbol("12")
]

const POSICOES = [:C, :V, :NA]

const CONTAS_ATIVO   = [:A10, :A20, :A30, :A40, :A50, :A90]
const CONTAS_PASSIVO = [:P10, :P20, :P30, :P40, :P50, :P90]
const CONTAS_DERIVATIVO = [
    :D10
    :D20
    :D30
    :D40
    :D41
    :D42
    :D43
    :D50
    :D51
    :D52
    :D53
    :D60
    :D90
]
const CONTAS_ATIVIDADE_FINANCEIRA = [:AFC]

const CONTAS = [CONTAS_ATIVO; CONTAS_PASSIVO; CONTAS_DERIVATIVO; CONTAS_ATIVIDADE_FINANCEIRA]

const TIPOS_ARQUIVO = [
    :I # inclusao
    :S # substituicao
]


text:
using Test
include(joinpath(@__DIR__, "../src/example_data.jl"))

@testset "rows" begin
    @test live_usb_devices !== []
end # rows

@testset "path" begin
    @test TableDirectories.path(partition) == "/sys/dev/block/8:1"
end # path

@testset "Cell.path" begin
    @test TableDirectories.path(file_value) == "/sys/dev/block/8:1/power"
end # Cell.path

@testset "cell_value" begin
    @test sum(ismissing.(data_frame[!, "partition"])) > 0
end # cell_value

@testset "subpath" begin
    #@test subpath(cell, basename(cell_path)) == "power/async"
end # subpath

@testset "column_name_from_path" begin
    #@test column_name_from_path(cell_path).subpath == "power"
end # column_name_from_path

@testset "find_all" begin
    # subsystems = filter(TableDirectories.find_all("/sys/devices", "subsystem")) do cell
    #     !isnothing(cell)
    # end
    # @test subsystems == [] broken=true
end # column_name_from_path

@testset "file_tree_roots()" begin
    #@test @subset(file_tree_roots(), :partition .!== missing) == @subset(file_tree_roots(), :bdi .== missing)
    @test TableDirectories.column_names(partition, "holders") == []
end

@testset "column_names" begin
    subpath = "power"
    cell_path = TableDirectories.path(file_value)
    #@test Main.cell.file_tree_root == Main.partition
    #@test cell.subpath == subpath
    @test !islink(cell_path)
    @test !isfile(cell_path)
    @test isdir(cell_path)
    @test basename(cell_path) == "power"
    #@test column_name_from_path(sub_column_path, parent_column.file_tree_root).subpath == "power/async"
    #@test column_names(partition, sub_column_name(cell, basename(cell_path)))

    path = joinpath(TableDirectories.path(partition), subpath)
    @test path == "/sys/dev/block/8:1/power"
    column_paths = readdir(path; join = true)
    #@test column_paths == ["/sys/dev/block/8:18/power/async", "/sys/dev/block/8:18/power/autosuspend_delay_ms", "/sys/dev/block/8:18/power/control", "/sys/dev/block/8:18/power/runtime_active_kids", "/sys/dev/block/8:18/power/runtime_active_time", "/sys/dev/block/8:18/power/runtime_enabled", "/sys/dev/block/8:18/power/runtime_status", "/sys/dev/block/8:18/power/runtime_suspended_time", "/sys/dev/block/8:18/power/runtime_usage"]
    #paths = paths(cell)
    #@test paths == ["/sys/dev/block/8:18/power/async", "/sys/dev/block/8:18/power/autosuspend_delay_ms", "/sys/dev/block/8:18/power/control", "/sys/dev/block/8:18/power/runtime_active_kids", "/sys/dev/block/8:18/power/runtime_active_time", "/sys/dev/block/8:18/power/runtime_enabled", "/sys/dev/block/8:18/power/runtime_status", "/sys/dev/block/8:18/power/runtime_suspended_time", "/sys/dev/block/8:18/power/runtime_usage"]
    #subpath * "/" .* basename.(cell_path)
    @test TableDirectories.column_names(partition, "power") == ["power/async", "power/autosuspend_delay_ms", "power/control", "power/runtime_active_kids", "power/runtime_active_time", "power/runtime_enabled", "power/runtime_status", "power/runtime_suspended_time", "power/runtime_usage"]
    #TableDirectories.column_names(partition, "")
    TableDirectories.column_names(partition)
end # column_names

@testset "dict" begin
end # dict

@testset "data_frame" begin
			#TableDirectories.data_frame()
end # data_frame

@testset "find" begin
    @test process_exited(run(pipeline(`find --help`; stdout=devnull)))
    #@test process_exited(run(`find --help`; wait=false))
    # subsystems = map(readlines(`find /sys -name subsystem -ls`)) do line
    #     link = line[60:end]
    #     split(link, " -> ")
    # end
end