text:
########################################
# Adapters
########################################

abstract type Postgres end
struct PostgresODBC <: Postgres end

struct SQLite end

struct Adapter{AT, CT}
    conn::CT
    adapter::Type{AT}
end

########################################
# Table building functions
########################################

struct DBField
    name::Symbol
    type::Type
    nullable::Bool
end

struct DeconstructClosure{PT}
    fields::Vector{DBField}
    prefix::Symbol
    parenttype::Type{PT}
end

DeconstructClosure(PT; prefix = Symbol()) = DeconstructClosure(DBField[], prefix, PT)

function (f::DeconstructClosure)(i, nm, TT; kw...)
    nametypeindex!(StructTypes.StructType(TT), TT, i, nm, f)
    return
end

function deconstruct(PT)
    c = DeconstructClosure(PT)
    StructTypes.foreachfield(c, PT)
    return c
end

function nametypeindex!(ST, T, i, nm, c)
    dbfield = DBField(Symbol(c.prefix, nm), T, false)
    push!(c.fields, dbfield)
    return
end

function nametypeindex!(::Union{StructTypes.Struct, StructTypes.Mutable}, ::Type{Union{Nothing, T}}, i, nm, c) where T
    dbfield = DBField(Symbol(c.prefix, nm), T, true)
    push!(c.fields, dbfield)
    return
end

function nametypeindex!(::Union{StructTypes.Struct, StructTypes.Mutable}, ::Type{Union{Missing, T}}, i, nm, c) where T
    dbfield = DBField(Symbol(c.prefix, nm), T, true)
    push!(c.fields, dbfield)
    return
end

function nametypeindex!(::Union{StructTypes.Struct, StructTypes.Mutable}, T, i, nm, c)
    prefix = Symbol(c.prefix, StructTypes.fieldprefix(c.parenttype, nm))
    c2 = DeconstructClosure(c.fields, prefix, T)
    StructTypes.foreachfield(c2, T)
    return
end

function build_autoincrement_column!(iob, adapter::Adapter{DBT}, T, field) where {DBT <: Postgres}
    print(iob, field.name, " ", dbautoincrement(adapter, field.type), idproperty(T) == field.name ? " PRIMARY KEY" : "")
end

function build_column!(iob, adapter, T, field)
    if autoincrement(T) == field.name
        build_autoincrement_column!(iob, adapter, T, field)
        return nothing
    end

    print(iob, field.name, " ", dbtype(adapter, field.type), field.nullable ? "" : " NOT NUll", idproperty(T) == field.name ? " PRIMARY KEY" : "")

    return nothing
end

function build_create_table(adapter, T)
    iob = IOBuffer()
    print(iob, "CREATE TABLE ", tablename(T), " (\n")
    c = deconstruct(T)
    field, rest = Iterators.peel(c.fields)
    build_column!(iob, adapter, T, field)
    for field in rest
        print(iob, ",\n")
        build_column!(iob, adapter, T, field)
    end
    print(iob, "\n)")

    query = String(take!(iob))
end

function create_table(adapter, T)
    query = build_create_table(adapter, T)
    execute(adapter, query)
end

########################################
# Auxiliary table functions
########################################
tablename(T) = lowercase(string(Symbol(T))) * "s"

dbtype(::Adapter{DBT}, ::Type{Int}) where {DBT <: Postgres} = "BIGINT"
dbtype(::Adapter{DBT}, ::Type{String}) where {DBT <: Postgres} = "TEXT"
dbautoincrement(::Adapter{DBT}, ::Type{Int}) where {DBT <: Postgres} = "BIGSERIAL"

idproperty(::Type{T}) where T = :_
autoincrement(::Type{T}) where T = :_

########################################
# Interaction utils
########################################

function execute(adapter::Adapter, query, vals = nothing)
    if vals === nothing || isempty(vals)
        res = DBInterface.execute(adapter.conn, query)
    else
        stmt = prepare(adapter, query)
        res = DBInterface.execute(stmt, vals)
        # TODO: I do not like that I have hanging statements, but this is the only way
        # to ensure select results...
        # DBInterface.close!(stmt)
    end

    res
end

function execute(f, adapter::Adapter, query, vals = nothing)
    if vals === nothing || isempty(vals)
        res = f(DBInterface.execute(adapter.conn, query))
    else
        stmt = prepare(adapter, query)
        res = f(DBInterface.execute(stmt, vals))
        # TODO: I do not like that I have hanging statements, but this is the only way
        # to ensure select results...
        DBInterface.close!(stmt)
    end

    res
end

prepare(adapter::Adapter, query) = DBInterface.prepare(adapter.conn, query)

function build_insert(adapter::Adapter, x::T) where T
    c = deconstruct(T)
    iob = IOBuffer()

    aicol = autoincrement(T)

    print(iob, "INSERT INTO ", tablename(T), " (")
    isfirst = true
    isaicol = false
    cnt = 0
    for field in c.fields
        if field.name == aicol
            isaicol = true
            continue
        end
        cnt += 1
        print(iob, isfirst ? "" : ", ", field.name)
        isfirst = false
    end

    print(iob, ") VALUES (")

    isfirst = true
    for i in 1:cnt
        print(iob, isfirst ? "?" : ",?")
        isfirst = false
    end
    print(iob, ")")

    if isaicol
        print(iob, " RETURNING ", aicol)
    end
    query = String(take!(iob))
end

function insert(adapter::Adapter, x::T) where T
    # TODO: All of this is so so so slow... Should be implemented better!
    # Currently it's just proof of concept
    query = build_insert(adapter, x)
    c = deconstruct(T)

    protovals = Strapping.deconstruct(x) |> only |> values
    vals = []
    aicol = autoincrement(T)
    isaicol = false
    for (i, v) in enumerate(protovals)
        if c.fields[i].name == aicol
            isaicol = true
            continue
        end
        push!(vals, v)
    end
    x = execute(adapter, query, vals) do res
        if isaicol
            l = Setfield.PropertyLens{aicol}()
            set(x, l, only(only(res.columns)))
        end
    end

    return x
end

function build_select(adapter::Adapter, ::Type{T}, filt = ()) where T
    iob = IOBuffer()
    print(iob, "SELECT * FROM ", tablename(T))

    if !isempty(filt)
        isfirst = true
        print(iob, " WHERE")
        for (k, _) in filt
            print(iob, isfirst ? " " : " AND ", k, " = ?")
        end
    end

    return String(take!(iob))
end

function select(adapter::Adapter, ::Type{T}, filt = ()) where T
    query = build_select(adapter, T, filt)
    vals = map(x -> x[2], filt)

    execute(adapter, query, vals) do dbres
        Strapping.construct(T, dbres)
    end
end

function select(adapter::Adapter, ::Type{Vector{T}}, filt = ()) where T
    query = build_select(adapter, T, filt)
    vals = map(x -> x[2], filt)

    execute(adapter, query, vals) do dbres
        Strapping.construct(Vector{T}, dbres)
    end
end

function build_delete(adapter::Adapter, ::Type{T}, filt) where T
    iob = IOBuffer()

    print(iob, "DELETE FROM ", tablename(T), " WHERE ")
    isfirst = true
    for (k, _) in filt
        print(iob, isfirst ? "" : " AND ", k, " = ?")
        isfirst = false
    end

    return String(take!(iob))
end

function delete(adapter::Adapter, ::Type{T}, filt) where T
    isempty(filt) && return nothing
    query = build_delete(adapter, T, filt)
    vals = map(x -> x[2], filt)
    execute(identity, adapter, query, vals)

    return nothing
end

function delete(adapter::Adapter, x::T, key = ()) where T
    if isempty(key)
        # Should delete over primary key
        id = getfield(x, idproperty(T))
        query = "DELETE FROM $(tablename(T)) WHERE $(idproperty(T)) = ?"
        execute(identity, adapter, query, (id, ))
    else
        # Not implemented
    end
end

function build_upsert(adapter::Adapter, x::T, onconflict = ()) where T
    c = deconstruct(T)
    iob = IOBuffer()

    pkcol = idproperty(T)

    print(iob, "INSERT INTO ", tablename(T), " (")
    isfirst = true
    cnt = 0
    for field in c.fields
        cnt += 1
        print(iob, isfirst ? "" : ", ", field.name)
        isfirst = false
    end

    print(iob, ") VALUES (")

    isfirst = true
    for i in 1:cnt
        print(iob, isfirst ? "?" : ",?")
        isfirst = false
    end
    print(iob, ")")

    print(iob, " ON CONFLICT")
    if isempty(onconflict)
        print(iob, " DO NOTHING")
    else
        print(iob, "(", pkcol, ")")
        print(iob, " DO UPDATE SET ")
        isfirst = true
        for key in onconflict
            print(iob, isfirst ? "" : ", ", key, " = EXCLUDED.", key)
            isfirst = false
        end
    end

    query = String(take!(iob))
end

function upsert(adapter::Adapter, x::T, onconflict = ()) where T
    # Yeah... I'll hardcode primary key upsert. Should be done better
    query = build_upsert(adapter, x, onconflict)
    c = deconstruct(T)

    vals = Strapping.deconstruct(x) |> only |> values
    res = execute(identity, adapter, query, vals)

    return nothing
end


text:
"""
Constructs a Differentiable Optimizer model from a MOI Optimizer.
Supports `forward` and `backward` methods for solving and differentiating the model respectectively.

## Note
Currently supports differentiating linear and quadratic programs only.
"""

Base.@kwdef struct QPCache
    lhs::SparseMatrixCSC{Float64, Int}
end

Base.@kwdef struct ConicCache
    M::SparseMatrixCSC{Float64, Int}
    vp::Vector{Float64}
    Dπv::BlockDiagonals.BlockDiagonal{Float64, Matrix{Float64}}
    A::SparseMatrixCSC{Float64, Int}
    b::Vector{Float64}
    c::Vector{Float64}
end

Base.@kwdef struct QPForwBackCache
    dz::Vector{Float64}
    dλ::Vector{Float64}
    dν::Vector{Float64}
end
Base.@kwdef struct ConicForwCache
    du::Vector{Float64}
    dv::Vector{Float64}
    dw::Vector{Float64}
end
Base.@kwdef struct ConicBackCache
    g::Vector{Float64}
    πz::Vector{Float64}
end

const MOIDD = MOI.Utilities.DoubleDicts

Base.@kwdef mutable struct DiffInputCache
    dx::Dict{VI, Float64} = Dict{VI, Float64}()# dz for QP
    # ds
    # dy #= [d\lambda, d\nu] for QP
    # FIXME Would it be possible to have a DoubleDict where the value depends
    #       on the function type ? Here, we need to create two dicts to have
    #       concrete value types.
    # `scalar_constraints` and `vector_constraints` includes `A` and `b` for CPs
    # or `G` and `h` for QPs
    scalar_constraints::MOIDD.DoubleDict{MOI.ScalarAffineFunction{Float64}} = MOIDD.DoubleDict{MOI.ScalarAffineFunction{Float64}}() # also includes G for QPs
    vector_constraints::MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}} = MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}}() # also includes G for QPs
    objective::Union{Nothing,MOI.AbstractScalarFunction} = nothing
end

function Base.empty!(cache::DiffInputCache)
    empty!(cache.dx)
    empty!(cache.scalar_constraints)
    empty!(cache.vector_constraints)
    cache.objective = nothing
    return
end

"""
    ForwardInObjective <: MOI.AbstractModelAttribute

A `MOI.AbstractModelAttribute` to set input data to forward differentiation, that
is, problem input data.
The possible values are any `MOI.AbstractScalarFunction`.
A `MOI.ScalarQuadraticFunction` can only be used in linearly constrained
quadratic models.

For instance, if the objective contains `θ * (x + 2y)`, for the purpose of
computing the derivative with respect to `θ`, the following should be set:
```julia
MOI.set(model, DiffOpt.ForwardInObjective(), 1.0 * x + 2.0 * y)
```
where `x` and `y` are the relevant `MOI.VariableIndex`.
"""
struct ForwardInObjective <: MOI.AbstractModelAttribute end

"""
    ForwardInConstraint <: MOI.AbstractConstraintAttribute

A `MOI.AbstractConstraintAttribute` to set input data to forward differentiation, that
is, problem input data.

For instance, if the scalar constraint of index `ci` contains `θ * (x + 2y) <= 5θ`,
for the purpose of computing the derivative with respect to `θ`, the following
should be set:
```julia
MOI.set(model, DiffOpt.ForwardInConstraint(), ci, 1.0 * x + 2.0 * y - 5.0)
```
Note that we use `-5` as the `ForwardInConstraint` sets the tangent of the
ConstraintFunction so we consider the expression `θ * (x + 2y - 5)`.
"""
struct ForwardInConstraint <: MOI.AbstractConstraintAttribute end

"""
    ForwardOutVariablePrimal <: MOI.AbstractVariableAttribute

A `MOI.AbstractVariableAttribute` to get output data from forward
differentiation, that is, problem solution.

For instance, to get the tangent of the variable of index `vi` corresponding to
the tangents given to `ForwardInObjective` and `ForwardInConstraint`, do the
following:
```julia
MOI.get(model, DiffOpt.ForwardOutVariablePrimal(), vi)
```
"""
struct ForwardOutVariablePrimal <: MOI.AbstractVariableAttribute end
MOI.is_set_by_optimize(::ForwardOutVariablePrimal) = true

"""
    BackwardInVariablePrimal <: MOI.AbstractVariableAttribute

A `MOI.AbstractVariableAttribute` to set input data to backward
differentiation, that is, problem solution.

For instance, to set the tangent of the variable of index `vi`, do the
following:
```julia
MOI.set(model, DiffOpt.BackwardInVariablePrimal(), x)
```
"""
struct BackwardInVariablePrimal <: MOI.AbstractVariableAttribute end

"""
    BackwardOutObjective <: MOI.AbstractModelAttribute

A `MOI.AbstractModelAttribute` to get output data to backward differentiation,
that is, problem input data.

For instance, to get the tangent of the objective function corresponding to
the tangent given to `BackwardInVariablePrimal`, do the
following:
```julia
func = MOI.get(model, DiffOpt.BackwardOutObjective)
```
Then, to get the sensitivity of the linear term with variable `x`, do
```julia
JuMP.coefficient(func, x)
```
To get the sensitivity with respect to the quadratic term with variables `x`
and `y`, do either
```julia
JuMP.coefficient(func, x, y)
```
or
```julia
DiffOpt.quad_sym_half(func, x, y)
```

!!! warning
    These two lines are **not** equivalent in case `x == y`, see
    [`quad_sym_half`](@ref) for the details on the difference between these two
    functions.
"""
struct BackwardOutObjective <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::BackwardOutObjective) = true

"""
    BackwardOutConstraint

An `MOI.AbstractConstraintAttribute` to get output data to backward differentiation, that
is, problem input data.

For instance, if the following returns `x + 2y + 5`, it means that the tangent
has coordinate `1` for the coefficient of `x`, coordinate `2` for the
coefficient of `y` and `5` for the function constant.
If the constraint is of the form `func == constant` or `func <= constant`,
the tangent for the constant on the right-hand side is `-5`.
```julia
MOI.get(model, DiffOpt.BackwardOutConstraint(), ci)
```
"""
struct BackwardOutConstraint <: MOI.AbstractConstraintAttribute end
MOI.is_set_by_optimize(::BackwardOutConstraint) = true

"""
    @enum ProgramClassCode QUADRATIC CONIC AUTOMATIC

Program class used by DiffOpt. DiffOpt implements differentiation of two
different program class:
1) Quadratic Program (QP): quadratic objective and linear constraints and
2) Conic Program (CP): linear objective and conic constraints.

`AUTOMATIC` which means that the class will be automatically selected given the
problem data: if any constraint is conic, CP is used and QP is used otherwise.
See [`ProgramClass`](@ref).
"""
@enum ProgramClassCode QUADRATIC CONIC AUTOMATIC

abstract type DiffModel <: MOI.ModelLike end

MOI.supports_incremental_interface(::DiffModel) = true

MOI.is_valid(model::DiffModel, idx::MOI.Index) = MOI.is_valid(model.model, idx)

function MOI.add_variable(model::DiffModel)
    return MOI.add_variable(model.model)
end

function MOI.add_variables(model::DiffModel, n)
    return MOI.add_variables(model.model, n)
end

function MOI.Utilities.pass_nonvariable_constraints(
    dest::DiffModel,
    src::MOI.ModelLike,
    idxmap::MOIU.IndexMap,
    constraint_types,
)
    MOI.Utilities.pass_nonvariable_constraints(dest.model, src, idxmap, constraint_types)
end

function MOI.Utilities.final_touch(model::DiffModel, index_map)
    MOI.Utilities.final_touch(model.model, index_map)
end

function MOI.supports_constraint(model::DiffModel, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    return MOI.supports_constraint(model.model, F, S)
end

function MOI.add_constraint(model::DiffModel, func::MOI.AbstractFunction, set::MOI.AbstractSet)
    return MOI.add_constraint(model.model, func, set)
end

function _enlarge_set(vec::Vector, idx, value)
    m = last(idx)
    if length(vec) < m
        n = length(vec)
        resize!(vec, m)
        fill!(view(vec, (n+1):m), NaN)
    end
    vec[idx] = value
    return
end

function MOI.set(model::DiffModel, ::MOI.VariablePrimalStart, vi::MOI.VariableIndex, value)
    MOI.throw_if_not_valid(model, vi)
    _enlarge_set(model.x, vi.value, value)
end

function MOI.set(model::DiffModel, ::ForwardInObjective, objective)
    model.input_cache.objective = objective
    return
end
function MOI.set(model::DiffModel, ::BackwardInVariablePrimal, vi::VI, val)
    model.input_cache.dx[vi] = val
    return
end
function MOI.set(model::DiffModel,
    ::ForwardInConstraint,
    ci::CI{MOI.ScalarAffineFunction{T},S},
    func::MOI.ScalarAffineFunction{T},
) where {T,S}
    model.input_cache.scalar_constraints[ci] = func
    return
end
function MOI.set(model::DiffModel,
    ::ForwardInConstraint,
    ci::CI{MOI.VectorAffineFunction{T},S},
    func::MOI.VectorAffineFunction{T},
) where {T,S}
    model.input_cache.vector_constraints[ci] = func
    return
end

function lazy_combination(op::F, α, a, β, b) where {F<:Function}
    return LazyArrays.ApplyArray(
        op,
        LazyArrays.@~(α .* a),
        LazyArrays.@~(β .* b),
    )
end
# Workaround for Julia v1.0
@static if VERSION < v"1.6"
    _view(x, I) = x[I]
else
    _view(x, I) = view(x, I)
end
function lazy_combination(op::F, α, a, β, b, I::UnitRange) where {F<:Function}
    return lazy_combination(op, α, _view(a, I), β, _view(b, I))
end
function lazy_combination(op::F, a, b, i::Integer, args::Vararg{Any,N}) where {F<:Function,N}
    return lazy_combination(op, a[i], b, b[i], a, args...)
end
function lazy_combination(op::F, a, b, i::UnitRange, I::UnitRange) where {F<:Function}
    return lazy_combination(op, _view(a, i), b, _view(b, i), a, I)
end

_lazy_affine(vector, constant::Number) = VectorScalarAffineFunction(vector, constant)
_lazy_affine(matrix, vector) = MatrixVectorAffineFunction(matrix, vector)
function MOI.get(model::DiffModel, ::BackwardOutConstraint, ci::CI)
    return _lazy_affine(_get_dA(model, ci), _get_db(model, ci))
end

"""
    π(v::Vector{Float64}, model::MOI.ModelLike, cones::ProductOfSets)

Given a `model`, its `cones`, find the projection of the vectors `v`
of length equal to the number of rows in the conic form onto the cartesian
product of the cones corresponding to these rows.
For more info, refer to https://github.com/matbesancon/MathOptSetDistances.jl
"""
function π(v::Vector{T}, model::MOI.ModelLike, cones::ProductOfSets) where T
    return map_rows(model, cones, Flattened{T}()) do ci, r
        MOSD.projection_on_set(
            MOSD.DefaultDistance(),
            v[r],
            MOI.dual_set(MOI.get(model, MOI.ConstraintSet(), ci))
        )
    end
end

"""
    Dπ(v::Vector{Float64}, model, cones::ProductOfSets)

Given a `model`, its `cones`, find the gradient of the projection of
the vectors `v` of length equal to the number of rows in the conic form onto the
cartesian product of the cones corresponding to these rows.
For more info, refer to https://github.com/matbesancon/MathOptSetDistances.jl
"""
function Dπ(v::Vector{T}, model::MOI.ModelLike, cones::ProductOfSets) where T
    return BlockDiagonals.BlockDiagonal(
        map_rows(model, cones, Nested{Matrix{T}}()) do ci, r
            MOSD.projection_gradient_on_set(
                MOSD.DefaultDistance(),
                v[r],
                MOI.dual_set(MOI.get(model, MOI.ConstraintSet(), ci)),
            )
        end
    )
end

# See the docstring of `map_rows`.
struct Nested{T} end
struct Flattened{T} end

# Store in `x` the values `y` corresponding to the rows `r` and the `k`th
# constraint.
function _assign_mapped!(x, y, r, k, ::Nested)
    x[k] = y
end
function _assign_mapped!(x, y, r, k, ::Flattened)
    x[r] = y
end

# Map the rows corresponding to `F`-in-`S` constraints and store it in `x`.
function _map_rows!(f::Function, x::Vector, model, cones::ProductOfSets, ::Type{F}, ::Type{S}, map_mode, k) where {F, S}
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        r = MOI.Utilities.rows(cones, ci)
        k += 1
        _assign_mapped!(x, f(ci, r), r, k, map_mode)
    end
    return k
end

# Allocate a vector for storing the output of `map_rows`.
_allocate_rows(cones, ::Nested{T}) where {T} = Vector{T}(undef, length(cones.dimension))
_allocate_rows(cones, ::Flattened{T}) where {T} = Vector{T}(undef, MOI.dimension(cones))

"""
    map_rows(f::Function, model, cones::ProductOfSets, map_mode::Union{Nested{T}, Flattened{T}})

Given a `model`, its `cones` and `map_mode` of type `Nested` (resp.
`Flattened`), return a `Vector{T}` of length equal to the number of cones (resp.
rows) in the conic form where the value for the index (resp. rows) corresponding
to each cone is equal to `f(ci, r)` where `ci` is the corresponding constraint
index in `model` and `r` is a `UnitRange` of the corresponding rows in the conic
form.
"""
function map_rows(f::Function, model, cones::ProductOfSets, map_mode::Union{Nested, Flattened})
    x = _allocate_rows(cones, map_mode)
    k = 0
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        # Function barrier for type unstability of `F` and `S`
        # `con_map` is a `MOIU.DoubleDicts.MainIndexDoubleDict`, we index it at `F, S`
        # which returns a `MOIU.DoubleDicts.IndexWithType{F, S}` which is type stable.
        # If we have a small number of different constraint types and many
        # constraint of each type, this mostly removes type unstabilities
        # as most the time is in `_map_rows!` which is type stable.
        k = _map_rows!(f, x, model, cones, F, S, map_mode, k)
    end
    return x
end

function _fill(neg::Function, gradient_cache, input_cache, cones, args...)
    _fill(S -> true, neg, gradient_cache, input_cache, cones, args...)
end
function _fill(filter::Function, neg::Function, gradient_cache, input_cache, cones, args...)
    for (F, S) in keys(input_cache.scalar_constraints.dict)
        filter(S) || continue
        _fill(args..., neg(S), cones, input_cache.scalar_constraints[F,S])
    end
    for (F, S) in keys(input_cache.vector_constraints.dict)
        _fill(args..., neg(S), cones, input_cache.vector_constraints[F,S])
    end
    return
end

function _fill(vector::Vector, neg::Bool, cones, dict)
    for (ci, func) in dict
        r = MOI.Utilities.rows(cones, ci)
        vector[r] = neg ? -MOI.constant(func) : MOI.constant(func)
    end
end
function _fill(I::Vector, J::Vector, V::Vector, neg::Bool, cones, dict)
    for (ci, func) in dict
        r = MOI.Utilities.rows(cones, ci)
        for term in func.terms
            _push_term(I, J, V, neg, r, term)
        end
    end
end
function _push_term(I::Vector, J::Vector, V::Vector, neg::Bool, r::Integer, term::MOI.ScalarAffineTerm)
    push!(I, r)
    push!(J, term.variable.value)
    push!(V, neg ? -term.coefficient : term.coefficient)
end
function _push_term(I::Vector, J::Vector, V::Vector, neg::Bool, r::UnitRange, term::MOI.VectorAffineTerm)
    _push_term(I, J, V, neg, r[term.output_index], term.scalar_term)
end

function MOI.supports(model::DiffModel, attr::MOI.AbstractModelAttribute)
    return MOI.supports(model.model, attr)
end

function MOI.set(model::DiffModel, attr::MOI.AbstractModelAttribute, value)
    MOI.set(model.model, attr, value)
end

function MOI.get(model::DiffModel, attr::MOI.AbstractModelAttribute)
    return MOI.get(model.model, attr)
end

ext: jl


text:
## ExoplanetsSysSim/src/transit_prob_geometric.jl
## (c) 2015 Eric B. Ford

# Transit probability expressions if treat planets separately
function calc_transit_prob_single_planet_approx(P::Float64, Rstar::Float64, Mstar::Float64)
  return min(rsol_in_au*Rstar/semimajor_axis(P,Mstar), 1.0)
end

function calc_transit_prob_single_planet_obs_ave(ps::PlanetarySystemAbstract, pl::Integer)
 ecc::Float64 = ps.orbit[pl].ecc
 a::Float64 = semimajor_axis(ps,pl)
 Rstar::Float64 = rsol_in_au*ps.star.radius
 return min(Rstar/(a*(1-ecc)*(1+ecc)), 1.0)
end
calc_transit_prob_single_planet_obs_ave(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_prob_single_planet_obs_ave(t.sys[s], p)

#=
# WARNING: This knows about e and w, but still returns a fraction rather than a 0 or 1.  Commented out for now, so no one uses it accidentally until we figure out why it was this way
function calc_transit_prob_single_planet_one_obs(ps::PlanetarySystemAbstract, pl::Integer)
 ecc::Float64 = ps.orbit[pl].ecc
 a::Float64 = semimajor_axis(ps,pl)
 Rstar::Float64 = rsol_in_au*ps.star.radius
 return min(Rstar*(1+ecc*sin(ps.orbit[pl].omega))/(a*(1-ecc)*(1+ecc)), 1.0)
end
calc_transit_prob_single_planet_one_obs(t::KeplerTarget, s::Integer, p::Integer) = calc_transit_prob_single_planet_one_obs(t.sys[s], p)
=#

# WARNING: Assumes that planets with b>1 won't be detected/pass vetting
function does_planet_transit(ps::PlanetarySystemAbstract, pl::Integer)
   ecc::Float64 = ps.orbit[pl].ecc
   incl::Float64 = ps.orbit[pl].incl
   a::Float64 = semimajor_axis(ps,pl)
   Rstar::Float64 = rsol_in_au*ps.star.radius
   if (Rstar >= (a*(1-ecc)*(1+ecc))/(1+ecc*sin(ps.orbit[pl].omega))*cos(incl))
     return true
   else
     return false
   end
end

function corbits_placeholder_obs_ave( ps::PlanetarySystemSingleStar, use_pl::Vector{Cint} )    # Might be useful for someone to test w/o CORBITS
  n = num_planets(ps)
  prob = 1.0
  for p in 1:n
     ptr = calc_transit_prob_single_planet_obs_ave(ps,p)
     prob *= (use_pl[p]==1) ? ptr : 1.0-ptr
     #=
     if(use_pl[p]==1)
       prob *= calc_transit_prob_single_planet_obs_ave(ps,p)
     else
       prob *= 1.0-calc_transit_prob_single_planet_obs_ave(ps,p)
     end
     =#
  end
  return prob
end

function calc_impact_parameter(ps::PlanetarySystemSingleStar, pl::Integer)
      one_minus_e2 = (1-ps.orbit[pl].ecc)*(1+ps.orbit[pl].ecc)
      a_semimajor_axis = semimajor_axis(ps,pl)
      b = a_semimajor_axis *cos(ps.orbit[pl].incl)/(ps.star.radius*rsol_in_au)
      b *= one_minus_e2/(1+ps.orbit[pl].ecc*sin(ps.orbit[pl].omega))
end

function prob_combo_transits_one_obs( ps::PlanetarySystemSingleStar, use_pl::Vector{Cint} )
  n = num_planets(ps)
  for p in 1:n
      #one_minus_e2 = (1-ps.orbit[p].ecc)*(1+ps.orbit[p].ecc)
      #a_semimajor_axis = semimajor_axis(ps,p)
      #b = a_semimajor_axis *cos(ps.orbit[p].incl)/ps.star.radius
      #b *= one_minus_e2/(1+ps.orbit[p].ecc*sin(ps.orbit[p].omega))
      b = calc_impact_parameter(ps, p)
      if ! ( (b<=1.0 && use_pl[p]==1) || (b> 1.0 && use_pl[p]!=1) )
        return 0.0
      end
  end
  return 1.0
end

struct prob_combo_transits_obs_ave_workspace_type
  a::Vector{Cdouble}
  r::Vector{Cdouble}
  ecc::Vector{Cdouble}
  Omega::Vector{Cdouble}
  omega::Vector{Cdouble}
  inc::Vector{Cdouble}
  function prob_combo_transits_obs_ave_workspace_type(n::Integer)
     @assert(1<=n<=100)
     new( Array{Cdouble}(n), Array{Cdouble}(n), Array{Cdouble}(n), Array{Cdouble}(n), Array{Cdouble}(n), Array{Cdouble}(n) )
  end
end

#=
# Attempt to reduce memory allocations.  It does that, but no noticable speed improvement, so I'm leaving it commented out for now.
const global corbits_max_num_planets_per_system = 20
prob_combo_transits_obs_ave_workspace =  prob_combo_transits_obs_ave_workspace_type(corbits_max_num_planets_per_system)

function prob_combo_transits_obs_ave( ps::PlanetarySystemSingleStar, use_pl::Vector{Cint}; print_orbit::Bool = false)
  n = num_planets(ps)
  @assert(n<=corbits_max_num_planets_per_system)
  for i in 1:n
       prob_combo_transits_obs_ave_workspace.a[i] = semimajor_axis(ps,i)
       prob_combo_transits_obs_ave_workspace.r[i] = ps.planet[i].radius * rsol_in_au
       prob_combo_transits_obs_ave_workspace.ecc[i] = ps.orbit[i].ecc
       prob_combo_transits_obs_ave_workspace.Omega[i] = ps.orbit[i].asc_node
       prob_combo_transits_obs_ave_workspace.omega[i] =ps.orbit[i].omega
       prob_combo_transits_obs_ave_workspace.inc[i] = ps.orbit[i].incl
  end

  r_star = convert(Cdouble,ps.star.radius *  rsol_in_au )
  prob = prob_of_transits_approx(prob_combo_transits_obs_ave_workspace.a, r_star, prob_combo_transits_obs_ave_workspace.r, prob_combo_transits_obs_ave_workspace.ecc, prob_combo_transits_obs_ave_workspace.Omega, prob_combo_transits_obs_ave_workspace.omega, prob_combo_transits_obs_ave_workspace.inc, use_pl)
  #prob = prob_of_transits_approx(a, r_star, r, ecc, Omega, omega, inc, use_pl)

  if print_orbit
  println("# a = ", prob_combo_transits_obs_ave_workspace.a)
  println("# r_star = ", r_star)
  println("# r = ", prob_combo_transits_obs_ave_workspace.r)
  println("# ecc = ", prob_combo_transits_obs_ave_workspace.ecc)
  println("# Omega = ", prob_combo_transits_obs_ave_workspace.Omega)
  println("# omega = ", prob_combo_transits_obs_ave_workspace.omega)
  println("# inc = ", prob_combo_transits_obs_ave_workspace.inc)
  println("# use_pl = ", use_pl)
  println("")
  end
  return prob
end
=#

function prob_combo_transits_obs_ave( ps::PlanetarySystemSingleStar, use_pl::Vector{Cint}; print_orbit::Bool = false)
  n = num_planets(ps)
  a =  Cdouble[ semimajor_axis(ps,i) for i in 1:n ]
  r_star = convert(Cdouble,ps.star.radius *  rsol_in_au )
  r = Cdouble[ ps.planet[i].radius * rsol_in_au for i in 1:n ]
  ecc = Cdouble[ ps.orbit[i].ecc for i in 1:n ]
  Omega = Cdouble[ ps.orbit[i].asc_node for i in 1:n ]
  omega = Cdouble[ ps.orbit[i].omega for i in 1:n ]
  inc = Cdouble[ ps.orbit[i].incl for i in 1:n ]
  #use_pl = Cint[0 for i in 1:n]
  #for i in 1:length(combo)
  #  use_pl[i] = 1
  #end
  prob = prob_of_transits_approx(a, r_star, r, ecc, Omega, omega, inc, use_pl)

  if print_orbit
  println("# a = ", a)
  println("# r_star = ", r_star)
  println("# r = ", r)
  println("# ecc = ", ecc)
  println("# Omega = ", Omega)
  println("# omega = ", omega)
  println("# inc = ", inc)
  println("# use_pl = ", use_pl)
  println("")
  end
  return prob
end

@compat abstract type SystemDetectionProbsAbstract  end
@compat abstract type SystemDetectionProbsTrait end
@compat abstract type SkyAveraged <: SystemDetectionProbsTrait end
@compat abstract type OneObserver <: SystemDetectionProbsTrait end
# Derived types will allow us to specialize depending on whether using sky-averaged values, values for actual geometry (both of which require the physical catalog), or estimates based on observed data

type SimulatedSystemDetectionProbs{T<:SystemDetectionProbsTrait} <: SystemDetectionProbsAbstract         # To be used for simulated systems where we can calculat everything
  # Inputs to CORBITS
  detect_planet_if_transits::Vector{Float64}       # Probability of detecting each planet, averaged over all observers for each planet individually, assumes b~U[0,1);  To be used in pass 1

  # Outputs from CORBITS
  pairwise::Matrix{Float64}                 # detection probability (incl geometry & detection probability) for each planet (diagonal) and each planet pair (off diagonal)
                                            # TODO OPT: Make matrix symmetric to save memory?
  n_planets::Vector{Float64}                # fraction of time would detect n planets (incl. geometry & detection probability)

  combo_detected::Vector{Vector{Int64}}     # List of combinations of planets drawn from full joint multi-transit probability
end
# Removed since not really need and previous had used alised typename and function
# Specialize so know whether the values are sky averaged or not
#typealias SkyAveragedSystemDetectionProbs SimulatedSystemDetectionProbs{SkyAveraged}
#typealias OneObserverSystemDetectionProbs SimulatedSystemDetectionProbs{OneObserver}
#SkyAveragedSystemDetectionProbs = SimulatedSystemDetectionProbs{SkyAveraged}
#OneObserverSystemDetectionProbs = SimulatedSystemDetectionProbs{OneObserver}

function SimulatedSystemDetectionProbs(traits::Type, p::Vector{Float64}; num_samples::Integer = 1)
  SimulatedSystemDetectionProbs{traits}( p, zeros(length(p),length(p)), zeros(length(p)),  fill(Array{Int64}(0), num_samples) )
end

function SimulatedSystemDetectionProbs(traits::Type, n::Integer; num_samples::Integer = 1)
  SimulatedSystemDetectionProbs{traits}( ones(n), zeros(n,n), zeros(n), fill(Array{Int64}(0), num_samples) )
end

SkyAveragedSystemDetectionProbs(p::Vector{Float64}; num_samples::Integer = 1) = SimulatedSystemDetectionProbs( SkyAveraged, p, num_samples=num_samples)
SkyAveragedSystemDetectionProbs(n::Integer; num_samples::Integer = 1) = SimulatedSystemDetectionProbs(SkyAveraged, n, num_samples=num_samples)
SkyAveragedSystemDetectionProbsEmpty() = SimulatedSystemDetectionProbs(SkyAveraged,0)

OneObserverSystemDetectionProbs(p::Vector{Float64}; num_samples::Integer = 1) = SimulatedSystemDetectionProbs( OneObserver, p, num_samples=num_samples)

OneObserverSystemDetectionProbs(n::Integer; num_samples::Integer = 1) = SimulatedSystemDetectionProbs(OneObserver, n, num_samples=num_samples)
OneObserverSystemDetectionProbsEmpty() = SimulatedSystemDetectionProbs(OneObserver,0)

# Functions common to various types of SystemDetectionProbs
num_planets{T<:SystemDetectionProbsTrait}(prob::SimulatedSystemDetectionProbs{T}) = length(prob.detect_planet_if_transits)
function prob_detect_if_transits{T<:SystemDetectionProbsTrait}(prob::SimulatedSystemDetectionProbs{T}, pl_id::Integer)
  @assert 1<=pl_id<=length(prob.detect_planet_if_transits)
  prob.detect_planet_if_transits[pl_id]
end
function prob_detect{T<:SystemDetectionProbsTrait}(prob::SimulatedSystemDetectionProbs{T}, pl_id::Integer)
  if ! (1<=pl_id<=size(prob.pairwise,1) )
    println("#ERROR:  pl_id =", pl_id, " prob.pairwise= ", prob.pairwise)
  end
  @assert 1<=pl_id<=size(prob.pairwise,1)
  return prob.pairwise[pl_id,pl_id]
end
function prob_detect_both_planets{T<:SystemDetectionProbsTrait}(prob::SimulatedSystemDetectionProbs{T}, pl_id::Integer, ql_id::Integer)
  @assert 1<=pl_id<=size(prob.pairwise,1)
  @assert 1<=ql_id<=size(prob.pairwise,1)
  prob.pairwise[pl_id,ql_id]
end
prob_detect_n_planets{T<:SystemDetectionProbsTrait}(prob::SimulatedSystemDetectionProbs{T}, n::Integer) = 1<=n<=length(prob.n_planets) ? prob.n_planets[n] : 0.0

# Compute sky-averaged transit probabilities from a planetary system with known physical properties, assuming a single host star
function calc_simulated_system_detection_probs(ps::PlanetarySystemSingleStar, prob_det_if_tr::Vector{Float64}; num_samples::Integer = 1, max_tranets_in_sys::Integer = 10, min_detect_prob_to_be_included::Float64 = 0.0, observer_trait::Type=SkyAveraged)
  @assert issubtype(observer_trait,SystemDetectionProbsTrait)
  @assert num_planets(ps) == length(prob_det_if_tr)
  idx_detectable = find(prob_det_if_tr)
  n = length(idx_detectable)
  @assert n <= max_tranets_in_sys       # Make sure memory to store these
  #if n==0
  #   println("# WARNING found no detectable planets based on ",prob_det_if_tr)
  #end
  invalid_prob_flag = false

  #ps_detectable = PlanetarySystemSingleStar(ps,idx_detectable)
  ps_detectable = PlanetarySystem(ps,idx_detectable)

  combo_sample_probs = rand(num_samples)
  combo_cum_probs = zeros(num_samples)
  #sdp = SimulatedSystemDetectionProbs{observer_trait}(prob_det_if_tr[idx_detectable])
  sdp = SimulatedSystemDetectionProbs(observer_trait,prob_det_if_tr[idx_detectable])

  planet_should_transit = zeros(Cint,n)
  for p in 1:num_planets(sdp)
      sdp.pairwise[p,p] = 0.0
  end
  for ntr in 1:min(n,max_tranets_in_sys)  # Loop over number of planets transiting
      sdp.n_planets[ntr] = 0.0
      for combo in combinations(1:n,ntr)  # Loop over specific combinations of detectable planets
        prob_det_this_combo = 1.0
        for p in combo  # Loop over each planet in this combination of detectable planets
            prob_det_this_combo *= prob_det_if_tr[idx_detectable[p]]
        end

        if prob_det_this_combo < min_detect_prob_to_be_included
           continue
        end

	fill!(planet_should_transit,zero(Cint))
	for i in 1:length(combo)
      	    planet_should_transit[combo[i]] = one(Cint)
	end
	local geo_factor::Float64
	if observer_trait == SkyAveraged
	   geo_factor = prob_combo_transits_obs_ave(ps_detectable,planet_should_transit)
	elseif observer_trait == OneObserver
	   geo_factor = prob_combo_transits_one_obs(ps_detectable,planet_should_transit)
	else
	   error(string("typeof(",observer_trait,") is not a valid trait."))
	end
        prob_det_this_combo *= geo_factor

	# Store samples of combinations of planets detected drawn from the full joint multi-planet density
	for i in 1:num_samples
	   if combo_cum_probs[i] < combo_sample_probs[i] <= combo_cum_probs[i]+prob_det_this_combo
	      sdp.combo_detected[i] = combo
	   end
	end
	combo_cum_probs += prob_det_this_combo

        sdp.n_planets[ntr] += prob_det_this_combo   # Accumulate the probability of detecting any n planets

        for p in combo                # Accumulate the probability of detecting each planet individually
            sdp.pairwise[p,p] += prob_det_this_combo
        end

        #=
        for pq in combinations(combo,2)                # Accumulate the probability of detecting each planet pair # TODO: OPT: replace with simply calculating integers for pairs to avoid allocations of small arrays
           sdp.pairwise[pq[1],pq[2]] = prob_det_this_combo
           sdp.pairwise[pq[2],pq[1]] = prob_det_this_combo   # TODO OPT: Remove if use symmetric matrix type.
        end
        =#
        if length(combo)>=2
          for pi in 2:length(combo)                # Accumulate the probability of detecting each planet pair # TODO: OPT: replace with simply calculating integers for pairs to avoid allocations of small arrays
            p = combo[pi]
            for qi in 1:(pi-1)
               q = combo[qi]
               sdp.pairwise[p,q] = prob_det_this_combo
               sdp.pairwise[q,p] = prob_det_this_combo   # TODO OPT: Remove if use symmetric matrix type.
            end # qi
          end # pi
        end # if
      end # combo
  end # ntr

  for p in 1:n
      if sdp.pairwise[p,p] > 1.0
          invalid_prob_flag = true
          println(string("Error! Invalid prob for planet ",p,": ", sdp.pairwise[p,p]))
      end
  end

  if invalid_prob_flag
      println("")
      for ntr in 1:min(n,max_tranets_in_sys)
          for combo in combinations(1:n,ntr)
              fill!(planet_should_transit,zero(Cint))
	      for i in 1:length(combo)
      	          planet_should_transit[combo[i]] = one(Cint)
	      end
              if length(combo) == n
                  geo_factor = prob_combo_transits_obs_ave(ps_detectable,planet_should_transit, print_orbit = true)
              else
                  geo_factor = prob_combo_transits_obs_ave(ps_detectable,planet_should_transit)
              end
              println(string("Geo. factor of ",combo," = ",geo_factor))
          end
      end
      println(string("Det. prob. = ", prob_det_if_tr[idx_detectable]))
      println("")
      #quit()
  end
  return sdp
end

if false # WARNING: Complicated and untested
function combine_system_detection_probs{T}(prob::Vector{SimulatedSystemDetectionProbs{T}}, s1::Integer, s2::Integer) # WARNING: Complicated and untested
    npl_s1 = min(num_planets(prob[s1]), max_tranets_in_sys)
    npl_s2 = min(num_planets(prob[s2]), max_tranets_in_sys)
    num_planets_across_systems = npl_s1 + npl_s2
    prob_merged = SimulatedSystemDetectionProbs{T}(num_planets_across_systems)
    # Copy probabilities for detecting planets and planet pairs within one system
    prob_merged.detect_planet_if_transits[1:npl_s1] = prob[s1].detect_planet_if_transits
    prob_merged.pairwise[1:npl_s1,1:npl_s1] = prob[s1].pairwise
    offset = npl_s1
    prob_merged.detect_planet_if_transits[offset+1:offset+npl_s2] = prob[s2].detect_planet_if_transits
    prob_merged.pairwise[offset+1:offset+npl_s2,offset+1:offset+npl_s2] = prob[s2].pairwise
    # Calculate probabilities of detecting pairs of planets in different systems, assuming uncorrelated orientations
    for p1 in 1:npl_s1
        for p2 in 1:npl_s2
	    prob_detect_both = prob_detect(prob[s1],p1)*prob_detect(prob[s2],p2)
	    idx1 = p1
	    idx2 = offset+p2
	    prob_merged.pairwise[idx1,idx2] = prob_detect_both
	    prob_merged.pairwise[idx2,idx1] = prob_detect_both
	end
    end
    # Merge probabilities of detecting n_planets, assuming uncorrelated orientations
    p_zero_pl_s1 = 1.0-sum(prob[s1].n_planets)
    p_zero_pl_s2 = 1.0-sum(prob[s2].n_planets)
    for n in 1:min(num_planets_across_systems,max_tranets_in_sys)
        prob_merged.n_planets[n] = p_zero_pl_s1*prob[s2].n_planets[n] + p_zero_pl_s2*prob[s1].n_planets[n]
	for i in 1:n-1
	   prob_merged.n_planets[n] += prob[s1].n_planets[i]*prob[s2].n_planets[n-i] + prob[s1].n_planets[n-i]*prob[s2].n_planets[i]
	end
    end
    # Combine samples of detected planet combinations
    prob_merged.combo_detected = fill(Array{Int64}(0), min(length(prob[s1].combo_detected), length(prob[s2].combo_detected) ) )
    for i in 1:length(prob_merged.combo_detected)
       prob_merged.combo_detected[i] = vcat( prob[s1].combo_detected, prob[s1].combo_detected+offset )
    end
    return prob_merged
end

function select_subset{T<:SystemDetectionProbsTrait}(prob::SimulatedSystemDetectionProbs{T}, idx::Vector{Int64}) # WARNING: Complicated and untested
    n = length(idx)
    subset = SimulatedSystemDetectionProbs{T}(n)
    subset.detect_planet_if_transits = prob.detect_planet_if_transits[idx]
    subset.pairwise = ones(n,n)
    subset.n_planets = zeros(max(n,length(prob.n_planets)))
    if 1<=n<=length(subset.n_planets)
        subset.n_planets[n] = 1.0
    end
    subset.combo_detected = Array{Int64,1}[collect(1:n)]
    return subset
end

# ASSUMING: Planetary systems for same target are uncorrelated
# Compute sky-averaged transit probabilities from a target with known physical properties
function calc_simulated_system_detection_probs(t::KeplerTarget, sim_param::SimParam ) # WARNING: Complicated and untested
  const max_tranets_in_sys = get_int(param,"max_tranets_in_sys",10)
  const min_detect_prob_to_be_included = get(param,"max_tranets_in_sys", 0.0)
  s1 = findfirst(num_planets,t.sys)
  if num_planets(t) == num_planets(t.sys[s1])
    # Target has only one system with planets
    prob_det_if_tr = Float64[calc_ave_prob_detect_if_transit(t, s1, p, sim_param) for p in 1:num_planets(t.sys[s1])]
    return calc_simulated_system_detection_probs(t.sys[s1], prob_det_if_tr, max_tranets_in_sys=max_tranets_in_sys, min_detect_prob_to_be_included=min_detect_prob_to_be_included )
  else
    # Target has multiple systems with planets
    sdp = SkyAveragedSystemDetectionProbs[ SkyAveragedSystemDetectionProbs( min(num_planets(t.sys[s]),max_tranets_in_sys) ) for s in 1:length(t.sys) ]
    num_planets_across_systems = 0
    # Calculate detection probabilities for each system separately
    for s in 1:length(t.sys)
        prob_det_if_tr = Float64[calc_ave_prob_detect_if_transit(t, s, p, sim_param) for p in 1:num_planets(t.sys[s])]
    	sdp[s] = calc_simulated_system_detection_probs(t.sys[s], prob_det_if_tr, max_tranets_in_sys=max_tranets_in_sys, min_detect_prob_to_be_included=min_detect_prob_to_be_included )
    	num_planets_across_systems += num_planets(t.sys[s])
    end
    @assert num_planets_across_systems <= max_tranets_in_sys # Make sure memory to store these  # QUERY: DETAIL: Should we relax?
    # Find system ids for first two systems with planets
    #s1, s2 = find_system_detection_probs_with_planets(::Vector{})
    more_than_two_systems_with_planets = false
    s1 = 0
    s2 = 0
    for s in 1:length(t.sys)
	if num_planets(t.sys[s])>=1
	   if s1==0
	      s1 = s
	   elseif s2==0
	      s2 = s
	   else
	      more_than_two_systems_with_planets = true
	   end
	end
    end
    @assert (s1!=0) && (s2!=0)
    @assert !more_than_two_systems_with_planets
    sdp_merged = combine_system_detection_probs(sdp,s1,s2)   # Merge SystemDetectionProbs across systems with common target
    return sdp_merged
  end
end
end

type ObservedSystemDetectionProbs <: SystemDetectionProbsAbstract          # TODO OPT:  For observed systems (or simulations of observed systems) were we can't know everything.  Is this even used?  Or should we just compute these on the fly, rather than storing them? Do we even want to keep this?
  planet_transits::Vector{Float64}                         # Probability that each planet transits individually for one observer based on actual i, e, and omega
  detect_planet_if_transits::Vector{Float64}               # Probability of detecting each planet given that it transits. Assumes one observer based on actual i, e and omega
  # snr::Vector{Float64}                      # Dimensionless SNR of detection for each planet QUERY: Should we store this here?
end
ObservedSystemDetectionProbs(p::Vector{Float64}) = ObservedSystemDetectionProbs( ones(length(p)), p )
ObservedSystemDetectionProbs(n::Integer) = ObservedSystemDetectionProbs( ones(n), zeros(n) )
ObservedSystemDetectionProbsEmpty() = ObservedSystemDetectionProbs(0)

# Functions common to various types of SystemDetectionProbs
num_planets(prob::ObservedSystemDetectionProbs) = length(prob.detect_planet_if_transits)
prob_detect_if_transits(prob::ObservedSystemDetectionProbs, pl_id::Integer) = prob.detect_planet_if_transits[pl_id]
prob_detect(prob::ObservedSystemDetectionProbs, pl_id::Integer) = prob.planet_transits[pl_id]*prob.detect_planet_if_transits[pl_id]
#prob_detect_both_planets(prob::ObservedSystemDetectionProbs, pl_id::Integer, ql_id::Integer) = prob_detect(prob,pl_id) * prob_detect(prob,ql_id)   # WARNING: Assumes independent.  Intent is for testing CORBITS.  Or should we delete?

#if false
function prob_detect_n_planets(prob::ObservedSystemDetectionProbs, n::Integer)   # WARNING: Assumes independent.  Intent is for testing CORBITS.  May wnat to comment out to prevent accidental use.
  if n<1 || n > num_planets(prob)  return 0.0 end
  sum_prob = 0.0
  for combo in combinations(1:num_planets(prob), n)
      prob_this_combo = 1.0
      for pl_id in combo
          prob_this_combo *= prob_detect(prob,pl_id)
      end
      sum_prob += prob_this_combo
  end
  return sum_prob
end

#=
# Compute transit probabilities for a single observer from a target with known physical properties
function calc_observed_system_detection_probs(targ::KeplerTarget, sim_param::SimParam)
  n = num_planets(targ)
  pdet = zeros(n)
  ptr = zeros(n)
  pl = 1
  for s in 1:length(targ.sys)
      for p in 1:length(targ.sys[s].planet)
         pdet[pl] = calc_prob_detect_if_transit_with_actual_b(targ, s, p, sim_param)
         ptr[pl]  = calc_transit_prob_single_planet_one_obs(targ, s, p)
         pl += 1
      end
  end
  ObservedSystemDetectionProbs( ptr, pdet )
end
=#

# Estimate transit probabilities for a single observer from a target with known physical properties.
if false    # Do we actually want this for anything?
function calc_observed_system_detection_probs(kto::KeplerTargetObs, sim_param::SimParam)
  n = num_planets(kto)
  pdet = ones(n)   # WARNING: We assume all observed objects were detected and we don't have enough info to calcualte a detection probability.  Do we want to do something different?
  ptr = zeros(n)
  for pl in 1:n
      ptr[pl]  = calc_transit_prob_single_planet_approx(kto.obs[pl].period, kto.star.radius, kto.star.mass )
  end
  ObservedSystemDetectionProbs( ptr, pdet )
end
end


text:
fname = joinpath(dirname(@__FILE__), "../data", "test_Hz19.5-testing.bdf")

a = read_SSR(fname)

a.processing["Carrier_Frequency"] = 500

a = subsample(a)
a = subsample(a, plot=true,  plot_channel=4)
a = subsample(a, valid_triggers = 1)

a = highpass_filter(a)

a = rereference(a, "Cz")

a = extract_epochs(a)

a = hotelling(a)

println(a.processing["statistics"])

# These values are not validated, but just to ensure it doesn't change
# TODO get true values to check with
# @test_approx_eq_eps a.processing["statistics"][:SNRdB] [NaN, -0.0474418, -1.42819, -0.334432, -2.56155, -1.39146] 0.002

a = read_SSR(fname)

a.processing["Carrier_Frequency"] = 500

a = blank(a, 0.0018, valid_triggers=1)


text:
using Codex
using Statistics
using Test

@testset "Stats" begin
    A = [1, 2, 3]
    B = [3, 4, 5]
    # Checked with effsize package in R.
    @test cohens_d(A, B) == -2
    @test cohens_d(length(A), mean(A), std(A), length(B), mean(B), std(B)) == -2
end


text:
cd("/Users/Frank/.julia/dev/Bridge")

using Bridge, StaticArrays, Distributions
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models
using DelimitedFiles
using RCall
using DataFrames
using CSV

@rlibrary ggplot2

T = 2.0
dt = 1/5000
tt = 0.:dt:T

aux_choice = ["linearised_end" "linearised_startend" "matching"][1]
endpoint = ["first", "extreme"][2]

# specify target process
struct FitzhughDiffusion <: ContinuousTimeProcess{ℝ{2}}
    ϵ::Float64
    s::Float64
    γ::Float64
    β::Float64
    σ::Float64
end

#Bridge.B(t, P::Union{FitzhughDiffusion,FitzhughDiffusionAux}) = @SMatrix [-1.0 1.0; 0.0 -1.0]

Bridge.b(t, x, P::FitzhughDiffusion) = ℝ{2}((x[1]-x[2]-x[1]^3+P.s)/P.ϵ, P.γ*x[1]-x[2] +P.β)
Bridge.σ(t, x, P::FitzhughDiffusion) = ℝ{2}(0.0, P.σ)
Bridge.constdiff(::FitzhughDiffusion) = true

# Generate Data
Random.seed!(3)
P = FitzhughDiffusion(0.1, 0.0, 1.5, 0.8, 0.3)
W = sample(tt, Wiener())
#x0 = ℝ{2}(-0.5, 0.0) # is a good choice, with either v=-0.5 or where the forward process ends
x0 = ℝ{2}(-0.5, -0.6)
X = solve(Euler(), x0, W, P)

XX = [X]
samples = 250
#draw 9 more path and save these
for j in 2:samples
    W = sample(tt, Wiener())
    X = solve(Euler(), x0, W, P)
    push!(XX, X)
end

# write forwards to csv file
f = open("/Users/Frank/Dropbox/DiffBridges/Rcode/integrated_diff/forwards.csv","w")
head = "iteration, time, component, value \n"
write(f, head)
iterates = [Any[i, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:length(X), i in 1:samples ][:]
writedlm(f, iterates, ',')
close(f)

# Choose truepath
if endpoint == "first"
    itrue = 1
elseif endpoint == "extreme"
    itrue = findmax([XX[i].yy[end][1] for i in 1:samples])[2]
else
    error("not implemented")
end

lt = size(X.tt)[1]
trueX = DataFrame(iteration=fill(-50000,2*lt),time=[X.tt; X.tt],
    component=[fill(1,lt);fill(2,lt)], value=[first.(XX[itrue].yy); last.(XX[itrue].yy)] )
ggplot(trueX,aes(x=:time,y=:value)) +
      geom_path()+facet_wrap(R"~component",ncol=1,scales="free_y")

# Specify observation scheme
L = @SMatrix [1. 0.]
Σ = @SMatrix [0.0]
#v = ℝ{1}(-0.5)

v = ℝ{1}(XX[itrue].yy[end][1] + √Σ[1,1]*randn())
#v = ℝ{2}(X.yy[end])

# L = @SMatrix [1. 0.; 0. 1.]
# Σ = @SMatrix [0.0 0.0 ; 0.0 0.0]
# v = ℝ{2}(-5.0, 2)    #ℝ{1}(X.yy[end][1] + √Σ[1,1]*randn())

# specify auxiliary process
struct FitzhughDiffusionAux <: ContinuousTimeProcess{ℝ{2}}
    ϵ::Float64
    s::Float64
    γ::Float64
    β::Float64
    σ::Float64
    t::Float64
    u::Float64
    T::Float64
    v::Float64
end

function uv(t, P::FitzhughDiffusionAux)
    λ = (t - P.t)/(P.T - P.t)
    P.v*λ + P.u*(1-λ)
end

if aux_choice=="linearised_end"
    Bridge.B(t, P::FitzhughDiffusionAux) = @SMatrix [1/P.ϵ-3*P.v^2/P.ϵ  -1/P.ϵ; P.γ -1.0]
    Bridge.β(t, P::FitzhughDiffusionAux) = ℝ{2}(P.s/P.ϵ+2*P.v^3/P.ϵ, P.β)
    ρ = endpoint=="extreme" ? 0.99 : 0.6
elseif aux_choice=="linearised_startend"
    Bridge.B(t, P::FitzhughDiffusionAux) = @SMatrix [1/P.ϵ-3*uv(t, P)^2/P.ϵ  -1/P.ϵ; P.γ -1.0]
    Bridge.β(t, P::FitzhughDiffusionAux) = ℝ{2}(P.s/P.ϵ+2*uv(t, P)^3/P.ϵ, P.β)
    ρ = endpoint=="extreme" ? 0.99 : 0.6
else
    Bridge.B(t, P::FitzhughDiffusionAux) = @SMatrix [1/P.ϵ  -1/P.ϵ; P.γ -1.0]
    Bridge.β(t, P::FitzhughDiffusionAux) = ℝ{2}(P.s/P.ϵ-(P.v^3)/P.ϵ, P.β)
    ρ = 0.99
end

Bridge.σ(t, x, P::FitzhughDiffusionAux) = ℝ{2}(0.0, P.σ)
Bridge.constdiff(::FitzhughDiffusionAux) = true

Bridge.b(t, x, P::FitzhughDiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P)
Bridge.a(t, P::FitzhughDiffusionAux) = Bridge.σ(t,0,P) * Bridge.σ(t, 0, P)'

#Pt = IntegratedDiffusionAux(0.7)
Pt = FitzhughDiffusionAux(0.1, 0.0, 1.5, 0.8, 0.3, tt[1], x0[1], tt[end], v[1])

# Solve Backward Recursion
Po = Bridge.PartialBridge(tt, P, Pt, L, v, Σ)

# initalisation
W = sample(tt, Wiener())
Xo = copy(X)
bridge!(Xo, x0, W, Po)
sample!(W, Wiener())
bridge!(X, x0, W, Po)
ll = llikelihood(Bridge.LeftRule(), X, Po)

# settings sampler
iterations = 10*10^4
subsamples = 0:1000:iterations

# further initialisation
Wo = copy(W)
W2 = copy(W)
XX = Any[]
if 0 in subsamples
    push!(XX, copy(X))
end

acc = 0

for iter in 1:iterations
    # Proposal
    sample!(W2, Wiener())
    Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy
    bridge!(Xo, x0, Wo, Po)

    llo = llikelihood(Bridge.LeftRule(), Xo, Po)
    print("ll $ll $llo, diff_ll: ",round(llo-ll,3))#, X[10], " ", Xo[10])

    if log(rand()) <= llo - ll
        X.yy .= Xo.yy
        W.yy .= Wo.yy
        ll = llo
        print("✓")
        acc +=1
    end
    println()
    if iter in subsamples
        push!(XX, copy(X))
    end
end

@info "Done."*"\x7"^6

# write mcmc iterates to csv file
f = open("/Users/Frank/Dropbox/DiffBridges/Rcode/integrated_diff/iterates.csv","w")
head = "iteration, time, component, value \n"
write(f, head)
iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:length(X), (i,s) in enumerate(subsamples) ][:]
writecsv(f,iterates)
close(f)

# write a true path to csv file
CSV.write("/Users/Frank/Dropbox/DiffBridges/Rcode/integrated_diff/truepath.csv",trueX)

println("Average acceptance percentage: ",100*round(acc/iterations,2),"\n")

# Plotting with ggplot
dd = DataFrame(iteration=map(x->x[1],iterates), time=map(x->x[2],iterates),
        component=map(x->x[3],iterates),value=map(x->x[4],iterates))

ind =  map(x-> x in [0:200:1000;  40000:2000:50000],dd[:iteration])
ddsub = dd[ind,:]

ggplot([trueX;ddsub],aes(x=:time,y=:value,group=:iteration,colour=:iteration)) +
        geom_path()+facet_wrap(R"~component",ncol=1,scales="free_y")+
         scale_colour_gradient(low="green",high="blue")            |>
 p -> ggsave("/Users/Frank/Dropbox/DiffBridges/Rcode/integrated_diff/p1.pdf", p, width=7,height=4)

ind2 =  map(x-> x in 50000:2000:100000,dd[:iteration])
ddsub2 = dd[ind2,:]

ggplot([trueX;ddsub2],aes(x=:time,y=:value,group=:iteration,colour=:iteration)) +
        geom_path()+facet_wrap(R"~component",ncol=1,scales="free_y")+
         scale_colour_gradient(low="green",high="blue")            |>
 p -> ggsave("/Users/Frank/Dropbox/DiffBridges/Rcode/integrated_diff/p2.pdf", p, width=7,height=4)


text:
"""
    ModelWithParams(args...; kwargs...)::JuMP.Model

This function is deprecated. Use `JuMP.Model` instead.
"""
function ModelWithParams(args...; kwargs...)
    @warn("ModelWithParams is deprecated and no longer needed. Use JuMP.Model instead.")
    return JuMP.Model(args...; kwargs...)
end

export ModelWithParams

function JuMP.fix(p::ParameterRef, val::Real)
    @warn(
        "JuMP.fix has been deprecated. Use `set_value(p, v)` instead.",
        maxlog = 1,
    )
    return set_value(p, val)
end

"""
    add_parameter(model::JuMP.Model, val::Real)::ParameterRef

Adds one parameter fixed at `val` to the model `model`.

add_parameter(model::JuMP.Model)::ParameterRef

Adds one parameter fixed at zero to the model `model`.
"""
function add_parameter(model::JuMP.Model, val::Real = 0.0)
    @warn(
        "This function is deprecated. Pass `Param` to `@variable` instead.",
        maxlog = 1,
    )
    return _add_parameter(model, val)
end

"""
    add_parameters(model::JuMP.Model, val::Vector{R})::Vector{ParameterRef}

Adds one parameter for each element of the vector `val`.

    add_parameters(model::JuMP.Model, N::Integer)::Vector{ParameterRef}

Adds `N` parameters to the model, all of them fixed in zero.
"""
function add_parameters(model::JuMP.Model, N::Integer)
    return add_parameters(model::JuMP.Model, zeros(N))
end

function add_parameters(model::JuMP.Model, val::AbstractArray)
    return add_parameter.(model, val)
end

export add_parameter, add_parameters


text:
using CosmoCentral
using Test
using QuadGK
using NumericalIntegration
using PyCall
using Interpolations

numpy = pyimport("numpy")

w0waCDMCosmology = CosmoCentral.Flatw0waCDMCosmology()
AnalitycalDensity = CosmoCentral.AnalitycalDensity()
InstrumentResponse = CosmoCentral.InstrumentResponse()
ConvolvedDensity = CosmoCentral.ConvolvedDensity()
CosmologicalGrid  = CosmoCentral.CosmologicalGrid(
ZArray=Array(LinRange(0.001, 4.0, 500)))
BackgroundQuantities = CosmoCentral.BackgroundQuantities(HZArray=
zeros(length(CosmologicalGrid.ZArray)),
rZArray=zeros(length(CosmologicalGrid.ZArray)))
GCWeightFunction = CosmoCentral.GCWeightFunction(WeightFunctionArray =
zeros(length(ConvolvedDensity.DensityNormalizationArray),
length(CosmologicalGrid.ZArray)))
WLWeightFunction = CosmoCentral.WLWeightFunction(WeightFunctionArray =
zeros(length(ConvolvedDensity.DensityNormalizationArray),
length(CosmologicalGrid.ZArray)),
LensingEfficiencyArray = zeros(length(
ConvolvedDensity.DensityNormalizationArray),
length(CosmologicalGrid.ZArray)))
AngularCoefficients = CosmoCentral.Cℓ(
CℓArray = zeros(length(CosmologicalGrid.ℓBinCenters),
length(GCWeightFunction.WeightFunctionArray[:, 1]),
length(GCWeightFunction.WeightFunctionArray[:, 1])))
classyParams = CosmoCentral.Initializeclassy(w0waCDMCosmology)
input_path_pmm = pwd()*"/p_mm"
input_path_Cℓ = pwd()*"/cl"
input_path_Forecast = pwd()
CosmoCentral.classy.Class()

@testset "Evaluation of background quantities" begin
    test_E_z = CosmoCentral.ComputeAdimensionalHubbleFactor(0., w0waCDMCosmology)
    @test test_E_z == 1.
    test_H_z = CosmoCentral.ComputeHubbleFactor(0., w0waCDMCosmology)
    @test test_H_z == w0waCDMCosmology.H0
    test_r_z = CosmoCentral.Computeχ(0., w0waCDMCosmology)
    @test test_r_z == 0.
    test_r_array = zeros(length(CosmologicalGrid.ZArray))
    test_H_array = zeros(length(CosmologicalGrid.ZArray))
    for (idxz, zvalue) in enumerate(CosmologicalGrid.ZArray)
        test_r_array[idxz] = CosmoCentral.Computeχ(zvalue,
        w0waCDMCosmology)
        test_H_array[idxz] = CosmoCentral.ComputeHubbleFactor(zvalue,
        w0waCDMCosmology)
    end
    CosmoCentral.ComputeBackgroundQuantitiesGrid!(CosmologicalGrid,
    BackgroundQuantities, w0waCDMCosmology)
    @test test_r_array ==BackgroundQuantities.rZArray
    @test test_H_array ==BackgroundQuantities.HZArray
end

@testset "Check the normalization of density function" begin
    CosmoCentral.NormalizeAnalitycalDensity!(AnalitycalDensity)
    int, err = QuadGK.quadgk(x -> CosmoCentral.ComputeDensity(x, AnalitycalDensity),
    AnalitycalDensity.ZMin, AnalitycalDensity.ZMax, rtol=1e-12)
    @test isapprox(int, AnalitycalDensity.SurfaceDensity, atol=1e-9)
end

@testset "Check the normalization of convolved density function" begin
    test_normalization = zeros(length(ConvolvedDensity.ZBinArray)-1)
    CosmoCentral.NormalizeConvolvedDensity!(ConvolvedDensity, AnalitycalDensity,
    InstrumentResponse, CosmologicalGrid)
    for idx in 1:length(test_normalization)
        int, err = QuadGK.quadgk(x ->
        CosmoCentral.ComputeConvolvedDensity(x, idx,
        ConvolvedDensity, AnalitycalDensity, InstrumentResponse),
        AnalitycalDensity.ZMin,
        AnalitycalDensity.ZMax, rtol=1e-12)
        test_normalization[idx] = int
    end
    @test isapprox(test_normalization, ones(length(test_normalization)), atol=1e-12)
end

@testset "Check the computation of the convolved density function on grid" begin
    test_array = zeros(Float64,length(ConvolvedDensity.ZBinArray)-1,
    length(CosmologicalGrid.ZArray))
    CosmoCentral.ComputeConvolvedDensityGrid!(CosmologicalGrid, ConvolvedDensity,
    AnalitycalDensity, InstrumentResponse)
    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
            test_array[idx_ZBinArray, idx_ZArray] =
            CosmoCentral.ComputeConvolvedDensity(
            CosmologicalGrid.ZArray[idx_ZArray],
            idx_ZBinArray,
            ConvolvedDensity,
            AnalitycalDensity,
            InstrumentResponse)
        end
    end
    @test isapprox(test_array, ConvolvedDensity.DensityGridArray, atol=1e-12)
end

@testset "Check the functions defined in MathUtils" begin
    minarray = 1e-5
    maxarray = 10.
    n = 100
    benchmark_numpy = numpy.logspace(log10(minarray), log10(maxarray), n)
    test_logspace = CosmoCentral.LogSpaced(minarray, maxarray, n)
    @test isapprox(benchmark_numpy, test_logspace, atol=1e-12)
    array = [1.,2.,3.]
    @test 1 == CosmoCentral.BinSearch(1.5, array)
    x = Array(LinRange(0., 10., 100))
    y = 11.2 .*x .+0.4
    c, m = CosmoCentral.CustomRegression(x, y)
    @test isapprox(c, 0.4, atol=1e-12)
    @test isapprox(m, 11.2, atol=1e-12)
end

@testset "Check the Piecewise bias evaluation" begin
    z = 0.4
    test_array = zeros(length(CosmologicalGrid.ZArray))
    bias = CosmoCentral.ComputeBias(z, GCWeightFunction.BiasKind,
    ConvolvedDensity)
    @test isapprox(1.0997727037892875, bias, atol=1e-12)
    for (idxz, zvalue) in enumerate(CosmologicalGrid.ZArray)
        test_array[idxz] = CosmoCentral.ComputeBias(zvalue,
        GCWeightFunction.BiasKind, ConvolvedDensity)
    end
    CosmoCentral.ComputeBiasGrid!(CosmologicalGrid,
        GCWeightFunction, ConvolvedDensity)
    @test isapprox(test_array, GCWeightFunction.BiasArray[1,:], atol=1e-10)
end

@testset "Check the Weight function evaluation" begin
    test_gc = zeros(length(ConvolvedDensity.ZBinArray)-1,
    length(CosmologicalGrid.ZArray))
    test_le = zeros(length(ConvolvedDensity.ZBinArray)-1,
    length(CosmologicalGrid.ZArray))
    test_wl = zeros(length(ConvolvedDensity.ZBinArray)-1,
    length(CosmologicalGrid.ZArray))
    for (idxz, zvalue) in enumerate(CosmologicalGrid.ZArray)
        for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
            test_gc[idx_ZBinArray, idxz] =
            CosmoCentral.ComputeWeightFunction(zvalue, idx_ZBinArray,
            ConvolvedDensity, AnalitycalDensity, InstrumentResponse,
            w0waCDMCosmology, GCWeightFunction)
            test_le[idx_ZBinArray, idxz] =
            CosmoCentral.ComputeLensingEfficiency(zvalue, idx_ZBinArray,
            ConvolvedDensity, AnalitycalDensity, InstrumentResponse,
            w0waCDMCosmology, CosmologicalGrid, WLWeightFunction)
            test_wl[idx_ZBinArray, idxz] =
            CosmoCentral.ComputeWeightFunction(zvalue, idx_ZBinArray,
                ConvolvedDensity, AnalitycalDensity,
                InstrumentResponse, w0waCDMCosmology, CosmologicalGrid,
                WLWeightFunction)
        end
    end
    CosmoCentral.ComputeWeightFunctionGrid!(GCWeightFunction,
        ConvolvedDensity,
        CosmologicalGrid,
        BackgroundQuantities,
        w0waCDMCosmology)
    CosmoCentral.ComputeLensingEfficiencyGrid!(
    WLWeightFunction, AnalitycalDensity,
    InstrumentResponse, ConvolvedDensity,
    CosmologicalGrid,
    BackgroundQuantities,
    w0waCDMCosmology,
    CosmoCentral.StandardLensingEfficiency())
    CosmoCentral.ComputeWeightFunctionGrid!(WLWeightFunction,
        ConvolvedDensity,
        CosmologicalGrid,
        BackgroundQuantities,
        w0waCDMCosmology)
    @test isapprox(test_gc, GCWeightFunction.WeightFunctionArray, rtol=1e-9)
    @test isapprox(test_wl, WLWeightFunction.WeightFunctionArray, rtol=1e-9)
    @test isapprox(test_le, WLWeightFunction.LensingEfficiencyArray, rtol=1e-9)
end

@testset "Check the Power Spectrum evaluated over the Limber Grid" begin
    ℓBinCenters = Array(LinRange(10.5, 2999.5, 2990))
    PowerSpectrum, BackgroundQuantitiesLoaded, CosmologicalGrid =
    CosmoCentral.ReadPowerSpectrumBackground(input_path_pmm, ℓBinCenters, Array(ones(2990)))
    CosmoCentral.InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid,
    BackgroundQuantitiesLoaded, PowerSpectrum, CosmoCentral.BSplineCubic())
    CosmoCentral.ComputeLimberArray!(CosmologicalGrid,
    BackgroundQuantitiesLoaded)
    test_k_limber = (CosmologicalGrid.ℓBinCenters[1]+0.5) /
    BackgroundQuantitiesLoaded.rZArray[1]
    @test test_k_limber == CosmologicalGrid.KLimberArray[1, 1]
    test_Omega_cdm = w0waCDMCosmology.ΩM-w0waCDMCosmology.ΩB-
    w0waCDMCosmology.Mν/(93.14*(w0waCDMCosmology.H0/100)^2)
    @test test_Omega_cdm == classyParams.classyParamsDict["Omega_cdm"]
    x = LinRange(log10(first(CosmologicalGrid.KArray)),
    log10(last(CosmologicalGrid.KArray)), length(CosmologicalGrid.KArray))
    y = LinRange(first(CosmologicalGrid.ZArray), last(CosmologicalGrid.ZArray),
    length(CosmologicalGrid.ZArray))
    InterpPmm = Interpolations.interpolate(
    log10.(PowerSpectrum.PowerSpectrumNonlinArray),
    BSpline(Cubic(Line(OnGrid()))))
    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
    test_power_spectrum = 10 .^InterpPmm(log10(
    CosmologicalGrid.KLimberArray[1, 1]), CosmologicalGrid.ZArray[1])
    CosmoCentral.InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid,
    BackgroundQuantitiesLoaded, PowerSpectrum, CosmoCentral.BSplineCubic())
    @test isapprox(test_power_spectrum,
    PowerSpectrum.InterpolatedPowerSpectrum[1, 1], rtol=1e-2)
    CosmoCentral.WritePowerSpectrumBackground(PowerSpectrum,
    BackgroundQuantitiesLoaded, CosmologicalGrid, "new_p_mm")
    NewPowerSpectrum, NewBackgroundQuantitiesLoaded, NewCosmologicalGrid =
    CosmoCentral.ReadPowerSpectrumBackground("new_p_mm", ℓBinCenters, Array(ones(2990)))
    @test isapprox(PowerSpectrum.PowerSpectrumNonlinArray,
    NewPowerSpectrum.PowerSpectrumNonlinArray, rtol=1e-6)
    @test isapprox(PowerSpectrum.PowerSpectrumLinArray,
    NewPowerSpectrum.PowerSpectrumLinArray, rtol=1e-6)
    @test isapprox(BackgroundQuantitiesLoaded.HZArray,
    NewBackgroundQuantitiesLoaded.HZArray, rtol=1e-6)
    @test isapprox(BackgroundQuantitiesLoaded.rZArray,
    NewBackgroundQuantitiesLoaded.rZArray, rtol=1e-6)
    PowerSpectrumDierckx, BackgroundQuantitiesLoaded, CosmologicalGrid =
    CosmoCentral.ReadPowerSpectrumBackground(input_path_pmm, ℓBinCenters, Array(ones(2990)))
    CosmoCentral.InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid,
    BackgroundQuantitiesLoaded, PowerSpectrumDierckx,
    CosmoCentral.RectBivSplineDierckx())
    @test isapprox(PowerSpectrumDierckx.InterpolatedPowerSpectrum[1, 1],
    PowerSpectrum.InterpolatedPowerSpectrum[1, 1], rtol=1e-2)
end

@testset "Test Angular coefficients evaluation" begin
    CℓLoaded = CosmoCentral.ReadCℓ(input_path_Cℓ, "PhotometricGalaxy_PhotometricGalaxy")
    ℓBinCenters = Array(LinRange(10.5, 2999.5, 2990))
    PowerSpectrum, BackgroundQuantities, CosmologicalGrid =
    CosmoCentral.ReadPowerSpectrumBackground(input_path_pmm, ℓBinCenters, Array(ones(2990)))
    CosmoCentral.ComputeLimberArray!(CosmologicalGrid, BackgroundQuantities)
    CosmoCentral.InterpolatePowerSpectrumLimberGrid!(CosmologicalGrid,
    BackgroundQuantities, PowerSpectrum, CosmoCentral.BSplineCubic())
    Cℓ = CosmoCentral.Cℓ(CℓArray = zeros(length(CosmologicalGrid.ℓBinCenters),
    length(GCWeightFunction.WeightFunctionArray[:, 1]),
    length(GCWeightFunction.WeightFunctionArray[:, 1])))
    CosmoCentral.ComputeCℓ!(Cℓ, GCWeightFunction, GCWeightFunction, BackgroundQuantities,
    w0waCDMCosmology, CosmologicalGrid, PowerSpectrum, CosmoCentral.CustomSimpson())
    @test isapprox(CℓLoaded.CℓArray,
    Cℓ.CℓArray, rtol=1e-6)
    CosmoCentral.WriteCℓ!("PhotometricGalaxy_PhotometricGalaxy", Cℓ, "new_cl")
    CℓReloaded = CosmoCentral.ReadCℓ("new_cl", "PhotometricGalaxy_PhotometricGalaxy")
    @test isapprox(CℓReloaded.CℓArray, Cℓ.CℓArray, rtol=1e-9)
end

@testset "Test Fisher Forecast: Integration Test" begin
    run(`wget https://zenodo.org/record/5270335/files/forecast_pmm.tar.xz\?download=1`);
    run(`mv forecast_pmm.tar.xz\?download\=1 forecast_pmm.tar.xz`);
    run(`tar xvf forecast_pmm.tar.xz`);
    DictCosmo = Dict{String,Array{Any,1}}()
    DictCosmo["w0"]  = [-1.,   "present"]
    DictCosmo["wa"]  = [0.,    "present"]
    DictCosmo["Mν"]  = [0.06,  "present"]
    DictCosmo["H0"]  = [67.,   "present"]
    DictCosmo["ΩM"]  = [0.32,  "present"]
    DictCosmo["ΩB"]  = [0.05,  "present"]
    DictCosmo["ns"]  = [0.96,  "present"]
    DictCosmo["σ8"]  = [0.816, "present"]
    CosmoModel = "Flatw0waCDMCosmology"

    DictIA = Dict{String,Array{Any,1}}()
    DictIA["𝓐IA"] = [1.72, "present"]
    DictIA["βIA"] = [2.17, "present"]
    DictIA["𝓒IA"] = [0.0134, "not_present"]
    DictIA["ηIA"] = [-0.41, "present"]
    IAModel = "ExtendedNLIA"

    DictBias = Dict{String,Array{Any,1}}()
    DictBias["A"] = [1.0, "present"]
    DictBias["B"] = [2.5, "present"]
    DictBias["C"] = [2.8, "present"]
    DictBias["D"] = [1.6, "present"]
    BiasModel = "EuclidBias"

    MultipolesArrayTemp = CosmoCentral.LogSpaced(10.,3000., 101)
    MultipolesArray = zeros(100)
    MultipolesWidths = CosmoCentral.Difference(MultipolesArrayTemp)
    for i in 1:100
        MultipolesArray[i] = (MultipolesArrayTemp[i+1]+MultipolesArrayTemp[i])/2
    end

    steps = Array([0.00625, 0.01250, 0.01875, 0.02500, 0.03750, 0.05000, 0.10000])
    Cosmologies, IntrinsAlignment, Bias = CosmoCentral.CreateCosmologies(DictCosmo, CosmoModel, DictIA, IAModel, DictBias, BiasModel, steps)
    CosmologicalGrid = CosmoCentral.CosmologicalGrid(ZArray = LinRange(0.001, 4., 500), KArray = CosmoCentral.LogSpaced(1e-5, 50., 1000), ℓBinCenters = MultipolesArray, ℓBinWidths = MultipolesWidths);
    CosmoCentral.CreateDirectoriesForecast!(Cosmologies, DictCosmo, IntrinsAlignment, DictIA, Bias, DictBias ,pwd()*"/test_forecast/")
    CosmoCentral.ForecastCℓ!(Cosmologies, IntrinsAlignment, Bias, pwd()*"/forecast_pmm/PowerSpectrum/",pwd()*"/test_forecast/Angular/", CosmologicalGrid, "../input_files/Angular.json", "../input_files/scaledmeanlum-E2Sa.txt")
    CosmoCentral.Forecast∂Cℓ!(DictCosmo, DictIA, DictBias, pwd()*"/test_forecast/", "../input_files/Angular.json", steps);
    PathCentralCℓ = pwd()*"/test_forecast/Angular/dvar_central_step_0/cl"
    Path∂Cℓ = pwd()*"/test_forecast/Derivative"
    InputList = [DictCosmo, DictIA, DictBias]
    Fisheraₗₘ = CosmoCentral.ForecastFisherαβ(PathCentralCℓ, Path∂Cℓ, InputList,
    CosmologicalGrid)
    FisherCℓ = CosmoCentral.ForecastFisherαβ(PathCentralCℓ, Path∂Cℓ, InputList,
    CosmologicalGrid, "test")
    CheckFisher = CosmoCentral.ReadFisher("CheckFisher", "Lensing_Lensing")
    CosmoCentral.SelectMatrixAndMarginalize!(CheckFisher.ParametersList, CheckFisher)
    @test isapprox(CheckFisher.FisherMatrix, Fisheraₗₘ.FisherMatrix, rtol=1e-6)
    @test isapprox(CheckFisher.FisherMatrix, FisherCℓ.FisherMatrix, rtol=1e-6)
    for key in Fisheraₗₘ.SelectedParametersList
        @test isapprox(CheckFisher.MarginalizedErrors[key],
        Fisheraₗₘ.MarginalizedErrors[key], rtol=1e-6)
        @test isapprox(CheckFisher.MarginalizedErrors[key],
        FisherCℓ.MarginalizedErrors[key], rtol=1e-6)
    end
end


text:
struct Params_SV{ComplexArray,Total,KF}
    # boundary conditions
    periodic::Bool
    # normalization
    total::Total

    # algorithm-specific
    # gradient kernel
    kernelfactor::KF

    # fft buffers
    complex_surface::ComplexArray
    complex_void::ComplexArray
end

function Params_SV(img; periodic::Bool=true, kernelfactor=KernelFactors.sobel)
    box = size(img)
    complex_box = periodic ? box : box .* 2

    # scale factor
    total = cnt_total(img, periodic)

    p = Params_SV(
        periodic,
        total,
        kernelfactor,
        similar(img, ComplexF32, complex_box),
        similar(img, ComplexF32, complex_box)
    )
    cf_type = :full
    p, cf_type
end

"""
    correllation_function!(res, img, params::Params_SV)

Compute surface-void correlation function in positive directions
"""
function correllation_function!(res, img, params::Params_SV)
    ix = CartesianIndices(img)

    f = params.complex_surface .= 0
    g = params.complex_void .= 0
    v_f = view(f, ix)
    v_g = view(g, ix)
    v_f .= gradient_norm(img, params.kernelfactor)
    v_g .= img .== 0

    cross_correlation!(f, g)

    res .= real.(v_f) ./ params.total
end

"""
    surfvoid(image; periodic = false)

Calculate `Fsv(x)` (surface-void) correlation function map
for the N-dimensional image and return a `CFMap` object.

The `image` contains the probability of the voxel being in the correct phase.

# Examples
```jldoctest
julia> surfvoid([1 0; 0 1]; periodic=true).result
3×3 Matrix{Float32}:
 0.176777  0.176777  0.176777
 0.176777  0.176777  0.176777
 0.176777  0.176777  0.176777
```
"""
function surfvoid(image; periodic::Bool=false, kernelfactor=KernelFactors.sobel)
    corr_function_map(image, Params_SV; periodic, kernelfactor)
end

function surfvoid(image, phase::Int; periodic::Bool=false, kernelfactor=KernelFactors.sobel)
    one_phase_img = image .== phase
    corr_function_map(one_phase_img, Params_SV; periodic, kernelfactor)
end


text:
module StringUtils

export is_palindrome
function is_palindrome(s)
    len = length(s)
    for i = 1:div(len, 2)
        if s[i] != s[len + 1 - i]
            return false
        end
    end
    return true
end

export conc
function conc(a::Array{String,1})
    buf = IOBuffer()
    for s in a
        print(buf, chomp(s))
    end
    takebuf_string(buf)
end

export grid_to_matrix
function grid_to_matrix(a::Array{String,1})
    result::Array{Array{Int64,1},1} = []
    for s in a
        row = map(x -> parse(Int64, x), split(chomp(s)))
        push!(result, row)
    end
    return result
end

end


text:
############################################################################################
"""
$(TYPEDEF)
Diagnostics fields for HMC sampler.

# Fields
$(TYPEDFIELDS)
"""
struct DiagnosticsHMC{T<:AbstractFloat} <: MCMCKernelDiagnostics
    "Discretization size"
    ϵ::T
    "Number of Leapfrog steps"
    steps::Int64
    function DiagnosticsHMC(ϵ::T, steps::Int64) where {T<:AbstractFloat}
        @argcheck steps > 0 "Minimum number of Leapfrog steps is 1."
        @argcheck ϵ > 0.0 "Discretization size has to be positive."
        return new{T}(ϵ, steps)
    end
end

function results(
    diagnostics::AbstractVector{DiagnosticsHMC{F}}, Ndigits::Integer, quantiles::Vector{T}
) where {F<:AbstractFloat,T<:Real}
    println(
        "HMC sampler has average (final) stepsize of ",
        round.(
            mean(diagnostics[iter].ϵ for iter in eachindex(diagnostics)); digits=Ndigits
        ),
        " (",
        round.(diagnostics[end].ϵ; digits=Ndigits),
        ") with ",
        round.(
            mean(diagnostics[iter].steps for iter in eachindex(diagnostics)); digits=Ndigits
        ),
        " (",
        diagnostics[end].steps,
        ") number of steps.",
    )
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Show relevant diagnostic results.

# Examples
```julia
```

"""
function generate_showvalues(diagnostics::D) where {D<:DiagnosticsHMC}
    return function showvalues()
        return (:HMC, "diagnostics"), (:ϵ, diagnostics.ϵ), (:steps, diagnostics.steps)
    end
end

############################################################################################
#export
export DiagnosticsHMC, generate_showvalues


text:
export MultiFunction, first

"""
$(TYPEDEF)

An immutable structure to hold an elementary function F: Rᵏ → R
that can be decomposed as the product of univariate basis,
where each basis is finite and contains constant and/or linear
and Hermite functions
F(x_1, x_2, ..., x_k) = f_1(x_1) × f_2(x_2) × ... × f_M(x_k)

## Fields

$(TYPEDFIELDS)

"""
struct MultiFunction
    m::Int64
    Nx::Int64
    B::MultiBasis
    α::Array{Int64,1}
    function MultiFunction(B::MultiBasis, α::Array{Int64,1})
        m = B.B.m
        Nx = B.Nx
        @assert Nx == size(α,1) "Dimension of the space doesn't match the size of α"
        for i=1:Nx
            @assert α[i]<=m "multi index α can't be greater than the size of the univariate basis "
        end
        return new(m, Nx, B, α)
    end
end

function MultiFunction(B::MultiBasis)
    return MultiFunction(B.B.m, B.Nx, B, ones(Int64, B.Nx))
end

function MultiFunction(B::Basis, Nx::Int64; scaled::Bool = true)
    return MultiFunction(B.B.m, Nx, MultiBasis(k, B), ones(Int64, Nx))
end

"""
$(TYPEDSIGNATURES)
Evaluates the `MultiFunction` `F` at `x`
"""
function (F::MultiFunction)(x::Array{T,1}) where {T <: Real}
    out = 1.0
    @inbounds for i=1:F.Nx
        out *= F.B.B[F.α[i]+1](x[i])
    end
    return out
end


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector

radius, height, offset1, offset2, height_fraction = 3, 2, 1, 1.5, .1
file_name = joinpath(@__DIR__, "current_carrying_wire_cylinder.geo")
geo_write!(file_name, characteristic_length_factor=.3,
    footer="""
    Cylinder(1) = {0, 0, 0, 0, 0, $height, $radius, 2*Pi};
    Point(3) = {0, 0, 0, 1.0};
    Point(4) = {0, 0, $height, 1.0};
    Point(5) = {$offset1, 0, $(height_fraction * height), 1.0};
    Point(6) = {$offset1, 0, $((1-height_fraction) * height), 1.0};
    Point(7) = {$offset2, 0, $(height_fraction * height), 1.0};
    Point(8) = {$offset2, 0, $((1-height_fraction) * height), 1.0};
    Point(9) = {$(-offset1), 0, 0, 1.0};
    Point(10) = {$(-offset1), 0, $height, 1.0};
    Line(4) = {3, 4};
    Line(5) = {5, 6};
    Line(6) = {6, 8};
    Line(7) = {8, 7};
    Line(8) = {7, 5};
    Line(9) = {9, 10};
    Point{3, 9} In Surface{3};
    Point{4, 10} In Surface{2};
    Point{6, 8, 7, 5} In Volume{1};
    Curve{4, 5, 6, 7, 8, 9} In Volume{1};
    Physical Curve("Wire center") = {4};
    Physical Curve("Loop") = {6, 7, 8, 5};
    Physical Curve("Wire offset") = {9};
    """)
initialize!()
gmsh_open!(file_name)
N, K = 3, 4
mesh!(K)

node_tags, points, tcomp = get_triangulated_complex(N, K)
group_dict = get_physical_groups(node_tags, points)
@test typeof(tcomp) <: TriangulatedComplex{N, K}
comp = tcomp.complex
orient!(comp)
sources = [get_current_source(comp, group_dict[k]) for k in
    ["Wire center", "Loop", "Wire offset"]]
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

μ⁻, Λ⁻ = 2, 0
μ⁻_form = get_material(comp, μ⁻, 3)
Λ⁻_form = get_material(comp, Λ⁻, 2)

bbox, null_basis = magnetostatics_blackbox(m, mesh, sources, boundary(comp),
    μ⁻_form, Λ⁻_form)
inductance = impedance_matrix(bbox, null_basis)
@testset "mutual inductance of wires and loops" begin
    @test isapprox(abs(inductance[1,3]),
        abs((log(radius) - log(offset1)) * height/(2π*μ⁻)), rtol=3e-2)
    @test isapprox(abs(inductance[1,2]),
        abs((log(offset1) - log(offset2)) * height * (1-2*height_fraction)/(2π*μ⁻)), rtol=3e-2)
end

if false # set to true to display plots with Gmsh
    A, _ = solve_statics(bbox, null_basis, [0.0, 1.0, 0.0])
    vec_A = sharp(m, comp, A)
    comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    add_field!("Vector potential", node_tags, vec_A[ordering])
    gui!()
end
rm(file_name)


text:
using Test
using SpiderMonkey: smo
using BlackBoxOptimizationBenchmarking

@testset "Optimize Sphere function" begin
    D = 100 # Dimensions
    P = 30 # Population
    fn = BlackBoxOptimizationBenchmarking.F1.f # Function to optimize (in example format)
    ɳ = 0.5 # Perturbation rate
    λ = 25 # Local leader Limit
    α = 20 # Global Leader Limit
    β = 8 # Max Group Count
    τ = 1000 # Iterations
    ub = fill(100,(D,)) # Upper bound for value
    lb = fill(-100,(D,)) # Lower bound for values
    verbose = false
    GL,GL_fitness = smo(D,P,fn,ɳ,λ,α,β,τ,ub,lb,verbose)
    @test GL_fitness ≈ BlackBoxOptimizationBenchmarking.F1.f_opt atol=0.01
    @test size(GL)[1] == D
end

@testset "Bounds Error" begin
    D = 100 # Dimensions
    P = 30 # Population
    fn = BlackBoxOptimizationBenchmarking.F1.f # Function to optimize (in example format)
    ɳ = 0.5 # Perturbation rate
    λ = 25 # Local leader Limit
    α = 20 # Global Leader Limit
    β = 8 # Max Group Count
    τ = 1000 # Iterations
    ub = fill(100,(90,)) # Upper bound for value
    lb = fill(-100,(90,)) # Lower bound for values
    verbose = false
    @test_throws BoundsError smo(D,P,fn,ɳ,λ,α,β,τ,ub,lb,verbose)
end


text:
# compare_coercive_kernel_discrepancies
#
# This experiment highlights the need to properly account
# for the possibility that distributions could diverge
# infinity and look suitable vis-a-vis the KSD.

include("experiment_utils.jl")

# dimension of the problem
@parseintcli d "d" "dimension" 8
# random number generator seed
@parseintcli seed "s" "seed" 7
# set the sample distribution {gaussian, studentt, epspacking, randomepspacking}
@parsestringcli samplesource "Q" "samplesource" "gaussian"
# set the discrepancy
@parsestringcli discrepancytype "k" "discrepancytype" "inversemultiquadric"
# should we use multiple cores?
@parseintcli numcores "n" "numcores" 1
# parallelize it
if (numcores > 1)
    if (numcores > CPU_CORES)
        error(@sprintf("Requested %d cores, only %d available.", numcores, CPU_CORES))
    end
    addprocs(numcores - 1)
    @assert (numcores == nprocs())
end

using Distributions: Exponential
using Iterators: product

using SteinDistributions:
    SteinGaussian,
    SteinScaleLocationStudentT,
    gradlogdensity
using SteinDiscrepancy:
    SteinGaussianKernel,
    SteinMaternRadialKernel,
    SteinGaussianPowerKernel,
    SteinInverseMultiquadricKernel,
    ksd

# specify the target distribution
target = SteinGaussian(d)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# set up student t distribution
df = 10.0
studentt = SteinScaleLocationStudentT(df, 0.0, sqrt((df-2.0)/df))
# set the kernel
if discrepancytype == "matern"
    kernel = SteinMaternRadialKernel()
elseif discrepancytype == "gaussian"
    kernel = SteinGaussianKernel()
elseif discrepancytype == "gaussianpower"
    kernel = SteinGaussianPowerKernel(2.0)
elseif discrepancytype == "inversemultiquadric"
    kernel = SteinInverseMultiquadricKernel(0.5)
end
# Sample sizes at which optimization problem will be solved
ns = vcat(100:100:1000, 2000:1000:10000, 20000:10000:100000)
# the radii to use for the lattace generation
rs = [0.2, 0.1, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03,
      0.02, 0.015, 0.012, 0.01, 0.0085, 0.007]

function packsphere3d(r::Float64)
    # we'll put the lattice in [-1, 1]^3 \cap B(0,1)
    I = ceil(1 / r)
    J = ceil(2 / (r * sqrt(3)))
    K = ceil(3 / (r * sqrt(6)))
    packing = Array(Float64, 0, 3)
    for (i,j,k) in product(0:I, 0:J, 0:K)
        pt = [2*i + mod(j + k, 2),
              sqrt(3)*(mod(k,2)/3.0 + j),
              2*sqrt(6)/3*k]
        pt = pt - ones(3)
        pt = pt .* r
        if norm(pt) <= 1
            packing = vcat(packing, pt')
        end
    end
    packing
end

function packsphererandom(n::Int, d::Int; deltan = n::Int -> 2*log(n))
    packing = Array(Float64, 0, d)
    @assert (d >= 3)

    delta = deltan(n)
    r = delta * n^(1/d)
    ii = 0
    depth = Exponential(1.0)

    while size(packing, 1) < n && (ii <= 100*n)
        ii += 1
        x0 = rand(target, 1)
        x0 = r * x0 / (norm(x0)^2 + rand(depth))^(0.5)
        shouldadd = true
        for jj in 1:size(packing,1)
            if norm(packing[jj:jj,:] - x0) <= delta
                shouldadd = false
                break
            end
        end
        if shouldadd
            packing = vcat(packing, x0)
        end
    end
    if (ii > 100*n)
        error(@sprintf("For n=%d the packing failed.", n))
    end
    packing
end

println("Beginning optimization")
continueloop = true
idx = 0
while continueloop
    idx += 1
    if samplesource == "gaussian"
        ii = ns[idx]
        X = @setseed seed rand(target, ii)
        continueloop = (idx < length(ns))
    elseif samplesource == "studentt"
        ii = ns[idx]
        Xvec = @setseed seed rand(studentt, ii*d)
        X = reshape(Xvec, d, ii)'
        continueloop = (idx < length(ns))
    elseif samplesource == "epspacking"
        @assert (d == 3)
        rr = rs[idx]
        X = packsphere3d(rr)
        ii = size(X,1)
        X = sqrt(ii / log(ii)) .* X
        continueloop = (idx < length(rs))
    elseif samplesource == "randomepspacking"
        ii = ns[idx]
        X = @setseed seed packsphererandom(ii, d)
        continueloop = (idx < length(ns))
    end
    kernelresult = ksd(points=X, gradlogdensity=gradlogp, kernel=kernel)
    kerneldiscrepancy = sqrt(kernelresult.discrepancy2)
    solvetime = kernelresult.solvetime
    # save data
    instance_data = Dict{Any, Any}(
        "samplesource" => samplesource,
        "discrepancytype" => discrepancytype,
        "n" => ii,
        "d" => d,
        "seed" => seed,
        "steindiscrepancy" => kerneldiscrepancy,
        "solvetime" => solvetime,
        "ncores" => nprocs(),
    )
    save_json(
        instance_data;
        dir="compare_coercive_kernel_discrepancies",
        samplesource=samplesource,
        discrepancytype=discrepancytype,
        n=ii,
        d=d,
        seed=seed,
    )
end
println("COMPLETE!")


text:
using GeneralizedSampling
import IntervalWavelets: unit
using Base.Test

println("Testing least squares solver...")

#=
Test the least squares solvers:
- Solution with Freq2Wave and the collected version should be identical.
- When sampling from a scaling function used in the Freq2Wave type, the solution should be a unit vector.
=#

# ------------------------------------------------------------
# Reconstruct unit vector

J = 3
M = 2^(J+1)
N = 2^J
xi = grid(M, 0.5)
wavename = "db2"
vm = van_moment(wavename)

T = Freq2Wave(xi, wavename, J)

for k = 1:N
	if k <= vm
		F = FourScalingFunc(xi, wavename, 'L', J)
		f = F[:,k]
	elseif k > N-vm
		F = FourScalingFunc(xi, wavename, 'R', J)
		f = F[:,k-N+vm]
	else
		trans = k - div(N,2) - 1
		f = FourDaubScaling(xi, vm, J, trans)
	end
	y = T \ f

	u = unit(Complex{Float64}, length(y), k)
	#@show k, norm( y - u )
	@test_approx_eq_eps y u 1e-4
end

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

begin
	J = 10
	M = 2^(J+2)
	N = 2^J

	# Uniform samples
	xi = grid(M, 0.5)
	TU = Freq2Wave(xi, "Haar", J)
	AU = collect(TU)

	b = rand(M)
	x0 = rand(N)

	x1 = AU \ b
	x2 = TU \ b

	#@show norm(x1 - x2, Inf)
	@test_approx_eq_eps x1 x2 1e-5

	# Non-uniform samples: Requires high M/N ratio
	K = N/2
	xi = 2*K*rand(M) - K
	sort!(xi)
	TNU = Freq2Wave(xi, "Haar", J, K)
	ANU = collect(TNU)
	b = rand(M) + rand(M)*im

	bb = b .* get(TNU.weights)
	z1 = ANU \ bb
	z2 = TNU \ b

	#@show norm(z1 - z2, Inf)
	@test_approx_eq_eps z1 z2 1e-4
end

# 2D
begin
	J = 5
	M = 2^(J+1)
	N = 2^J

	# Uniform samples
	xi = grid((M,M), 0.5)
	TU = Freq2Wave(xi, "Haar", J)
	AU = collect(TU)
	b = rand(size(TU,1))

	x1 = AU \ b
	x2 = TU \ b

	#@show norm(x1 - vec(x2), Inf)
	@test_approx_eq_eps x1 x2 1e-5

	# Non-uniform samples: Requires high M/N ratio
	K = N/2
	xi = 2*K*rand(M^2,2) - K
	TNU = Freq2Wave(xi, "Haar", J, K)
	ANU = collect(TNU)
	b = rand(M^2) + rand(M^2)*im

	bb = b .* get(TNU.weights)
	z1 = ANU \ bb
	z2 = TNU \ b

	#@show norm(z1 - z2, Inf)
	@test_approx_eq_eps z1 z2 1e-4
end


text:
module LocalCoverage

using CoverageTools
using DocStringExtensions
using Printf
using PrettyTables
import Pkg

export generate_coverage, open_coverage, clean_coverage, coverage_summary

"Directory for coverage results."
const COVDIR = "coverage"

"Coverage tracefile."
const LCOVINFO = "lcov.info"

const PYTHON = get!(ENV, "PYTHON", isnothing(Sys.which("python3")) ? "python" : "python3")

"""
$(SIGNATURES)

Get the root directory of a package.
"""
function pkgdir(pkgstr::AbstractString)
    joinpath(dirname(Base.locate_package(Base.PkgId(pkgstr))), "..")
end
pkgdir(m::Module) = joinpath(dirname(pathof(m)), "..")

"""
$(SIGNATURES)

Open the HTML coverage results in a browser for `pkg` if they exist.

See [`generate_coverage`](@ref).
"""
function open_coverage(pkg;
                       coverage_file::AbstractString=joinpath(COVDIR, "index.html"))
    htmlfile = joinpath(pkgdir(pkg), coverage_file)
    if !isfile(htmlfile)
        @warn("Not found, run generate_coverage(pkg) first.")
        return nothing
    end
    try
        if Sys.isapple()
            run(`open $htmlfile`)
        elseif Sys.islinux() || Sys.isbsd()
            run(`xdg-open $htmlfile`)
        elseif Sys.iswindows()
            run(`start $htmlfile`)
        end
    catch e
        error("Failed to open the generated $(htmlfile)\n",
              "Error: ", sprint(Base.showerror, e))
    end
    nothing
end

"""
$(SIGNATURES)

Returns a table giving coverage details by file in human readable form:
- column 1: source file name
- column 2: tuple (lines_hit, lines_tracked)
- column 3: coverage fraction
"""
function coverage_summary_table(coverage)
    n = length(coverage)
    tab = Array{Any, 2}(undef, 1+n, 3)
    total_hit = 0
    total_tracked    = 0

    for (i, f) in enumerate(coverage)
        hit     = count(x->!isnothing(x) && x>0, f.coverage)
        tracked = count(x->!isnothing(x),        f.coverage)

        total_hit     += hit
        total_tracked += tracked

        tab[i,:] .= (f.filename, (hit, tracked),
                     100 * hit / tracked)
    end

    tab[n+1, :] .= ("TOTAL", (total_hit, total_tracked),
                    100 * total_hit / total_tracked)
    tab
end

"""
$(SIGNATURES)

Pretty-prints a table giving coverage details by file.
"""
function coverage_summary(coverage)
    tab = coverage_summary_table(coverage)

    highlighters = (
        Highlighter((data,i,j)->j==3 && data[i,j] <= 50,
                    bold       = true,
                    foreground = :red),
        Highlighter((data,i,j)->j==3 && data[i,j] <= 70,
                    foreground = :yellow),
        Highlighter((data,i,j)->j==3 && data[i,j] >= 90,
                    foreground = :green),
    )

    formatter(value, i, j) = if j==3
        isnan(value) ? "-" : @sprintf("%3.0f%%", value)
    elseif j==2
        hit, tracked = value
        @sprintf("%3d / %3d", hit, tracked)
    else
        value
    end

    pretty_table(tab,
                 ["File name", "Lines hit", "Coverage"],
                 alignment = [:l, :r, :r],
                 highlighters = highlighters,
                 body_hlines = [size(tab, 1)-1],
                 formatters = formatter)
end

"""
    generate_xml(pkg, filename="cov.xml")

Generate a coverage Cobertura XML in the package `coverage` directory.

This requires the Python package `lcov_cobertura`, available in PyPl via
`pip install lcov_cobertura`.
"""
function generate_xml(pkg, filename="cov.xml")
    run(Cmd(Cmd([PYTHON, "-m", "lcov_cobertura", "lcov.info", "-o", filename]),
            dir=joinpath(pkgdir(pkg),COVDIR)))
    @info("generated cobertura XML $filename")
end

"""
$(SIGNATURES)

Generate a coverage report for package `pkg`.

When `genhtml`, the corresponding external command will be called to generate a
HTML report. This can be found in eg the package `lcov` on Debian/Ubuntu.

If `genxml` is true, will generate a Cobertura XML in the `coverage` directory
(requires Python package `lcov_cobertura`, see `generate_xml`).

If `show_summary` is true, a summary will be printed to `stdout`.

`*.cov` files are near the source files as generated by Julia, everything else
is placed in `Pkg.dir(pkg, \"$(COVDIR)\")`. The summary is in
`Pkg.dir(pkg, \"$(COVDIR)\", \"$(LCOVINFO)\")`.

Use [`clean_coverage`](@ref) for cleaning.
"""
function generate_coverage(pkg; genhtml=true, show_summary=true, genxml=false)
    Pkg.test(pkg; coverage = true)
    coverage = cd(pkgdir(pkg)) do
        coverage = CoverageTools.process_folder()
        isdir(COVDIR) || mkdir(COVDIR)
        tracefile = "$(COVDIR)/lcov.info"
        CoverageTools.LCOV.writefile(tracefile, coverage)
        if genhtml
            branch =
                try
                    strip(read(`git rev-parse --abbrev-ref HEAD`, String))
                catch
                    @warn "git branch could not be detected"
                end
            title = "on branch $(branch)"
            try
                run(`genhtml -t $(title) -o $(COVDIR) $(tracefile)`)
            catch e
                error("Failed to run genhtml. Check that lcov is installed (see the README).",
                      "\nError message: ", sprint(Base.showerror, e))
            end
            @info("generated coverage HTML")
        end
        coverage
    end
    genxml && generate_xml(pkg)
    show_summary && coverage_summary(coverage)
    coverage
end

"""
$(SIGNATURES)

Clean up after [`generate_coverage`](@ref).

If `rm_directory`, will delete the coverage directory, otherwise only deletes
`*.cov` coverage output.
"""
function clean_coverage(pkg;
                        coverage_directory::AbstractString=COVDIR,
                        rm_directory::Bool=true)
    CoverageTools.clean_folder(pkgdir(pkg))
    rm_directory && rm(joinpath(pkgdir(pkg), coverage_directory); force = true, recursive = true)
end

end # module


text:
module JuliaCon

import Base: string
using Preferences
using Distributed
using Dates: Dates, Date, DateTime, Time, Hour, Minute
using TimeZones
using JSON
using Downloads: download
using PrettyTables
using TimerOutputs

include("schedule_structs.jl")
include("preferences.jl")
include("countries.jl")
include("tshirtcode.jl")
include("schedule.jl")

function __init__()
    if isdefined(Main, :Distributed)
        env = Base.active_project()
        @everywhere env = $env

        # activate local env of the master process on all the workers
        @eval Main @everywhere begin
            if myid() != 1
                # suppress REPL warnings and "activating ..." messages
                using Pkg, Logging
                with_logger(NullLogger()) do
                    Pkg.activate(env; io=devnull)
                end
            end
        end
        # load JuliaCon.jl on all workers
        @eval Main @everywhere using JuliaCon
    end
end

export juliacon2021

end


text:
export SGD,DebugOptAlg,ADAM

# ------------------------------------------------------------------------------
# SGD
"""
  SGD(;η=t -> 1/(1+t), λ=2)

Stochastic Gradient Descent algorithm (default)

# Fields:
- `η`: Learning rate, as a function of the current epoch [def: t -> 1/(1+t)]
- `λ`: Multiplicative constant to the learning rate [def: 2]
"""
struct SGD <: OptimisationAlgorithm
    η::Function
    λ::Float64
    function SGD(;η=t -> 1/(1+t), λ=2)
        return new(η,λ)
    end
end

function singleUpdate!(θ,▽,optAlg::SGD;nEpoch,nBatch,nBatches,xbatch,ybatch)
    η    = optAlg.η(nEpoch)*optAlg.λ
    #newθ = gradSub.(θ,gradMul.(▽,η))
    θ =  θ - ▽ .* η
    #newθ = gradientDescentSingleUpdate(θ,▽,η)
    return (θ=θ,stop=false)
end

#gradientDescentSingleUpdate(θ::Number,▽::Number,η) = θ .- (η .* ▽)
#gradientDescentSingleUpdate(θ::AbstractArray,▽::AbstractArray,η) = gradientDescentSingleUpdate.(θ,▽,Ref(η))
#gradientDescentSingleUpdate(θ::Tuple,▽::Tuple,η) = gradientDescentSingleUpdate.(θ,▽,Ref(η))

#maxEpochs=1000, η=t -> 1/(1+t), λ=1, rShuffle=true, nMsgs=10, tol=0

# ------------------------------------------------------------------------------
# ADAM
#

"""
  ADAM(;η, λ, β₁, β₂, ϵ)

The [ADAM](https://arxiv.org/pdf/1412.6980.pdf) algorithm, an adaptive moment estimation optimiser.

# Fields:
- `η`:  Learning rate (stepsize, α in the paper), as a function of the current epoch [def: t -> 0.001 (i.e. fixed)]
- `λ`:  Multiplicative constant to the learning rate [def: 1]
- `β₁`: Exponential decay rate for the first moment estimate [range: ∈ [0,1], def: 0.9]
- `β₂`: Exponential decay rate for the second moment estimate [range: ∈ [0,1], def: 0.999]
- `ϵ`:  Epsilon value to avoid division by zero [def: 10^-8]
"""
mutable struct ADAM <: OptimisationAlgorithm
    η::Function
    λ::Float64
    β₁::Float64
    β₂::Float64
    ϵ::Float64
    m::Vector{Learnable}
    v::Vector{Learnable}
    function ADAM(;η=t -> 0.001, λ=1.0, β₁=0.9, β₂=0.999, ϵ=1e-8)
        return new(η,λ,β₁,β₂,ϵ,[],[])
    end
end

"""
   initOptAlg!(optAlg::ADAM;θ,batchSize,x,y,rng)

Initialize the ADAM algorithm with the parameters m and v as zeros and check parameter bounds
"""
function initOptAlg!(optAlg::ADAM;θ,batchSize,x,y,rng = Random.GLOBAL_RNG)
    optAlg.m = θ .- θ # setting to zeros
    optAlg.v = θ .- θ # setting to zeros
    if optAlg.β₁ <= 0 || optAlg.β₁ >= 1 @error "The parameter β₁ must be ∈ [0,1]" end
    if optAlg.β₂ <= 0 || optAlg.β₂ >= 1 @error "The parameter β₂ must be ∈ [0,1]" end
end

function singleUpdate!(θ,▽,optAlg::ADAM;nEpoch,nBatch,nBatches,xbatch,ybatch)
    β₁,β₂,ϵ  = optAlg.β₁, optAlg.β₂, optAlg.ϵ
    η        = optAlg.η(nEpoch)*optAlg.λ
    t        = (nEpoch-1)*nBatches+nBatch
    optAlg.m = @. β₁ * optAlg.m + (1-β₁) * ▽
    optAlg.v = @. β₂ * optAlg.v + (1-β₂) * (▽*▽)
    #optAlg.v = [β₂ .* optAlg.v.data[i] .+ (1-β₂) .* (▽.data[i] .* ▽.data[i]) for i in 1:size(optAlg.v.data)]
    m̂        = @. optAlg.m /(1-β₁^t)
    v̂        = @. optAlg.v /(1-β₂^t)
    θ        = @. θ - (η * m̂) /(sqrt(v̂)+ϵ)
    return     (θ=θ,stop=false)
end

# ------------------------------------------------------------------------------
# DebugOptAlg

struct DebugOptAlg <: OptimisationAlgorithm
    dString::String
    function DebugOptAlg(;dString="Hello World, I am a Debugging Algorithm. I done nothing to your Net.")
        return new(dString)
    end
end

function singleUpdate!(θ,▽,optAlg::DebugOptAlg;nEpoch,nBatch,batchSize,ϵ_epoch,ϵ_epoch_l)
    println(optAlg.dString)
    return (θ=θ,stop=false)
end


text:
using PyPlot

using SubmodularMaximization

pygui(false)
fig_path = "./fig/probabilistic_coverage_test"
mkpath(fig_path)

########
# Params
########

num_events = 50
num_agents = 10
num_sensors = 5
nominal_area = 0.6

max_success_probability = 1.0
sensor_radius = sqrt(nominal_area / (num_agents * pi))
station_radius = 2*sensor_radius

#############
# environment
#############

figure()
events = generate_events(num_events)

########
# agents
########

agent_specification = ProbabilisticAgentSpecification(max_success_probability,
                                                      sensor_radius,
                                                      station_radius,
                                                      num_sensors)

agents = generate_agents(agent_specification, num_agents)
f(x) = mean_detection_probability(x, events)
problem = ExplicitPartitionProblem(f, agents)

####################
# visualize scenario
####################

figure()
visualize_pdf(standard_mixture())
visualize_events(events)
colors = generate_colors(agents)
xlim([0, 1])
ylim([0, 1])
savefig("$(fig_path)/scenario.png", pad_inches=0.00, bbox_inches="tight")
title("Scenario")

############
# evaluation
############
function evaluate_solver(solver, name)
  println("$name solver running")
  @time solution = solver(problem)

  figure()
  xlim([0, 1])
  ylim([0, 1])
  colors = generate_colors(agents)
  visualize_agents(agents, colors)
  visualize_solution(problem, solution, colors)
  visualize_events(events)

  coverage = solution.value

  savefig("$(fig_path)/$(to_file(name)).png", pad_inches=0.00, bbox_inches="tight")

  title("$name Solver ($(round(coverage, digits=3)))")

  @show coverage
end

# evaluate basic solvers
#evaluate_solver(solve_myopic, "Myopic")
#evaluate_solver(solve_random, "Random")
evaluate_solver(solve_sequential, "Sequential")

# evaluate partition solvers
for num_partitions in [2, 4, 8]
  solve_n(p) = solve_n_partitions(num_partitions, p)
  evaluate_solver(solve_n, "Partition-$num_partitions")
end

#@show mean_weight(problem)
#@show total_weight(problem)

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing