text:
"""
    NETSETTINGS

Generate settings by specifying values for non-default fields. Parsed with
`Parameters.jl`.

# Arguments
- `modelname::String`: name of model (see `src/weights.jl`).
- `test_iter::Int`: Number of independent consecutive Carleo/Troyer runs.
- `repetitions::Int`: Number of wavefunction optimizations.
- `n::Int`: Number of visible neurons.
- `α::Int`: Hidden unit density (No. hidden / No. visible).
- `γ_init::Float`: Learning parameter (default: `0.05`).
- `γ_decay::Float`: Decay of learning parameter (default: `1.0`).
- `mc_trials::Int`: Number of variational derivatives PER repetition PER CPU.
- `mc_steps::Int`: Number of MH steps per sample (default: `n`).
- `regulator::Float`: Regularization parameter of the S-matrix.
- `init_therm_steps::Int`: Burn-in steps on first repetition.
- `therm_steps::Int`: Burn-in on other repetitions.
- `nflips::Int`: Number of flips to be performed per mc_trial (default: `1`).
- `pbc::Bool`: Periodic boundary conditions for lattice (default: `true`).
- `mag0::Bool`: Sample only in sector ∑s_z = 0 (default: `false`).
- `mfree::Bool`: Matrix-free inversion (default: `true`). Overridden to `false`.
    when `iterative_inverse = false`.
- `iterative_inverse::Bool`: Invert S-matrix iteratively (default: `true`).
    Assumed true when `mfree = true`.
- `calc_stat::Bool`: Calculate single-site statistics and correlations after
    last repetition.
- `stat_samples::Int`: No. of samples PER CPU for statistics calculation.
- `save_figures::Bool `: Save statistics in figure (default: `true`).
- `writetofile::Bool`: Write optimized NQS parameters, energy[, statistics] to a
    file.
- `use_meter::Bool`: Use a ProgressMeter bar to visualize progress (default:
    `true`).

# Examples
```julia-repl
julia> NetSettings = NETSETTINGS(
         modelname = \"U_afh\",
         repetitions =1000,
         n = 6,
         α = 3,
         mc_trials = 500,
         init_therm_steps = 100,
         therm_steps = 50,
         stat_samples = 2000)
NETSETTINGS{String,Int64,Float64,Bool}
modelname: String \"U_afh\"
test_iter: Int64 1
repetitions: Int64 1000
n: Int64 6
α: Int64 3
m: Int64 18
dim_rbm: Int64 132
γ_init: Float64 0.05
γ_decay: Float64 1.0
mc_trials: Int64 500
mc_steps: Int64 6
regulator: Float64 0.1
init_therm_steps: Int64 100
therm_steps: Int64 50
nflips: Int64 1
pbc: Bool true
mag0: Bool false
mfree: Bool true
iterative_inverse: Bool true
calc_stat: Bool true
stat_samples: Int64 5000
save_figures: Bool true
writetofile: Bool true
use_meter: Bool true
```
"""
@with_kw struct NETSETTINGS{S<:String, T<:Int, U<:Float64, V<:Bool}
    modelname::S
    test_iter::T = 1            # Number of independent consecutive runs.
    repetitions::T              # Number of wavefunction optimizations.
    n::T                        # Number of visible neurons.
    α::T                        # No. hidden / No. visible.
    m::T = α * n                # Number of hidden neurons.
    dim_rbm::T = m * n + m + n  # Total amount of free parameters.
    γ_init::U = 0.05            # Learning parameter.
    γ_decay::U = 1.0            # Learning parameter at rep. t: γ*(γ_decay)^t.
    mc_trials::T                # Number of var derivatives per repetition.
    mc_steps::T = n             # Number of MH steps per sample.
    regulator::U = 0.1          # Regularization for S-matrix inversion.
    init_therm_steps::T         # Burn-in / thermalization first iter.
    therm_steps::T              # Burn-in / thermalization.
    nflips::T = 1               # Number of flips to be performed per MH step
    pbc::V = true               # Periodic boundary conditions. See also `weights.jl`.
    mag0::V = false             # Sample only in sector ∑s_z = 0.
    mfree::V = true             # Assumed false when iterative_inverse = false.
    iterative_inverse::V = true # Assumed true when mfree = true. true/false when mfree = false.
    calc_stat::V = true         # Calculate statistics
    stat_samples::T             # No. of samples for statistics calculation.
    save_figures::V = true      # Save statistics as figure
    writetofile::V = true       # Write network parameters / energy to a file.
    use_meter::V = true         # Use ProgressMeter.
end

mutable struct WEIGHTS{S<:Array{Float64, 1}, T<:Array{Float64, 2}}
    w_x::S
    w_y::S
    w_z::S
    w_xx::T
    w_yy::T
    w_zz::T
end

mutable struct STATISTICS{S<:Array{Float64, 1},
    T<:Array{Float64, 2}}
    s_x::S
    s_y::S
    s_z::S
    s_xx::T
    s_yy::T
    s_zz::T
end

mutable struct NETPARAMS{S<:SharedArray{Complex{Float64}, 1},
    T<:SharedArray{Complex{Float64}, 2}}
    a::S           # Bias on visible neurons
    b::S           # Hidden neurons
    w::T           # Connection matrix
end

mutable struct NETSTATE{S<:Array{Int, 1},
        T<:Array{Complex{Float64}, 1},
        U<:Complex{Float64}}
    state::S        # State of visible neurons.
    theta::T        # Look-up table.
    theta_upd::T    # Updated look-up table for a MH step.
    logcoshsum::U   # Logarithmic cosh product lookup.
    logcoshsum_upd::U
end

mutable struct VARDERIVS{S<:SharedArray{Complex{Float64}, 3},
    T<:Float64}
    o_tot::S
    reg::T
end

mutable struct OPTPARAMS{S<:Array{Complex{Float64}, 2},
    T<:Array{Complex{Float64}, 1}}
    s::S           # S-matrix
    f::T           # Forces
    dw_tot::T
end

mutable struct OPTPARAMSMFREE{S<:Array{Complex{Float64}, 1}}
    f::S
    dw_tot::S
end

mutable struct CGPARAMS{S<:Array{Complex{Float64}, 1},
        T<:SharedArray{Complex{Float64}, 1},
        U<:SharedArray{Complex{Float64}, 2},
        V<:Float64}
    r::S
    dir::T
    dw::S
    Sdir_w::U
    Sdir::S
    res::V
    res_prev::V
end


text:
function ChainRulesCore.rrule(::typeof(Base.:(/)), A::AbstractVecOrMat, B::Cholesky)
    Y, back = Zygote.pullback((A, U) -> (A / U) / U', A, B.U)
    function rdiv_callback(Ȳ)
        A̅, B̅_factor = back(Ȳ)
        return (
            NoTangent(), A̅, (uplo=NoTangent(), info=NoTangent(), factors=B̅_factor)
        )
    end
    return Y, rdiv_callback
end

function ChainRulesCore.rrule(::typeof(StatsFuns.softmax), x)
    y = StatsFuns.softmax(x)
    function softmax_pullback(Δ)
        out = Δ .* y
        return (NoTangent(), out .-= y .* sum(out))
    end
    return y, softmax_pullback
end

ext: jl


text:
# QuEST_jl/src/base/gates.jl
#

function  collapseToOutcome(qureg        ::QuEST_Types.Qureg,
                            measureQubit ::Integer,
                            outcome      ::Integer)   ::Qreal

    return ccall(:collapseToOutcome, Qreal,
                 (QuEST_Types.Qureg, Cint,         Cint),
                 qureg,              measureQubit, outcome)
end

function measure(qureg        ::QuEST_Types.Qureg,
                 measureQubit ::Integer)            ::Int
    return ccall(:measure, Cint,
                 (QuEST_Types.Qureg, Cint),
                 qureg,              measureQubit)
end

function measureWithStats(qureg             ::QuEST_Types.Qureg,
                          measureQubit      ::Integer)            :: Tuple{Int,Qreal}

    outcomeProb = Ref{Qreal}(-1)

    outcome = ccall(:measureWithStats, Cint,
                    (QuEST_Types.Qureg, Cint,            Ref{Qreal}),
                    qureg,              measureQubit,    outcomeProb)

    return outcome, outcomeProb[]

end

#EOF


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the API function to compose rotations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/compose_rotations.jl
# ================================

# Functions: compose_rotations
# ----------------------------

@testset "Compose rotations (Float64)" begin
    T = Float64

    # DCMs
    # ==========================================================================

    # Sample 4 DCMs.
    for sym in [:D1, :D2, :D3, :D4]
        @eval $sym = angle_to_dcm(
            _rand_ang($T),
            _rand_ang($T),
            _rand_ang($T),
            rand($valid_rot_seqs)
        )
    end

    # Test the function `compose_rotation`.
    @test                D1 === compose_rotation(D1)
    @test           D2 * D1  ≈  compose_rotation(D1,D2)
    @test      D3 * D2 * D1  ≈  compose_rotation(D1,D2,D3)
    @test D4 * D3 * D2 * D1  ≈  compose_rotation(D1,D2,D3,D4)

    # Euler angle and axis
    # ==========================================================================

    # Sample 4 Euler angle and axis.
    for sym in [:ea1, :ea2, :ea3, :ea4]
        α   = _rand_ang(T)
        v   = @SVector randn(T, 3)
        v   = v/norm(v)

        @eval $sym = EulerAngleAxis($α, $v)
    end

    ear1 = ea2 * ea1
    eac1 = compose_rotation(ea1, ea2)

    ear2 = ea3 * ea2 * ea1
    eac2 = compose_rotation(ea1, ea2, ea3)

    ear3 = ea4 * ea3 * ea2 * ea1
    eac3 = compose_rotation(ea1, ea2, ea3, ea4)

    # Test the function `compose_rotation`.
    @test  ea1 === compose_rotation(ea1)
    @test ear1  ≈  eac1
    @test ear2  ≈  eac2
    @test ear3  ≈  eac3

    # Euler angles
    # ==========================================================================

    # Sample 4 Euler angles.
    for sym in [:Θ1, :Θ2, :Θ3, :Θ4]
        θx      = _rand_ang(T)
        θy      = _rand_ang(T)
        θz      = _rand_ang(T)
        rot_seq = Meta.quot(rand(valid_rot_seqs))

        @eval $sym = EulerAngles($θx, $θy, $θz, $rot_seq)
    end

    Θr1 = Θ2 * Θ1
    Θc1 = compose_rotation(Θ1, Θ2)

    Θr2 = Θ3 * Θ2 * Θ1
    Θc2 = compose_rotation(Θ1, Θ2, Θ3)

    Θr3 = Θ4 * Θ3 * Θ2 * Θ1
    Θc3 = compose_rotation(Θ1, Θ2, Θ3, Θ4)

    # Test the function `compose_rotation`.
    @test  Θ1 === compose_rotation(Θ1)
    @test Θr1  ≈  Θc1
    @test Θr2  ≈  Θc2
    @test Θr3  ≈  Θc3

    # Quaternions
    # ==========================================================================

    # Sample 4 quaternions.
    for sym in [:q1, :q2, :q3, :q4]
        @eval $sym = angle_to_quat(
            _rand_ang($T),
            _rand_ang($T),
            _rand_ang($T),
            rand($valid_rot_seqs)
        )
    end

    # Test the function `compose_rotation`.
    @test               q1 === compose_rotation(q1)
    @test       (q1*q2)[:]  ≈  compose_rotation(q1,q2)[:]
    @test    (q1*q2*q3)[:]  ≈  compose_rotation(q1,q2,q3)[:]
    @test (q1*q2*q3*q4)[:]  ≈  compose_rotation(q1,q2,q3,q4)[:]
end

@testset "Compose rotations (Float32)" begin
    T = Float32

    # DCMs
    # ==========================================================================

    # Sample 4 DCMs.
    for sym in [:D1, :D2, :D3, :D4]
        @eval $sym = angle_to_dcm(
            _rand_ang($T),
            _rand_ang($T),
            _rand_ang($T),
            rand($valid_rot_seqs)
        )
    end

    # Test the function `compose_rotation`.
    @test                D1 === compose_rotation(D1)
    @test           D2 * D1  ≈  compose_rotation(D1,D2)
    @test      D3 * D2 * D1  ≈  compose_rotation(D1,D2,D3)
    @test D4 * D3 * D2 * D1  ≈  compose_rotation(D1,D2,D3,D4)

    # Euler angle and axis
    # ==========================================================================

    # Sample 4 Euler angle and axis.
    for sym in [:ea1, :ea2, :ea3, :ea4]
        α   = _rand_ang(T)
        v   = @SVector randn(T, 3)
        v   = v/norm(v)

        @eval $sym = EulerAngleAxis($α, $v)
    end

    ear1 = ea2 * ea1
    eac1 = compose_rotation(ea1, ea2)

    ear2 = ea3 * ea2 * ea1
    eac2 = compose_rotation(ea1, ea2, ea3)

    ear3 = ea4 * ea3 * ea2 * ea1
    eac3 = compose_rotation(ea1, ea2, ea3, ea4)

    # Test the function `compose_rotation`.
    @test  ea1 === compose_rotation(ea1)
    @test ear1  ≈  eac1
    @test ear2  ≈  eac2
    @test ear3  ≈  eac3

    # Euler angles
    # ==========================================================================

    # Sample 4 Euler angles.
    for sym in [:Θ1, :Θ2, :Θ3, :Θ4]
        θx      = _rand_ang(T)
        θy      = _rand_ang(T)
        θz      = _rand_ang(T)
        rot_seq = Meta.quot(rand(valid_rot_seqs))

        @eval $sym = EulerAngles($θx, $θy, $θz, $rot_seq)
    end

    Θr1 = Θ2 * Θ1
    Θc1 = compose_rotation(Θ1, Θ2)

    Θr2 = Θ3 * Θ2 * Θ1
    Θc2 = compose_rotation(Θ1, Θ2, Θ3)

    Θr3 = Θ4 * Θ3 * Θ2 * Θ1
    Θc3 = compose_rotation(Θ1, Θ2, Θ3, Θ4)

    # Test the function `compose_rotation`.
    @test  Θ1 === compose_rotation(Θ1)
    @test Θr1  ≈  Θc1
    @test Θr2  ≈  Θc2
    @test Θr3  ≈  Θc3

    # Quaternions
    # ==========================================================================

    # Sample 4 quaternions.
    for sym in [:q1, :q2, :q3, :q4]
        @eval $sym = angle_to_quat(
            _rand_ang($T),
            _rand_ang($T),
            _rand_ang($T),
            rand($valid_rot_seqs)
        )
    end

    # Test the function `compose_rotation`.
    @test               q1 === compose_rotation(q1)
    @test       (q1*q2)[:]  ≈  compose_rotation(q1,q2)[:]
    @test    (q1*q2*q3)[:]  ≈  compose_rotation(q1,q2,q3)[:]
    @test (q1*q2*q3*q4)[:]  ≈  compose_rotation(q1,q2,q3,q4)[:]
end


text:
# Separable sum, using slices of an array as variables

export SlicedSeparableSum

"""
**Sliced separable sum of functions**

    SlicedSeparableSum((f₁,…,fₖ), (J₁,…,Jₖ))

Returns the function
```math
g(x) = ∑_{i=1}^k f_i(x_{J_i}).
```

    SlicedSeparableSum(f, (J₁,…,Jₖ))

Analogous to the previous one, but applies the same function `f` to all slices of the variable `x`:
```math
g(x) = ∑_{i=1}^k f(x_{J_i}).
```
"""
struct SlicedSeparableSum{S <: Tuple, T <: AbstractArray, N} <: ProximableFunction
  fs::S    # Tuple, where each element is a Vector with elements of the same type; the functions to prox on
  # Example: S = Tuple{Array{ProximalOperators.NormL1{Float64},1}, Array{ProximalOperators.NormL2{Float64},1}}
  idxs::T  # Vector, where each element is a Vector contining the indices to prox on
  # Example: T = Array{Array{Tuple{Colon,UnitRange{Int64}},1},1}
end

function SlicedSeparableSum(fs::S, idxs::T) where {N,
						   S <: Tuple{Vararg{<:ProximableFunction,N}},
						   M,
						   I <: Integer,
						   T1 <: NTuple{M,Union{I,
								       AbstractArray{I},
								       Colon,
								       AbstractRange
								       }
							       },
						   T <:NTuple{N,T1}
						   }
    ftypes = DataType[]
    fsarr = Array{Any,1}[]
    indarr = Array{eltype(idxs),1}[]
    for (i,f) in enumerate(fs)
      t = typeof(f)
      fi = findfirst(isequal(t), ftypes)
      if fi === nothing
        push!(ftypes, t)
        push!(fsarr, Any[f])
        push!(indarr, eltype(idxs)[idxs[i]])
      else
        push!(fsarr[fi], f)
        push!(indarr[fi], idxs[i])
      end
    end
    fsnew = ((Array{typeof(fs[1]),1}(fs) for fs in fsarr)...,)
    @assert typeof(fsnew) == Tuple{(Array{ft,1} for ft in ftypes)...}
    SlicedSeparableSum{typeof(fsnew),typeof(indarr),length(fsnew)}(fsnew, indarr)
end

# Constructor for the case where the same function is applied to all slices
SlicedSeparableSum(f::F, idxs::T) where {F <: ProximableFunction, T <: Tuple} =
SlicedSeparableSum(([f for k in eachindex(idxs)]...,), idxs)

# Unroll the loop over the different types of functions to evaluate
@generated function (f::SlicedSeparableSum{A, B, N})(x::T) where {A, B, N, T <: AbstractArray}
	ex = :(v = 0.0)
  for i = 1:N # For each function type
    ex = quote $ex;
      for k in eachindex(f.fs[$i]) # For each function of that type
				v += f.fs[$i][k](view(x,f.idxs[$i][k]...))
      end
    end
  end
  ex = :($ex; return v)
end

# Unroll the loop over the different types of functions to prox on
@generated function prox!(y::AbstractArray{T}, f::SlicedSeparableSum{A, B, N}, x::AbstractArray{T}, gamma::Real=1.0) where {T <: RealOrComplex, A, B, N}
  ex = :(v = 0.0)
  for i = 1:N # For each function type
    ex = quote $ex;
      for k in eachindex(f.fs[$i]) # For each function of that type
        g = prox!(view(y, f.idxs[$i][k]...), f.fs[$i][k], view(x,f.idxs[$i][k]...), gamma)
        v += g
      end
    end
  end
  ex = :($ex; return v)
end

is_prox_accurate(f::SlicedSeparableSum) = all([all(is_prox_accurate.(f.fs[k])) for k in eachindex(f.fs)])
is_convex(f::SlicedSeparableSum) = all([all(is_convex.(f.fs[k])) for k in eachindex(f.fs)])
is_set(f::SlicedSeparableSum) = all([all(is_set.(f.fs[k])) for k in eachindex(f.fs)])
is_cone(f::SlicedSeparableSum) = all([all(is_cone.(f.fs[k])) for k in eachindex(f.fs)])

fun_name(f::SlicedSeparableSum) = "sliced separable sum"
fun_dom(f::SlicedSeparableSum) = "n/a" # for now
fun_expr(f::SlicedSeparableSum) = "hard to explain"
fun_params(f::SlicedSeparableSum) = "n/a" # for now

function prox_naive(f::SlicedSeparableSum, x::AbstractArray, gamma)
	fy = 0;
	y = similar(x);
	for t in eachindex(f.fs)
		for k in eachindex(f.fs[t])
			y[f.idxs[t][k]...], fy1 = prox_naive(f.fs[t][k], x[f.idxs[t][k]...], gamma)
			fy += fy1
		end
	end
	return y, fy
end


text:
module Trapdoor

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Trapdoor" => Ahorn.EntityPlacement(
        Maple.Trapdoor
    )
)

# Doesn't seem to have a sprite?
doorColor = (22, 27, 48, 1) ./ (255, 255, 255, 1)

function Ahorn.selection(entity::Maple.Trapdoor)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x, y + 5, 24, 4)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Trapdoor, room::Maple.Room) = Ahorn.drawRectangle(ctx, 0, 6, 24, 2, doorColor)

end


text:
_stand_medium_proc_file() = procdir(RathData, "rath2017___42_MAX_UB_standard_medium.tsv")

_a1at_aa_rel_abundance_proc_file() = procdir(RathData, "a1at_aa_rel_abundance.csv")

function _cont_culture_proc_file(expid)
    expid = string(expid)
    procdir(RathData, "rath2017___cont_exp_$(expid).tsv")
end

_max_invitro_fluxs_proc_file() = procdir(RathData, "rath2017___max_invitro_fluxs.tsv")


text:
using JuMP
using ParameterJuMP
using Polyhedra
mutable struct BalancedRealPolytope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim} <: Polyhedra.VPolytope{T}
    d::D
    points::Vector{VT}
    factory::Union{Nothing, JuMP.OptimizerFactory}
    model::Union{Nothing, JuMP.Model}
    z::Union{Nothing, Vector{ParameterJuMP.ParameterRef}}
    t_0::Union{Nothing, JuMP.VariableRef}
    function BalancedRealPolytope{T, VT, D}(
        d::Polyhedra.FullDim, points::Polyhedra.PointIt,
        factory::Union{Nothing, JuMP.OptimizerFactory}=nothing) where {T, VT, D}
        new{T, VT, D}(Polyhedra.FullDim_convert(D, d), Polyhedra.lazy_collect(points), factory, nothing, nothing, nothing)
    end
end
function BalancedRealPolytope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
    return BalancedRealPolytope{Polyhedra.coefficient_type(eltype(points)), eltype(points),
                                typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::BalancedRealPolytope) = v.d
Polyhedra.vvectortype(::Type{<:BalancedRealPolytope{T, AT}}) where {T, AT} = AT
function Polyhedra.similar_type(PT::Type{<:BalancedRealPolytope}, d::Polyhedra.FullDim, ::Type{T}) where {T}
    BalancedRealPolytope{T, Polyhedra.similar_type(Polyhedra.vvectortype(PT), d, T), typeof(d)}
end

Polyhedra.vreptype(::Type{BalancedRealPolytope{T, AT, D}}) where {T, AT, D} = Polyhedra.Hull{T, AT, D}
Polyhedra.fulltype(::Type{BalancedRealPolytope{T, AT, D}}) where {T, AT, D} = Polyhedra.Hull{T, AT, D}

Base.length(idxs::Polyhedra.PointIndices{T, <:BalancedRealPolytope{T}}) where {T} = 2length(idxs.rep.points)
Base.isempty(idxs::Polyhedra.PointIndices{T, <:BalancedRealPolytope{T}}) where {T} = isempty(idxs.rep.points)
function Polyhedra.startindex(idxs::Polyhedra.PointIndices{T, <:BalancedRealPolytope{T}}) where {T}
    if isempty(idxs.rep.points)
        return nothing
    else
        return eltype(idxs)(1)
    end
end
function Base.get(rep::BalancedRealPolytope{T}, idx::Polyhedra.PointIndex{T}) where {T}
    return sign(idx.value) * rep.points[abs(idx.value)]
end
function Polyhedra.nextindex(rep::BalancedRealPolytope{T}, idx::Polyhedra.PointIndex{T}) where {T}
    if idx.value < 0
        if -idx.value >= length(rep.points)
            return nothing
        else
            return typeof(idx)(idx.value - 1)
        end
    else
        if idx.value >= length(rep.points)
            return typeof(idx)(-1)
        else
            return typeof(idx)(idx.value + 1)
        end
    end
end

function _build_model(brp::BalancedRealPolytope)
    n = length(brp.points)
    brp.model = ParameterJuMP.ModelWithParams(brp.factory)
    t = @variable(brp.model, [1:n])
    q = @variable(brp.model, [1:n])
    # |t| ≤ q
    @constraint(brp.model,  t .≤ q)
    @constraint(brp.model, -t .≤ q)
    # sum |t| ≤ sum q ≤ 1
    @constraint(brp.model, sum(q) ≤ 1)
    return sum(t[i] .* brp.points[i] for i in 1:n)
end

function _parametrized_model(brp::BalancedRealPolytope, v)
    hull = _build_model(brp)
    brp.z = ParameterJuMP.add_parameters(brp.model, collect(v))
    @constraint(brp.model, brp.z .== hull)
end

function _ratio_model(brp::BalancedRealPolytope, v)
    hull = _build_model(brp)
    brp.t_0 = @variable(brp.model)
    @constraint(brp.model, brp.t_0 .* v .== hull)
end


text:
autostyle{T <: Integer}(tp::Type{T}) = generate_formatter("%i")
autostyle{T <: FloatingPoint}(tp::Type{T}) = generate_formatter("%.2f")
autostyle(tp::Type{String}) = generate_formatter("%s")
# autostyle(tp::Any) = generate_formatter("%s")

function autostyle{T}(tp::Type{Nullable{T}})
  g(content) = autostyle(T)(content)
  function f(content)
    if isnull(content)
      return ""
    else
      return g(get(content))
    end
  end
  return f
end

function joinkeyvalpairs(dict)
  buf = IOBuffer()
  for (k,v) in dict
    print(buf, " $k=\"$v\"")
  end
  return takebuf_string(buf)
end

function htmlclasses(s, default)
  if haskey(s.d, :html)
    htmlkey = s.d[:html]
    return joinkeyvalpairs(htmlkey)
  end
  return default
end

function style(::Type{HTMLStyleElement}, data, s::UniformStyleElement)
  classes = htmlclasses(s, "datael")
  f = autostyle(eltype(data))
  HTMLStyleElement(TextProp(f), StaticSpanProp(), classes)
end

function style(::Type{HTMLStyleElement}, data, s::ExpandStyleElement)
  classes = htmlclasses(s, "datael")
  f = generate_formatter("%s")
  HTMLStyleElement(TextProp(f), ExpandSpanProp(), classes)
end

function style(::Type{HTMLStyleElement}, data, s::CartesianStyleElement)
  classes = htmlclasses(s, "datael")
  f = generate_formatter("%s")
  HTMLStyleElement(TextProp(f), CartesianSpanProp(), classes)
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

# To be implemented by user.

"""
    BlackBox <: AbstractSimulation

Parent type for user simulation in which environment is not exposed.
"""
abstract type BlackBox <: AbstractSimulation end

"""
    reset!(sim::BlackBox)

Reset simulation.
"""
reset!(sim::BlackBox) = unimplemented()

"""
    observe(sim::BlackBox)

Return quasi-normalized observation of simulation. Observation should copy, not reference,
simulation state. If simulation is unobservable, leave unimplemented.
"""
observe(sim::BlackBox)::AbstractVector{<:Real} = unimplemented()

"""
    step([rng::AbstractRNG], sim::BlackBox)

Step simulation and return log probability of transition. If simulation uses global RNG,
the one-argument step! should be implemented. If simulation requires non-global RNG, the
two-argument function should be implemented instead.
"""
step!(sim::BlackBox)::Real = unimplemented()
step!(rng::AbstractRNG, sim::BlackBox)::Real = step!(sim)

"""
    isterminal(sim::BlackBox)

Check whether simulation has finished due to time limit or terminal state, independent of
failure status.
"""
isterminal(sim::BlackBox)::Bool = unimplemented()

"""
   isevent(sim::BlackBox)

Check whether simulation is in event (failure) state. For episodic simulation, check whether
event occurred during episode.
"""
isevent(sim::BlackBox)::Bool = unimplemented()

"""
    distance(sim::BlackBox)

Return custom metric of distance to event. For best results, metric should depend only on
current state. For episodic simulation, return minimum distance to event across entire
episode (miss distance).
"""
distance(sim::BlackBox)::Real = unimplemented()

"""
    reward(sim::BlackBox)

Advanced option: additional reward or heuristic. Returns scalar value or partially-applied
function, allowing efficient calculation of R(s,s') = R(s)(s'). Should not alter simulation.
See documentation for more details.
"""
reward(sim::BlackBox) = 0.0


text:
function mean_extrema_pad( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return mean_extrema_pad!( img, zeros(Float64, padsize .+ 1),
                              rad, th_hard,
                              zeros(Bool,size(img)), zeros(Bool,size(img)),
                              fmin=fmin, fmax=fmax, ovp=ovp );
end

function mean_extrema_pad!( img::Array{T,2}, intA::Array{<:AbstractFloat,2},
                            rad::Tuple{Int64,Int64}, th_hard,
                            minima::Array{Bool,2}, maxima::Array{Bool,2};
                            fmin=1, fmax=1, ovp=0 ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    integralArray_pad!( img, intA, 3 .* rad .+ 1 )

    # convenient quantities
    lows  = (  1,  1  ) .+ ( 3 .* rad .+ 1 );
    highs = size( img ) .+ ( 3 .* rad .+ 1 );
    npix  = prod(  (  2 .* rad .+ 1  )  );
    vert  = ( 2*rad[1] + 1 - ovp[1], 0  );
    horz  = ( 0, 2*rad[2] + 1  - ovp[2] );

    # these will accomodate the means and signs of the 9x9 mean neighbourhood
    sums9 = zeros( Float64, 3, 3 );
    signs = zeros(   Int64, 3, 3 );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl = ( y,x ) .- rad .- 1;
        br = ( y,x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );

        ( s0/npix <= th_hard ) && ( continue; ) # s0/npix == mean of the central square

        idx = 1;
        for xoff in -1:1, yoff in -1:1
            off = vert .* yoff .+ horz .* xoff;
            sums9[idx] = ImageAnalysis.integralArea( intA, tl .+ off, br .+ off );
            signs[idx] = sign( s0 - sums9[idx] );
            idx += 1;
        end

        #=
        sums9[1,1] = ImageAnalysis.integralArea( intA, tl .- vert .- horz, br .- vert .- horz );
        sums9[2,1] = ImageAnalysis.integralArea( intA, tl         .- horz, br         .- horz );
        sums9[3,1] = ImageAnalysis.integralArea( intA, tl .+ vert .- horz, br .+ vert .- horz );
        sums9[1,2] = ImageAnalysis.integralArea( intA, tl .- vert        , br .- vert         );
        sums9[2,2] = ImageAnalysis.integralArea( intA, tl                , br                 );
        sums9[3,2] = ImageAnalysis.integralArea( intA, tl .+ vert        , br .+ vert         );
        sums9[1,3] = ImageAnalysis.integralArea( intA, tl .- vert .+ horz, br .- vert .+ horz );
        sums9[2,3] = ImageAnalysis.integralArea( intA, tl         .+ horz, br         .+ horz );
        sums9[3,3] = ImageAnalysis.integralArea( intA, tl .+ vert .+ horz, br .+ vert .+ horz );
        =#

        score_min = 0
        score_max = 0
        for idx in 1:4
            score_min += ( signs[idx] == -1 && signs[end-idx+1] == -1 )
            score_max += ( signs[idx] ==  1 && signs[end-idx+1] ==  1 )
        end

        minima[y-3*rad[1]-1,x-3*rad[2]-1] = score_min > fmin;
        maxima[y-3*rad[1]-1,x-3*rad[2]-1] = score_max > fmax;
    end

    return minima, maxima
end

# sum ( x - mean )^2 = sum( x^2 + mean^2 - 2*x*mean ) = sum x^2 + N*mean^2 - 2*mean*sum x
#                                                     = sum x^2 + mean*sum x - 2*mean*sum x
#                                                     = sum x^2 - mean*sum x

function mean_std_pad( img::Array{<:Real,2}, rad::Tuple{Int64,Int64}; typ=Float64 )

    padsize = size( img ) .+ 2 .* rad;
    return mean_std_pad!( img, zeros( typ, padsize .+ 1), zeros( typ, padsize .+ 1),
                          rad, zeros( typ, size( img ) ) );
end

function mean_std_pad!( img::Array{<:Real,2}, intA::Array{T,2}, intA2::Array{T,2},
                        rad::Tuple{Int64,Int64}, stds::Array{T,2} ) where {T<:AbstractFloat}

    # integral array
    @inbounds @simd for idx in 1:length(intA)
         intA[idx] = 0.0
        intA2[idx] = 0.0
    end
    integralArray_pad!( img,  intA, rad, fun=(x)->(T(x)) );
    integralArray_pad!( img, intA2, rad, fun=(x)->(T(x)^2) );

    # convenience variables
    lows  = (  1,  1  ) .+ rad;
    highs = size( img ) .+ rad;
    npix  = prod( ( 2 .* rad .+ 1 ) );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl  = ( y, x ) .- rad .- 1;
        br  = ( y, x ) .+ rad;
        s0  = ImageAnalysis.integralArea(  intA, tl, br );
        s02 = ImageAnalysis.integralArea( intA2, tl, br );
        num = s02 - s0*s0/npix;

        stds[y-rad[1],x-rad[2]] = sqrt( num / npix )
    end

    return stds
end

function mean_thresh_pad( img::Array{T,2}, rad::Tuple{Int64,Int64}; fun=(x,y)->(x>y) ) where {T<:Real}

    padsize = size( img ) .+ 2 .*rad;
    return mean_thresh_pad!( img, zeros(Float64, padsize .+ 1),
                             rad, zeros(Bool,size(img)), fun=fun );
end

function mean_thresh_pad!(  img::Array{T,2},
                            intA::Array{<:AbstractFloat,2},
                            rad::Tuple{Int64,Int64},
                            filt::Array{Bool,2}; fun=(x,y)->(x>y) ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    integralArray_pad!( img, intA, rad )

    # convenient quantities
    lows  = (  1,  1  ) .+ rad;
    highs = size( img ) .+ rad;
    npix  = prod( ( 2 .* rad .+ 1 ) );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl = ( y,x ) .- rad .- 1;
        br = ( y,x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );
        mean = s0/npix;

        filt[y-rad[1],x-rad[2]] = fun( img[ y-rad[1],x-rad[2] ], mean ) ;
    end

    return filt
end

function mean_extrema_pad2( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}

    return mean_extrema_pad2!( img, zeros(Float64, size(img) .+ 1), rad, th_hard, zeros(Bool,size(img)), zeros(Bool,size(img)),
                              fmin=fmin, fmax=fmax, ovp=ovp );
end

function mean_extrema_pad2!( img::Array{T,2}, intA::Array{<:AbstractFloat,2},
                             rad::Tuple{Int64,Int64}, th_hard,
                             minima::Array{Bool,2}, maxima::Array{Bool,2};
                             fmin=1, fmax=1, ovp=0 ) where {T<:Real}

    # integral array
    @inbounds @simd for idx in 1:length(intA)
    intA[idx] = 0.0
    end
    integralArray!( img, intA )

    if th_hard == nothing
    th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    npix  = prod( ( 2 .* rad .+ 1 ) );
    vert  = ( 2*rad[1] + 1 - ovp[1], 0 );
    horz  = ( 0, 2*rad[2] + 1  - ovp[2]);
    offs  = ( -1 .* (vert .+ horz), -1 .* horz, vert .- horz, -1 .* vert, (0,0), vert, horz .- vert, horz, horz .+ vert )

    # these will accomodate the means and signs of the 9x9 mean neighbourhood
    sums9 = zeros( Float64, 3, 3 );
    signs = zeros(   Int64, 3, 3 );

    @inbounds for x in 1:size(img,2), y in 1:size(img,1)

    tl = max.( 1, min.( size(img), ( y,x ) .- rad .- 1 ) );
    br = max.( 1, min.( size(img), ( y,x ) .+ rad ) );

    s0 = ImageAnalysis.integralArea( intA, tl, br );
    ( s0/npix <= th_hard ) && ( continue; ) # s0/npix == mean of the central square

    cont = 1;
    for off in offs
        tl_ = max.( 1, min.( size(img), tl .+ off ) );
        br_ = max.( 1, min.( size(img), br .+ off ) );
        sums9[cont] = ImageAnalysis.integralArea( intA, tl_, br_ );
        cont += 1;
    end

    @inbounds @simd for idx in 1:9
        signs[idx] = sign( s0 .- sums9[idx] );
    end

    score_min = 0
    score_max = 0
    for idx in 1:4
        score_min += ( signs[idx] == -1 && signs[end-idx+1] == -1 )
        score_max += ( signs[idx] ==  1 && signs[end-idx+1] ==  1 )
    end

    minima[y,x] = score_min > fmin;
    maxima[y,x] = score_max > fmax;
    end

    return minima, maxima
end

function std_extrema_pad_filt( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fmin=1, fmax=1, m=1, fun=(x)->(x) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return std_extrema_pad_filt!( img,
                                  zeros(Float64, padsize .+ 1), zeros(Float64, padsize .+ 1),
                                  zeros(Float64, padsize .+ 1), zeros(Float64, padsize .+ 1), zeros(Float64, padsize .+ 1),
                                  rad, th_hard, zeros(Float32,size(img)),
                                  fmin=fmin, fmax=fmax, m=m, fun=fun );
end

# sum( ( x - mean )^2 - ( xf - mean )^2 ) = sum( x^2 + mean^2 - 2*x*mean ) - sum( xf^2 + mean^2 - 2*xf*mean )

function std_extrema_pad_filt!( img::Array{T,2},
                                intA::Array{<:AbstractFloat,2}, intA2::Array{<:AbstractFloat,2},
                                intF::Array{<:AbstractFloat,2}, intF2::Array{<:AbstractFloat,2}, Ns::Array{<:AbstractFloat,2},
                                rad::Tuple{Int64,Int64}, th_hard, stds::Array{Float32,2};
                                fmin=1, fmax=1, m=1, fun=(x)->(x) ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
         intA[idx] = 0.0
        intA2[idx] = 0.0
    end

    integralArray_pad!( img,  intA, 3 .* rad .+ 1 )
    integralArray_pad!( img, intA2, 3 .* rad .+ 1, fun=(x)->(x*x) )
    integralArray_pad!( img,  intF, 3 .* rad .+ 1, fun=fun )
    integralArray_pad!( img, intF2, 3 .* rad .+ 1, fun=(x)->(fun(x)^2) )
    integralArray_pad!( img,    Ns, 3 .* rad .+ 1, fun=(x)->(fun(x)>0) );

    if th_hard == nothing
        th_hard = ( intA[ end, end ] - intF[ end, end ] )/ length(img); # == mean of the image
    end

    # convenient quantities
    lows  = (  1,  1  ) .+ ( 3 .* rad .+ 1 );
    highs = size( img ) .+ ( 3 .* rad .+ 1 );
    npix  = prod( ( 2 .* rad .+ 1 ) );
    vert  = ( 2*rad[1] + 1, 0 );
    horz  = ( 0, 2*rad[2] + 1 );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl   = ( y,x ) .- rad .- 1;
        br   = ( y,x ) .+ rad;

        s0   = ImageAnalysis.integralArea(  intA, tl, br );
        s02  = ImageAnalysis.integralArea( intA2, tl, br );
        sF   = ImageAnalysis.integralArea(  intF, tl, br );
        sF2  = ImageAnalysis.integralArea( intF2, tl, br );
        nF   = ImageAnalysis.integralArea(    Ns, tl, br );

        mean = ( s0 - sF )/( npix - nF );
        std  = ( s02 + npix*mean*mean - 2*s0*mean ) - ( sF2 + nF*mean*mean - 2*sF*mean );

        stds[y-3*rad[1]-1,x-3*rad[2]-1] = sqrt( abs(std) / npix ) # img[y-3*rad[1]-1,x-3*rad[2]-1] - mean # sqrt( std / npix );
    end

    return stds
end

function mean_thresh_pad_filt( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fun=(x)->(x) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return mean_thresh_pad_filt!( img,
                                  zeros(Float64, padsize .+ 1),
                                  zeros(Float64, padsize .+ 1), zeros(Float64, padsize .+ 1),
                                  rad, th_hard, zeros(Bool,size(img)),
                                  fun=fun );
end

function mean_thresh_pad_filt!(  img::Array{T,2},
                                 intA::Array{<:AbstractFloat,2},
                                 intF::Array{<:AbstractFloat,2}, numF::Array{<:AbstractFloat,2},
                                 rad::Tuple{Int64,Int64}, th_hard,
                                 filt::Array{Bool,2};
                                 fun=(x)->(x) ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    integralArray_pad!( img, intA, 3 .* rad .+ 1 )
    integralArray_pad!( img, intF, 3 .* rad .+ 1, fun=fun )  # fun is a condition, x = cond * x
    integralArray_pad!( img, numF, 3 .* rad .+ 1, fun=(x)->(fun(x)>0) ); # counts the num of pixels that fulfill the condition

    if th_hard == nothing
        th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    lows  = (  1,  1  ) .+ ( 3 .* rad .+ 1 );
    highs = size( img ) .+ ( 3 .* rad .+ 1 );
    npix  = prod( ( 2 .* rad .+ 1 ) );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl = ( y,x ) .- rad .- 1;
        br = ( y,x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );
        sF = ImageAnalysis.integralArea( intF, tl, br );
        nF = ImageAnalysis.integralArea( numF, tl, br );
        mean = ( s0 - sF )/( npix - nF );

        filt[y-3*rad[1]-1,x-3*rad[2]-1] = img[y-3*rad[1]-1,x-3*rad[2]-1] > mean;
    end

    return filt
end

function mean_slope_pad( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fun=(x)->(x>1), ovp=(0,0) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return mean_slope_pad!( img, zeros(Float64, padsize .+ 1), rad, th_hard, zeros(Bool,size(img)),
                            fun=fun, ovp=ovp );
end

function mean_slope_pad!( img::Array{T,2}, intA::Array{<:AbstractFloat,2},
                          rad::Tuple{Int64,Int64}, th_hard,
                          slopes::Array{Bool,2};
                          fun=(x)->(x>1), ovp=0 ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    integralArray_pad!( img, intA, 3 .* rad .+ 1 )

    if th_hard == nothing
        th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    lows  = (  1,  1  ) .+ ( 3 .* rad .+ 1 );
    highs = size( img ) .+ ( 3 .* rad .+ 1 );
    npix  = prod( ( 2 .* rad .+ 1 ) );
    vert  = ( 2*rad[1] + 1 - ovp[1], 0 );
    horz  = ( 0, 2*rad[2] + 1  - ovp[2]);

    # these will accomodate the means and signs of the 9x9 mean neighbourhood
    sums9 = zeros( Float64, 3, 3 );
    signs = zeros(   Int64, 3, 3 );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl = ( y,x ) .- rad .- 1;
        br = ( y,x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );
        ( s0/npix <= th_hard ) && ( continue; ) # s0/npix == mean of the central square

        sums9[1,1] = ImageAnalysis.integralArea( intA, tl .- vert .- horz, br .- vert .- horz );
        sums9[2,1] = ImageAnalysis.integralArea( intA, tl         .- horz, br         .- horz );
        sums9[3,1] = ImageAnalysis.integralArea( intA, tl .+ vert .- horz, br .+ vert .- horz );
        sums9[1,2] = ImageAnalysis.integralArea( intA, tl .- vert        , br .- vert         );
        sums9[2,2] = ImageAnalysis.integralArea( intA, tl                , br                 );
        sums9[3,2] = ImageAnalysis.integralArea( intA, tl .+ vert        , br .+ vert         );
        sums9[1,3] = ImageAnalysis.integralArea( intA, tl .- vert .+ horz, br .- vert .+ horz );
        sums9[2,3] = ImageAnalysis.integralArea( intA, tl         .+ horz, br         .+ horz );
        sums9[3,3] = ImageAnalysis.integralArea( intA, tl .+ vert .+ horz, br .+ vert .+ horz );

        @inbounds @simd for idx in 1:9
            signs[idx] = sign( s0 .- sums9[idx] );
        end

        score_slope = 0
        for idx in 1:4
            score_slope += ( signs[idx]*signs[end-idx+1] == -1 )
        end

        slopes[y-3*rad[1]-1,x-3*rad[2]-1] = fun( score_slope );
    end

    return slopes
end

function mean_thresh_count( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fun=(x,y)->(x<y), f=1.0 ) where {T<:Real}

    padsize = size( img ) .+ 2 .* rad; #( 3 .* rad .+ 1 );
    return mean_thresh_count!( img,
                               zeros(Float64, padsize .+ 1),
                               rad, th_hard, zeros(Int64,size(img)), fun=fun, f=f );
end

function mean_thresh_count!(  img::Array{T,2},
                              intA::Array{<:AbstractFloat,2},
                              rad::Tuple{Int64,Int64}, th_hard,
                              counts::Array{Int64,2}; fun=(x,y)->(x<y), f=1.0 ) where {T<:Real}
    # Building integral array
    # 1-. Setting intA to zero (redudant when calling out-of-place).
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    # 2-. Populating integral array
    integralArray_pad!( img, intA, rad ); #3 .* rad .+ 1 )

    # Checking the threshold value
    if th_hard == nothing
        th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # Convenient quantities
    npix   = prod( ( 2 .* rad .+ 1 ) );
    y0, x0 = (  1,  1  ) .+ rad; #( 3 .* rad .+ 1 );
    y1, x1 = size( img ) .+ rad; #( 3 .* rad .+ 1 );
    ry, rx = rad;

    for x in x0:x1, y in y0:y1

        # top left (tl) and bottom right (br) of the central patch
        tl = ( y, x ) .- rad .- 1;
        br = ( y, x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );
        local_mean = s0*f/npix;

        if local_mean < th_hard
            counts[ y-y0+1, x-x0+1 ] = 1000000;
        end

        # For each pixel in the central patch that falls inside the image (we added padding to the integral array)
        # add +1 if the pixel is above the local mean.

        for col in max(tl[2]+1,x0):min(br[2],x1)
            for row in max(tl[1]+1,y0):min(br[1],y1)

                counts[ row-y0+1,col-x0+1 ] += fun( img[ row-y0+1,col-x0+1 ], local_mean )
        end end

        #filt[y-3*rad[1]-1,x-3*rad[2]-1] = img[ y-3*rad[1]-1,x-3*rad[2]-1 ] > local_mean*0.95 ;
    end

    for idx in 1:length(counts)
        if counts[idx] > 1000000
            counts[idx] = -1
        end
    end

    return counts
end

function mean_laplacian_pad( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return mean_laplacian_pad!( img, zeros(Float64, padsize .+ 1), rad, th_hard, zeros(Float32,size(img)),
                              fmin=fmin, fmax=fmax, ovp=ovp );
end

function mean_laplacian_pad!( img::Array{T,2}, intA::Array{<:AbstractFloat,2},
                              rad::Tuple{Int64,Int64}, th_hard,
                              laplacian::Array{Float32,2};
                              fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    integralArray_pad!( img, intA, 3 .* rad .+ 1 )

    if th_hard == nothing
        th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    lows  = (  1,  1  ) .+ ( 3 .* rad .+ 1 );
    highs = size( img ) .+ ( 3 .* rad .+ 1 );
    npix  = prod( ( 2 .* rad .+ 1 ) );
    vert  = ( 2*rad[1] + 1 - ovp[1], 0 );
    horz  = ( 0, 2*rad[2] + 1  - ovp[2]);

    sums9 = zeros( Float32, 3, 3 );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl = ( y,x ) .- rad .- 1;
        br = ( y,x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );
        ( s0/npix <= th_hard ) && ( continue; ) # s0/npix == mean of the central square

        sums9[1,1] = ImageAnalysis.integralArea( intA, tl .- vert .- horz, br .- vert .- horz );
        sums9[2,1] = ImageAnalysis.integralArea( intA, tl         .- horz, br         .- horz );
        sums9[3,1] = ImageAnalysis.integralArea( intA, tl .+ vert .- horz, br .+ vert .- horz );
        sums9[1,2] = ImageAnalysis.integralArea( intA, tl .- vert        , br .- vert         );
        sums9[2,2] = ImageAnalysis.integralArea( intA, tl                , br                 );
        sums9[3,2] = ImageAnalysis.integralArea( intA, tl .+ vert        , br .+ vert         );
        sums9[1,3] = ImageAnalysis.integralArea( intA, tl .- vert .+ horz, br .- vert .+ horz );
        sums9[2,3] = ImageAnalysis.integralArea( intA, tl         .+ horz, br         .+ horz );
        sums9[3,3] = ImageAnalysis.integralArea( intA, tl .+ vert .+ horz, br .+ vert .+ horz );

        laplacian[y-3*rad[1]-1,x-3*rad[2]-1] = abs( 4*s0 - sums9[1,2] - sums9[2,1] - sums9[2,3] - sums9[3,2] );
    end

    return laplacian
end

function mean_laplacian_pad_filt( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fun=(x)->(x), ovp=(0,0) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return mean_laplacian_pad_filt!( img,
                                  zeros(Float64, padsize .+ 1),
                                  zeros(Float64, padsize .+ 1), zeros(Float64, padsize .+ 1),
                                  rad, th_hard, zeros(Float32,size(img)),
                                  fun=fun, ovp=ovp );
end

function mean_laplacian_pad_filt!(  img::Array{T,2},
                                 intA::Array{<:AbstractFloat,2},
                                 intF::Array{<:AbstractFloat,2}, numF::Array{<:AbstractFloat,2},
                                 rad::Tuple{Int64,Int64}, th_hard,
                                 laplacian::Array{Float32,2};
                                 fun=(x)->(x), ovp=(0,0) ) where {T<:Real}
    # integral array
    @inbounds @simd for idx in 1:length(intA)
        intA[idx] = 0.0
    end
    integralArray_pad!( img, intA, 3 .* rad .+ 1 )
    integralArray_pad!( img, intF, 3 .* rad .+ 1, fun=fun )
    integralArray_pad!( img, numF, 3 .* rad .+ 1, fun=(x)->(fun(x)>0) );

    if th_hard == nothing
        th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    lows  = (  1,  1  ) .+ ( 3 .* rad .+ 1 );
    highs = size( img ) .+ ( 3 .* rad .+ 1 );
    npix  = prod( ( 2 .* rad .+ 1 ) );
    vert  = ( 2*rad[1] + 1 - ovp[1], 0 );
    horz  = ( 0, 2*rad[2] + 1  - ovp[2]);

    sums9 = zeros( Float32, 3, 3 );

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl = ( y,x ) .- rad .- 1;
        br = ( y,x ) .+ rad;
        s0 = ImageAnalysis.integralArea( intA, tl, br );
        sF = ImageAnalysis.integralArea( intF, tl, br );
        nF = ImageAnalysis.integralArea( numF, tl, br );
        mean = ( s0 - sF )/( 1 + npix - nF );

        sums9[1] = ( ImageAnalysis.integralArea( intA, tl .- vert .- horz, br .- vert .- horz ) - ImageAnalysis.integralArea( intF, tl .- vert .- horz, br .- vert .- horz ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl .- vert .- horz, br .- vert .- horz ) );
        sums9[2] = ( ImageAnalysis.integralArea( intA, tl         .- horz, br         .- horz ) - ImageAnalysis.integralArea( intF, tl         .- horz, br         .- horz ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl         .- horz, br         .- horz ) );
        sums9[3] = ( ImageAnalysis.integralArea( intA, tl .+ vert .- horz, br .+ vert .- horz ) - ImageAnalysis.integralArea( intF, tl .+ vert .- horz, br .+ vert .- horz ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl .+ vert .- horz, br .+ vert .- horz ) );
        sums9[4] = ( ImageAnalysis.integralArea( intA, tl .+ vert        , br .+ vert         ) - ImageAnalysis.integralArea( intF, tl .+ vert        , br .+ vert         ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl .+ vert        , br .+ vert         ) );
        sums9[5] = ( ImageAnalysis.integralArea( intA, tl .+ vert .+ horz, br .+ vert .+ horz ) - ImageAnalysis.integralArea( intF, tl .+ vert .+ horz, br .+ vert .+ horz ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl .+ vert .+ horz, br .+ vert .+ horz ) );
        sums9[6] = ( ImageAnalysis.integralArea( intA, tl         .+ horz, br         .+ horz ) - ImageAnalysis.integralArea( intF, tl         .+ horz, br         .+ horz ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl         .+ horz, br         .+ horz ) );
        sums9[7] = ( ImageAnalysis.integralArea( intA, tl .- vert .+ horz, br .- vert .+ horz ) - ImageAnalysis.integralArea( intF, tl .- vert .+ horz, br .- vert .+ horz ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl .- vert .+ horz, br .- vert .+ horz ) );
        sums9[8] = ( ImageAnalysis.integralArea( intA, tl .- vert        , br .- vert         ) - ImageAnalysis.integralArea( intF, tl .- vert        , br .- vert         ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl .- vert        , br .- vert         ) );
        sums9[9] = sums9[1]; #( ImageAnalysis.integralArea( intA, tl                , br                 ) - ImageAnalysis.integralArea( intF, tl                , br                 ) )/( 1 + npix - ImageAnalysis.integralArea( numF, tl                , br                 ) );

        difs = 0.0
        for idx in 1:8
            difs = max( abs( sums9[idx] - sums9[idx+1] ), difs );
        end
        laplacian[y-3*rad[1]-1,x-3*rad[2]-1] = difs;
    end

    return laplacian
end

function mean_corr_pad( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}

    padsize = size( img ) .+ 2 .* ( 3 .* rad .+ 1 );
    return mean_corr_pad!( img, zeros(Float64, padsize .+ 1), zeros(Float64, padsize .+ 1),
                           rad, th_hard, zeros(Float32,size(img)), fmin=fmin, fmax=fmax, ovp=ovp );
end

function mean_corr_pad!( img::Array{T,2},
                         intA::Array{<:AbstractFloat,2}, intA2::Array{<:AbstractFloat,2},
                         rad::Tuple{Int64,Int64}, th_hard,
                         corrs::Array{Float32,2};
                         fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}

    # integral array
    @inbounds @simd for idx in 1:length(intA)
         intA[idx] = 0.0;
        intA2[idx] = 0.0;
    end
    integralArray!( img, intA );
    integralArray!( img, intA2, fun=(x)->(x^2) );

    if th_hard == nothing
        th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    lows  = (  1,  1  ) .+ rad;
    highs = size( img ) .- rad;
    npix  = prod( ( 2 .* rad .+ 1 ) );
    dists = zeros( Float32, 2 .* rad .+ 1 );

    # x = image intensity
    # y = distance from minimum value
    # corr( x, y ) = sum( ( x - meanX )( y - meanY ) )/ ( stdx * stdy );
    #              = sum( x*y - x*meanY - y*meanX + meanX*meanY )
    #              = sum( x*y ) - meanY*sum( x ) - meanX*sum( y ) + N*meanX*meanY
    #              =   loop       loop * IA          IA * loop         IA * loop

    # stdx = sum( ( x - meanX )^2 ) = sum( x^2 ) + N*meanX^2 - meanX*sum( x )  -> IA
    # stdy = sum( ( y - meanY )^2 ) = sum( y^2 ) + N*meanY^2 - meanY*sum( y )  -> loop

    @inbounds for x in lows[2]:highs[2], y in lows[1]:highs[1]

        tl    = ( y,x ) .- rad .- 1;
        br    = ( y,x ) .+ rad;
        sumI  = ImageAnalysis.integralArea( intA , tl, br );
        sumI2 = ImageAnalysis.integralArea( intA2, tl, br );
        meanI = sumI/npix;
        stdI  = sqrt( ( sumI2 + npix*meanI^2 - meanI*sumI )/npix );

        ( meanI <= th_hard ) && ( continue; )

        #println( tlimg, "  ", brimg )

        sumD, sumD2, meanD, sumID = maxdist!( img, dists, tl[1] + 1, br[1], tl[2] + 1, br[2] );
        stdD = sqrt( ( sumD2 + npix*meanD^2 - meanD*sumD )/npix );

        corr = ( sumID - meanD*sumI - meanI*sumD + npix*meanI*meanD )/( stdI * stdD )

        corrs[y,x] = corr;
    end

    return corrs
end

function mindist!( a, dists, y0, y1, x0, x1 )
	minpos = ( y0, x0 );
	minval = a[minpos...];
	N = 1;
	for col in x0:x1, row in y0:y1

        if     ( a[row,col] <  minval )
            minval = a[row,col];
            minpos = (row,col);
            N = 1;
        elseif ( a[row,col] == minval )
            minpos = minpos .+ (row,col);
            N += 1;
        end
	end

    minpos = minpos ./ N

    sumdist2 = 0.0
    sumdist  = 0.0
    sumXY    = 0.0
	@inbounds for col in 1:size( dists, 2 )
		d2 = ( x0 + col - 1 - minpos[2] )^2;
		@simd for row in 1:size( dists, 1 )
			d1 = ( y0 + row - 1 - minpos[1] )^2;
			dists[row,col] = sqrt( d1 + d2 )
            sumdist  += dists[row,col]
            sumdist2 += d1 + d2
            sumXY    += a[ y0+row-1, x0+col-1 ] * dists[row,col];
		end
	end

	return sumdist, sumdist2, sumdist/length(dists ), sumXY # mean
end

function maxdist!( a, dists, y0, y1, x0, x1 )
	maxpos = ( y0, x0 );
	maxval = a[maxpos...];
	N = 1;
	for col in x0:x1, row in y0:y1

        if     ( a[row,col] > maxval )
            maxval = a[row,col];
            maxpos = (row,col);
            N = 1;
        elseif ( a[row,col] == maxval )
            maxpos = maxpos .+ (row,col);
            N += 1;
        end
	end

    maxpos = maxpos ./ N;

    sumdist2 = 0.0
    sumdist  = 0.0
    sumXY    = 0.0
	@inbounds for col in 1:size( dists, 2 )
		d2 = ( x0 + col - 1 - maxpos[2] )^2;
		@simd for row in 1:size( dists, 1 )
			d1 = ( y0 + row - 1 - maxpos[1] )^2;
			dists[row,col] = sqrt( d1 + d2 )
            sumdist  += dists[row,col]
            sumdist2 += d1 + d2
            sumXY    += a[ y0+row-1, x0+col-1 ] * dists[row,col];
		end
	end

	return sumdist, sumdist2, sumdist/length(dists ), sumXY # mean
end

"""  3D """

function mean_extrema_pad3( img::Array{T,2}, rad::Tuple{Int64,Int64}, th_hard; fmin=1, fmax=1, ovp=(0,0) ) where {T<:Real}

    return mean_extrema_pad3!( img, zeros(Float64, size(img) .+ 1), rad, th_hard, zeros(Bool,size(img)), zeros(Bool,size(img)),
                              fmin=fmin, fmax=fmax, ovp=ovp );
end

function mean_extrema_pad3!( img::Array{T,2}, intA::Array{<:AbstractFloat,2},
                             rad::Tuple{Int64,Int64}, th_hard,
                             minima::Array{Bool,2}, maxima::Array{Bool,2};
                             fmin=1, fmax=1, ovp=0 ) where {T<:Real}

    # integral array
    @inbounds @simd for idx in 1:length(intA)
    intA[idx] = 0.0
    end
    integralArray!( img, intA )

    if th_hard == nothing
    th_hard = intA[ end, end ] / length(img); # == mean of the image
    end

    # convenient quantities
    npix  = prod( ( 2 .* rad .+ 1 ) );
    vert  = ( 2*rad[1] + 1 - ovp[1], 0 );
    horz  = ( 0, 2*rad[2] + 1  - ovp[2]);

    # these will accomodate the means and signs of the 9x9 mean neighbourhood
    sums9 = zeros( Float64, 3, 3 );
    signs = zeros(   Int64, 3, 3 );

    @inbounds for x in 1:size(img,2), y in 1:size(img,1)

    tl = max.( 1, min.( size(img), ( y,x ) .- rad .- 1 ) );
    br = max.( 1, min.( size(img), ( y,x ) .+ rad ) );

    s0 = ImageAnalysis.integralArea( intA, tl, br );
    ( s0/npix <= th_hard ) && ( continue; ) # s0/npix == mean of the central square

    cont = 1;
    for off in ( -1 .* (vert .+ horz), -1 .* horz, vert .- horz, -1 .* vert, (0,0), vert, horz .- vert, horz, horz .+ vert )
        tl_ = max.( 1, min.( size(img), tl .+ off ) );
        br_ = max.( 1, min.( size(img), br .+ off ) );
        sums9[cont] = ImageAnalysis.integralArea( intA, tl_, br_ );
        cont += 1;
    end

    @inbounds @simd for idx in 1:9
        signs[idx] = sign( s0 .- sums9[idx] );
    end

    score_min = 0
    score_max = 0
    for idx in 1:4
        score_min += ( signs[idx] == -1 && signs[end-idx+1] == -1 )
        score_max += ( signs[idx] ==  1 && signs[end-idx+1] ==  1 )
    end

    minima[y,x] = score_min > fmin;
    maxima[y,x] = score_max > fmax;
    end

    return minima, maxima
end

function int_maxima_pad!( vol::Array{T,3},
                          intA::Array{Float32,3},
                          maxima::Array{Bool,3}, minima::Array{Bool,3},
                          rad::Tuple{Int64,Int64,Int64}, th_hard; f=12, m=1 ) where {T<:Real}

    h, w, d = size( vol );
    npix = prod( ( 2 .* rad .+ 1 ) );

    # integral array of padvolume

    intA = ImageAnalysis.integralArray( vol, padvol, type=Float32 );

    # Convenience vectors to move 1 side in each direction

    vert = (2*rad[1] + 1,0,0);
    horz = (0,2*rad[2] + 1,0);
    depz = (0,0,2*rad[3] + 1);
    offs = ( -1 .* (vert .+ horz .+ depz ), -1 .* ( horz .+ depz ), vert .- horz .- depz, -1 .* ( vert .+ depz ), -1 .* depz, vert .- depz, horz .- vert .- depz, horz .- depz, vert .+ horz .- depz,
             -1 .* (vert .+ horz ), -1 .* ( horz ), vert .- horz, -1 .* ( vert ), (0,0,0), vert, horz .- vert, horz, vert .+ horz,
             -1 .* (vert .+ horz .- depz ), -1 .* ( horz .- depz ), vert .- horz .+ depz, -1 .* ( vert .- depz ), depz, vert .+ depz, horz .- vert .+ depz, horz .+ depz, vert .+ horz .+ depz,
           )

    # array to store the 27-neighbour sums

    sums27 = zeros( Float64, 3, 3, 3 );
    signs  = zeros(   Int64, 3, 3, 3 );

    for z in 1:d, x in 1:w, y in 1:h

        # top-left and bottom-right corners
        tl = min.( 1, max.( size(vol), ( y,x,z ) .- rad .- 1 ) );
        br = min.( 1, max.( size(vol), ( y,x,z ) .+ rad ) );

        s0   = ImageAnalysis.integralArea( intA, tl, br );
        mean = s0/npix;

        if mean <= th_hard
            continue;
        end

        cont = 1;
        for off in offs
            tl_ = min.( 1, max.( size(vol), tl .+ off ) );
            br_ = min.( 1, max.( size(vol), br .+ off ) );

            sums27[cont] = ImageAnalysis.integralArea( intA, tl_, br_ );
            cont += 1;
        end

        signs = sign.( s0*m .- sums27 );
        score = 0
        for idx in 1:13
            score += signs[idx] == 1 && signs[end-idx+1] == 1
        end

        maxima[y-3*rad[1]-1,x-3*rad[2]-1,z-3*rad[3]-1] = score > f;
    end

    return maxima;
end

#=
    sums27[1,1,1] = ImageAnalysis.integralArea( intA, tl .- vert .- horz .- depz, br .- vert .- horz .- depz );
    sums27[2,1,1] = ImageAnalysis.integralArea( intA, tl         .- horz .- depz, br         .- horz .- depz );
    sums27[3,1,1] = ImageAnalysis.integralArea( intA, tl .+ vert .- horz .- depz, br .+ vert .- horz .- depz );
    sums27[1,2,1] = ImageAnalysis.integralArea( intA, tl .- vert         .- depz, br .- vert         .- depz );
    sums27[2,2,1] = ImageAnalysis.integralArea( intA, tl                 .- depz, br                 .- depz );
    sums27[3,2,1] = ImageAnalysis.integralArea( intA, tl .+ vert         .- depz, br .+ vert         .- depz );
    sums27[1,3,1] = ImageAnalysis.integralArea( intA, tl .- vert .+ horz .- depz, br .- vert .+ horz .- depz );
    sums27[2,3,1] = ImageAnalysis.integralArea( intA, tl         .+ horz .- depz, br         .+ horz .- depz );
    sums27[3,3,1] = ImageAnalysis.integralArea( intA, tl .+ vert .+ horz .- depz, br .+ vert .+ horz .- depz );

    sums27[1,1,2] = ImageAnalysis.integralArea( intA, tl .- vert .- horz        , br .- vert .- horz         );
    sums27[2,1,2] = ImageAnalysis.integralArea( intA, tl         .- horz        , br         .- horz         );
    sums27[3,1,2] = ImageAnalysis.integralArea( intA, tl .+ vert .- horz        , br .+ vert .- horz         );
    sums27[1,2,2] = ImageAnalysis.integralArea( intA, tl .- vert                , br .- vert                 );
    sums27[2,2,2] = s0 #
    sums27[3,2,2] = ImageAnalysis.integralArea( intA, tl .+ vert                , br .+ vert                 );
    sums27[1,3,2] = ImageAnalysis.integralArea( intA, tl .- vert .+ horz        , br .- vert .+ horz         );
    sums27[2,3,2] = ImageAnalysis.integralArea( intA, tl         .+ horz        , br         .+ horz         );
    sums27[3,3,2] = ImageAnalysis.integralArea( intA, tl .+ vert .+ horz        , br .+ vert .+ horz         );

    sums27[1,1,3] = ImageAnalysis.integralArea( intA, tl .- vert .- horz .+ depz, br .- vert .- horz .+ depz );
    sums27[2,1,3] = ImageAnalysis.integralArea( intA, tl         .- horz .+ depz, br         .- horz .+ depz );
    sums27[3,1,3] = ImageAnalysis.integralArea( intA, tl .+ vert .- horz .+ depz, br .+ vert .- horz .+ depz );
    sums27[1,2,3] = ImageAnalysis.integralArea( intA, tl .- vert         .+ depz, br .- vert         .+ depz );
    sums27[2,2,3] = ImageAnalysis.integralArea( intA, tl                 .+ depz, br                 .+ depz );
    sums27[3,2,3] = ImageAnalysis.integralArea( intA, tl .+ vert         .+ depz, br .+ vert         .+ depz );
    sums27[1,3,3] = ImageAnalysis.integralArea( intA, tl .- vert .+ horz .+ depz, br .- vert .+ horz .+ depz );
    sums27[2,3,3] = ImageAnalysis.integralArea( intA, tl         .+ horz .+ depz, br         .+ horz .+ depz );
    sums27[3,3,3] = ImageAnalysis.integralArea( intA, tl .+ vert .+ horz .+ depz, br .+ vert .+ horz .+ depz );
=#

# gr = grid radius; sr = square radius; big_radius = (2 .* sr .+  1) .* gr .+ sr
# padvol = sizevol + 2.* big_radius
# you can avoid padding by testing wether coordinates are within range, but that is somewhat slower

function local_grid_pad( vol, gr, sr )
    big_radius = ( 2 .* sr .+  1 ) .* gr .+ sr;
    padsize    = size( vol ) .+ 2 .* big_radius;
    local_grid_pad!( vol, zeros( padsize .+ 1 ), zeros( size(vol) ), gr, sr )
end

function local_grid_pad!( vol::Array{T,3},
                          intA::Array{Float32,3},
                          out::Array{Bool,3},
                          gr::Tuple{Int64,Int64,Int64},
                          sr::Tuple{Int64,Int64,Int64} ) where {T<:Real}

    # integral array of padvolume
    big_rad = ( 2 .* sr .+  1 ) .* gr .+ sr;
    intA = ImageAnalysis.integralArray_pad!( vol, intA, ones(3).*big_rad, type=Float32 );

    h, w, d = size( vol );
    n = prod( ( 2 .* square_rad .+ 1 ) ); # number of pi/voxels in each grid element

    # Convenience vectors to move 1 side in each direction

    vert = (2*rad[1] + 1,0,0);
    horz = (0,2*rad[2] + 1,0);
    depz = (0,0,2*rad[3] + 1);

    offs = ( -1 .* (vert .+ horz .+ depz ), -1 .* ( horz .+ depz ), vert .- horz .- depz, -1 .* ( vert .+ depz ), -1 .* depz, vert .- depz, horz .- vert .- depz, horz .- depz, vert .+ horz .- depz,
    -1 .* (vert .+ horz ), -1 .* ( horz ), vert .- horz, -1 .* ( vert ), (0,0,0), vert, horz .- vert, horz, vert .+ horz,
    -1 .* (vert .+ horz .- depz ), -1 .* ( horz .- depz ), vert .- horz .+ depz, -1 .* ( vert .- depz ), depz, vert .+ depz, horz .- vert .+ depz, horz .+ depz, vert .+ horz .+ depz,
    )

    # array to store the 27-neighbour sums

    sums27 = zeros( Float64, 3, 3, 3 );
    signs  = zeros(   Int64, 3, 3, 3 );

    for z in 1:d, x in 1:w, y in 1:h

        # top-left and bottom-right corners
        tl = min.( 1, max.( size(vol), ( y,x,z ) .- rad .- 1 ) );
        br = min.( 1, max.( size(vol), ( y,x,z ) .+ rad ) );

        s0   = ImageAnalysis.integralArea( intA, tl, br );
        mean = s0/npix;

        if mean <= th_hard
        continue;
        end

        cont = 1;
        for off in offs
            tl_ = min.( 1, max.( size(vol), tl .+ off ) );
            br_ = min.( 1, max.( size(vol), br .+ off ) );

            sums27[cont] = ImageAnalysis.integralArea( intA, tl_, br_ );
            cont += 1;
        end

        signs = sign.( s0*m .- sums27 );
        score = 0
        for idx in 1:13
            score += signs[idx] == 1 && signs[end-idx+1] == 1
        end

        maxima[y-3*rad[1]-1,x-3*rad[2]-1,z-3*rad[3]-1] = score > f;
    end

    return maxima;
end

max_stars_repo_name: BrucePBon/ImageAnalysis

max_issues_repo_name: BrucePBon/ImageAnalysis

max_forks_repo_name: BrucePBon/ImageAnalysis


text:
module Examples

  using ..AlphaZero

  include("../games/tictactoe/main.jl")
  export Tictactoe

  include("../games/connect-four/main.jl")
  export ConnectFour

  include("../games/grid-world/main.jl")
  export GridWorld

  include("../games/mancala/main.jl")
  export Mancala

  const games = Dict(
    "grid-world" => GridWorld.GameSpec(),
    "tictactoe" => Tictactoe.GameSpec(),
    "connect-four" => ConnectFour.GameSpec(),
    "mancala" => Mancala.GameSpec())
    # "ospiel_ttt" => OSpielTictactoe.GameSpec()
  # ospiel_ttt is added from openspiel_example.jl when OpenSpiel.jl is imported

  const experiments = Dict(
    "grid-world" => GridWorld.Training.experiment,
    "tictactoe" => Tictactoe.Training.experiment,
    "connect-four" => ConnectFour.Training.experiment,
    "mancala" => Mancala.Training.experiment)
    # "ospiel_ttt" => OSpielTictactoe.Training.experiment

end


text:
using Test
push!(LOAD_PATH, joinpath(@__DIR__, "../src"))
using Gravityfield_unitless

let
    @test V(1.) == 4π / 3
    @test ρ(m = 4π / 3, r = 1) == 1.
    @test ρ(r = 1, m = 4π / 3) == 1.
    @test ρ_earth |> typeof == Float64
end

let
    "Example position"
    p1 = [0.,0.,0.]
    "Example position"
    p2 = [0.,1.,0.]
    "Example mass"
    m1 = 316227. #kg
    "Example mass"
    m2 = 129099. #kg
end
@testset "Gravity acceleration" begin
    @test g(1, 0) == 0.
    @test g(1., 1.) == γ
    @test g(m_earth, r_earth) ≈ 10.0391664
    ω_earth_sun = 2π / yr
    a_earth_sun = ω_earth_sun^2 * d_sun_earth
    @test round(g(m_sun, d_sun_earth), digits = 4) == round(a_earth_sun, digits = 4)
end

@testset "Force vector on 1 due to 2" begin
    p1 = [0.,0.,0.] #m
    p2 = [0.,1.,0.] #m
    m1 = 316227. #kg
    m2 = 129099. #kg
    f1_2 = m1* uvec(p1, p2) * g(m2, len(p2-p1))
    @test f1_2 |> length == 3
    @test f1_2 |> typeof == Vector{Float64}
    @test f1_2 ≈ [0.0, 2.7247474, 0.]
end
@testset "Euclidean vector length" begin
    @test len(1) == 1.
    @test len([1]) == 1.
    @test len([2,0]) == 2.
    @test len([1,1]) == sqrt(2)
end
@testset "Vector between bodies or points" begin
    p1 = [0.,0.,0.] #m
    p2 = [0.,1.,0.] #m
    s = Sphere(p1, m_earth, r_earth)
    @test rvec(s,s) == [0., 0., 0.]
    @test len(s,s) == 0.
    @test uvec(s,s) == [0., 0., 0.]
    @test rvec(p2,s) == [0., -1.0, 0.]
    @test len(p2,s) == 1.
    @test uvec(p2,s) == [0., -1.0, 0.]
    s2 = Sphere(p2, m_earth, r_earth)
    @test rvec(s,s2) == [0., 1.0, 0.]
    @test len(s,s2) == 1
    @test uvec(s,s2) == [0., 1.0, 0.]
    s3 = Sphere([-d_sun_earth], 1, 1)
    s4 = Sphere([0.], 1, 1)
    @test rvec(s3,s4) == [d_sun_earth]
    @test len(s4,s3) == d_sun_earth
    @test uvec(s3,s4) == [1]
    @test uvec(s4,s3) == [-1]
end

@testset "Gravity acceleration vector btw. bodies or point and body" begin
    p1 = [0.,0.,0.] #m
    s = Sphere(p1, 1., 1.)
    @test gvec(s,s) == [0., 0., 0.]
    @test gvec([-1, 0., 0.], s) == [γ, 0., 0.]
    s2 = Sphere([0.0], m_earth, r_earth)
    @test gvec([r_earth], s2) == [-10.039166439909296]
    @test gvec([r_earth/2], s2) == [-10.039166439909296]/2
end

@testset "Volume" begin
    @test V(1) == 4π/3
    @test round(V(r_sun) * ρ_sun /m_sun, digits=3) == 0.999
    @test typeof(V(r_sun)) == Float64
end

@testset "Sphere constructors" begin
    s1 = Sphere(p = [0], r = r_earth)
    s2 = Sphere([0], m_earth, r_earth)
    @test s1.m == s2.m
    @test s1.r == s2.r
    s3 = Sphere()
    @test s3.m ≈ s2.m
    @test s3.r ≈ s2.r
    s4 = earthlikeSphere_r_known(r = 10.E6)
    s5 = earthlikeSphere_m_known(m = ρ_earth * V(10.E6))
    @test s4.m ≈ s5.m
    @test s4.r ≈ s5.r
end

@testset "Position vectors" begin
    @test [1, 2, 3] isa Pos
    @test !isa([1 2; 3 4], Pos)
    @test !isa([1, "2", 3], Pos)
    @test !isa([1, sin, 3], Pos)
    @test !isa([[1, 2, 3, 4,5,6]], Pos)
    @test Pos([1, 2, 3, 4,5,6]) isa Pos
    @test !isa([[1, 2, 3, 4,5,6]], Pos)
end

@testset "Gravity acceleration vector in system of bodies" begin
    earth = Sphere()
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0., 0.], m = m_moon)
    sun = Sphere([-d_sun_earth, 0., 0.], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    bodies = [sun, earth, moon]
    @test gvec(earth.p, bodies)[1] < 0 # Sun stronger than moon in centre of earth
    @test gvec([r_earth/2, 0., 0.], bodies)[1] < -5 # Halfway to centre of earth, earth dominates
    xgrid, ygrid = xygrids( xmin = -d_moon_earth/2.,
    xmax = d_moon_earth * 3/2,
    ymin =-d_moon_earth/2,
    ymax = d_moon_earth * 3/2, n = 100)
    gmagn = gfieldvals(nest(xgrid, ygrid, moon.r),  bodies)
    @test foldr(+, gmagn)  ≈  222.5649149449058
    @test gvec(earth, bodies) == gvec(earth.p, bodies)
end

@testset "Nested position vectors" begin
    earth = Sphere()
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0., 0.], m = m_moon)
    sun = Sphere([-d_sun_earth, 0., 0.], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    bodies = [sun, earth, moon]
    p1 = Pos([1., 2.])
    p2 = Pos([-1., -2.])
    @test Npos([p1, p2]) == [[1., 2.], [-1., -2.]]
    @test nest(1:3, 2, 3) == [[1, 2, 3], [2, 2, 3], [3, 2, 3]]
    @test nest(1,2:3, 4) == [[1, 2, 4], [1, 3, 4]]
    @test nest(1,2, 3:4) == [[1, 2, 3], [1, 2, 4]]
    @test nest(1:2, 1:2) == [[1,1], [2,2]]
    @test nest(1,2:3) == [[1,2], [1,3]]
    @test nest(1:2,3) == [[1,3], [2,3]]
    @test nest(range(1.5, 3.5, length = 2)) == [[1.5], [3.5]]
    @test gfieldvals(nest(1,2:3,4), bodies)[1] ≈ 0.0059000661748
    @test nest([1, 2, 3], [4, 5, 6], [7, 8, 9]) ==  [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
    np = nest([1, 10], [2, 20], [3, 30], [4, 40], [5, 50], [6, 60])
    @test np ==  [[1, 2, 3, 4, 5, 6], [10, 20, 30, 40, 50, 60]]
    @test slice_1(np) == [1, 10]
    @test slice_2(np) == [2, 20]
    @test slice_3(np) == [3, 30]
    @test slice_4(np) == [4, 40]
    @test slice_5(np) == [5, 50]
    @test slice_6(np) == [6, 60]
    @test gvec([[0.,0.,0]], bodies)[1][1] ≈ -0.005898468352763
end
nothing

max_stars_repo_name: hustf/Gravityfield

max_issues_repo_name: hustf/Gravityfield

max_forks_repo_name: hustf/Gravityfield


text:
pprint(it...;sep=" ",ed="\n")=length(it)>0 ? (@inbounds for i in 1:length(it)-1;print(it[i],sep);end;print(last(it),ed)) : print(ed)

n=parse(Int64, readline())
pprint((n)*(n-1)÷2)

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
function fit(::Type{Exponential}, m::AbstractMoments)
    # https://en.wikipedia.org/wiki/Exponential_distribution
    n_moments(m) >= 1 || error("Need mean to estimate exponential")
    return Exponential(mean(m))
end

function fit(::Type{Exponential}, lower::QuantilePoint, upper::QuantilePoint)
    # return average for the two quantiles
    θ_lower = -lower.q/log(1-lower.p)
    θ_upper = -upper.q/log(1-upper.p)
    θ = (θ_lower + θ_upper)/2
    Exponential(θ)
end

function fit_mean_quantile(::Type{Exponential}, mean::Real, qp::QuantilePoint)
    # only fit to mean
    fit(Type{Exponential}, AbstractMoments(mean))
end

function fit_mode_quantile(::Type{Exponential}, mode::Real, qp::QuantilePoint)
    # ignore mode (its always at 0)
    θ = -qp.q/log(1-qp.p)
    Exponential(θ)
end


text:
using Distributed
@everywhere using DelimitedFiles
@everywhere using DifferentialEquations
@everywhere using ProgressMeter
include("where_did_it_come_from_parallel.jl")
@everywhere using PyCall
pygui(:qt)
@everywhere using PyPlot
pygui(true)

@everywhere max_hit=3
@everywhere Energy=.25

@everywhere H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
@everywhere Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
@everywhere ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
@everywhere ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
@everywhere function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
return
end
@everywhere condition(u,t,integrator)= u[5]==max_hit
@everywhere function condition2(u,t,integrator) # Event when event_f(u,t) == 0
u[1]
end
@everywhere condition3b(u,t,integrator)=  u[5]==2
@everywhere condition3f(u,t,integrator)=  u[5]==2

@everywhere affect!(integrator) = terminate!(integrator)
@everywhere function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
end
# @everywhere cb2 = ContinuousCallback(condition2,affect2!, save_positions = (false,true) )
@everywhere cb2 = ContinuousCallback(condition2, affect2!, rootfind=true, save_positions = (true,false))
@everywhere cb2f = ContinuousCallback(condition2, affect2!, rootfind=true,  save_positions = (true,false))
@everywhere cb2b = ContinuousCallback(condition2, affect2!,rootfind=true, save_positions = (true,false))

@everywhere cb1 = DiscreteCallback(condition,affect!)
@everywhere cb3b = DiscreteCallback(condition3b,affect!)
@everywhere cb3f = DiscreteCallback(condition3f,affect!)

@everywhere cb=CallbackSet(cb2,cb1)

@everywhere cb_backward=CallbackSet(cb2b,cb3b)
@everywhere cb_forward=CallbackSet(cb2f,cb3f)

@everywhere X_check=zeros(0)
@everywhere Y_check=zeros(0)
@everywhere Q_check=zeros(0)
@everywhere P_check=zeros(0)
@everywhere hit_check1=zeros(0)
@everywhere hit_check2=zeros(0)

@everywhere Q_init=zeros(0)
@everywhere P_init=zeros(0)
@everywhere Y_initS=zeros(0)
@everywhere Y_initB=zeros(0)
@everywhere test=zeros(0)

@everywhere Y_from_nowhere=zeros(0)
@everywhere Q_from_nowhere=zeros(0)
@everywhere P_from_nowhere=zeros(0)

@everywhere X_1before=zeros(0)
@everywhere Y_1before=zeros(0)
@everywhere Q_1before=zeros(0)
@everywhere P_1before=zeros(0)
# ArrQ[k]=.1*rand()
# ArrP[j]=.01*rand()
@everywhere Q_end=1.5
@everywhere Q_start=-Q_end
@everywhere P_end=3
@everywhere P_start=-P_end
@everywhere n_iter_P=8001
@everywhere n_iter_Q=8001
@everywhere N=n_iter_P*n_iter_Q

@everywhere delta_P=(P_end-P_start)/(n_iter_P)
@everywhere delta_Q=(Q_end-Q_start)/(n_iter_Q)

@everywhere ArrP=P_start: delta_P:P_end
@everywhere ArrQ=Q_start: delta_Q: Q_end

T=@showprogress pmap(where_did_it_come_from_parallel,1:N)

for i=1:N
    Q=ArrQ[Int(ceil(i/n_iter_Q))]
    P=ArrP[mod(i,n_iter_P)+1]
    if T[i]!=nothing
        push!(test,0)
        # if T[i][5]<2
        #     push!(Q_from_nowhere,Q)
        #     push!(P_from_nowhere,P)
        #     push!(Y_from_nowhere, T[i][6])
        # elseif T[i][5]==2
            push!(Q_init,Q)
            push!(P_init,P)
            push!(Y_initS,T[i][6])
            # push!(Y_initB,T[i][7])

            q=T[i][1]
            p=T[i][2]
            y=Yfind(q,p,.25)
            # if ~isempty(y) && ~isempty(yB)
            # y_min=min(abs(yS),abs(yB))
            if ~isempty(y) && T[i][5]==2 && abs(T[i][4]-y[1])<1e-2
            push!(Q_1before,T[i][1])
            push!(P_1before,T[i][2])
            push!(X_1before,T[i][3])
            push!(Y_1before,T[i][4])
            # push!(Q_check,T[i][6])
            # push!(P_check,T[i][7])
            # push!(X_check,T[i][8])
            # push!(Y_check,T[i][9])
            push!(hit_check1,T[i][5])
             # push!(hit_check2,T[i][6])
         end
    end
end

figure()
 axis([ -3, 3,-1.5, 1.5])
Q_exit=vcat(Q_init,Q_init,-Q_init,-Q_init)
P_exit=vcat(P_init,-P_init,P_init,-P_init)
# plot(P_exit,Q_exit ,color="r",".",markersize=2, markeredgewidth=.1)
plot(P_init,Q_init ,color="r",".",markersize=2, markeredgewidth=.1)

figure()
axis([ -3, 3,-1.5, 1.5])
Q_previous=vcat(Q_1before,Q_1before,-Q_1before,-Q_1before)
P_previous=vcat(P_1before,-P_1before,P_1before,-P_1before)
# plot(P_previous,Q_previous ,color="b",".",markersize=2, markeredgewidth=.1)
plot(P_1before,Q_1before ,color="b",".",markersize=2, markeredgewidth=.1)

# figure()
# axis([ -3, 3,-1.5, 1.5])
# Q_nowhere=vcat(Q_from_nowhere,Q_from_nowhere,-Q_from_nowhere,-Q_from_nowhere)
# P_nowhere=vcat(P_from_nowhere,-P_from_nowhere,P_from_nowhere,-P_from_nowhere)
# plot(P_nowhere,Q_nowhere ,color="r",".",markersize=2, markeredgewidth=.1)
#
# figure()
# axis([ -3, 3,-1.5, 1.5])
# Qcheck=vcat(Q_check,Q_check,-Q_check,-Q_check)
# Pcheck=vcat(P_check,-P_check,P_check,-P_check)
# plot(Pcheck,Qcheck ,color="k",".",markersize=2, markeredgewidth=.1)

figure()
axis([ -3, 3,-1.5, 1.5])
# plot(P_exit,Q_exit ,color="r",".",markersize=2, markeredgewidth=.1)
# plot(P_previous,Q_previous ,color="b",".",markersize=2, markeredgewidth=.1)
#plot(P_nowhere,Q_nowhere ,color="r",".",markersize=2, markeredgewidth=.1)
# plot(Pcheck,Qcheck ,color="k",".",markersize=2, markeredgewidth=.1)
plot(P_init,Q_init ,color="r",".",markersize=2, markeredgewidth=.1)
plot(P_1before,Q_1before ,color="b",".",markersize=2, markeredgewidth=.1)

#
# outfile_0 = "ParVern9_outfile0_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
# outfile_1 = "ParVern9_outfile1_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
# outfile_2 = "ParVern9_outfile2_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
# outfile_3 = "ParVern9_outfile3_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
# outfile_4 = "ParVern9_outfile4_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
# outfile_5 = "ParVern9_outfile5_"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
# outfile_bound = "ParVern9_outfile1_bound"*string(Energy)*"_nP_"*string(n_iter_P)*"_nQ_"*string(n_iter_Q)*".dat"
#
#
# a = open(outfile_0, "w")
# b = open(outfile_1, "w")
# c = open(outfile_2, "w")
# d = open(outfile_3, "w")
# e = open(outfile_4, "w")
# f = open(outfile_5, "w")
# g = open(outfile_bound, "w")
#
# writedlm(outfile_0, T0)
# writedlm(outfile_1, T1)
# writedlm(outfile_2, T2)
# writedlm(outfile_3, T3)
# writedlm(outfile_4, T4)
# writedlm(outfile_5, T5)
# writedlm(outfile_bound, Tbound)
#
#
# close(a)
# close(b)
# close(c)
# close(d)
# close(e)
# close(f)
# close(g)


text:
using Test, Turing, Distributions
using Turing: VarInfo, gradient_forward, gradient_reverse

@model foo_ad() = begin
    x ~ Normal(3, 1)
    y ~ Normal(x, 1)
end

# Check that gradient_forward doesn't change the RV values or logp of a VarInfo.
let
    # Construct model.
    model, sampler, vi = foo_ad(), nothing, VarInfo()
    model(vi, sampler)

    # Record initial values.
    θ, logp = deepcopy(vi.vals), deepcopy(vi.logp)

    # Compute gradient using ForwardDiff.
    gradient_forward(deepcopy(θ), vi, model, sampler)

    # Check that θ and logp haven't changed.
    @test θ == vi.vals
    @test logp == vi.logp
end

# Check that gradient_reverse doesn't change the RV values or logp of a VarInfo.
let
    # Construct model
    model, sampler, vi = foo_ad(), nothing, VarInfo()
    model(vi, sampler)

    # Record initial values.
    θ, logp = deepcopy(vi.vals), deepcopy(vi.logp)

    # Compute gradient using ForwardDiff.
    gradient_reverse(deepcopy(θ), vi, model, sampler)

    # Check that θ and logp haven't changed.
    @test θ == vi.vals
    @test logp == vi.logp
end