text:
#=
Part of this work is derived from astropy/photutils. The relevant derivations
are considered under a BSD 3-clause license. =#

"""
    EllipticalAperture(x, y, a, b, θ)
    EllipticalAperture(position, a, b, θ)

An elliptical aperture with semi-major axis `a`, semi-minor axis `b`, and position angle `θ`. `a` and `b` must be ≥ 0, `θ` is measured in degrees counter-clockwise the standard x-axis.

# Examples
```jldoctest
julia> ap = EllipticalAperture(0, 0, 4, 2, 35)
EllipticalAperture(0, 0, a=4, b=2, θ=35°)
```
"""
struct EllipticalAperture{T <: Number} <: AbstractAperture
    x::T
    y::T
    a::T
    b::T
    theta::T

    function EllipticalAperture(x::T, y::T, a::T, b::T, theta::T) where T <: Number
        a < 0 && error("Invalid axis a=$a. a must be greater than or equal to 0")
        b < 0 && error("Invalid axis b=$b. b must be greater than or equal to 0")
        new{T}(x, y, a, b, mod(theta, 360))
    end
end

EllipticalAperture(x, y, a, b, theta) = EllipticalAperture(promote(x, y, a, b, theta)...)
EllipticalAperture(center::AbstractVector, a, b, theta) = EllipticalAperture(center..., a, b, theta)

function Base.show(io::IO, e::EllipticalAperture)
    print(io, "EllipticalAperture($(e.x), $(e.y), a=$(e.a), b=$(e.b), θ=$(e.theta)°)")
end

function oblique_coefficients(a, b, theta)
    sintheta, costheta = sincos(deg2rad(theta))
    a2 = a^2
    b2 = b^2
    cxx = costheta^2 / a2 + sintheta^2 / b2
    cyy = sintheta^2 / a2 + costheta^2 / b2
    cxy = 2 * costheta * sintheta * (1 / a2 - 1 / b2)
    return cxx, cyy, cxy
end

function bbox(e::EllipticalAperture)

    sintheta, costheta = sincos(deg2rad(e.theta))

    t = atan((-e.b * tand(e.theta)) / e.a)

    sint, cost = sincos(t)
    xmin = e.x + e.a * cost * costheta - e.b * sint * sintheta
    xmax = xmin

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        xmin = min(xmin, e.x + e.a * cost * costheta - e.b * sint * sintheta)
    end

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        xmax = max(xmax, e.x + e.a * cost * costheta - e.b * sint * sintheta)
    end

    t = atan((e.b * cotd(e.theta)) / e.a)

    sint, cost = sincos(t)
    ymin = e.y + e.b * sint * costheta + e.a * cost * sintheta
    ymax = ymin

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        ymin = min(ymin, e.y + e.b * sint * costheta + e.a * cost * sintheta)
    end

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        ymax = max(ymax, e.y + e.b * sint * costheta + e.a * cost * sintheta)
    end

    xmin = floor(Int, xmin)
    xmax = ceil(Int, xmax)
    ymin = floor(Int, ymin)
    ymax = ceil(Int, ymax)

    return xmin, xmax, ymin, ymax
end

function mask(e::EllipticalAperture; method = :center)
    bounds = edges(e)
    ny, nx = size(e)
    return elliptical_overlap(bounds..., nx, ny, e.a, e.b, e.theta, method = method)
end

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

"""
    EllipticalAnnulus(x, y, a_in, a_out, b_out, θ)
    EllipticalAnnulus(position, a_in, a_out, b_out, θ)
An elliptical annulus with inner semi-major axis `a_in`, outer semi-major axis `a_out`, outer semi-minor axis `b_out`, and position angle `θ`.
`a_out` ≥ `a_in` ≥ 0 and `b_out` must be ≥ 0, `θ` is measured in degrees counter-clockwise the standard x-axis.

`b_in` will automatically be calculated from `(a_in / a_out) * b_out`. Note this may cause a type instability.
# Examples
```jldoctest
julia> ap = EllipticalAnnulus(0, 0, 4, 10, 5, 45)
EllipticalAnnulus(0.0, 0.0, a_in=4.0, a_out=10.0, b_in=2.0, b_out=5.0, θ=45.0°)
```
"""
struct EllipticalAnnulus{T <: Number} <: AbstractAperture
    x::T
    y::T
    a_in::T
    b_in::T
    a_out::T
    b_out::T
    theta::T

    function EllipticalAnnulus(x::T, y::T, a_in::T, b_in::T, a_out::T, b_out::T, theta::T) where T <: Number
        0 ≤ a_in ≤ a_out || error("Invalid axis a_in=$a_in or a_out=$a_out. a_out must be greater than a_in which must be greater than or equal to 0")
        0 ≤ b_in ≤ b_out || error("Invalid axis b_in=$b_in or b_out=$b_out. b_out must be greater than b_in which must be greater than or equal to 0")
        new{T}(x, y, a_in, b_in, a_out, b_out, mod(theta, 360))
    end
end

EllipticalAnnulus(x, y, a_in, a_out, b_out, theta) = EllipticalAnnulus(promote(x, y, a_in, a_in / a_out * b_out, a_out, b_out, theta)...)
EllipticalAnnulus(center::AbstractVector, a_in, a_out, b_out, theta) = EllipticalAnnulus(center..., a_in, a_out, b_out, theta)

function Base.show(io::IO, e::EllipticalAnnulus)
    print(io, "EllipticalAnnulus($(e.x), $(e.y), a_in=$(e.a_in), a_out=$(e.a_out), b_in=$(e.b_in), b_out=$(e.b_out), θ=$(e.theta)°)")
end

function bbox(e::EllipticalAnnulus)

    sintheta, costheta = sincos(deg2rad(e.theta))

    t = atan((-e.b_out * tand(e.theta)) / e.a_out)

    sint, cost = sincos(t)
    xmin = e.x + e.a_out * cost * costheta - e.b_out * sint * sintheta
    xmax = xmin

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        xmin = min(xmin, e.x + e.a_out * cost * costheta - e.b_out * sint * sintheta)
    end

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        xmax = max(xmax, e.x + e.a_out * cost * costheta - e.b_out * sint * sintheta)
    end

    t = atan((e.b_out * cotd(e.theta)) / e.a_out)

    sint, cost = sincos(t)
    ymin = e.y + e.b_out * sint * costheta + e.a_out * cost * sintheta
    ymax = ymin

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        ymin = min(ymin, e.y + e.b_out * sint * costheta + e.a_out * cost * sintheta)
    end

    for n in -2:2
        sint, cost = sincos(t + n * pi)
        ymax = max(ymax, e.y + e.b_out * sint * costheta + e.a_out * cost * sintheta)
    end

    xmin = floor(Int, xmin)
    xmax = ceil(Int, xmax)
    ymin = floor(Int, ymin)
    ymax = ceil(Int, ymax)

    return xmin, xmax, ymin, ymax
end

function mask(e::EllipticalAnnulus; method = :exact)
    bounds = edges(e)
    ny, nx = size(e)
    out = elliptical_overlap(bounds..., nx, ny, e.a_out, e.b_out, e.theta, method = method)
    out .-= elliptical_overlap(bounds..., nx, ny, e.a_in, e.b_in, e.theta, method = method)

    return out
end


text:
"""
    _make_jvp_call(fdm, f, y, xs, ẋs, ignores)

Call `FiniteDifferences.jvp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `y`: The primal output `y=f(xs...)` or at least something of the right type
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ẋs`: The directional derivatives of `xs` w.r.t. some real number `t`.
- `ignores`: Collection of `Bool`s, the same length as `xs` and `ẋs`.
   If `ignores[i] === true`, then `ẋs[i]` is ignored for derivative estimation.

# Returns
- `Ω̇`: Derivative of output w.r.t. `t` estimated by finite differencing.
"""
function _make_jvp_call(fdm, f, y, xs, ẋs, ignores)
    f2 = _wrap_function(f, xs, ignores)

    ignores = collect(ignores)
    all(ignores) && return ntuple(_ -> NoTangent(), length(xs))
    sigargs = zip(xs[.!ignores], ẋs[.!ignores])
    return ProjectTo(y)(jvp(fdm, f2, sigargs...))
end

"""
    _make_j′vp_call(fdm, f, ȳ, xs, ignores) -> Tuple

Call `FiniteDifferences.j′vp`, with the option to ignore certain `xs`.

# Arguments
- `fdm::FiniteDifferenceMethod`: How to numerically differentiate `f`.
- `f`: The function to differentiate.
- `ȳ`: The adjoint w.r.t. output of `f`.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === NoTangent()`.

# Returns
- `∂xs::Tuple`: Derivatives estimated by finite differencing.
"""
function _make_j′vp_call(fdm, f, ȳ, xs, ignores)
    f2 = _wrap_function(f, xs, ignores)

    ignores = collect(ignores)
    args = Any[NoTangent() for _ in 1:length(xs)]
    all(ignores) && return (args...,)
    sigargs = xs[.!ignores]
    arginds = (1:length(xs))[.!ignores]
    fd = j′vp(fdm, f2, ȳ, sigargs...)
    @assert length(fd) == length(arginds)

    for (dx, ind) in zip(fd, arginds)
        args[ind] = ProjectTo(xs[ind])(dx)
    end
    return (args...,)
end

"""
    _wrap_function(f, xs, ignores)

Return a new version of `f`, `fnew`, that ignores some of the arguments `xs`.

# Arguments
- `f`: The function to be wrapped.
- `xs`: Inputs to `f`, such that `y = f(xs...)`.
- `ignores`: Collection of `Bool`s, the same length as `xs`.
  If `ignores[i] === true`, then `xs[i]` is ignored; `∂xs[i] === NoTangent()`.
"""
function _wrap_function(f, xs, ignores)
    function fnew(sigargs...)
        callargs = Any[]
        j = 1

        for (i, (x, ignore)) in enumerate(zip(xs, ignores))
            if ignore
                push!(callargs, x)
            else
                push!(callargs, sigargs[j])
                j += 1
            end
        end
        @assert j == length(sigargs) + 1
        @assert length(callargs) == length(xs)
        return f(callargs...)
    end
    return fnew
end

ext: jl


text:
using ..Math
using ..Internal

import Iterators

tensorise_weights{T}(w::Vector{Vector{T}}) = kron(w...)
tensorise_weights(w::Vector{Any}) = begin; assert(isempty(w)); Vector(0); end

tensorise_nodes{T}(x::Vector{Vector{T}}) = tensorise(x...)
tensorise_nodes(x::Vector{Any}) = begin; assert(isempty(x)); []; end

tensorise_nodes_as_tuples{T}(x::Vector{Vector{T}}) = collect(Iterators.product(x...))
tensorise_nodes_as_tuples(x::Vector{Any}) = begin; assert(isempty(x)); []; end

tensorise_rules{T,S}(rules::Vector{Tuple{T,S}}) =
  (tensorise_nodes([r[1] for r in rules]), tensorise_weights([r[2] for r in rules]))
tensorise_rulesX(rules::Vector{Tuple{Any,Any}}) =
  ([r[1] for r in rules], [r[2] for r in rules])

smolyak_grid(d::Integer, n::Integer) = begin
  I0 = multiindex(d, n-1) + 1
  ordI0 = multiindex_order(I0)
  ind = ordI0.>=n
  I = I0[ind, :]
  ordI = ordI0[ind]
  factor = ord_α -> ((iseven(d+n-1-ord_α) ? 1 : -1) * binomial(d-1,ord_α-n))
  factors = map(factor, ordI)
  (I, factors)
end

collect_rules(I, factors, rules) = begin
  #TODO Make faster by preallocating space for x and w
  d = size(I, 2)
  n = size(I, 1)
  #@assert d==length(rules)
  @assert n==length(factors)
  x = Matrix{Float64}(d, 0)
  w = Vector{Float64}(0)
  for i=1:n
    rule_i = rules[vec(I[i,:])]
    (x_i, w_i) = tensorise_rules(rule_i)
    x = hcat(x, x_i)
    w = vcat(w, factors[i]*w_i)
  end
  (x, w)
end

export
  tensorise_rules,
  smolyak_grid,
  collect_rules


text:
# Parameters
β = 1 / 1.05
P = [0.8 0.2 0.0
     0.0 0.5 0.5
     0.0 0.0 1.0]

# Possible states of the world
# Each column is a state of the world. The rows are [g d b s 1]
x_vals = [0.5 0.5 0.25;
          0.0 0.0  0.0;
          2.2 2.2  2.2;
          0.0 0.0  0.0;
          1.0 1.0  1.0]
Sg = [1.0 0.0 0.0 0.0 0.0]
Sd = [0.0 1.0 0.0 0.0 0.0]
Sb = [0.0 0.0 1.0 0.0 0.0]
Ss = [0.0 0.0 0.0 1.0 0.0]
discrete = true
proc = DiscreteStochProcess(P, x_vals)

econ = Economy(β, Sg, Sd, Sb, Ss, discrete, proc)
T = 15
path = compute_paths(econ, T)

gen_fig_1(path)


text:
@testset "eval_strain_on_gauss_pts1" begin
    m = 20
    n = 10
    h = 0.1
    u3 = zeros((m+1)*(n+1))
    a = zeros(m, n)
    b = zeros(m, n)
    for i = 1:m+1
        for j = 1:n+1
            x = (i-1)*h
            y = (j-1)*h
            u3[ i+(j-1)*(m+1) ] = x^2+y^2
            if i<=m && j<=n
                a[i, j] = 2x
                b[i,j] = 2y
            end
        end
    end
    s = eval_strain_on_gauss_pts1(constant(u3), m, n, h)

    S = run(sess, s)
    A = reshape(S[1:4:end,1], m, n)
    B = reshape(S[1:4:end,2], m, n)

    close("all")
    mesh(A)
    mesh(a, color="orange")

    close("all")
    mesh(B)
    mesh(b, color="orange")

end

@testset "compute_fem_stiffness_matrix1" begin
using Revise
using AdFem
using ADCME
    m = 10
    n = 5
    h = 0.1
    U = zeros((m+1)*(n+1))
    for i = 1:m+1
        for j = 1:n+1
            x = (i-1)*h
            y = (j-1)*h
            U[i+(j-1)*(m+1)] = x^2 + y^2
        end
    end

    hmat = zeros(4m*n, 2, 2)
    for i = 1:4*m*n
        hmat[i,:,:] = [2. 1.;3. 4.]
    end
    A = compute_fem_stiffness_matrix1(constant(hmat), m, n, h)
    A2 = compute_fem_stiffness_matrix1([2. 1.;3. 4.], m, n, h)
    sess = Session()
    run(sess, A) - A2

end


text:
function test_intset_validated_setdiff()
    s = IntSet(1:9)

    @test validated_setdiff!(s, [1,3,5,7,9]) == true
    @test s == IntSet([2,4,6,8])

    @test validated_setdiff!(s, [11,13]) == false
    @test s == IntSet([2,4,6,8])

    @test validated_setdiff!(s, [2,4,6,10]) == false
    @test s == IntSet([2,4,6,8])

    @test validated_setdiff!(s, [2,4,6,8]) == true
    @test s == IntSet([])

    @test validated_setdiff!(s, [1]) == false
    @test s == IntSet([])

    @test validated_setdiff!(s, []) == true
    @test s == IntSet([])

    @test validated_setdiff!(IntSet(1:3), IntSet(2)) == true
end

function test_intset_all()
    print(rpad("DataStructure.IntSet...", 50, ' '))

    test_intset_validated_setdiff()

    println("PASS")
end


text:
# # Intro to observations
#
# This example demonstrates the features of `SyntheticObservations`
# when constructed from  "synthetic observations" generated by an Oceananigans `Simulation`.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add ParameterEstimocean, Oceananigans, CairoMakie"
# ```

# First we load few things

using ParameterEstimocean
using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity
using CairoMakie

# # Generating synthetic observations
#
# We define a utility function for constructing synthetic observations,

default_closure = ConvectiveAdjustmentVerticalDiffusivity(; convective_κz = 1.0,
                                                            convective_νz = 0.9,
                                                            background_κz = 1e-4,
                                                            background_νz = 1e-5)

function generate_synthetic_observations(name = "convective_adjustment"; Nz = 32, Lz = 64,
                                         Qᵇ = +1e-8, Qᵘ = -1e-5, f₀ = 1e-4, N² = 1e-6,
                                         Δt = 10.0, stop_time = 12hours, overwrite=false,
                                         tracers = :b, closure = default_closure)

    data_path = name * ".jld2"

    if isfile(data_path)
        @warn("Using existing data at $data_path. " *
              "Please delete this file if you wish to generate new data.")

        return data_path
    end

    grid = RectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded))
    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(; grid, tracers, closure,
                                          buoyancy = BuoyancyTracer(),
                                          boundary_conditions = (u=u_bcs, b=b_bcs),
                                          coriolis = FPlane(f=f₀))

    set!(model, b = (x, y, z) -> N² * z)
    simulation = Simulation(model; Δt, stop_time)
    init_with_parameters(file, model) = file["parameters"] = (; Qᵇ, Qᵘ, Δt, N², tracers=keys(model.tracers))

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(stop_time/3),
                                                          filename = name,
                                                          array_type = Array{Float64},
                                                          with_halos = true,
                                                          init = init_with_parameters,
                                                          overwrite_existing = true)

    run!(simulation)

    return data_path
end

# and invoke it:

data_path = generate_synthetic_observations()

# # Specifying observations
#
# When synthetic observations are constructed from simulation data, we
# can select
#
# * The fields to include via `field_names`
#
# * Which data in the time-series to include via the `times` keyword.
#   This can be used to change the initial condition for a calibration run.
#
# For example, to build observations with a single field we write,

single_field_observations = SyntheticObservations(data_path, field_names=:b, transformation=ZScore())

# To build observations with two fields we write

two_field_observations = SyntheticObservations(data_path, field_names=(:u, :b), transformation=ZScore())

# And to build observations with specified times we write

times = single_field_observations.times[2:end]
specified_times_observations = SyntheticObservations(data_path, field_names=(:u, :b), transformation=ZScore(), times=times)

# Notice that in the last case, `specified_times_observations.times` is missing `0.0`.

# # Visualizing observations

# For this we include the initial condition and ``v`` velocity component,

observations = SyntheticObservations(data_path, field_names=(:u, :v, :b), transformation=ZScore())

fig = Figure()

ax_b = Axis(fig[1, 1], xlabel = "Buoyancy [m s⁻²]", ylabel = "z [m]")
ax_u = Axis(fig[1, 2], xlabel = "Velocities [m s⁻¹]", ylabel = "z [m]")

z = znodes(Center, observations.grid)

colorcycle = [:black, :red, :blue, :orange, :pink]

for i = 1:length(observations.times)
    b_ = observations.field_time_serieses.b[i]
    u_ = observations.field_time_serieses.u[i]
    v_ = observations.field_time_serieses.v[i]
    t_ = observations.times[i]

    label = "t = " * prettytime(t_)
    u_label = i == 1 ? "u, " * label : label
    v_label = i == 1 ? "v, " * label : label

    lines!(ax_b, 1e4 * interior(b_)[1, 1, :], z; label, color=colorcycle[i]) # convert units from m s⁻² to 10⁻⁴ m s⁻²
    lines!(ax_u, interior(u_)[1, 1, :], z; linestyle=:solid, color=colorcycle[i], label=u_label)
    lines!(ax_u, interior(v_)[1, 1, :], z; linestyle=:dash, color=colorcycle[i], label=v_label)
end

axislegend(ax_b, position=:rb)
axislegend(ax_u, position=:lb, merge=true)

save("intro_to_observations.svg", fig)

# ![](intro_to_observations.svg)

# Hint: if using a REPL or notebook, try
# `using Pkg; Pkg.add("ElectronDisplay"); using ElectronDisplay; display(fig)`
# To see the figure in a window.


text:
module FastBEAST

include("tree/tree.jl")

include("aca.jl")
include("skeletons.jl")
include("utils.jl")
include("hmatrix.jl")
include("beast.jl")

export BoundingBox
export getboxframe
export getchildbox
export whichchildbox

export BoxTreeNode
export create_tree
export BoxTreeOptions
export KMeansTreeOptions
export KMeansTreeNode

export aca_compression

export FullMatrixView, FullMatrixView2, MatrixView, LowRankMatrixView

export HMatrix
export buildhmatrix
export adjoint
export estimate_norm
export estimate_reldifference
export nnz
export compressionrate

export hassemble
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
using Revise
using PyPlot

using CovSel
using ROCKET

function generatePrecisionXia(p)
    D = diagm(0.5 + 2*rand(p))
    sqrtD = sqrt.(D)

    Ωstar = eye(p) + diagm(2*0.6*ones(p-1),1) + diagm(2*0.3*ones(p-2),2)
    Ωstar = (Ωstar + Ωstar')/2
    Ω = sqrtD*Ωstar*sqrtD

    return (Ω + Ω')/2
end

p = 500
n = 400
numRep = 500
T = zeros(numRep)

Ω = generatePrecisionXia(p)

for rep=1:numRep
  if mod(rep, 10) == 0
    @show rep
  end
  data = CovSel.generateData(Ω, 2*n)
  X = data[1:n, :]
  Y = data[n+1:end, :]

  methodType = 4
  covarianceType = 4
  ex, vx = teInference(X, 1, 2, methodType, covarianceType)
  ey, vy = teInference(Y, 1, 2, methodType, covarianceType)

  eθ = ex - ey
  stdθ = sqrt( vx + vy)
  T[rep] = eθ / stdθ
end

plt[:hist](T, 50)
plt[:axvline](x=0, linewidth=2, color="red")


text:
export @q_cmd
KDB_HANDLE[] = 0

function kerror()
    e = asarray(ee(K_(C_NULL)))
    unsafe_string(e[])
end
function apply(f, args...)
    x = K(args)
    r = dot_(kpointer(f), kpointer(x))
    if r == C_NULL
        throw(KdbException(kerror()))
    end
    K(r)
end

for T in (K_Lambda, K_Other, K_symbol)
    @eval (f::$T)() = f(nothing)
    @eval (f::$T)(args...) = apply(f, args...)
end

const q_parse = K(k(0, "parse"))
const q_eval = K(k(0, "eval"))

macro q_cmd(s) q_parse(s) end
Base.run(x::K_List, args...) = (r = q_eval(x); length(args) == 0 ? r : r(args...))

const b100 = [  # in alphabetical order
    "aj",
    "aj0",
    "asc",
    "asof",
    "avgs",
    "cols",
    "cross",
    "cut",
    "desc",
    "each",
    "ej",
    "ema",
    "except",
    "fby",
    "fkeys",
    "gtime",
    "iasc",
    "idesc",
    "ij",
    "ijf",
    "inter",
    "keys",
    "lj",
    "ljf",
    "lower",
    "ltrim",
    "mavg",
    "mcount",
    "mdev",
    "med",
    "meta",
    "mmax",
    "mmin",
    "mod",
    "msum",
    "next",
    "over",
    "parse",
    "peach",
    "pj",
    "prior",
    "rand",
    "rank",
    "rotate",
    "rtrim",
    "scan",
    "scov",
    "sdev",
    "set",
    "show",
    "signum",
    "ssr",
    "sublist",
    "sv",
    "svar",
    "tables",
    "til",
    "trim",
    "uj",
    "ujf",
    "ungroup",
    "upper",
    "view",
    "views",
    "vs",
    "wj",
    "wj1",
    "ww",
    "xasc",
    "xbar",
    "xcol",
    "xcols",
    "xdesc",
    "xgroup",
    "xkey",
    "xlog",
    "xprev",
    "xrank",
]

const b101 = [
    # "::",       #  0
    "flip",       #  1 +:
    "neg",        #  2 -:
    "first",      #  3 *:
    "reciprocal", #  4 %:
    "where",      #  5 &:
    "reverse",    #  6 |:
    "null",       #  7 ^:
    "group",      #  8 =:
    "hopen",      #  9 <:
    "hclose",     # 10 >:
    "string",     # 11 $:
    #             # 12 ,:
    "count",      # 13 #:
    "floor",      # 14 _:
    "not",        # 15 ~:
    "key",        # 16 !:
    "distinct",   # 17 ?:
    "type_",      # 18 @:
    "value",      # 19 .:
    "read0",      # 20 0::
    "read1",      # 21 1::
    #             # 22 2::
    "avg",        # 23
    "last",       # 24
    "sum",        # 25
    "prd",        # 26
    "min",        # 27
    "max",        # 28
    "exit",       # 29
    "getenv",     # 30
    "abs",        # 31
    "sqrt",       # 32
    "log",        # 33
    "exp",        # 34
    "sin",        # 35
    "asin",       # 36
    "cos",        # 37
    "acos",       # 38
    "tan",        # 39
    "atan",       # 40
    "enlist",     # 41
    "var",        # 42
    "dev",        # 43
]

const b102 = [
    "and",        #  5 &
    "or",         #  6 |
    # ...
    "mmu",        # 11 $
    # ...
    "lsq",        # 16 !
    "in",         # 23
    "within",     # 24
    "like",       # 25
    "bin",        # 26
    "ss",         # 27
    "insert",     # 28
    "wsum",       # 29
    "wavg",       # 30
    "div",        # 31
    "xexp",       # 32
    "setenv",     # 33
    "binr",       # 34
    "cov",        # 35
    "cor",        # 36
]

const b104 = [
    "md5",
    "attr",
    "upsert",
    "hcount",
    "eval",
    "reval",
]

const res = [b100; b101; b102; b104]

@eval struct _Q
    $([Symbol(x) for x in res]...)
    _Q() = new($([K(k(0, rstrip(x, ['_']))) for x in res]...))
end

# function Base.show(io::IO, x::Union{K_Other,K_Lambda})
#     s = k(0, "{` sv .Q.S[40 80;0;x]}", K_new(x))
#     try
#         write(io, strip(transcode(String, kG(s))))
#     finally
#         r0(s)
#     end
#     nothing
# end
function __init__()
    f = dl(_eval_string_c, 1)
    r0(k(0, "{.J.e::x}", f))
end


text:
"""
    _power(d0::UnivariateDistribution, d1::UnivariateDistribution, alpha::Real, tail::Tail)

Return `power` for which `quantile(d0, alpha) == quantile(d1, beta)` and `beta = 1 - power`.
"""
function _power(d0::UnivariateDistribution, d1::UnivariateDistribution, alpha::Real, tail::Tail)
    right_tail = tail == one_tail ? alpha : alpha / 2
    except_right_tail = 1 - right_tail
    critical_value = quantile(d0, except_right_tail)
    beta = cdf(d1, critical_value)
    return 1 - beta
end

"""
    _alpha(d0::UnivariateDistribution, d1::UnivariateDistribution, power::Real, tail::Tail)

Return `alpha` for which `quantile(d0, alpha) == quantile(d1, beta)` and `beta = 1 - power`.
"""
function _alpha(d0::UnivariateDistribution, d1::UnivariateDistribution, power::Real, tail::Tail)
    beta = 1 - power
    critical_value = quantile(d1, beta)
    except_right_tail = cdf(d0, critical_value)
    right_tail = 1 - except_right_tail
    return tail == one_tail ? right_tail : 2 * right_tail
end

_distribution_parameters(T::IndependentSamplesTTest, n) = n - 2 # n1 + n2 - 2
_distribution_parameters(T::PointBiserialTTest, n) = n - 2
_distribution_parameters(T::TTest, n) = n - 1
function _distribution_parameters(T::DeviationFromZeroMultipleRegression, n)
    return (T.n_predictors, n - T.n_predictors - 1)
end
function _distribution_parameters(T::IncreaseMultipleRegression, n)
    return (T.n_tested_predictors, n - T.n_predictors - 1)
end
_distribution_parameters(T::OneWayANOVA, n) = (T.n_groups - 1, n - T.n_groups)
function _distribution_parameters(T::ConstantVectorHotellingTsqTest, n)
    return (T.n_response_variables, n - T.n_response_variables)
end
function _distribution_parameters(T::TwoVectorsHotellingTsqTest, n::AbstractVector)
    return (T.n_response_variables, sum(n) - T.n_response_variables - 1)
end
_distribution_parameters(T::MultifactorFixedEffectsANOVA, n) = (T.df, n - T.n_groups)
_distribution_parameters(T::ConstantVarianceChisqTest, n) = n - 1
_distribution_parameters(T::ChisqTest, n) = T.df

_noncentrality_parameter(T::IndependentSamplesTTest, es, n) = sqrt(n / 2) * es
_noncentrality_parameter(T::PointBiserialTTest, es, n) = sqrt(es^2 / (1 - es^2)) * sqrt(n)
_noncentrality_parameter(T::TTest, es, n) = sqrt(n) * es
_noncentrality_parameter(T::MultifactorFixedEffectsANOVA, es, n) = es^2 * n
function _noncentrality_parameter(T::TwoVectorsHotellingTsqTest, es, n::AbstractVector)
    return es^2 * ((n[1] * n[2]) / (n[1] + n[2]))
end
_noncentrality_parameter(T::FTest, es, n::AbstractVector) = sum(n) * es^2
_noncentrality_parameter(T::FTest, es, n) = n * es^2
_noncentrality_parameter(T::ChisqTest, es, n) = n * es^2

_distribution(T::TTest) = NoncentralT
_distribution(T::ChisqTest) = NoncentralChisq
_distribution(T::FTest) = NoncentralF

_tail(T::TTest) = T.tail
_tail(T::StatisticalTest) = one_tail

"""
    _null_distribution(T::StatisticalTest, v)

Return the distribution for the null hypothesis of test `T` with distribution parameters `v`.
"""
function _null_distribution(T::StatisticalTest, v)
    D = _distribution(T)
    d0 = D(v..., 0)
    return d0
end

"""
    _alternative_distribution(T::StatisticalTest, v; es, n)

Return the distribution for the null hypothesis of test `T` with distribution parameters `v`.
"""
function _alternative_distribution(T::StatisticalTest, v, es, n)
    D = _distribution(T)
    λ = _noncentrality_parameter(T, es, n)
    d1 = D(v..., λ)
    return d1
end

function _alternative_distribution(T::ConstantVarianceChisqTest, v, es, n)
    D = _distribution(T)
    d1 = D(v..., 0)
    d2 = v + es * d1
    return d2
end

"""
    get_power(T::StatisticalTest; es::Real, alpha::Real, n)

Return the power for some test `T` with effect size `es`, required significance level `alpha` and sample size `n`.
"""
function get_power(T::StatisticalTest; es::Real, alpha::Real, n)
    v = _distribution_parameters(T, n)
    d0 = _null_distribution(T, v)
    d1 = _alternative_distribution(T, v, es, n)
    return _power(d0, d1, alpha, _tail(T))
end

"""
    get_alpha(T::StatisticalTest; es::Real, power::Real, n)

Return the significance level for some test `T` with effect size `es`, power `power` and sample size `n`.
"""
function get_alpha(T::StatisticalTest; es::Real, power::Real, n)
    v = _distribution_parameters(T, n)
    d0 = _null_distribution(T, v)
    d1 = _alternative_distribution(T, v, es, n)
    return _alpha(d0, d1, power, _tail(T))
end

"""
    get_es(T::StatisticalTest; alpha::Real, power::Real, n)

Return the minimum effect size for some test `T` with significance level `alpha`, power `power` and sample size `n`.
"""
function get_es(T::StatisticalTest; alpha::Real, power::Real, n)
    f(es) = get_alpha(T; es, power, n) - alpha
    initial_value = 0.5
    return find_zero(f, initial_value)
end

"""
    get_n(T::StatisticalTest; alpha::Real, power::Real, es::Real)

Return minimum sample size `n` for some test `T` with significance level `alpha`, power `power` and effect size `es`.
"""
function get_n(T::StatisticalTest; alpha::Real, power::Real, es::Real)
    f(n) = get_alpha(T; es, power, n) - alpha
    initial_value = 50
    return find_zero(f, initial_value)
end


text:
module ActorsWorld

import Base.==
import Base.copy

export
    Location,
    Orientation,
    orientation_rotate,
    Actor,
    World,
    Actor_Definition,
    actor_create!,
    actor_delete!,
    actor_rotate!,
    actor_move!,
    actor_moveto!,
    actor_putdown!,
    actor_pickup!,
    get_actors_at_location,
    location_move,
    location_fix_ooBound,
    location_within_world,
    wrap_actor_move,
    wrap_actor_rotate,
    wrap_actor_putdown,
    wrap_actor_pickup,
    actor_definition_at_location,
    is_actor_definition_left,
    is_actor_definition_right,
    is_actor_definition_front,
    is_actor_definition_here,
    wrap_actor_definition_left,
    wrap_actor_definition_right,
    wrap_actor_definition_front,
    wrap_actor_definition_here,
    orientation_to_rad,
    AbstractActorsWorldException,
    InvalidDirectionError,
    LocationFullError,
    ActorOnWrongLayerError,
    LocationOutsideError,
    ActorNotFound,
    ActorInvalidRotationError,
    ActorInvalidMovementError,
    ActorInvalidMultipleMovementError,
    ActorInvalidGrabError,
    ActorGrabNotFoundError,
    ActorNotPlaceableError,
    WorldCorrupError,
    copy,
    world_state_save,
    reset!,
    World_State

abstract type AbstractActorsWorldException <: Exception end
struct InvalidDirectionError <: AbstractActorsWorldException
    direction::Symbol
end
function Base.showerror(io::IO,e::InvalidDirectionError)
    print(io,"direction $(e.direction) is invalid")
end
struct LocationFullError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::LocationFullError)
    print(io,"can't place more than two actors at one location")
end
struct ActorOnWrongLayerError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorOnWrongLayerError)
    print(io,"two actors on the same lable can't share a cell")
end
struct ActorNotPlaceableError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorNotPlaceableError)
    print(io,"to place an actor it must be exactly one layer below")
end
struct LocationOutsideError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::LocationOutsideError)
    print(io,"location is outside of this world")
end
struct ActorNotFound <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorNotFound)
    print(io,"actor is not in this world")
end
struct ActorInvalidRotationError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorInvalidRotationError)
    print(io,"actor is not turnable")
end
struct ActorInvalidMovementError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorInvalidMovementError)
    print(io,"actor is not moveable")
end
struct ActorInvalidMultipleMovementError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorInvalidMultipleMovementError)
    print(io,"too many actors to move")
end
struct ActorInvalidGrabError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorInvalidGrabError)
    print(io,"actor is not grabable")
end
struct ActorGrabNotFoundError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::ActorGrabNotFoundError)
    print(io,"there is not actor to grab")
end
struct WorldCorrupError <: AbstractActorsWorldException end
function Base.showerror(io::IO,e::WorldCorrupError)
    print(io,"this world is corrupted. This should not happen")
end

const DIRECTIONS = [
    :NORTH,
    :EAST,
    :SOUTH,
    :WEST
]

"""
    Orientation(value::Symbol)

Defines a orientation. Possible values for `value` are
defined in DIRECTIONS.
"""
struct Orientation
    value::Symbol
    Orientation(value) = begin
        value in DIRECTIONS || throw(InvalidDirectionError(value))
        new(value)
    end
end

"""
    orientation_rotate(or::Orientation,::Type{Val{bool}})

Rotates a `Orientation` counter-clockwise for Val{false} and clockwise
for Val{true}. Basically jumps to the next enty in `DIRECTIONS`. The
last jumps to the first and the first to the last.
"""
function orientation_rotate(or::Orientation,::Type{Val{false}})
    if or.value == DIRECTIONS[1]
        return Orientation(DIRECTIONS[4])
    elseif or.value == DIRECTIONS[2]
        return Orientation(DIRECTIONS[1])
    elseif or.value == DIRECTIONS[3]
        return Orientation(DIRECTIONS[2])
    elseif or.value == DIRECTIONS[4]
        return Orientation(DIRECTIONS[3])
    else
        throw(InvalidDirectionError(or.value))
    end
end

function orientation_rotate(or::Orientation,::Type{Val{true}})
    if or.value == DIRECTIONS[1]
        return Orientation(DIRECTIONS[2])
    elseif or.value == DIRECTIONS[2]
        return Orientation(DIRECTIONS[3])
    elseif or.value == DIRECTIONS[3]
        return Orientation(DIRECTIONS[4])
    elseif or.value == DIRECTIONS[4]
        return Orientation(DIRECTIONS[1])
    else
        throw(InvalidDirectionError(or.value))
    end
end

function orientation_to_rad(or::Orientation)
    if or.value == DIRECTIONS[1]
        return π/2
    elseif or.value == DIRECTIONS[2]
        return 0
    elseif or.value == DIRECTIONS[3]
        return 3π/2
    elseif or.value == DIRECTIONS[4]
        return π
    else
        throw(InvalidDirectionError(or.value))
    end
end

"""
    Location(x::Int,y::Int)

Stores a location defined by x and y on a gridded space.
"""
struct Location
    x::Int
    y::Int
end

"""
    location_move(lo::Location,or::Orientation)

Moves one step into the direction defined by the Orientation `or`.
"""
function location_move(lo::Location,or::Orientation)
    if or.value == DIRECTIONS[1]
        return Location(lo.x,lo.y+1)
    elseif or.value == DIRECTIONS[2]
        return Location(lo.x+1,lo.y)
    elseif or.value == DIRECTIONS[3]
        return Location(lo.x,lo.y-1)
    else
        return Location(lo.x-1,lo.y)
    end
end

==(a::Location,b::Location) = a.x == b.x && a.y == b.y

"""
    Size(width::Int,height::Int)

Stores the size of a grid.
"""
struct Size
    width::Int
    height::Int
end

"""
    Actor_Definition(;<keyword arguments>)

Defines the behavior and the constraints of an actor.

# Argmuments
- `moveable::Bool`: Defines the movement of this actor.
- `turnable::Bool`: Defines the rotation of this actor.
- `grabable::Bool`: Defines if the actor can be picked-up and put-down
- `layer::Int`    : Defines the leayer the actor moves on

"""
struct Actor_Definition
    moveable::Bool
    turnable::Bool
    grabable::Bool
    layer::Int
    Actor_Definition(;moveable::Bool=false,
                     turnable::Bool=false,
                     grabable::Bool=false,
                     layer::Int=1
                     ) = new(moveable,turnable,grabable,layer)
end

"""
    Actor(actor_definition::Actor_Definition,location::Location,orientation::Orientation)

Defines the actual actor which is placed on the world.

# Examples
The following creates an actor which can be moved and turned. It is placed at
(0,0) in the world and looks north.

```julia-repl
julia> Actor(
    Actor_Definition(moveable=true,turnable=true),
    Location(0,0),
    Orientation(:NORTH)
)
```
"""
mutable struct Actor
    actor_definition::Actor_Definition
    location::Location
    orientation::Orientation
end

function copy(ac::Actor)
    Actor(
        ac.actor_definition,
        ac.location,
        ac.orientation
    )
end

"""
    World(width::Int,height::Int)

Creates a new world with a given height and a given width.
"""
struct World
    size::Size
    actors::Vector{Actor}
    World(width::Int,height::Int) = new(Size(width,height),Actor[])
end

function copy(wo::World)
    c_wo = World(wo.size.width,wo.size.height)
    for ac in wo.actors
        push!(c_wo.actors,copy(ac))
    end
    return c_wo
end

"""
    actor_create!(wo::World,a_def::Actor_Definition,lo::Location,or::Orientation)

Creates a new actor defined by the actor definition `a_def` at Location `lo`,
oriented in `or`. The actor is added to the world `wo`.

The functions returns the newly generated actor, thus to enable interaction it
should be stored.

# Examples
```julia-repl
julia> wo = World(10,10)
julia> adef = Actor_Definition(
    moveable=true,
    turnable=true
)
julia> ac_new = actor_create(
    wo,adef,
    Location(1,1),Orientation(:NORTH)
)
julia> actor_move!(wo,ac_new,:NORTH)
```
"""
function actor_create!(wo::World,a_def::Actor_Definition,lo::Location,or::Orientation)
    actor_validate_location_move(wo,a_def,lo)
    ac = Actor(a_def,lo,or)
    push!(wo.actors,ac)
    return ac
end

"""
    actor_moveto!(wo::World,ac::Actor,lo::Location)

Moves `ac` to `lo` after validating.
"""
function actor_moveto!(wo::World,ac::Actor,lo::Location)
    actor_validate_location_move(wo,ac.actor_definition,lo)
    ac.location = lo
    return nothing
end

"""
    actor_validate_location_move(wo::World,a_def::Actor_Definition,lo::Location)

Validate if it's possible to place an actor of type `a_def` at `lo`.
"""
function actor_validate_location_move(wo::World,a_def::Actor_Definition,lo::Location)
    # Check if actors already exist at this location
    ac_at_lo = get_actors_at_location_on_layer(wo,lo,a_def.layer)
    if length(ac_at_lo) > 0
        throw(LocationFullError())
    end
    # Check if the position is within the world
    if !location_within_world(wo,lo)
        throw(LocationOutsideError())
    end
    return true
end

"""
    actor_delete!(wo::World,ac::Actor)

Delete the actor `ac` from the World `wo`.
"""
function actor_delete!(wo::World,ac::Actor)
    i = something(findnext(isequal(ac),wo.actors,1),0)
    i != 0 || throw(ActorNotFound())
    deleteat!(wo.actors,i)
end

"""
    get_actors_at_location(wo::World,lo::Location)

Return a list of actors at `lo`. If no actor is at `lo` return `[]`.
"""
function get_actors_at_location(wo::World,lo::Location)
    filter(a->a.location == lo,wo.actors)
end

"""
    get_actors_at_location_on_layer(wo::World,lo::Location,layer::Int)

Return a list of actors at `lo` on `layer`. If no actor is at `lo` return `[]`.
"""
function get_actors_at_location_on_layer(wo::World,lo::Location,layer::Int)
    filter(a->a.location == lo && a.actor_definition.layer == layer,wo.actors)
end

"""
    actor_definition_at_location(wo::World,lo::Location,acd::Actor_Definition)

Checks if an actor of type `acd` is at `lo` in `wo`.
"""
function actor_definition_at_location(wo::World,lo::Location,acd::Actor_Definition)
    for ac in get_actors_at_location(wo,lo)
        if ac.actor_definition == acd
            return true
        end
    end
    return false
end

"""
    actor_rotate!(ac::Actor,direction::Bool)

Rotate an actor `ac` by 1 step counter-clockwise for `false` and clockwise
for `true`.
"""
function actor_rotate!(ac::Actor,direction::Bool)
    ac.actor_definition.turnable || throw(ActorInvalidRotationError())
    ac.orientation = orientation_rotate(ac.orientation,Val{direction})
end
"""
    location_within_world(wo::World,lo::Location)

Check if `lo` is within the bounds of the worlds size.
"""
function location_within_world(wo::World,lo::Location)
    lo.x > 0 && lo.x <= wo.size.width && lo.y > 0 && lo.y <= wo.size.height
end
"""
    location_fix_ooBound(wo::World,lo::Location)

Fix a location `lo` which is out of bounds of the worlds size.
The fix is made such that when leaving the world at one end
the world is entered at the opposit end.
"""
function location_fix_ooBound(wo::World,lo::Location)
    if location_within_world(wo,lo)
        return lo
    else
        x = lo.x; y=lo.y
        if lo.x <= 0
            x = wo.size.width + lo.x
        elseif lo.x > wo.size.width
            x = lo.x - wo.size.width
        end
        if lo.y <= 0
            y = wo.size.height + lo.y
        elseif lo.y > wo.size.height
            y = lo.y - wo.size.height
        end
        return Location(x,y)
    end
end

"""
    actor_move!(wo::World,ac::Actor,direction::Symbol[,parent::Bool])

Move the actor `ac` one step in the direction `direction` with the world `wo`.
The optional attribute `parent` should never be used directly as its purpos
is to only allow the movemnt of one consecutive moveable actor.
It actually stops the movement recursion by switching from true to false, which
only allows one nested layer of recursion.
"""
function actor_move!(wo::World,ac::Actor,direction::Symbol,parent::Bool=true)
    ac.actor_definition.moveable || throw(ActorInvalidMovementError())
    new_lo = location_move(ac.location,Orientation(direction))
    new_lo = location_fix_ooBound(wo,new_lo)

    for a in get_actors_at_location_on_layer(wo,new_lo,ac.actor_definition.layer)
        if a.actor_definition.moveable && !parent
            throw(ActorInvalidMultipleMovementError())
        end
        actor_move!(wo,a,direction,false)
    end
    ac.location = new_lo
end

"""
    actor_pickup!(wo::World,ac::Actor)

Remove an `grabable` actor from the same location `ac` is at.
Only elements one layer beneath the actors can be picked up.
"""
function actor_pickup!(wo::World,ac::Actor)
    actrs = get_actors_at_location_on_layer(wo,ac.location,ac.actor_definition.layer-1)
    if length(actrs) == 0
        throw(ActorGrabNotFoundError())
    end
    if length(actrs) > 1
        throw(WorldCorruptError())
    end
    if !actrs[1].actor_definition.grabable
        throw(ActorInvalidGrabError())
    end
    actor_delete!(wo,actrs[1])
end

"""
    actor_putdown!(wo::Word,ac::Actor,acd_put::Actor_Definition)

Create an actor of type `acd_put` at `ac`'s location with `ac`'s orientation.
Only works if `acd_put` has `grabable=true`.
"""
function actor_putdown!(wo::World,ac::Actor,acd_put::Actor_Definition)
    !acd_put.grabable && throw(ActorInvalidGrabError())
    !(acd_put.layer + 1 == ac.actor_definition.layer) && throw(ActorNotPlaceableError())
    actor_create!(wo,acd_put,ac.location,Orientation(DIRECTIONS[1]))
end

"""
    is_actor_definition_left(wo::World,ac::Actor,acd::Actor_Definition)

Checks if an actor of type `acd` is left of `ac`.
"""
function is_actor_definition_left(wo::World,ac::Actor,acd::Actor_Definition)
    lo_left = location_move(ac.location,orientation_rotate(ac.orientation,Val{false}))
    lo_left = location_fix_ooBound(wo,lo_left)
    return actor_definition_at_location(wo,lo_left,acd)
end

"""
    is_actor_definition_right(wo::World,ac::Actor,acd::Actor_Definition)

Checks if an actor of type `acd` is right of `ac`.
"""
function is_actor_definition_right(wo::World,ac::Actor,acd::Actor_Definition)
    lo_right = location_move(ac.location,orientation_rotate(ac.orientation,Val{true}))
    lo_right = location_fix_ooBound(wo,lo_right)
    return actor_definition_at_location(wo,lo_right,acd)
end

"""
    is_actor_definition_front(wo::World,ac::Actor,acd::Actor_Definition)

Checks if an actor of type `acd` is front of `ac`.
"""
function is_actor_definition_front(wo::World,ac::Actor,acd::Actor_Definition)
    lo_front = location_move(ac.location,ac.orientation)
    lo_front = location_fix_ooBound(wo,lo_front)
    return actor_definition_at_location(wo,lo_front,acd)
end

"""
    is_actor_definition_here(wo::World,ac::Actor,acd::Actor_Definition)

Checks if an actor of type `acd` is here of `ac`.
"""
function is_actor_definition_here(wo::World,ac::Actor,acd::Actor_Definition)
    return actor_definition_at_location(wo,ac.location,acd)
end

struct World_State_Entity
    location::Location
    orientation::Orientation
end

struct World_State
    actors::Dict{Actor,World_State_Entity}
end

"""
    world_state_save(wo::World)

Conerts a world `wo` to a structure that holds all actors with their
current location and orientation.  Can be used with `reset!` to revert
a world to a certain `World_State`.

# Examples
```julia-repl
julia> ws = world_state_save(some_world)
julia> # Do something in some_world
julia> reset!(some_world,ws)` # World is back to the saved state
```
"""
function world_state_save(wo::World)
    World_State(Dict([ac => World_State_Entity(
        ac.location,
        ac.orientation
    ) for ac in wo.actors]))
end

function actor_reset!(ac::Actor,wse::World_State_Entity)
    ac.location = wse.location
    ac.orientation = wse.orientation
    return nothing
end

"""
    reset!(wo::World,state::World_State)

Resets a wo back to a given `state`. The state is obtained from `world_state_save`.

# Examples
```julia-repl
julia> ws = world_state_save(some_world)
julia> # Do something in some_world
julia> reset!(some_world,ws)` # World is back to the saved state
```
"""
function reset!(wo::World,state::World_State)
    # Empty world
    while length(wo.actors) > 0
        pop!(wo.actors)
    end
    # Refill

    for ast in pairs(state.actors)
        actor_reset!(ast.first,ast.second)
        push!(wo.actors,ast.first)
    end
    return nothing
end

end


text:
# ---
# notebook: False
# ---


text:
using RiskAdjustedLinearizations, Test, LinearAlgebra
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "wachter_disaster_risk", "wachter.jl"))

# Solve model
m_wachter = WachterDisasterRisk()
m = inplace_wachter_disaster_risk(m_wachter)
try
    solve!(m, m.z, m.y; verbose = :none)
catch e
    local sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference", "iterative_sss_output.jld2"), "r")
    update!(m, sssout["z"], sssout["y"], sssout["Psi"])
end

# Verify impulse responses with a zero shock is the same as simulate with no shocks
horizon = 100
@testset "Calculate impulse responses for Wachter (2013)" begin

    # No shocks and start from steady state
    state1, jump1 = simulate(m, horizon)
    state2, jump2 = impulse_responses(m, horizon, 1, 0.; deviations = false)
    state3, jump3 = impulse_responses(m, horizon, 2, 0.; deviations = false)
    state4, jump4 = impulse_responses(m, horizon, 3, 0.; deviations = false)

    @test state1 ≈ state2
    @test state1 ≈ state3
    @test state1 ≈ state4
    @test jump1 ≈ jump2
    @test jump1 ≈ jump3
    @test jump1 ≈ jump4
    @test state1 ≈ repeat(m.z, 1, horizon)
    @test jump1 ≈ repeat(m.y, 1, horizon)

    # No shocks but perturb away from steady state
    state1, jump1 = simulate(m, horizon, 1.01 * m.z)
    state2, jump2 = impulse_responses(m, horizon, 1, 0., 1.01 * m.z; deviations = false)
    state3, jump3 = impulse_responses(m, horizon, 2, 0., 1.01 * m.z; deviations = false)
    state4, jump4 = impulse_responses(m, horizon, 3, 0., 1.01 * m.z; deviations = false)

    @test state1 ≈ state2
    @test state1 ≈ state3
    @test state1 ≈ state4
    @test jump1 ≈ jump2
    @test jump1 ≈ jump3
    @test jump1 ≈ jump4
    @test !(state1[:, 2] ≈ m.z)
    @test !(jump1[:, 2] ≈ m.y)

    # Now with shocks, from steady state
    shocks = zeros(3, horizon)
    shocks[2] = -3.

    state1, jump1 = impulse_responses(m, horizon, 2, -3.; deviations = false)
    state2, jump2 = impulse_responses(m, horizon, 2, -3., m.z; deviations = false)
    state3, jump3 = simulate(m, horizon, shocks)
    state4, jump4 = impulse_responses(m, horizon, 1, -3; deviations = false)

    @test state1 ≈ state2
    @test state1 ≈ state3
    @test !(state1 ≈ state4)
    @test jump1 ≈ jump2
    @test jump1 ≈ jump3
    @test !(jump1 ≈ jump4)
end


text:
export radius, midpoint, lbound, ubound, abs_lbound, abs_ubound, getinterval, getball

"""
    radius([T, ] x::ArbOrRef)

Returns the radius of `x` as a `Mag`. If `T` is given convert to this
type, supports `Mag`, `Arf` and `Arb`.
"""
radius(::Type{Mag}, x::ArbOrRef) = Mag(radref(x))
radius(::Type{Arf}, x::ArbOrRef) = Arf(radref(x), prec = precision(x))
radius(::Type{Arb}, x::ArbOrRef) = Arb(radref(x), prec = precision(x))
radius(x::ArbOrRef) = radius(Mag, x)

"""
    midpoint([T, ] x::ArbOrRef)

Returns the midpoint of `x` as an `Arf`. If `T` is given convert to this
type, supports `Arf` and `Arb`.
"""
midpoint(::Type{Arf}, x::ArbOrRef) = Arf(midref(x))
midpoint(::Type{Arb}, x::ArbOrRef) = Arb(midref(x))
midpoint(x::ArbOrRef) = midpoint(Arf, x)

"""
    lbound([T, ] x::ArbOrRef)

Returns a lower bound of `x` as an `Arf`. If `T` is given convert to
this type, supports `Arf` and `Arb`.

If `x` contains `NaN` it returns `NaN`.
"""
lbound(::Type{Arf}, x::ArbOrRef) = get_lbound!(Arf(prec = precision(x)), x)
function lbound(::Type{Arb}, x::ArbOrRef)
    res = zero(x)
    get_lbound!(midref(res), x)
    return res
end
lbound(x::ArbOrRef) = lbound(Arf, x)

"""
    ubound([T, ] x::ArbOrRef)

Returns an upper bound of `x` as an `Arf`. If `T` is given convert to
this type, supports `Arf` and `Arb`.

If `x` contains `NaN` it returns `NaN`.
"""
ubound(::Type{Arf}, x::ArbOrRef) = get_ubound!(Arf(prec = precision(x)), x)
function ubound(::Type{Arb}, x::ArbOrRef)
    res = zero(x)
    get_ubound!(midref(res), x)
    return res
end
ubound(x::ArbOrRef) = ubound(Arf, x)

"""
    abs_lbound([T, ] x::Union{ArbOrRef,AcbOrRef})

Returns a lower bound of `abs(x)` as an `Arf`. If `T` is given convert
to this type, supports `Arf` and `Arb`.

If `x` contains `NaN` it returns `NaN`.
"""
abs_lbound(::Type{Arf}, x::Union{ArbOrRef,AcbOrRef}) =
    get_abs_lbound!(Arf(prec = precision(x)), x)
function abs_lbound(::Type{Arb}, x::Union{ArbOrRef,AcbOrRef})
    res = Arb(prec = precision(x))
    get_abs_lbound!(midref(res), x)
    return res
end
abs_lbound(x::Union{ArbOrRef,AcbOrRef}) = abs_lbound(Arf, x)

"""
    abs_ubound([T, ] x::Union{ArbOrRef,AcbOrRef})

Returns an upper bound of `abs(x)` as an `Arf`. If `T` is given
convert to this type, supports `Arf` and `Arb`.

If `x` contains `NaN` it returns `NaN`.
"""
abs_ubound(::Type{Arf}, x::Union{ArbOrRef,AcbOrRef}) =
    get_abs_ubound!(Arf(prec = precision(x)), x)
function abs_ubound(::Type{Arb}, x::Union{ArbOrRef,AcbOrRef})
    res = Arb(prec = precision(x))
    get_abs_ubound!(midref(res), x)
    return res
end
abs_ubound(x::Union{ArbOrRef,AcbOrRef}) = abs_ubound(Arf, x)

"""
    getinterval([T, ] x::ArbOrRef)

Returns a tuple `(l, u)` representing an interval `[l, u]` enclosing
the ball `x`, both of them are of type `Arf`. If `T` is given convert
to this type, supports `Arf`, `BigFloat` and `Arb`.

If `x` contains `NaN` both `l` and `u` will be `NaN`.
"""
function getinterval(::Type{Arf}, x::ArbOrRef)
    l, u = Arf(prec = precision(x)), Arf(prec = precision(x))
    get_interval!(l, u, x)
    return (l, u)
end
function getinterval(::Type{BigFloat}, x::ArbOrRef)
    l, u = BigFloat(precision = precision(x)), BigFloat(precision = precision(x))
    get_interval!(l, u, x)
    return (l, u)
end
function getinterval(::Type{Arb}, x::ArbOrRef)
    l, u = zero(x), zero(x)
    get_interval!(midref(l), midref(u), x)
    return (l, u)
end
getinterval(x::ArbOrRef) = getinterval(Arf, x)

"""
    getball([T, ] x::ArbOrRef)

Returns a tuple `(m::Arf, r::Mag)` where `m` is the midpoint of the
ball and `r` is the radius. If `T` is given convert both `m` and `r`
to this type, supports `Arb`.
"""
getball(x::ArbOrRef) = (Arf(midref(x)), Mag(radref(x)))
getball(::Type{Arb}, x::ArbOrRef) = (Arb(midref(x)), Arb(radref(x), prec = precision(x)))

"""
    union(x::ArbOrRef, y::ArbOrRef)
    union(x::AcbOrRef, y::AcbOrRef)
    union(x, y, z...)

`union(x, y)` returns a ball containing the union of `x` and `y`.

`union(x, y, z...)` returns a ball containing the union of all given
balls.
"""
Base.union(x::ArbOrRef, y::ArbOrRef) = union!(Arb(prec = _precision((x, y))), x, y)
Base.union(x::AcbOrRef, y::AcbOrRef) = union!(Acb(prec = _precision((x, y))), x, y)
# TODO: Could be optimized, both for performance and enclosure
Base.union(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef, xs...) =
    foldl(union, xs, init = union(union(x, y), z))
Base.union(x::AcbOrRef, y::AcbOrRef, z::AcbOrRef, xs...) =
    foldl(union, xs, init = union(union(x, y), z))

"""
    intersect(x::ArbOrRef, y::ArbOrRef)
    intersect(x::AcbOrRef, y::AcbOrRef)
    intersect(x, y, z...)

`intersect(x, y)` returns a ball containing the intersection of `x`
and `y`. If `x` and `y` do not overlap (as given by `overlaps(a, b)`)
throws an `ArgumentError`.

`intersect(x, y, z...)` returns a ball containing the intersection of
all given balls. If all the balls do not overlap throws an
`ArgumentError`.
"""
function Base.intersect(x::ArbOrRef, y::ArbOrRef)
    overlaps(x, y) ||
        throw(ArgumentError("intersection of non-intersecting balls not allowed"))
    res = Arb(prec = _precision((x, y)))
    intersection!(res, x, y)
    return res
end
# TODO: Could be optimized, both for performance and enclosure
Base.intersect(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef, xs...) =
    foldl(intersect, xs, init = intersect(intersect(x, y), z))


text:
using StackOverflow
using Test

qHolder = getrecentquestionsfortag()
@test_nowarn qHolder[1].answer_count

@test_nowarn getrecentbadgerecipients()


text:
module ZernikePolynomials

export mn2OSA, mn2Noll, OSA2mn, Noll2mn, Noll2OSA, OSA2Noll, Zernike, Zernikecoefficients, evaluateZernike, normalization

# Zernike polynomial-based functions
# Some formulas from Wikipedia, some from: "Standards for Reporting the Optical Aberrations of Eyes", Journal of Refractive Surgery Volume 18 September/October 2002
# Reinier Doelman, 23-12-2018

"""
    mn2OSA(m::Int,n::Int)

Convert the integer pair (n,m) that defines the Zernike polynomial Z_n^m(ρ,θ) to the sequential OSA/ANSI stardard index number.

Returns NaN for invalid integer pairs.

See also: [`mn2Noll`], [`OSA2mn`]

Source: "Standards for Reporting the Optical Aberrations of Eyes", Journal of Refractive Surgery Volume 18 September/October 2002

# Example:
```julia-repl
julia> [mn2OSA(n,m) for n in 0:4, m in -5:5]
```
"""
function mn2OSA(m::Int,n::Int)
  if n < abs(m) || isodd(n-m)
    return NaN
  else
    return Int((1//2)*(n*(n+2)+m))
  end
end

"""
    mn2Noll(m::Int,n::Int)

Convert the integer pair (n,m) that defines the Zernike polynomial Z_n^m(ρ,θ) to the sequential Noll index number.

Returns NaN for invalid integer pairs.

See also: [`mn2OSA`], [`Noll2mn`]

Source: (https://en.wikipedia.org/wiki/Zernike_polynomials)

# Example:
```julia-repl
julia> [mn2Noll(n,m) for n in 0:4, m in -5:5]
```
"""
function mn2Noll(m::Int,n::Int)
  if n < abs(m) || isodd(n-m)
    return NaN
  else
    if m > 0 && (mod(n,4) ∈ (0,1))
      p = 0
    elseif m < 0 && (mod(n,4) ∈ (2,3))
      p = 0
    elseif m ≥ 0 && (mod(n,4) ∈ (2,3))
      p = 1
    elseif m ≤ 0 && (mod(n,4) ∈ (0,1))
      p = 1
    else
      error("Invalid combination of (m,n) in Noll indexing")
    end
    return Int((1//2)*n*(n+1) + abs(m) + p)
  end
end

"""
    OSA2mn(j::Int)

Convert the sequential OSA/ANSI stardard index number j to the integer pair (n,m) that defines the Zernike polynomial Z_n^m(ρ,θ).

See also: [`Noll2mn`], [`mn2OSA`], [`OSA2Noll`]

Source: "Standards for Reporting the Optical Aberrations of Eyes", Journal of Refractive Surgery Volume 18 September/October 2002

# Example:
```julia-repl
julia> [OSA2mn(j) for j in 0:10]
```
"""
function OSA2mn(j::Int)
  n = Int(ceil((-3 + sqrt(9+8j))/2))
  m = 2j-n*(n+2)
  return (Int(m),Int(n))
end

"""
    Noll2mn(j::Int)

Convert the Noll index number j to the integer pair (n,m) that defines the Zernike polynomial Z_n^m(ρ,θ).

See also: [`OSA2mn`], [`mn2Noll`], [`Noll2OSA`]

Source: (https://en.wikipedia.org/wiki/Zernike_polynomials)

# Example:
```julia-repl
julia> [OSA2mn(j) for j in 0:10]
```
"""
function Noll2mn(j::Int)
  n = Int(ceil((-3 + sqrt(1+8j))/2))
  jr = j - Int(n*(n+1)/2)
  if mod(n,4) ∈ (0,1)
    m1 = jr
    m2 = -(jr-1)
    if iseven(n-m1)
      m = m1
    else
      m = m2
    end
  else # mod(n,4) ∈ (2,3)
    m1 = jr-1
    m2 = -(jr)
    if iseven(n-m1)
      m = m1
    else
      m = m2
    end
  end
  return (m,n)
end

"""
    Noll2OSA(j::Int)

Convert the Noll index number j to the OSA/ANSI stardard index number for Zernike polynomials.

See also: [`OSA2Noll`], [`Noll2mn`]

# Example:
```julia-repl
julia> [Noll2OSA(OSA2Noll(j)) for j = 0:10]
```
"""
function Noll2OSA(j::Int)
  return mn2OSA(Noll2mn(j)...)
end

"""
    OSA2Noll(j::Int)

Convert OSA/ANSI stardard index number to the Noll index number for Zernike polynomials.

See also: [`OSA2Noll`], [`Noll2mn`]

# Example:
```julia-repl
julia> [Noll2OSA(OSA2Noll(j)...) for j = 0:10]
```
"""
function OSA2Noll(j::Int)
  return mn2Noll(OSA2mn(j)...)
end

"""
    R(m::Int,n::Int)

Obtain the function ρ -> R_n^|m|(ρ), where R is the radial polynomial in Zernike polynomials

# Example:
```julia-repl
julia> R(1,1)
```
"""
function R(m::Int,n::Int) # radial polynomial
  p = s -> ((-1)^s * factorial(n-s))/ (factorial(s)*factorial(Int(0.5*(n+abs(m))-s))*factorial(Int(0.5*(n-abs(m))-s)))
  f = x -> sum([p(s)*x.^(n-2s) for s in 0:Int((n-abs(m))/2)])
  return f
end

function normalization(m::Int,n::Int) # normalization constant of the zernike polynomial
  δ = (x,y) -> Int(x == y)
  c = sqrt(2(n+1)/(1+δ(m,0)))
  return c
end

"""
    Zernike(m::Int,n::Int;coord=:polar)

Obtain the function (ρ,θ) -> Z_n^m(ρ,θ), where Z is the Zernike polynomial with coefficients n and m. ρ is the radius and θ the angle.

If coord=:cartesian the function is (x,y) -> Z_n^m(x,y) in Cartesian coordinates.

# Example:
```julia-repl
julia> Zernike(1,1)
julia> Zernike(1,1,coord=:polar)
julia> Z = Zernike(1,1,:cartesian)
julia> Z(0.5,0.2)
```
"""
function Zernike(m::Int,n::Int;coord=:polar)
  δ = ρ -> (abs(ρ) <= 1.)
  if m ≥ 0
    Z = (ρ,θ) -> normalization(m,n)*R(m,n)(ρ) * cos(m*θ) * δ(ρ)
  else
    Z = (ρ,θ) -> -normalization(m,n)*R(m,n)(ρ) * sin(m*θ) * δ(ρ)
  end

  if coord == :cartesian
    g = (x,y) -> (sqrt(x.^2 + y.^2), atan(y,x))
    return (x,y) -> Z(g(x,y)...)
  elseif coord == :polar
    return Z
  else
    error("Unrecognized coordinate system")
  end
end

"""
    Zernike(j::Int;index=:OSA,coord=:polar)

Obtain the function (ρ,θ) -> Z_j(ρ,θ), where Z is the Zernike polynomial with sequential index j, according to the indexing by :OSA or :Noll.

If coord=:cartesian the function is (x,y) -> Z_j(x,y) in Cartesian coordinates.

# Example:
```julia-repl
julia> Zernike(1)
julia> Zernike(1,coord=:polar)
julia> Z = Zernike(5,index=:Noll,coord=:cartesian)
julia> Z(0.5,0.2)
```
"""
function Zernike(j::Int;index=:OSA,coord=:polar)
  if index == :OSA
    return Zernike(OSA2mn(j)...,coord=coord)
  elseif index == :Noll
    return Zernike(Noll2mn(j)...,coord=coord)
  else
    error("Unknown Zernike sequential index")
  end
end

"""
    Zernikecoefficients(phase::AbstractArray{Float64,2}, J::Vector{Int}; index=:OSA)

Compute the Zernike coefficients (OSA normalization) that in a least-squares sense optimally describe the phase.
The Zernike polynomials used are specified with an index vector J, according to OSA indexing.

"""
function Zernikecoefficients(phase::AbstractArray{Float64,2}, J::Vector{Int}; index=:OSA)
  s = size(phase)
  X = range(-1.,stop=1,length=s[1])
  Y = range(-1.,stop=1,length=s[2])

  D = [[Zernike(j,coord=:cartesian,index=index)(x,y) for x in X, y in Y] for j in J ]
  G = cat([D[i][:] for i in 1:length(J)]...; dims=2)
  return G\phase[:]
end

"""
    Zernikecoefficients(x::AbstractArray{<: AbstractFloat,1}, phase::AbstractArray{Float64,2}, J::Vector{Int}; index=:OSA)

Compute the Zernike coefficients (OSA normalization) that in a least-squares sense optimally describe the phase.
The Zernike polynomials used are specified with an index vector J, according to OSA indexing.

"""
function Zernikecoefficients(X::AbstractArray{<: AbstractFloat,1}, phase::AbstractArray{Float64,2}, J::Vector{Int}; index=:OSA)
  s = size(phase)
  if !(s[1] == s[2] && s[1] == length(X))
    error("Non-matching size")
  end

  D = [[Zernike(j,coord=:cartesian,index=index)(x,y) for x in X, y in X] for j in J ]
  G = cat([D[i][:] for i in 1:length(J)]...; dims=2)
  return G\phase[:]
end

"""
    evaluateZernike(N::Int, J::Vector{Int}, coefficients::Vector{Float64}; index=:OSA)

Evaluate the Zernike polynomials on an N-by-N grid as specified by the Zernike coefficients of the polynomials J

# Example:
```julia-repl
julia> W = evaluateZernike(64,[5, 6],[0.3, 4.1])
```
"""
function evaluateZernike(N::Int, J::Vector{Int}, coefficients::AbstractArray{Float64,1}; index=:OSA)
  X = range(-1.,stop=1,length=N)
  Y = range(-1.,stop=1,length=N)

  D = [[Zernike(j,coord=:cartesian,index=index)(x,y) for x in X, y in Y] for j in J ]
  return reduce(+,map(*,D,coefficients))
end

"""
    evaluateZernike(x::AbstractArray{<: AbstractFloat,1}, J::Vector{Int}, coefficients::Vector{Float64}; index=:OSA)

Evaluate the Zernike polynomials on a grid as specified by the Zernike coefficients of the polynomials J, on a grid with points in x

# Example:
```julia-repl
julia> W = evaluateZernike(64,[5, 6],[0.3, 4.1])
```
"""
function evaluateZernike(X::AbstractArray{<: AbstractFloat,1}, J::Vector{Int}, coefficients::Vector{Float64}; index=:OSA)
  D = [[Zernike(j,coord=:cartesian,index=index)(x,y) for x in X, y in X] for j in J ]
  return reduce(+,map(*,D,coefficients))
end

function evaluateZernike(n::Int, J::Int, coefficients::Float64; index=:OSA)
  return evaluateZernike(n, [J], [coefficients]; index=index)
end

function evaluateZernike(X::AbstractArray{<: AbstractFloat,1}, J::Int, coefficients::Float64; index=:OSA)
  return evaluateZernike(X, [J], [coefficients]; index=index)
end

end # module


text:
export PIBasis

# ---------------------- Implementation of the PIBasisSpec

"""
`struct PIBasisSpec`
"""
struct PIBasisSpec
   orders::Vector{Int}     # order (length) of ith basis function
   iAA2iA::Matrix{Int}     # where in A can we find the ith basis function
end

==(B1::PIBasisSpec, B2::PIBasisSpec) = _allfieldsequal(B1, B2)

Base.length(spec::PIBasisSpec) = length(spec.orders)

maxcorrorder(spec::PIBasisSpec) = size(spec.iAA2iA, 2)

function _get_pibfcn(spec0, Aspec, vv)
   vv1 = vv[2:end]
   vvnz = vv1[findall(vv1 .!= 0)]
   return (spec0[vv[1]], Aspec[vvnz])
end

function _get_pibfcn(Aspec, vv)
   vvnz = vv[findall(vv .!= 0)]
   return Aspec[vvnz]
end

# TODO: maybe instead of property == nothing, there should be a
#       generic property with no symmetry attached to it.

function PIBasisSpec( basis1p::OneParticleBasis,
                      symgrp::SymmetryGroup,
                      Bsel::DownsetBasisSelector;
                      property = nothing,
                      filterfun = _->true,
                      init1pbasis = true )

   # we initialize the 1p-basis here; to prevent this it must be manually
   # avoided by passing in init1pbasis = false
   if init1pbasis
      init1pspec!(basis1p, Bsel)
   end

   # get the basis spec of the one-particle basis
   #  Aspec[i] described the basis function that will get written into A[i]
   Aspec = get_spec(basis1p)

   # we assume that `Aspec` is sorted by degree, but best to double-check this
   # since the notion of degree used to construct `Aspec` might be different
   # from the one used to construct AAspec.
   if !issorted(Aspec; by = b -> level(b, Bsel, basis1p))
      error("""PIBasisSpec : AAspec construction failed because Aspec is not
               sorted by degree. This could e.g. happen if an incompatible
               notion of degree was used to construct the 1-p basis spec.""")
   end
   # An AA basis function is given by a tuple 𝒗 = vv. Each index 𝒗ᵢ = vv[i]
   # corresponds to the basis function Aspec[𝒗ᵢ] and the tuple
   # 𝒗 = (𝒗₁, ...) to a product basis function
   #   ∏ A_{vₐ}
   tup2b = vv -> _get_pibfcn(Aspec, vv)

   #  degree or level of a basis function ↦ is it admissible?
   admissible = bb -> (level(bb, Bsel, basis1p) <= maxlevel(bb, Bsel, basis1p))

   if property != nothing
      filter1 = bb -> filterfun(bb) && filter(bb, Bsel, basis1p) && filter(property, symgrp, bb)
   else
      filter1 = bb -> filterfun(bb) && filter(bb, Bsel, basis1p)
   end

   # we can now construct the basis specification; the `ordered = true`
   # keyword signifies that this is a permutation-invariant basis
   maxord = maxorder(Bsel)
   AAspec = gensparse(; NU = maxorder(Bsel),
                        tup2b = tup2b,
                        admissible = admissible,
                        ordered = true,
                        maxvv = [length(Aspec) for _=1:maxord],
                        filter = filter1)

   return PIBasisSpec(AAspec)
end

function PIBasisSpec(AAspec)
   orders = zeros(Int, length(AAspec))
   iAA2iA = zeros(Int, (length(AAspec), length(AAspec[1])))
   for (iAA, vv) in enumerate(AAspec)
      # we use reverse because gensparse constructs the indices in
      # ascending order, but we want descending here.
      # (I don't remember why though)
      iAA2iA[iAA, :] .= reverse(vv)
      orders[iAA] = length( findall( vv .!= 0 ) )
   end
   return PIBasisSpec(orders, iAA2iA)
end

get_spec(AAspec::PIBasisSpec, i::Integer) = AAspec.iAA2iA[i, 1:AAspec.orders[i]]

# ------------------ PISpec sparsification

"""
returns a new `PIBasisSpec` constructed from the old one, but keeping only
the basis indices `Ikeep`. This maintains the order of the basis functions.
"""
sparsify(spec::PIBasisSpec, Ikeep::AbstractVector{<: Integer}) =
      PIBasisSpec(spec.orders[Ikeep], spec.iAA2iA[Ikeep, :])

function _fix_A_indices!(spec::PIBasisSpec, new_inds::AbstractVector{<: Integer})
   for iAA = 1:size(spec.iAA2iA, 1)
      for α = 1:spec.orders[iAA]
         vα = spec.iAA2iA[iAA, α]
         new_vα = new_inds[vα]
         @assert new_vα > 0
         spec.iAA2iA[iAA, α] = new_vα
      end
   end
   return nothing
end

# --------------------------------- PIBasis implementation

"""
`mutable struct PIBasis:` implementation of a permutation-invariant
basis based on the density projection trick.

The standard constructor is
```
PIBasis(basis1p, N, D, maxdeg)
```
* `basis1p` : a one-particle basis
* `N` : maximum interaction order
* `D` : an abstract degee specification, e.g., SparsePSHDegree
* `maxdeg` : the maximum polynomial degree as measured by `D`
"""
mutable struct PIBasis{BOP, REAL} <: ACEBasis
   basis1p::BOP             # a one-particle basis
   spec::PIBasisSpec
   real::REAL     # could be `real` or `identity` to keep AA complex
end               # TODO -> chain??

cutoff(basis::PIBasis) = cutoff(basis.basis1p)

==(B1::PIBasis, B2::PIBasis) =
      ( (B1.basis1p == B2.basis1p) &&
        (B1.spec == B2.spec) &&
        (B1.real == B2.real) )

Base.length(basis::PIBasis) = length(basis.spec)

# default symmetry group
PIBasis(basis1p, Bsel::AbstractBasisSelector; kwargs...) =
   PIBasis(basis1p, O3(), Bsel; kwargs...)

PIBasis(basis1p, symgrp, Bsel::AbstractBasisSelector;
        isreal = false, kwargs...) =
   PIBasis(basis1p,
           PIBasisSpec(basis1p, symgrp, Bsel; kwargs...),
           isreal ? Base.real : Base.identity )

get_spec(pibasis::PIBasis) =
   [ get_spec(pibasis, i) for i = 1:length(pibasis) ]

get_spec(pibasis::PIBasis, i::Integer) =
      get_spec.( Ref(pibasis.basis1p), get_spec(pibasis.spec, i) )

setreal(basis::PIBasis, isreal::Bool) =
   PIBasis(basis.basis1p, basis.spec, isreal)

maxcorrorder(basis::PIBasis) = maxcorrorder(basis.spec)

# ------------------ sparsification

function sparsify!(basis::PIBasis, Ikeep::AbstractVector{<: Integer})
   basis.spec = sparsify(basis.spec, Ikeep)
   return basis
end

"""
This should allow the 1p basis to sparsify itself, then feed back to the
pibasis what the correct indices are.
"""
function clean_1pbasis!(basis::PIBasis)
   spec = get_spec(basis)
   B1p = basis.basis1p
   spec1p = NamedTuple[]
   for bb in spec
      append!(spec1p, bb)
   end
   identity.(unique!(spec1p))
   # sparsify the product 1p basis
   _, new_inds = sparsify!(basis.basis1p, spec1p)
   # now fix the indexing of the PIBasis specification
   _fix_A_indices!(basis.spec, new_inds)
   return basis
end

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

_scaling_absvalue(x::Number) = abs(x)
_scaling_absvalue(x::Symbol) = 0

# TODO: this is a hack; cf. #68
function scaling(pibasis::PIBasis, p)
   ww = zeros(Float64, length(pibasis))
   bspec = get_spec(pibasis)
   for i = 1:length(pibasis)
      for b in bspec[i]
         # TODO: revisit how this should be implemented for a general basis
         ww[i] += sum(x -> _scaling_absvalue(x)^p, b)  #  abs.(values(b)).^p
      end
   end
   return ww
end

# -------------------------------------------------
# FIO codes

write_dict(basis::PIBasis) =
   Dict(  "__id__" => "ACE_PIBasis",
         "basis1p" => write_dict(basis.basis1p),
            "spec" => write_dict(basis.spec),
            "real" => basis.real == Base.real ? true : false )

read_dict(::Val{:ACE_PIBasis}, D::Dict) =
   PIBasis( read_dict(D["basis1p"]),
            read_dict(D["spec"]),
            D["real"] ? Base.real : Base.identity )

write_dict(spec::PIBasisSpec) =
   Dict( "__id__" => "ACE_PIBasisSpec",
         "orders" => spec.orders,
         "iAA2iA" => write_dict(spec.iAA2iA) )

read_dict(::Val{:ACE_PIBasisSpec}, D::Dict) =
   PIBasisSpec( D["orders"], read_dict(D["iAA2iA"]) )

# -------------------------------------------------
# Evaluation codes

function evaluate!(AA, basis::PIBasis, config::UConfig)
   A = evaluate(basis.basis1p, config)
   evaluate!(AA, basis, A)
   release!(A)
   return AA
end

function evaluate!(AA, basis::PIBasis, A::AbstractVector{<: Number})
   fill!(AA, 1)
   for iAA = 1:length(basis)
      aa = one(eltype(A))
      for t = 1:basis.spec.orders[iAA]
         aa *= A[ basis.spec.iAA2iA[ iAA, t ] ]
      end
      AA[iAA] = basis.real(aa)
   end
   return AA
end

_valtype(basis::PIBasis, A::AbstractVector{<: Number}) =
      basis.real(eltype(A))

# draft of defining bases via chains
function evaluate(basis::PIBasis, A::AbstractVector{<: Number})
   VT = _valtype(basis, A)
   AA = Vector{VT}(undef, length(basis))
   evaluate!(AA, basis, A)
   return AA
end

# draft of defining bases via chains
function evaluate(basis::PIBasis, config::UConfig)
   A = evaluate(basis.basis1p, config)
   AA = evaluate(basis, A)
   release!(A)
   return AA
end

# -------------------------------------------------
# gradients

evaluate_d(basis::PIBasis, cfg::UConfig, args...) =
         evaluate_ed(basis, cfg, args...)[2]

function evaluate_ed(basis::PIBasis, cfg::UConfig, args...)
   A, dA = evaluate_ed(basis.basis1p, cfg, args...)
   AA, dAA = evaluate_ed(basis, A, dA)
   release!(A)
   release!(dA)
   return AA, dAA
end

_gradtype(basis::PIBasis, A, dA) =
      basis.real( promote_type(eltype(A), eltype(dA)) )

# this is really an frule I think
function evaluate_ed(basis::PIBasis, A, dA)
   VT = _valtype(basis, A)
   GT = _gradtype(basis, A, dA)
   AA = Vector{VT}(undef, length(basis))
   dAA = Matrix{GT}(undef, length(basis), size(dA, 2))
   evaluate_ed!(AA, dAA, basis, A, dA)
   release!(A)
   release!(dA)
   return AA, dAA
end

function evaluate_ed!(AA, dAA, basis::PIBasis,
                      cfg::UConfig, args...)
   A, dA = evaluate_ed(basis.basis1p, cfg, args...)
   evaluate_ed!(AA, dAA, basis, A, dA)
   release!(dA)
   release!(A)
   return AA, dAA
end

function _AA_local_adjoints!(dAAdA, A, iAA2iA, iAA, ord, _real)
   if ord == 1
      return _AA_local_adjoints_1!(dAAdA, A, iAA2iA, iAA, ord, _real)
   elseif ord == 2
      return _AA_local_adjoints_2!(dAAdA, A, iAA2iA, iAA, ord, _real)
   else
      return _AA_local_adjoints_x!(dAAdA, A, iAA2iA, iAA, ord, _real)
   end
end

function _AA_local_adjoints_1!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @inbounds dAAdA[1] = 1
   @inbounds A1 = A[iAA2iA[iAA, 1]]
   return _real(A1)
end

function _AA_local_adjoints_2!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @inbounds A1 = A[iAA2iA[iAA, 1]]
   @inbounds A2 = A[iAA2iA[iAA, 2]]
   @inbounds dAAdA[1] = A2
   @inbounds dAAdA[2] = A1
   return _real(A1 * A2)
end

function _AA_local_adjoints_x!(dAAdA, A, iAA2iA, iAA, ord, _real)
   @assert length(dAAdA) >= ord
   @assert ord >= 2
   # TODO - optimize a bit more? can move one operation out of the loop
   # Forward pass:
   @inbounds A1 = A[iAA2iA[iAA, 1]]
   @inbounds A2 = A[iAA2iA[iAA, 2]]
   @inbounds dAAdA[1] = 1
   @inbounds dAAdA[2] = A1
   @inbounds AAfwd = A1 * A2
   @inbounds for a = 3:ord-1
      dAAdA[a] = AAfwd
      AAfwd *= A[iAA2iA[iAA, a]]
   end
   @inbounds dAAdA[ord] = AAfwd
   @inbounds Aend = A[iAA2iA[iAA, ord]]
   aa = _real(AAfwd * Aend)
   # backward pass
   @inbounds AAbwd = Aend
   @inbounds for a = ord-1:-1:3
      dAAdA[a] *= AAbwd
      AAbwd *= A[iAA2iA[iAA, a]]
   end
   dAAdA[2] *= AAbwd
   AAbwd *= A2
   dAAdA[1] *= AAbwd

   return aa
end

_acquire_dAAdA!(basis::PIBasis, A) = Vector{eltype(A)}(undef, maxcorrorder(basis))

function evaluate_ed!(AA, dAA, basis::PIBasis,
                      A::AbstractVector, dA::AbstractMatrix)
   dAAdA = _acquire_dAAdA!(basis, A)
   _evaluate_ed!(AA, dAA, basis, A, dA, dAAdA)
end

function _evaluate_ed!(AA, dAA, basis::PIBasis,
                       A::AbstractVector, dA::AbstractMatrix, dAAdA)
   orders = basis.spec.orders
   iAA2iA = basis.spec.iAA2iA

   # Must treat the constants separately. This is not so elegant and could
   # maybe be improved?
   if orders[1] == 0  # SHOULD BE THE ONLY ONE with ord=0!!
      iAAinit = 2
      AA[1] = 1.0
      dAA[1, :] .= Ref(zero(eltype(dAA)))
   else
      iAAinit = 1
   end

   for iAA = iAAinit:length(basis)
      ord = orders[iAA]

      # ----- compute the local adjoints dAA / dA
      # dAAdA[a] ← ∏_{t ≂̸ a} A_{v_t}
      AA[iAA] = _AA_local_adjoints!(dAAdA, A, iAA2iA, iAA, orders[iAA], basis.real)

      # ----- now convert them into dAA / dX
      for j = 1:size(dA, 2)
         dAA[iAA, j] = sum(dAAdA[a] * dA[iAA2iA[iAA, a], j]
                           for a = 1:ord) |> basis.real
      end
   end

   return AA, dAA
end


text:
##
path = @__DIR__()
using DelimitedFiles

##
data = readdlm(joinpath(path,"input.txt"), '\n', Int) |> vec

##
function compute_fuel_requirement(data::Vector)
    sum(data) do x
        compute_fuel_requirement(x)
    end
end
function compute_fuel_requirement(x::Real)
    amount = x ÷ 3 - 2
    amount <= 0 && return 0
    amount + compute_fuel_requirement(amount)
end

##
compute_fuel_requirement(data)


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__
cd(ProjDir) #do

N = 100
b_AU = 0.5
b_AC = 3
b_UX = -1
b_UB = 2
b_CB = -1.5
b_CY = 1
b_XY = 5

df = DataFrame(
  :a => rand(Normal(), N)
);
df[!, :u] = rand(Normal(1, 2), N) + b_AU * df[:, :a]
df[!, :c] = rand(Normal(0, 1), N) + b_AC * df[:, :a]
df[!, :b] = rand(Normal(1, 1), N) + b_UB * df[:, :u]
df[!, :x] = rand(Normal(-2, 1), N) + b_UX * df[:, :u]
df[!, :y] = rand(Normal(1, 2), N) + b_XY * df[:, :x] + b_CY * df[:, :c]

d = OrderedDict(
  :y => [:x],
  [:c, :u] => :a,
  [:y, :b] => :c,
  [:x, :b] => :u
);
u = [:u]

dag = DAG("sr6_4_2a", d; df=df);
show(dag)

fn = joinpath(mktempdir(), "sr6_4_2a.dot")
to_graphviz(dag, fn)
Sys.isapple() && run(`open -a GraphViz.app $(fn)`)

adjustmentsets = adjustment_sets(dag, :x, :y)
println("Adjustment sets:")
adjustmentsets |> display


text:
struct VoxelBinaryLatentLikelihood <: Distribution{Matrix{Float64}} end

const voxel_binary_latent_likelihood = VoxelBinaryLatentLikelihood()

function Gen.logpdf(
        ::VoxelBinaryLatentLikelihood, Y::Matrix{Float64}, X::Matrix{Float64},
        p_outlier::Float64, radius::Float64, bounds::Geometry.Bounds)
    size(X, 1) == 3 || error("X must have size 3 × something, each column represents a point.  Got size(X) = $(size(X))")
    size(Y, 1) == 3 || error("Y must have size 3 × something, each column represents a point.  Got size(Y) = $(size(Y))")
    all((Y .≥ [bounds.xmin, bounds.ymin, bounds.zmin]) .&
        (Y .≤ [bounds.xmax, bounds.ymax, bounds.zmax])
       ) || @warn("Some points in Y were out of bounds")

    tree = get_tree_from_cloud(Y)

    m = size(X, 2)
    n = size(Y, 2)

    # all_idxs is an array of arrays, where the array all_idxs[i] contains all indices j such that
    # the points X[:,i] and Y[:,j] are within a distance of r units of each other.
    all_idxs = NearestNeighbors.inrange(tree, X, r)

    latent_has_corresponding_point_in_obs = fill(false, m)
    obs_has_no_corresponding_latent = fill(true, n)

    for (i,idxs) in enumerate(all_idxs)
        if length(idx) > 0
            latent_has_corresponding_point_in_obs[i] = true
        end
        obs_has_no_corresponding_latent[idx] .= false
    end

    num_latent_cells = 50000

    latent_active = m

    latent_active_fired = sum(latent_has_corresponding_point_in_obs)
    latent_active_not_fired = latent_active - latent_active_fired

    latent_inactive = num_latent_cells - latent_active

    latent_inactive_fired = sum(obs_has_no_corresponding_latent)
    latent_inactive_not_fired = latent_inactive - latent_inactive_fired

    @assert all(map(x -> x>= 0.0, [latent_active_fired, latent_active_not_fired, latent_inactive_fired, latent_inactive_not_fired,
        latent_active, latent_inactive]))
    @assert sum([latent_active_fired, latent_active_not_fired, latent_inactive_fired, latent_inactive_not_fired]) == num_latent_cells

    ACTIVE_FIRE_PROB = 0.999
    INACTIVE_NOT_FIRE_PROB = 0.999

    (
        log(ACTIVE_FIRE_PROB) * latent_active_fired + log(1.0 - ACTIVE_FIRE_PROB) * latent_active_not_fired +
        log(INACTIVE_NOT_FIRE_PROB) * latent_inactive_not_fired + log(1.0 - INACTIVE_NOT_FIRE_PROB) * latent_inactive_fired
    )
end

struct VoxelBinaryLatentLikelihoodMultiCloud <: Distribution{Matrix{Float64}} end
const voxel_binary_latent_likelihood_multi_cloud = VoxelBinaryLatentLikelihoodMultiCloud()

function Gen.logpdf(
        ::VoxelBinaryLatentLikelihoodMultiCloud, Y::Matrix{Float64}, X::Array{Matrix{Float64}},
        p_outlier::Float64, radius::Float64, bounds::Geometry.Bounds)
   log_pdfs = [logpdf(voxel_binary_latent_likelihood, Y, X[i], p_outlier, radius, bounds) for i in 1:length(X)]
   sum(log_pdfs)
end

export uniform_mixture_from_template, uniform_mixture_from_template_multi_cloud


text:
module _Inplace

using Brainstorm._DataStructure: PreAllocatedStack
using Brainstorm._Math._GMP: set!, add!, mul!, div!, neg!, pow!, isqrt!

include("binarysplitting.jl")

const c3_over_24 = Int64(640320)^3 ÷ 24
const digits_per_term = log10(c3_over_24 / (6 * 2 * 6))
const guard_terms = 5

preallocatedstack(size) = PreAllocatedStack(0, [(big(0), big(0), big(0)) for i in 1:size])

# Compute int(pi * 10^digits)
# This is done using Chudnovsky's series with binary splitting
pi(digits::Integer) = begin
    # how many terms to compute
    n = ceil(Int64, digits / digits_per_term) + guard_terms

    depth = 1
    while (1 << depth) <= n; depth += 1; end
    s = preallocatedstack(depth + 1)

    # calculate P(0, N), Q(0, N) and T(0, N)
    binarysplitting!(s, zero(n), n)
    p, q, t = s[1]

    sqrtC = p
    set!(sqrtC, 10)
    pow!(sqrtC, 2*digits) # one_squared
    mul!(sqrtC, 10005)
    isqrt!(sqrtC)

    mul!(q, 426880)
    mul!(q, sqrtC)
    div!(q, t)
end

#
# Source: http://www.craig-wood.com/nick/articles/pi-chudnovsky/
#
# Computes the terms for binary splitting the Chudnovsky infinite series
#
# a(a) = +/- (13591409 + 545140134*a)
# p(a) = (6*a-5)*(2*a-1)*(6*a-1)
# b(a) = 1
# q(a) = a*a*a*C3_OVER_24
#
compute!(s, a, b) = begin
    p, q, t = s[1]

    if a == 0
        set!(p, 1)
        set!(q, 1)
    else
        # p = (6a - 5) * (2a - 1) * (6a - 1)
        set!(p, 6*a - 5)
        mul!(p, 2*a - 1)
        mul!(p, 6*a - 1)

        # q = a^3 * c3_over_24
        set!(q, a)
        pow!(q, 3)
        mul!(q, c3_over_24)
    end

    # t(a) = p(a) * a(a)
    # t = p * (13591409 + 545140134 * a)
    set!(t, a)
    mul!(t, 545140134)
    add!(t, 13591409)
    mul!(t, p)
    isodd(a) && neg!(t)
end

combine!(s) = begin
    p₁, q₁, t₁ = s[1]
    p₂, q₂, t₂ = s[2]

    # t = q₂ * t₁ + p₁ * t₂
    mul!(t₁, q₂)
    mul!(t₂, p₁)
    add!(t₁, t₂)

    # p = p₁ * p₂
    mul!(p₁, p₂)

    # q = q₁ * q₂
    mul!(q₁, q₂)
end

end


text:
function ExpandableODE(arg0::OrdinaryDifferentialEquation)
    return ExpandableODE((OrdinaryDifferentialEquation,), arg0)
end

function add_secondary_equations(obj::ExpandableODE, arg0::SecondaryODE)
    return jcall(obj, "addSecondaryEquations", jint, (SecondaryODE,), arg0)
end

function compute_derivatives(obj::ExpandableODE, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "computeDerivatives", Vector{jdouble}, (jdouble, Vector{jdouble}), arg0, arg1)
end

function get_mapper(obj::ExpandableODE)
    return jcall(obj, "getMapper", EquationsMapper, ())
end

function get_primary(obj::ExpandableODE)
    return jcall(obj, "getPrimary", OrdinaryDifferentialEquation, ())
end

function init(obj::ExpandableODE, arg0::ODEState, arg1::jdouble)
    return jcall(obj, "init", void, (ODEState, jdouble), arg0, arg1)
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using StatsFuns
gkl_gen(Q, P) = sum(xlogy(p, p/q) for (p, q) in zip(P, Q)) + sum(Q) - sum(P)

# %%
n = 10^7
P = rand(n)
Q = abs.(randn(n))
gkl_gen(Q, P), gkl_gen(P, Q)

# %%
gkl_gen(Q, Q)

# %%
gkl_gen(reverse(Q), Q)

# %%
using BenchmarkTools
@benchmark gkl_gen(Q, P)

# %%
using LoopVectorization

unsafe_xlogy(x, y) = x*log(y)
safe_xlogy(x, y) = iszero(x) ? x : x*log(y)

function gkl_lv(Q, P)
    s = sum(Q) - sum(P)
    @tturbo for i in eachindex(P, Q)
        s += unsafe_xlogy(P[i], P[i]/Q[i])
    end
    s
end

gkl_lv(Q, P)

# %%
@benchmark gkl_lv(Q, P)

# %%
Threads.nthreads()

# %%

max_stars_repo_path: 0009/generalized Kullback-Leibler.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0009/generalized Kullback-Leibler.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0009/generalized Kullback-Leibler.jl

max_forks_repo_name: genkuroki/public


text:
using Xpress

function print_help()
    println("""
    Usage
        benchmark.jl [arg] [name]
    [arg]
        --new       Begin a new benchmark comparison
        --compare   Run another benchmark and compare to existing
    [name]          A name for the benchmark test
    Examples
        git checkout master
        julia benchmark.jl --new master
        git checkout approach_1
        julia benchmark.jl --new approach_1
        git checkout approach_2
        julia benchmark.jl --compare master
        julia benchmark.jl --compare approach_1
    """)
end

if length(ARGS) != 2
    print_help()
else
    const Benchmarks = Xpress.MOI.Benchmarks
    const suite = Benchmarks.suite(() -> Xpress.Optimizer())
    if ARGS[1] == "--new"
        Benchmarks.create_baseline(
            suite, ARGS[2]; directory = @__DIR__, verbose = true
        )
    elseif ARGS[1] == "--compare"
        Benchmarks.compare_against_baseline(
            suite, ARGS[2]; directory = @__DIR__, verbose = true
        )
    else
        print_help()
    end
end


text:
abstract type Basis{T} <: LazyQuasiMatrix{T} end
abstract type Weight{T} <: LazyQuasiVector{T} end

const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}

struct WeightLayout <: MemoryLayout end
abstract type AbstractBasisLayout <: MemoryLayout end
struct BasisLayout <: AbstractBasisLayout end
struct SubBasisLayout <: AbstractBasisLayout end
struct MappedBasisLayout <: AbstractBasisLayout end
struct WeightedBasisLayout <: AbstractBasisLayout end

abstract type AbstractAdjointBasisLayout <: MemoryLayout end
struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
struct AdjointSubBasisLayout <: AbstractAdjointBasisLayout end
struct AdjointMappedBasisLayout <: AbstractAdjointBasisLayout end

MemoryLayout(::Type{<:Basis}) = BasisLayout()
MemoryLayout(::Type{<:Weight}) = WeightLayout()

adjointlayout(::Type, ::BasisLayout) = AdjointBasisLayout()
adjointlayout(::Type, ::SubBasisLayout) = AdjointSubBasisLayout()
adjointlayout(::Type, ::MappedBasisLayout) = AdjointMappedBasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::BasisLayout) = WeightedBasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = WeightedBasisLayout()

combine_mul_styles(::AbstractBasisLayout) = LazyQuasiArrayApplyStyle()
combine_mul_styles(::AbstractAdjointBasisLayout) = LazyQuasiArrayApplyStyle()

ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
pinv(J::Basis) = apply(pinv,J)

function ==(A::Basis, B::Basis)
    axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
    false
end

@inline quasildivapplystyle(::AbstractBasisLayout, ::AbstractBasisLayout) = LdivApplyStyle()
@inline quasildivapplystyle(::AbstractBasisLayout, _) = LdivApplyStyle()
@inline quasildivapplystyle(_, ::AbstractBasisLayout) = LdivApplyStyle()

@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

@inline function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)}})
    a,b = arguments(L.B)
    (L.A\a)-(L.A\b)
end

@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

function copy(P::Ldiv{BasisLayout,BasisLayout})
    A, B = P.A, P.B
    A == B || throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
    SquareEye{eltype(P)}(size(A,2))
end
function copy(P::Ldiv{SubBasisLayout,SubBasisLayout})
    A, B = P.A, P.B
    (parent(A) == parent(B) && parentindices(A) == parentindices(B)) ||
        throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
    SquareEye{eltype(P)}(size(A,2))
end

@inline function copy(P::Ldiv{MappedBasisLayout,MappedBasisLayout})
    A, B = P.A, P.B
    demap(A)\demap(B)
end

@inline copy(L::Ldiv{BasisLayout,SubBasisLayout}) = apply(\, L.A, ApplyQuasiArray(L.B))
@inline function copy(L::Ldiv{SubBasisLayout,BasisLayout})
    P = parent(L.A)
    kr, jr = parentindices(L.A)
    lazy_getindex(apply(\, P, L.B), jr, :) # avoid sparse arrays
end

for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
    @eval ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2}) =
        all(parentindices(A) == parentindices(B)) && parent(A) == parent(B)
end

# expansion
_grid(_, P) = error("Overload Grid")
_grid(::MappedBasisLayout, P) = igetindex.(Ref(parentindices(P)[1]), grid(demap(P)))
_grid(::SubBasisLayout, P) = grid(parent(P))
_grid(::WeightedBasisLayout, P) = grid(last(P.args))
grid(P) = _grid(MemoryLayout(typeof(P)), P)

struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
    grid::Grid
    plan::Plan
    iplan::IPlan
end

TransformFactorization(grid, plan) =
    TransformFactorization{promote_type(eltype(grid),eltype(plan)),typeof(grid),typeof(plan),Nothing}(grid, plan, nothing)

TransformFactorization(grid, ::Nothing, iplan) =
    TransformFactorization{promote_type(eltype(grid),eltype(iplan)),typeof(grid),Nothing,typeof(iplan)}(grid, nothing, iplan)

grid(T::TransformFactorization) = T.grid

\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiVector) = a.iplan \  convert(Array, b[a.grid])
\(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])

\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVector) = a.iplan \  b
\(a::TransformFactorization, b::AbstractVector) = a.plan * b

function _factorize(::AbstractBasisLayout, L)
    p = grid(L)
    TransformFactorization(p, nothing, factorize(L[p,:]))
end

struct ProjectionFactorization{T, FAC<:Factorization{T}, INDS} <: Factorization{T}
    F::FAC
    inds::INDS
end

\(a::ProjectionFactorization, b::AbstractQuasiVector) = (a.F \ b)[a.inds]
\(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]

_factorize(::SubBasisLayout, L) = ProjectionFactorization(factorize(parent(L)), parentindices(L)[2])

transform_ldiv(A, B, _) = factorize(A) \ B
transform_ldiv(A, B) = transform_ldiv(A, B, axes(A))

copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

function copy(L::Ldiv{ApplyLayout{typeof(*)},<:AbstractBasisLayout})
    args = arguments(L.A)
    @assert length(args) == 2 # temporary
    apply(\, last(args), apply(\, first(args), L.B))
end

function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix,<:AbstractQuasiVector})
    p,T = factorize(L.A)
    T \ L.B[p]
end

## materialize views

# materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =
#     *(arguments(S)...)

# mass matrix
# y = p(x), dy = p'(x) * dx
# \int_a^b f(y) g(y) dy = \int_{-1}^1 f(p(x))*g(p(x)) * p'(x) dx

_sub_getindex(A, kr, jr) = A[kr, jr]
_sub_getindex(A, ::Slice, ::Slice) = A

function copy(M::QMul2{<:QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}},
                        <:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}})
    Ac, B = M.args
    A = Ac'
    PA,PB = parent(A),parent(B)
    kr,jr = parentindices(B)
    _sub_getindex((PA'PB)/kr.A,parentindices(A)[2],jr)
end

# Differentiation of sub-arrays
function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}}})
    A, B = M.args
    P = parent(B)
    (Derivative(axes(P,1))*P)[parentindices(B)...]
end

function copy(M::QMul2{<:Derivative,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}})
    A, B = M.args
    P = parent(B)
    kr,jr = parentindices(B)
    (Derivative(axes(P,1))*P*kr.A)[kr,jr]
end

function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:AbstractQuasiMatrix})
    args = arguments(L.B)
    # this is a temporary hack
    if args isa Tuple{AbstractQuasiMatrix,Number}
        (L.A \  first(args))*last(args)
    elseif args isa Tuple{Number,AbstractQuasiMatrix}
        first(args)*(L.A \ last(args))
    else
        error("Not implemented")
    end
end

# we represent as a Mul with a banded matrix
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()

@inline sub_materialize(::AbstractBasisLayout, V::AbstractQuasiArray) = V
@inline sub_materialize(::AbstractBasisLayout, V::AbstractArray) = V

demap(x) = x
demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) = parent(V)
function demap(V::SubQuasiArray{<:Any,2})
    kr, jr = parentindices(V)
    demap(parent(V)[kr,:])[:,jr]
end

##
# SubLayout behaves like ApplyLayout{typeof(*)}

combine_mul_styles(::SubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
_arguments(::SubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
call(::SubBasisLayout, ::SubQuasiArray) = *

combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
_arguments(::AdjointSubBasisLayout, A) = _arguments(ApplyLayout{typeof(*)}(), A)
arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
call(::AdjointSubBasisLayout, ::SubQuasiArray) = *

function arguments(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
    A = parent(V)
    _,jr = parentindices(V)
    first(jr) ≥ 1 || throw(BoundsError())
    P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
    A,P
end

####
# sum
####

_sum(V::AbstractQuasiArray, dims) = __sum(MemoryLayout(typeof(V)), V, dims)
_sum(V::AbstractQuasiArray, ::Colon) = __sum(MemoryLayout(typeof(V)), V, :)
sum(V::AbstractQuasiArray; dims=:) = _sum(V, dims)

__sum(L, Vm, _) = error("Override for $L")
function __sum(::SubBasisLayout, Vm, dims)
    @assert dims == 1
    sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
end
function __sum(::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
    a = arguments(V)
    first(apply(*, sum(a[1]; dims=1), tail(a)...))
end

function __sum(::MappedBasisLayout, V::AbstractQuasiArray, dims)
    kr, jr = parentindices(V)
    @assert kr isa AbstractAffineQuasiVector
    sum(demap(V); dims=dims)/kr.A
end

include("splines.jl")


text:
"""
There are extremely high values in the data due to fires, gas flare etc. You may find some values even greater than the aggregate radiance of large cities. Such pixels also have high standard deviation. These pixels may not be of importantance from the point of view of measureming prosperity. The ```outlier_mask``` function generates a mask of pixels with standard deviation less that the 99.9th percentile. Essentially, this function can be used to removed top 1 percent of pixels by standard deviation. A mask can be provided to the function, so that it calculates the percentile based on the lit pixel of the mask. For example, if the datacube is a box around India and the mask is the polygon mask of India, the outlier_mask function will calculate the 99th percentile of the standard deviation of the pixels inside India's boundary.

```julia
outlier_mask(datacube, mask)
```
"""
function outlier_mask(datacube, mask)
    stds = zeros(size(datacube)[1], size(datacube)[2])
    @showprogress for i in 1:size(datacube)[1]
        for j in 1:size(datacube)[2]
            if mask[i, j]==0
                continue
            end
            if count(i->(ismissing(i)),datacube[i, j, :])/length(datacube[i, j, :]) > 0.50  # Don't do anything if there are too many missings
                continue
            end
            stds[i, j] = std(detrend_ts(filter(x -> !ismissing(x), datacube[i, j, :])))
        end
    end
    function std_mask(std, threshold)
        if std < threshold
            return 1
        else
            return 0
        end
    end
    threshold   = quantile(mask_vals(stds, mask), 0.999)
    outlierMask = std_mask.(stds, threshold)
    outlierMask = outlierMask.*mask
    outlierMask = convert(Array{Int8, 2}, outlierMask)
    return outlierMask
end

"""
The time series of a pixel may show a few outliers, but as a whole the pixel may be of importantance in measuring economic activity. The outlier_ts function uses replaces the outlier observations with interpolated values. This is done using the tsclean function the forecast package of R.

```julia
sample_timeseries = datacube[1, 2, :] # The time series of pixel [1, 2]
outlier_ts(sample_timeseries)
```
"""

function outlier_ts(timeseries, window_size = 5, n_sigmas = 3)
    # Credit: https://gist.github.com/erykml/d15525855f2ef455bd7969240f6f4073#file-hampel_filter_forloop-py
    missings = findall(ismissing, timeseries)
    new_series = Array{Union{Float64, Missing}}(filter!(!ismissing, copy(timeseries)))
    n = length(new_series)
    k = 1.4826 # scale factor for Gaussian distribution
    indices = []
    for i in (window_size):(n - window_size)
        x0 = median(new_series[(i - window_size+1):(i + window_size)]) # Use this for the time series outlier package
        # x0 = missing
        S0 = k * median(abs.(new_series[(i - window_size +1):(i + window_size)] .- x0))
        if (abs(new_series[i] - x0) > n_sigmas * S0)
            # new_series[i] = missing # Can use x0
            push!(indices, i)
        end
    end
    new_series[indices] .= missing
    for i in missings #putting back the missing values
        new_series = insert!(new_series, i, missing)
    end
    return new_series
    # Think about how this would behave in the presence of missing.
end


text:
export FixedParameterSystem, fix_parameters

"""
    FixedParameterSystem(F:AbstractSystem, parameters)

Construct a system from the given [`AbstractSystem`](@ref) `F` with the given `parameters`
fixed.
"""
struct FixedParameterSystem{S<:AbstractSystem,T} <: AbstractSystem
    system::S
    parameters::Vector{T}
end
FixedParameterSystem(F::System, p) = FixedParameterSystem(ModelKitSystem(F), p)
Base.size(F::FixedParameterSystem) = size(F.system)

ModelKit.variables(F::FixedParameterSystem) = variables(F.system)
ModelKit.parameters(F::FixedParameterSystem) = Variable[]
ModelKit.variable_groups(F::FixedParameterSystem) = variable_groups(F.system)

(F::FixedParameterSystem)(x, p = nothing) = F.system(x, F.parameters)

evaluate!(u, F::FixedParameterSystem, x, p = nothing) =
    evaluate!(u, F.system, x, F.parameters)
evaluate_and_jacobian!(u, U, F::FixedParameterSystem, x, p = nothing) =
    evaluate_and_jacobian!(u, U, F.system, x, F.parameters)
taylor!(u, v::Val, F::FixedParameterSystem, tx, p = nothing) =
    taylor!(u, v, F.system, tx, F.parameters)

"""
    fix_parameters(F::AbstractSystem, p)

Fix the parameters of the given system `F`. Returns a [`FixedParameterSystem`](@ref).
"""
fix_parameters(F::AbstractSystem, p) = FixedParameterSystem(F, p)


text:
using JuMP
using CPLEX

include("parser.jl")

start = time()

function myModel(MyFileName::String, lambada::Float64)
  n, l, cout, d, Amax, Amin, B = read_instance(MyFileName)

  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m, x[1:(n*l)], Bin)
  @variable(m, y[1:(n*l),1:(n*l)], Bin)

  ## Constraints
  @constraint(m, sum(cout[i]*x[i] for i in 1:(n*l)) <= B) #couts
  @constraint(m, sum(x[i] for i in 1:(n*l)) <= Amax) # aire max
  @constraint(m, sum(x[i] for i in 1:(n*l)) >= Amin) # aire min
  @constraint(m, [i in 1:(n*l), j in 1:(n*l)], y[i,j] <= x[i])
  @constraint(m, [i in 1:(n*l), j in 1:(n*l)], y[i,j] <= x[j])
  @constraint(m, [i in 1:(n*l)], sum(y[i,j] for j in 1:(n*l)) == x[i])
  @constraint(m, [i in 1:(n*l)], y[i,i] == 0)

  ## Objective
  @objective(m, Min, sum(d[i,j]*y[i,j] for i in 1:(n*l), j in 1:(n*l))-lambada*sum(x[i] for i in 1:(n*l)))

  #resolution
  optimize!(m)

  objective = JuMP.objective_value.(m)
  x_opt = JuMP.getvalue.( m[:x] )
  y_opt = JuMP.getvalue.( m[:y] )

  return objective, x_opt, y_opt
end

function Dinkelbach(MyFileName::String)
  n, l, cout, d, Amax, Amin, B = read_instance(MyFileName)
  lambada = 0.5 # step 1
  objective, x_opt, y_opt = myModel(MyFileName, lambada) # step 2 & 3
  while objective > 0 + 1e-5
    lambada = sum(d[i,j]*y_opt[i,j] for i in 1:(n*l), j in 1:(n*l)) / sum(x_opt[i] for i in 1:(n*l))
    objective, x_opt, y_opt = myModel(MyFileName, lambada)
  end
  println("#####################", lambada)
end

Dinkelbach("Min_fragmentation_ampl.dat")

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD