text:
expr_removebackslash = GeneralTools.expr_removebackslash;

@testset "Check if regular expression pattern `$expr_removebackslash` works fine:" begin
	@test replace("\\\\using_gadfly_10_1.png", expr_removebackslash => "") == "using_gadfly_10_1.png"
	@test replace("\\blas\\(0)cu2\\using_gadfly_10_1.png", expr_removebackslash => "") == "blas\\(0)cu2\\using_gadfly_10_1.png"
	@test replace("/using_gadfly\\asdf_fd/_10_1.png", expr_removebackslash => "") == "using_gadfly\\asdf_fd/_10_1.png"
	@test replace("/using_gadfly/asdf_fd/_10_1.png", expr_removebackslash => "") == "using_gadfly/asdf_fd/_10_1.png"
	@test replace("using_gadfly/asdf_fd/_10_1.png", expr_removebackslash => "") == "using_gadfly/asdf_fd/_10_1.png"
	@test replace("using_gadfly_10_1.png", expr_removebackslash => "") == "using_gadfly_10_1.png"
end


text:
#=function parsestring(T::Type, s::AbstractString)
  v::Union{Nothing,T} = tryparse(T,s)

  return v
end

parsestring(Int, "3")
parsestring(Float64, "3.4")

function parsestringFloat64(s::AbstractString)
  v::Union{Nothing,Float64} = tryparse(Float64,s)

  return v
end

function parsestringInt(s::AbstractString)
  v::Union{Nothing,Int} = tryparse(Int,s)

  return v
end

parsestringInt("3")
parsestringFloat64("3.4")

function parsestring(s::AbstractString, ::Val{T}) where T
  v::Union{Nothing,T} = tryparse(T,s)

  return v
end

parsestring("3", Val(Int))
parsestring("3.4", Val(Float64))=#

#=function testreadcomp()
  df = GzipDecompressorStream(open("$(pwd())\\data\\CRSP-M.csv.gz")) |> CSV.read |> DataFrame
  println("compressed-size: $(size(df))")
end

function testread()
  df = CSV.read("$(pwd())\\data\\CRSP-M.csv") |> DataFrame
  #ostream = GzipCompressorStream(open("$(pwd())\\data\\CRSP-M.jls.gz", "w"))
  #serialize(ostream,df)
  #close(ostream)
  println("direct-size: $(size(df))")
end

function testreadbin()
  df = deserialize("$(pwd())\\data\\CRSP-M.jls")
  println("direct-size: $(size(df))")
end

function testreadbincomp()

b = open(GzipDecompressorStream, "$(pwd())\\data\\CRSP-M.jls.gz") do iStream
  b = IOBuffer(read(iStream))
  (b)
end

df = deserialize(b)

  println("direct-size: $(size(df))")
end

@time testreadcomp()
@time testread()
@time testreadbin()
@time testreadbincomp()=#
#@time gz2lz4("data\\CS-A.csv.gz")

#=
using CodecLz4, CodecZlib

function gz2lz4mwe(infile::String, outfile::String)

  @assert isfile(infile)

  indata = open(infile, "r")
  outdata = open(outfile, "w")

  s = GzipDecompressorStream(LZ4CompressorStream(outdata))

  write(s, indata)
  close(s)
end

gz2lz4mwe("data\\CS-A.csv.gz", "data\\CS-A.csv.lz4")=#

using Revise
using DataFrames, CodecLz4, Serialization

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

function lz4mwe(infile::String, outfile::String)
  @assert isfile(infile)
  df = deserialize(infile)
  println("size after read: $(size(df))")
  outdata = open(outfile, "w")

  s = LZ4CompressorStream(open(outfile, "w"))

  serialize(s, df)
  close(s)

  s = LZ4DecompressorStream(open(outfile))

  df = deserialize(s)

  close(s)
  println("size after decompression: $(size(df))")

end

#lz4mwe("data\\CRSP-D.jls", "data\\CRSP-D.csv.lz4")

ext: jl

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Test
using ThreeBodyDecay

@testset "Line-shape functions" begin
    lsh = BreitWigner(0.77,0.15)
    @test real(amp(0.77^2, lsh)) ≈ 0.0
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
using Distributions
using Turing
using Test
using Random

Random.seed!(125)

x = [1.5 2.0]

@model gibbstest(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  for i in 1:length(x)
    x[i] ~ Normal(m, sqrt(s))
  end
  s, m
end

alg = Gibbs(3000, CSMC(15, 1, :s), HMC(1, 0.2, 4, :m))
chain = sample(gibbstest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.1
println("mean(chain[:m]) = $(mean(chain[:m]))")
@test mean(chain[:m]) ≈ 7/6 atol=0.1

alg = CSMC(15, 5000)
chain = sample(gibbstest(x), alg)
@test mean(chain[:s]) ≈ 49/24 atol=0.1
@test mean(chain[:m]) ≈ 7/6 atol=0.1


text:
#=
Solves the "Offer Distribution Unknown" Model by value function
iteration and a second faster method discussed in the corresponding
quantecon lecture.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-14

References
----------

http://quant-econ.net/jl/odu.html

=#

using Distributions

"""
Unemployment/search problem where offer distribution is unknown

##### Fields

- `bet::Real` : Discount factor on (0, 1)
- `c::Real` : Unemployment compensation
- `F::Distribution` : Offer distribution `F`
- `G::Distribution` : Offer distribution `G`
- `f::Function` : The pdf of `F`
- `g::Function` : The pdf of `G`
- `n_w::Int` : Number of points on the grid for w
- `w_max::Real` : Maximum wage offer
- `w_grid::AbstractVector` : Grid of wage offers w
- `n_pi::Int` : Number of points on grid for pi
- `pi_min::Real` : Minimum of pi grid
- `pi_max::Real` : Maximum of pi grid
- `pi_grid::AbstractVector` : Grid of probabilities pi
- `quad_nodes::Vector` : Notes for quadrature ofer offers
- `quad_weights::Vector` : Weights for quadrature ofer offers

"""
type SearchProblem
    bet::Real
    c::Real
    F::Distribution
    G::Distribution
    f::Function
    g::Function
    n_w::Int
    w_max::Real
    w_grid::AbstractVector
    n_pi::Int
    pi_min::Real
    pi_max::Real
    pi_grid::AbstractVector
    quad_nodes::Vector
    quad_weights::Vector
end

"""
Constructor for `SearchProblem` with default values

##### Arguments

- `bet::Real(0.95)` : Discount factor in (0, 1)
- `c::Real(0.6)` : Unemployment compensation
- `F_a::Real(1), F_b::Real(1)` : Parameters of `Beta` distribution for `F`
- `G_a::Real(3), G_b::Real(1.2)` : Parameters of `Beta` distribution for `G`
- `w_max::Real(2)` : Maximum of wage offer grid
- `w_grid_size::Int(40)` : Number of points in wage offer grid
- `pi_grid_size::Int(40)` : Number of points in probability grid

##### Notes

There is also a version of this function that accepts keyword arguments for
each parameter

"""
function SearchProblem(bet=0.95, c=0.6, F_a=1, F_b=1, G_a=3, G_b=1.2,
                       w_max=2, w_grid_size=40, pi_grid_size=40)

    F = Beta(F_a, F_b)
    G = Beta(G_a, G_b)

    # NOTE: the x./w_max)./w_max in these functions makes our dist match
    #       the scipy one with scale=w_max given
    f(x) = pdf(F, x./w_max)./w_max
    g(x) = pdf(G, x./w_max)./w_max

    pi_min = 1e-3  # avoids instability
    pi_max = 1 - pi_min

    w_grid = linspace_range(0, w_max, w_grid_size)
    pi_grid = linspace_range(pi_min, pi_max, pi_grid_size)

    nodes, weights = qnwlege(21, 0.0, w_max)

    SearchProblem(bet, c, F, G, f, g,
                  w_grid_size, w_max, w_grid,
                  pi_grid_size, pi_min, pi_max, pi_grid, nodes, weights)
end

# make kwarg version
function SearchProblem(;bet=0.95, c=0.6, F_a=1, F_b=1, G_a=3, G_b=1.2,
                       w_max=2, w_grid_size=40, pi_grid_size=40)
    SearchProblem(bet, c, F_a, F_b, G_a, G_b, w_max, w_grid_size,
                  pi_grid_size)
end

function q(sp::SearchProblem, w, pi_val)
    new_pi = 1.0 ./ (1 + ((1 - pi_val) .* sp.g(w)) ./ (pi_val .* sp.f(w)))

    # Return new_pi when in [pi_min, pi_max] and else end points
    return clamp(new_pi, sp.pi_min, sp.pi_max)
end

"""
Apply the Bellman operator for a given model and initial value.

##### Arguments

- `sp::SearchProblem` : Instance of `SearchProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output.
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions

##### Returns

None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.

"""
function bellman_operator!(sp::SearchProblem, v::Matrix, out::Matrix;
                          ret_policy::Bool=false)
    # Simplify names
    f, g, bet, c = sp.f, sp.g, sp.bet, sp.c
    nodes, weights = sp.quad_nodes, sp.quad_weights

    vf = CoordInterpGrid((sp.w_grid, sp.pi_grid), v, BCnan, InterpLinear)

    # set up quadrature nodes/weights
    # q_nodes, q_weights = qnwlege(21, 0.0, sp.w_max)

    for w_i=1:sp.n_w
        w = sp.w_grid[w_i]

        # calculate v1
        v1 = w / (1 - bet)

        for pi_j=1:sp.n_pi
            _pi = sp.pi_grid[pi_j]

            # calculate v2
            function integrand(m)
                quad_out = similar(m)
                for i=1:length(m)
                    mm = m[i]
                    quad_out[i] = vf[mm, q(sp, mm, _pi)] * (_pi*f(mm) +
                                                            (1-_pi)*g(mm))
                end
                return quad_out
            end
            integral = do_quad(integrand, nodes, weights)
            # integral = do_quad(integrand, q_nodes, q_weights)
            v2 = c + bet * integral

            # return policy if asked for, otherwise return max of values
            out[w_i, pi_j] = ret_policy ? v1 > v2 : max(v1, v2)
        end
    end
    return out
end

function bellman_operator(sp::SearchProblem, v::Matrix;
                          ret_policy::Bool=false)
    out_type = ret_policy ? Bool : Float64
    out = Array(out_type, sp.n_w, sp.n_pi)
    bellman_operator!(sp, v, out, ret_policy=ret_policy)
end

"""
Extract the greedy policy (policy function) of the model.

##### Arguments

- `sp::SearchProblem` : Instance of `SearchProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
function get_greedy!(sp::SearchProblem, v::Matrix, out::Matrix)
    bellman_operator!(sp, v, out, ret_policy=true)
end

get_greedy(sp::SearchProblem, v::Matrix) = bellman_operator(sp, v,
                                                            ret_policy=true)

"""
Updates the reservation wage function guess phi via the operator Q.

##### Arguments

- `sp::SearchProblem` : Instance of `SearchProblem`
- `phi::Vector`: Current guess for phi
- `out::Vector` : Storage for output

##### Returns

None, `out` is updated in place to hold the updated levels of phi
"""
function res_wage_operator!(sp::SearchProblem, phi::Vector, out::Vector)
    # Simplify name
    f, g, bet, c = sp.f, sp.g, sp.bet, sp.c

    # Construct interpolator over pi_grid, given phi
    phi_f = CoordInterpGrid(sp.pi_grid, phi, BCnearest, InterpLinear)

    # set up quadrature nodes/weights
    q_nodes, q_weights = qnwlege(7, 0.0, sp.w_max)

    for (i, _pi) in enumerate(sp.pi_grid)
        integrand(x) = max(x, phi_f[q(sp, x, _pi)]).*(_pi*f(x) + (1-_pi)*g(x))
        integral = do_quad(integrand, q_nodes, q_weights)
        out[i] = (1 - bet)*c + bet*integral
    end
end

"""
Updates the reservation wage function guess phi via the operator Q.

See the documentation for the mutating method of this function for more details
on arguments
"""
function res_wage_operator(sp::SearchProblem, phi::Vector)
    out = similar(phi)
    res_wage_operator!(sp, phi, out)
    return out
end


text:
export PrioritizedDQNLearner

using Random
using Flux
using Zygote
using StatsBase: mean

mutable struct PrioritizedDQNLearner{
    Tq<:AbstractApproximator,
    Tt<:AbstractApproximator,
    Tf,
    S<:Union{Int,Nothing},
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    stack_size::S
    γ::Float32
    batch_size::Int
    update_horizon::Int
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    default_priority::Float32
    rng::R
end

function PrioritizedDQNLearner(;
    approximator::Tq,
    target_approximator::Tt,
    loss_func::Tf,
    stack_size::Union{Int,Nothing} = 4,
    γ::Float32 = 0.99f0,
    batch_size::Int = 32,
    update_horizon::Int = 1,
    min_replay_history::Int = 32,
    update_freq::Int = 1,
    target_update_freq::Int = 100,
    update_step::Int = 0,
    default_priority::Float32 = 100f0,
    seed = nothing,
) where {Tq,Tt,Tf}
    copyto!(approximator, target_approximator)
    rng = MersenneTwister(seed)
    PrioritizedDQNLearner(
        approximator,
        target_approximator,
        loss_func,
        stack_size,
        γ,
        batch_size,
        update_horizon,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        default_priority,
        rng,
    )
end

"""

!!! note
    The state of the observation is assumed to have been stacked,
    if `!isnothing(stack_size)`.
"""
(learner::PrioritizedDQNLearner)(obs) =
    obs |> get_state |>
    x ->
        send_to_device(device(learner.approximator), x) |> learner.approximator |>
        send_to_host

function RLBase.update!(learner::PrioritizedDQNLearner, batch)
    learner.update_step += 1
    learner.update_step % learner.update_freq == 0 || return nothing

    Q, Qₜ, γ, loss_func, update_horizon, batch_size = learner.approximator,
    learner.target_approximator,
    learner.γ,
    learner.loss_func,
    learner.update_horizon,
    learner.batch_size
    states, rewards, terminals, next_states = map(
        x -> send_to_device(device(Q), x),
        (batch.states, batch.rewards, batch.terminals, batch.next_states),
    )
    actions = CartesianIndex.(batch.actions, 1:batch_size)

    priorities = send_to_device(device(Q), Vector{Float32}())

    gs = gradient(params(Q)) do
        q = batch_estimate(Q, states)[actions]
        q′ = dropdims(maximum(batch_estimate(Qₜ, next_states); dims = 1), dims = 1)
        G = rewards .+ γ^update_horizon .* (1 .- terminals) .* q′

        batch_losses = loss_func(G, q)
        priorities = (Zygote.dropgrad(batch_losses) .+ 1f-10)
        mean(batch_losses)
    end

    update!(Q, gs)

    if learner.update_step % learner.target_update_freq == 0
        copyto!(Qₜ, Q)
    end

    send_to_host(priorities)
end

function RLBase.extract_experience(t::AbstractTrajectory, learner::PrioritizedDQNLearner)
    s = learner.stack_size
    h = learner.update_horizon
    n = learner.batch_size
    γ = learner.γ
    if length(t) > learner.min_replay_history
        # 1. sample indices based on priority
        inds = Vector{Int}(undef, n)
        valid_ind_range = isnothing(s) ? (1:(length(t)-h)) : (s:(1:(length(t)-h)))
        for i in 1:n
            ind, p = sample(learner.rng, get_trace(t, :priority))
            while ind ∉ valid_ind_range
                ind, p = sample(learner.rng, get_trace(t, :priority))
            end
            inds[i] = ind
        end

        # 2. extract SARTS
        states = consecutive_view(get_trace(t, :state), inds; n_stack = s)
        actions = consecutive_view(get_trace(t, :action), inds)
        next_states = consecutive_view(get_trace(t, :state), inds .+ h; n_stack = s)
        consecutive_rewards = consecutive_view(get_trace(t, :reward), inds; n_horizon = h)
        consecutive_terminals =
            consecutive_view(get_trace(t, :terminal), inds; n_horizon = h)
        rewards, terminals = zeros(Float32, n), fill(false, n)

        # 3. make sure that we only consider experiences in current episode for rewards and terminals
        for i in 1:n
            m = findfirst(view(consecutive_terminals, :, i))
            if isnothing(m)
                terminals[i] = false
                rewards[i] = discount_rewards_reduced(view(consecutive_rewards, :, i), γ)
            else
                terminals[i] = true
                rewards[i] = discount_rewards_reduced(view(consecutive_rewards, 1:m, i), γ)
            end
        end
        inds,
        (
            states = states,
            actions = actions,
            rewards = rewards,
            terminals = terminals,
            next_states = next_states,
        )
    else
        nothing
    end
end

function RLBase.update!(p::QBasedPolicy{<:PrioritizedDQNLearner}, t::AbstractTrajectory)
    indexed_experience = extract_experience(t, p)
    if !isnothing(indexed_experience)
        inds, experience = indexed_experience
        priorities = update!(p.learner, experience)
        if !isnothing(priorities)
            get_trace(t, :priority)[inds] .= priorities
        end
    end
end

function (
    agent::Agent{
        <:QBasedPolicy{<:PrioritizedDQNLearner},
        <:CircularCompactPSARTSATrajectory,
    }
)(
    ::PostActStage,
    obs,
)
    push!(
        agent.trajectory;
        reward = get_reward(obs),
        terminal = get_terminal(obs),
        priority = agent.policy.learner.default_priority,
    )
    nothing
end


text:
using LinearAlgebra: AbstractTriangular, AbstractRotation

# Copied from symmetric.jl
const AdjTransVec = Union{Transpose{<:Any,<:AbstractVector},Adjoint{<:Any,<:AbstractVector}}
const RealHermSym{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}}
const RealHermSymComplexHerm{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S}}
const RealHermSymComplexSym{T<:Real,S} = Union{Hermitian{T,S}, Symmetric{T,S}, Symmetric{Complex{T},S}}

# Ambiguities
for (a, b) in (
    (AbstractDimVector, AbstractDimMatrix),
    (AbstractDimMatrix, AbstractDimVector),
    (AbstractDimMatrix, AbstractDimMatrix),
    (AbstractDimMatrix, AbstractVector),
    (AbstractDimVector, AbstractMatrix),
    (AbstractDimMatrix, AbstractMatrix),
    (AbstractMatrix, AbstractDimVector),
    (AbstractVector, AbstractDimMatrix),
    (AbstractMatrix, AbstractDimMatrix),
    (AbstractDimVector, Adjoint{<:Any,<:AbstractMatrix}),
    (AbstractDimVector, AdjTransVec),
    (AbstractDimVector, Transpose{<:Any,<:AbstractMatrix}),
    (AbstractDimMatrix, Diagonal),
    (AbstractDimMatrix, Adjoint{<:Any,<:RealHermSymComplexHerm}),
    (AbstractDimMatrix, Adjoint{<:Any,<:AbstractTriangular}),
    (AbstractDimMatrix, Transpose{<:Any,<:AbstractTriangular}),
    (AbstractDimMatrix, Transpose{<:Any,<:RealHermSymComplexSym}),
    (AbstractDimMatrix, AbstractTriangular),
    (Diagonal, AbstractDimVector),
    (Diagonal, AbstractDimMatrix),
    (Transpose{<:Any,<:AbstractTriangular}, AbstractDimVector),
    (Transpose{<:Any,<:AbstractTriangular}, AbstractDimMatrix),
    (Transpose{<:Any,<:AbstractVector}, AbstractDimVector),
    (Transpose{<:Real,<:AbstractVector}, AbstractDimVector),
    (Transpose{<:Any,<:AbstractVector}, AbstractDimMatrix),
    (Transpose{<:Any,<:RealHermSymComplexSym}, AbstractDimMatrix),
    (Transpose{<:Any,<:RealHermSymComplexSym}, AbstractDimVector),
    (AbstractTriangular, AbstractDimVector),
    (AbstractTriangular, AbstractDimMatrix),
    (Adjoint{<:Any,<:AbstractTriangular}, AbstractDimVector),
    (Adjoint{<:Any,<:AbstractVector}, AbstractDimMatrix),
    (Adjoint{<:Any,<:RealHermSymComplexHerm}, AbstractDimMatrix),
    (Adjoint{<:Any,<:AbstractTriangular}, AbstractDimMatrix),
    (Adjoint{<:Number,<:AbstractVector}, AbstractDimVector{<:Number}),
    (AdjTransVec, AbstractDimVector),
    (Adjoint{<:Any,<:RealHermSymComplexHerm}, AbstractDimVector),
)
    @eval Base.:*(A::$a, B::$b) = _rebuildmul(A, B)
end

Base.:*(A::AbstractDimVector, B::Adjoint{T,<:AbstractRotation}) where T = _rebuildmul(A, B)
Base.:*(A::Adjoint{T,<:AbstractRotation}, B::AbstractDimMatrix) where T = _rebuildmul(A, B)
Base.:*(A::Transpose{<:Any,<:AbstractMatrix{T}}, B::AbstractDimArray{S,1}) where {T,S} = _rebuildmul(A, B)
Base.:*(A::Adjoint{<:Any,<:AbstractMatrix{T}}, B::AbstractDimArray{S,1}) where {T,S} = _rebuildmul(A, B)

function _rebuildmul(A::AbstractDimVector, B::AbstractDimMatrix)
    # Vector has no dim 2 to compare
    rebuild(A, parent(A) * parent(B), (first(dims(A)), last(dims(B)),))
end
function _rebuildmul(A::AbstractDimMatrix, B::AbstractDimVector)
    comparedims(last(dims(A)), first(dims(B)); val=true)
    rebuild(A, parent(A) * parent(B), (first(dims(A)),))
end
function _rebuildmul(A::AbstractDimMatrix, B::AbstractDimMatrix)
    comparedims(last(dims(A)), first(dims(B)); val=true)
    rebuild(A, parent(A) * parent(B), (first(dims(A)), last(dims(B))))
end
function _rebuildmul(A::AbstractDimVector, B::AbstractMatrix)
    rebuild(A, parent(A) * B, (first(dims(A)), AnonDim(Base.OneTo(size(B, 2)))))
end
function _rebuildmul(A::AbstractDimMatrix, B::AbstractVector)
    newdata = parent(A) * B
    if newdata isa AbstractArray
        rebuild(A, parent(A) * B, (first(dims(A)),))
    else
        newdata
    end
end
function _rebuildmul(A::AbstractDimMatrix, B::AbstractMatrix)
    rebuild(A, parent(A) * B, (first(dims(A)), AnonDim(Base.OneTo(size(B, 2)))))
end
function _rebuildmul(A::AbstractVector, B::AbstractDimMatrix)
    rebuild(B, A * parent(B), (AnonDim(Base.OneTo(size(A, 1))), last(dims(B))))
end
function _rebuildmul(A::AbstractMatrix, B::AbstractDimVector)
    newdata = A * parent(B)
    if newdata isa AbstractArray
        rebuild(B, A * parent(B), (AnonDim(Base.OneTo(1)),))
    else
        newdata
    end
end
function _rebuildmul(A::AbstractMatrix, B::AbstractDimMatrix)
    rebuild(B, A * parent(B), (AnonDim(Base.OneTo(size(A, 1))), last(dims(B))))
end


text:
function TV.as(d::ProductMeasure{F,S,A}) where {F,S,A<:AbstractArray}
    d1 = marginals(d).f(first(marginals(d).data))
    as(Array, as(d1), size(marginals(d))...)
end

###############################################################################
# I <: Base.Generator

function TV.as(d::ProductMeasure{F,S,I}) where {F,S,I<:Base.Generator}
    d1 = marginals(d).f(first(marginals(d).iter))
    as(Array, as(d1), size(marginals(d))...)
end

function TV.as(d::ProductMeasure{Returns{T},F,A}) where {T, F, A <: AbstractArray}
    as(Array, as(d.f.f.value), size(d.pars))
end

function Base.rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure, d1::Dists.Distribution) where {T}
    mar = marginals(d)

    # Distributions doens't (yet) have the three-argument form
    elT = typeof(rand(rng, first(mar)))

    sz = size(mar)
    x = Array{elT, length(sz)}(undef, sz)
    for (j,parj) in enumerate(d.pars)
        x[j] = rand(rng, d.f(parj))
    end
    x
end

# e.g. set(Normal(μ=2)^5, params, randn(5))
function Accessors.set(d::ProductMeasure{F,S,A}, ::typeof(params), p::AbstractArray) where {F,S,A<:AbstractArray}
    set.(marginals(d), params, p)
end

function Accessors.set(d::ProductMeasure{F,S,A}, ::typeof(params), p) where {F,S,A<:AbstractArray}
    par = typeof(d.pars[1])(p)
    ProductMeasure(d.f, Fill(par, size(d.pars)))
end


text:
using RowEchelon
n = 2
m = 6
M = [1.0000 1.5868  2.2477  2.7603  3.6690  5.2387;
     1.5868 2.7603  3.6690  5.1073  6.5115  8.8245;
     2.2477 3.6690  5.2387  6.5115  8.8245 12.7072;
     2.7603 5.1073  6.5115  9.8013 12.1965 15.9960;
     3.6690 6.5115  8.8245 12.1965 15.9960 22.1084;
     5.2387 8.8245 12.7072 15.9960 22.1084 32.1036] + 1.44e-5 * eye(m)

#d, U = eigs(M, nev=3)
#d, U = eigs(M)
#function getU(U)
#    V = [U[:,1]*sqrt(d[1]) U[:,2]*sqrt(d[2]) U[:,3]*sqrt(d[3])]
#    Matrix{Int}(round.(rref(V')'))
#end
r = rank(M, 1e-4)
V = chol(M)[1:r,:]
U = rref(V)'
Ns = [[ 0 1 0;
       -2 3 0;
       -4 2 2],
      [ 0 0 1;
       -4 2 2;
       -6 0 5]]
#N = 0.6906Ns[1] + 0.3091Ns[2]
for i in 1:10
    λ = rand()
    N = λ * Ns[1] + (1 - λ) * Ns[2]
    Z = schurfact(N)[:Z]
    x = [Vector{Float64}(n) for j in 1:r]
    for j in 1:r
        qj = Z[:, j]
        for i in 1:n
            x[j][i] = dot(qj, Ns[i] * qj)
        end
    end
    @show x
end


text:
# import JuliaRL

# import ArcadeLearningEnvironment
# const ALE = ArcadeLearningEnvironment

# """
#     Atari

# An interface adapted from https://github.com/JuliaML/AtariAlgos.jl/blob/master/src/AtariAlgos.jl with a backend implemented by https://github.com/JuliaReinforcementLearning/ArcadeLearningEnvironment.jl
# """
# mutable struct Atari <: JuliaRL.AbstractEnvironment
#     ale::ALE.ALEPtr
#     lives::Int
#     died::Bool
#     reward::Float64
#     score::Float64
#     frame_skip::Int64
#     nframes::Int
#     width::Int
#     height::Int
#     rawscreen::Vector{Cuchar}  # raw screen data from the most recent frame
#     state::Vector{Float64}  # the game state... raw screen data converted to Float64
#     screen::Matrix{RGB{Float64}}

#     function Atari(gamename::AbstractString; frame_skip=4)
#         @assert frame_skip >= 1
#         ale = ALE.ALE_new()
#         ALE.loadROM(ale, gamename)
#         w = ALE.getScreenWidth(ale)
#         h = ALE.getScreenHeight(ale)
#         rawscreen = Array{Cuchar}(undef, w * h * 3)
#         state = similar(rawscreen, Float64)
#         screen = fill(RGB{Float64}(0,0,0), h, w)
#         new(ale, 0, false, 0., 0., frame_skip, 0, w, h, rawscreen, state, screen)
#     end
# end

# function Base.close(game::Atari)
#     game.state = Closed
#     ALE.ALE_del(game.ale)
# end

# JuliaRL.get_actions(env::Atari) = ALE.getLegalActionSet(env.ale)
# valid_action(env::Atari, action) = action in JuliaRL.get_actions(env)

# function update_screen(game::Atari)
#     idx = 1
#     for i in 1:game.height, j in 1:game.width
#         game.screen[i,j] = RGB{Float64}(game.state[idx], game.state[idx+1], game.state[idx+2])
#         idx += 3
#     end
#     game.screen
# end

# @recipe function f(game::Atari)
#     ticks := nothing
#     foreground_color_border := nothing
#     grid := false
#     legend := false
#     aspect_ratio := 1

#     # convert to Image
#     update_screen(game)
# end

# function update_state(game::Atari)
#     # get the raw screen data
#     ALE.getScreenRGB!(game.ale, game.rawscreen)
#     for i in eachindex(game.rawscreen)
#         game.state[i] = game.rawscreen[i] / 256
#     end
#     game.lives = ALE.lives(game.ale)
#     game.state
# end

# function set_seed!(env::Atari, seed=0)
#     ALE.setInt(env.ale, seed)
# end

# # Set seed default to 0
# function JuliaRL.reset!(env::Atari; seed::Int64=0, kwargs...)
#     ALE.reset_game(env.ale)
#     ALE.setInt(env.ale, seed)
#     env.lives = 0
#     env.died = false
#     env.reward = 0
#     env.score = 0
#     env.nframes = 0
#     update_state(env)
# end

# function JuliaRL.environment_step!(env::Atari,
#                                    action;
#                                    rng=Random.GLOBAL_RNG, kwargs...)
#     # act and get the reward and new state
#     for 1:env.frame_skip
#         game.reward = ALE.act(game.ale, a)
#         game.score += game.reward
#         if JuliaRL.is_terminal(env)
#             break
#         end
#     end
#     update_state(game)
#     return
# end

# JuliaRL.get_reward(env::Atari) = env.reward
# JuliaRL.is_terminal(env::Atari) = ALE.game_over(env.ale)
# JuliaRL.get_state(env::Atari) = game.state


text:
using Effects
using Test

@testset "Effects.jl" begin

@testset "List" begin
  include("list.jl")
end

end


text:
data = CSV.read(abspath(joinpath(@__DIR__, "..", "data", "word_freq_probs.csv")),DataFrame)
const words, freqs, depth_score, safety_score  = data.word, data.frequency, data.depth_score, data.safety_score
const word_freq_dict = Dict(zip(words,freqs))


text:
struct SpectralDense{W,B,S,F}
    weight::W
    bias::B
    spectralnormalizer::S
    act::F
end

function SpectralDense(
    in::Integer,
    out::Integer,
    σ = identity;
    bias = true,
    iterations::Integer = 1,
    norm_multiplier = 0.95,
    init = glorot_uniform,
    device = cpu,
)
    weight = init(out, in)
    b = Flux.create_bias(weight, bias, size(weight, 1))
    norm_multiplier = convert(Float32, norm_multiplier)

    kernel_shape = size(weight)
    spectralnormalizer = SpectralNormalizer(
        kernel_shape,
        iterations = iterations,
        norm_multiplier = norm_multiplier,
        init = init,
        device = device,
    )

    return SpectralDense(weight, b, spectralnormalizer, σ)
end

@functor SpectralDense

function (sd::SpectralDense)(x; eps = 1e-12)
    # Spectral regularize the weight
    bounded_kernel = sd.spectralnormalizer(sd.weight)
    logits = bounded_kernel * x .+ sd.bias
    return sd.act.(logits)
end


text:
#========================================================================================#
#	Laboratory 3
#
# Ranges, Arrays and Vectors
#
# Author: Niall Palfreyman, 10/01/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we look at a very important structure in Julia: arrays. Once you
		understand how arrays work, and how amazingly useful they are, you will (almost)
		never need to write a for-loop again!

		... But first we start with Ranges. A Range is an interval of values between start
		and stop boundaries. Here is a range of values between 1 and 5:

			r = 1:5

		What is the type of r?
		""",
		"typeof(1:5)",
		(x -> x <: UnitRange)
	),
	Activity(
		"""
		We can collect() the individual values of a Range into a collection:

			collect(1:5)

		What is the type of this data structure?
		""",
		"typeof(collect(1:5))",
		(x -> x <: Vector)
	),
	Activity(
		"""
		We can also construct ranges for other types, and we can also specify our own
		steplength between the values:

			0.0:0.2:1.0

		Use collect() to list the values in this range and find out how many different
		values the range contains:
		""",
		"length(collect(0.0:0.2:1.0))",
		x -> x==6
	),
	Activity(
		"""
		collect() turns a Range into an Array. Arrays are the daily workhorse of data science:
		they collect elements into tables of 1 or more dimensions. A 1-dimensional Array is
		called a Vector. Here is a COLUMN vector:

			v = [5,4,7,6,3]

		What is its size()?
		""",
		"size(v)",
		x -> x==(5,)
	),
	Activity(
		"""
		What is the size of this ROW vector: w = [5 4 7 6 3 2 1] ?
		""",
		"size(w)",
		x -> x==(1,7)
	),
	Activity(
		"""
		As you can see, the FIRST dimension of an array runs down its columns, and
		the SECOND dimension runs across its rows.

		If you want to write a really super-duper-speedy program, you can create arrays
		using the values that are accidentally hanging around in memory:

			A = Array{Int}(undef,2,3)

		However, more usually we want to INITIALISE the elements of a new array. What
		is the value of the elements in this array:

			A = ones(2,3)
		""",
		"A[1]",
		x -> x == 1.0
	),
	Activity(
		"""
		What about this array: zeros(3,4). What is its size?
		""",
		"size(zeros(3,4))",
		x -> x==(3,4)
	),
	Activity(
		"""
		A MATRIX is an Array with just two dimensions. We often want to fill a matrix
		with a particular value. For example:

			B = Matrix{Int}(undef,3,2)
			fill!( B, 42)

		What is the size of this array?
		""",
		"size(B)",
		x -> x==(3,2)
	),
	Activity(
		"""
		We can create an array literal by explicitly naming the array's elements
		in square brackets:

			C = [2 4 7;3 9 1.2]

		Which value determines the type of the array elements here?
		""",
		"The element from the least restrictive number set",
		x -> x==1.2
	),
	Activity(
		"""
		We can also use square brackets to concatenate arrays, for example:

			[ones(2,3) zeros(2,4)]

		Find the symbol for using square brackets to concatenate two arrays vertically:
		""",
		"",
		x -> occursin( ";", x)
	),
	Activity(
		"""
		A very powerful way of constructing arrays is by COMPREHENSION. That is, we
		write a formula in square brackets that specifies (comprehends) the entire
		contents of the array. For example, here is a vector of square values:

			[x^2 for x in 1:12]

		What do you get if you add the condition "if isodd(x)" to the end of this
		comprehension?
		""",
		"[x^2 for x in 1:12 if isodd(x)]",
		x -> x == [x^2 for x in 1:12 if isodd(x)]
	),
	Activity(
		"""
		We can use several comprehension variables:

			[x*y for x in 1:5, y in 1:5]

		Use array comprehension to generate a matrix of values for the function
		exp(-(x^2 + y^2)) with x and y in the range from -1 to +1 in steps of 0.2. This
		function is a 'hill' with a peak value of 1.0 at the centre where x = y = 0.
		""",
		"[exp(-(x^2 + y^2)) for x in -1:0.2:1, y in -1:0.2:1]",
		x -> abs(sum(x) - 61.02) < 0.01
	),
	Activity(
		"""
		Once we have created an array, we usually want to inspect it. First create
		this array:

			my_matrix = ((1:3)*(1:4)')//3

		Now use eltype() to find the type of the elements in the array my_matrix:
		""",
		"eltype(my_matrix)",
		x -> x <: Rational
	),
	Activity(
		"""
		The length() of a matrix is the number of elements in it. What is the length
		of my_matrix?
		""",
		"length(my_matrix)",
		x -> x==12
	),
	Activity(
		"""
		What is the ndims() of my_matrix?
		""",
		"ndims(my_matrix)",
		x -> x==2
	),
	Activity(
		"""
		You have already used size() to return the size of a matrix, but you can also
		add a second argument to specify which dimension of the matrix you wish to know
		the size of. What is the size of the second dimension of my_matrix?
		""",
		"size(my_matrix,2)",
		x -> x==4
	),
	Activity(
		"""
		Matrix-based languages like Julia offer very many different ways of inspecting
		individual parts of an array. To experiment with this various ways, first define
		the following two arrays:

			my_vector = collect(1:5)
			my_matrix = [1 2 3;4 5 6;7 8 9]

		We can access individual elements of these arrays easily using square brackets
		and subscripts. How do we use subscripts to view the element 8 in my_matrix?
		""",
		"my_matrix[r,c]",
		x -> occursin("my_matrix[3,2]",replace(x," "=>""))
	),
	Activity(
		"""
		The keyword 'end' denotes the last subscript in an array dimension. For example,
		what value is denoted by my_matrix[1,end-1]?
		""",
		"my_matrix[2,end-1]",
		x -> x==2
	),
	Activity(
		"""
		Use subscripting with '=' to change the first entry in my_vector from 1 to 0, then
		give me the new vector that results from this change:
		""",
		"reply(my_vector)",
		x -> x==[0,2,3,4,5]
	),
	Activity(
		"""
		Often, we wish to work not just with individual array elements, but with a whole
		slice of the array. What is the result of the expression my_vector[2:4]?
		""",
		"my_vector[2:4]",
		x -> x==Main.my_vector[2:4]
	),
	Activity(
		"""
		The colon also stands for a complete row or column of an array. What is the
		result of the expression my_matrix[:,3]?
		""",
		"my_matrix[:,3]",
		x -> x==Main.my_matrix[:,3]
	),
	Activity(
		"""
		What is the result of the expression my_matrix[begin+1:end,end]?
		""",
		"my_matrix[begin+1:end,end]",
		x -> x==Main.my_matrix[begin+1:end,end]
	),
	Activity(
		"""
		We can also use slices to change the values in an array. What is the new
		value of the array my_matrix if we perform the following operation:

			my_matrix[3,:] = [17,18,19]
		""",
		"reply(my_matrix)",
		x -> x==[1 2 3;4 5 6;17 18 19]
	),
	Activity(
		"""
		What is the value of the array new_matrix if we construct it like this:

			new_matrix = reshape(my_matrix,1,9)
		""",
		"reshape(my_matrix,1,9)",
		x -> x==reshape(Main.my_matrix,1,9)
	),
	Activity(
		"""
		Notice the strange order of elements in new_matrix. This order arises
		because Julia uses "column-major" ordering of matrix elements: that is,
		these elements have a linear ordering that runs first down the columns
		and then across the rows. You can see this if you list all elements of
		my_matrix using ':'. Show me the result of doing this:
		""",
		"my_matrix[:]",
		x -> x==Main.my_matrix[:]
	),
	Activity(
		"""
		We can perform any function on all elements of an array by using the
		broadcast operator "." like this:

			log.(my_matrix)

		What is the sin() of all elements of my_matrix?
		""",
		"sin.(my_matrix)",
		x -> x==sin.(Main.my_matrix)
	),
	Activity(
		"""
		We can also use broadcasting with infix operators, for example:

			my_matrix .+ 5

		What is the result of taking the reciprocal (1/x) of each element
		in my_matrix?
		""",
		"1 ./ my_matrix",
		x -> x == 1 ./ Main.my_matrix
	),
	Activity(
		"""
		Another way of operating on all elements in an array is to map() the
		operation over the array:

			map(log,my_matrix)

		What is the sin() of my_matrix?
		""",
		"map(sin,my_matrix)",
		x -> x==map(sin,Main.my_matrix)
	),
	Activity(
		"""
		The map() function is particularly useful when we are using anonymous
		functions:

			map(x->5x,my_matrix)

		What do we obtain if we apply the function (x->5sin(x+3)) to my_matrix?
		""",
		"map(x->5sin(x+3),my_matrix)",
		x -> x==map(x->5sin(x+3),Main.my_matrix)
	),
	Activity(
		"""
		We can combine mapping with slicing:

			map(x->5x,my_matrix[:,3])

		What is the result of applying the function (x->5sin(x+3)) to the central
		column of my_matrix?
		""",
		"map(x->5sin(x+3),my_matrix[:,2])",
		x -> x==map(x->5sin(x+3),Main.my_matrix[:,2])
	),
	Activity(
		"""
		We can also iterate (i.e.: loop) over the elements in any array:

			for m in my_array
				print( m, ", ")
			end

		Write a loop to calculate the sum of all elements in my_matrix:
		""",
		"reply(ans)",
		x -> x==sum(Main.my_matrix)
	),
]


text:
using Documenter # , NamedTupleTools

makedocs(
    # modules = [NamedTupleTools],
    sitename = "NamedTupleTools",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages    = Any[
        "Overview"   => "index.md",
        "Concepts"   => "concepts.md",
        "NamedTuples" => "namedtuples.md",
        "Operations" => "operations.md",
        "Functions" =>  "functions.md",
        "Index"  => "generalindex.md"
        ]
    )

deploydocs(
    repo = "github.com/JeffreySarnoff/NamedTupleTools.jl.git",
    target = "build"
)


text:
# Probably this should not be used in the periodic system
function init_pot_Hps_HGH_G( Gvec::GvectorsT, center )

    Ns = Gvec.Ns
    Ngvec = prod(Ns)
    Npoints = Ngvec
    Ω = Gvec.Ω
    G2 = Gvec.G2

    Vpot = zeros( Float64,Npoints )

    const zval = 1
    const rloc = 0.2
    const c1 = -4.0663326
    const c2 = 0.6678322

    Vg = zeros(Complex128,Npoints)
    pre1 = -4*pi*zval/Ω
    pre2 = sqrt(8*pi^3)*rloc^3/Ω
    #
    for ig=2:Ngvec
        Gr = sqrt(G2[ig])*rloc
        expGr2 = exp(-0.5*Gr^2)
        Vg[ig] = pre1/G2[ig]*expGr2 + pre2*expGr2 * (c1 + c2*(3-Gr^2))
    end
    # limiting value, with minus sign ?
    Vg[1] = 2*pi*zval*rloc^2 + (2*pi)^1.5 * rloc^3 * (c1 + 3.0*c2)

    V_ionic = real( G_to_R(Ns,Vg) ) * Npoints

    return V_ionic
end

function init_pot_Hps( LF, center )
    Npoints = LF.Nx * LF.Ny * LF.Nz
    Vpot = zeros( Float64,Npoints )

    const rc1 = 0.25
    const rc2 = 0.284
    const aa  = -1.9287
    const bb  = 0.3374

    # TODO Add journal reference for this pseudopotential
    for ip=1:Npoints
        r = norm( LF.lingrid[:,ip] - center[:] )
        Vpot[ip] = -1.0/r * erf( r/rc1 ) + (aa + bb*r^2)*exp(-(r/rc2)^2)
    end
    return Vpot
end

# Probably this should not be used in the periodic system
function init_pot_Hps_HGH( LF, center )
    Npoints = LF.Nx * LF.Ny * LF.Nz
    Vpot = zeros( Float64,Npoints )

    const Zion = 1
    const rloc = 0.2
    const C1 = -4.0663326
    const C2 = 0.6678322

    # TODO Add journal reference
    for ip = 1:Npoints
        r = norm( LF.lingrid[:,ip] - center[:] )
        if r < 1e-7
            println("WARNING: small r:", r);
            println("Using limiting value:")
            Vpot[ip] = -2*Zion/(sqrt(2*pi)*rloc) + C1
        else
            rrloc = r/rloc
            Vpot[ip] = -Zion/r * erf( r/(sqrt(2.0)*rloc) ) +
                     (C1 + C2*rrloc^2)*exp(-0.5*(rrloc)^2)
        end
    end
    return Vpot
end


text:
using DataFrames, CSV, Dates, Statistics, CategoricalArrays
data = dropmissing(CSV.read("data/purchaseData.csv", DataFrame))

data[!,:Date] = Date.(data[!,:Date], "d/m/y")
println(first(sort(data, :Date), 3),"\n")

println(first(filter(row -> row[:Price] > 50000, data),3 ),"\n")

#categorical!(data, :Grade)
data[!, :Grade] = categorical(data[!, :Grade])
println(first(data, 3), "\n")

println(
#    by(data, :Grade, :Price =>
#        x -> ( NumSold=length(x), AvgPrice=mean(x)) )
    combine(
        groupby(data, :Grade),
        :Price=> (x -> [length(x) mean(x)]) => [:NumSold, :AvgPrice]
        )
)

max_stars_repo_name: newptcai/StatsWithJuliaBook

max_issues_repo_name: newptcai/StatsWithJuliaBook

max_forks_repo_name: newptcai/StatsWithJuliaBook


text:
include("test_matrixcore.jl")
include("test_fresnelltools.jl")


text:
#######################################################################
# Coverage.jl
# Input: Code coverage and memory allocations
# Output: Useful things
# https://github.com/IainNZ/Coverage.jl
#######################################################################
module Coverage

    # The code coverage results produced by Julia itself report some
    # lines as "null" (cannot be run), when they could have been run
    # but were not (should be 0). We use JuliaParser to augment the
    # coverage results by identifying this code.

    using Compat

    export process_folder, process_file
    export clean_folder, clean_file
    export process_cov, amend_coverage_from_src!
    export get_summary
    export analyze_malloc, merge_coverage_counts

    # The unit for line counts. Counts can be >= 0 or nothing, where
    # the nothing means it doesn't make sense to have a count for this
    # line (e.g. a comment), but 0 means it could have run but didn't.
    const CovCount = Union{Void,Int}

    export FileCoverage
    """
    FileCoverage

    Represents coverage info about a file, including the filename, the source
    code itself, and a `Vector` of run counts for each line. If the
    line was expected to be run the count will be an `Int` >= 0. Other lines
    such as comments will have a count of `nothing`.
    """
    type FileCoverage
        filename::AbstractString
        source::AbstractString
        coverage::Vector{CovCount}
    end

    """
        get_summary(fcs)

    Summarize results from a single `FileCoverage` instance or a `Vector` of
    them, returning a 2-tuple with the covered lines and total lines.
    """
    function get_summary end

    function get_summary(fc::FileCoverage)
        if length(fc.coverage) > 0
            cov_lines = sum(x -> x!=nothing && x > 0, fc.coverage)
            tot_lines = sum(x -> x!=nothing, fc.coverage)
            return cov_lines, tot_lines
        else
            return 0, 0
        end
    end
    function get_summary(fcs::Vector{FileCoverage})
        cov_lines, tot_lines = 0, 0
        for fc in fcs
            c, t = get_summary(fc)
            cov_lines += c
            tot_lines += t
        end
        return cov_lines, tot_lines
    end

    """
        merge_coverage_counts(a1::Vector{CovCount}, a2::Vector{CovCount})

    Given two vectors of line coverage counts, take the pairwise
    maximum of both vectors, preseving null counts if both are null.
    """
    function merge_coverage_counts(a1::Vector{CovCount},
                                   a2::Vector{CovCount})
        n = max(length(a1),length(a2))
        a = Vector{CovCount}(n)
        for i in 1:n
            a1v = isassigned(a1, i) ? a1[i] : nothing
            a2v = isassigned(a2, i) ? a2[i] : nothing
            a[i] = a1v == nothing ? a2v :
                   a2v == nothing ? a1v : max(a1v, a2v)
        end
        return a
    end

    """
        process_cov(filename, folder)

    Given a filename for a Julia source file, produce an array of
    line coverage counts by reading in all matching .{pid}.cov files.
    On Julia 0.3 there was just a .cov file, but this code works fine.
    """
    function process_cov(filename, folder)
        # Find all coverage files in the folder that match the file we
        # are currently working on
        files = readdir(folder)
        files = map!(file -> joinpath(folder,file), files, files)
        filter!(file -> contains(file,filename) && contains(file,".cov"), files)
        # If there are no coverage files...
        if isempty(files)
            # ... we will assume that, as there is a .jl file, it was
            # just never run. We'll report the coverage as all null.
            println( """Coverage.process_cov: Coverage file(s) for $filename do not exist.
                                              Assuming file has no coverage.""")
            lines = open(filename) do fp
                readlines(fp)
            end
            coverage = Vector{CovCount}(length(lines))
            return fill!(coverage, nothing)
        end
        # Keep track of the combined coverage
        full_coverage = Vector{CovCount}(0)
        for file in files
            lines = open(file, "r") do fp
                readlines(fp)
            end
            num_lines = length(lines)
            coverage = Vector{CovCount}(num_lines)
            for i in 1:num_lines
                # Columns 1:9 contain the coverage count
                cov_segment = lines[i][1:9]
                # If coverage is NA, there will be a dash
                coverage[i] = cov_segment[9] == '-' ? nothing : parse(Int, cov_segment)
            end
            full_coverage = merge_coverage_counts(full_coverage, coverage)
        end
        return full_coverage
    end

    """
        amend_coverage_from_src!(coverage::Vector{CovCount}, srcname)

    The code coverage functionality in Julia can miss code lines, which
    will be incorrectly recorded as `nothing` but should instead be 0
    This function takes a coverage count vector and a the filename for
    a Julia code file, and updates the coverage vector in place.
    """
    function amend_coverage_from_src!(coverage::Vector{CovCount}, srcname)
        # To make sure things stay in sync, parse the file position
        # corresonding to each new line
        linepos = Int[]
        open(srcname) do io
            while !eof(io)
                push!(linepos, position(io))
                readline(io)
            end
            push!(linepos, position(io))
        end
        open(srcname) do io
            while !eof(io)
                pos = position(io)
                linestart = minimum(searchsorted(linepos, pos))
                ast = Base.parse(io)
                isa(ast, Expr) || continue
                flines = function_body_lines(ast)
                if !isempty(flines)
                    flines += linestart-1
                    for l in flines
                        if l > length(coverage)
                            error("source file is longer than .cov file; source might have changed")
                        end
                        if coverage[l] == nothing
                            coverage[l] = 0
                        end
                    end
                end
            end
        end
        nothing
    end
    # function_body_lines is located in parser.jl
    include("parser.jl")

    """
        process_file(filename[, folder]) -> FileCoverage

    Given a .jl file and its containing folder, produce a corresponding
    `FileCoverage` instance from the source and matching coverage files. If the
    folder is not given it is extracted from the filename.
    """
    function process_file end

    function process_file(filename, folder)
        println("Coverage.process_file: Detecting coverage for $filename")
        coverage = process_cov(filename,folder)
        amend_coverage_from_src!(coverage, filename)
        return FileCoverage(filename, readstring(filename), coverage)
    end
    process_file(filename) = process_file(filename,splitdir(filename)[1])

    """
        process_folder(folder="src") -> Vector{FileCoverage}

    Process the contents of a folder of Julia source code to collect coverage
    statistics for all the files contained within. Will recursively traverse
    child folders. Default folder is "src", which is useful for the primary case
    where Coverage is called from the root directory of a package.
    """
    function process_folder(folder="src")
        println("""Coverage.process_folder: Searching $folder for .jl files...""")
        source_files = FileCoverage[]
        files = readdir(folder)
        for file in files
            fullfile = joinpath(folder,file)
            if isfile(fullfile)
                # Is it a Julia file?
                if splitext(fullfile)[2] == ".jl"
                    push!(source_files, process_file(fullfile,folder))
                else
                    println("Coverage.process_folder: Skipping $file, not a .jl file")
                end
            else isdir(fullfile)
                # If it is a folder, recursively traverse
                append!(source_files, process_folder(fullfile))
            end
        end
        return source_files
    end

    # matches julia coverage files with and without the PID
    iscovfile(filename) = ismatch(r"\.jl\.?[0-9]*\.cov$", filename)
    # matches a coverage file for the given sourcefile. They can be full paths
    # with directories, but the directories must match
    function iscovfile(filename, sourcefile)
        startswith(filename, sourcefile) || return false
        ismatch(r"\.jl\.?[0-9]*\.cov$", filename)
    end

    """
        clean_folder(folder::AbstractString)

    Cleans up all the `.cov` files in the given directory and subdirectories.
    Unlike `process_folder` this does not include a default value
    for the root folder, requiring the calling code to be more explicit about
    which files will be deleted.
    """
    function clean_folder(folder::AbstractString)
        files = readdir(folder)
        for file in files
            fullfile = joinpath(folder, file)
            if isfile(fullfile) && iscovfile(file)
                # we have ourselves a coverage file. eliminate it
                println("Removing $fullfile")
                rm(fullfile)
            elseif isdir(fullfile)
                clean_folder(fullfile)
            end
        end
        nothing
    end

    """
        clean_file(filename::AbstractString)

    Cleans up all `.cov` files associated with a given source file. This only
    looks in the directory of the given file, i.e. the `.cov` files should be
    siblings of the source file.
    """
    function clean_file(filename::AbstractString)
        folder = splitdir(filename)[1]
        files = readdir(folder)
        for file in files
            fullfile = joinpath(folder, file)
            if isfile(fullfile) && iscovfile(fullfile, filename)
                println("Removing $fullfile")
                rm(fullfile)
            end
        end
    end

    include("coveralls.jl")
    include("codecovio.jl")
    include("lcov.jl")
    include("memalloc.jl")
end


text:
@testset "nash_conv" begin
    env = OpenSpielEnv("kuhn_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = TabularRandomPolicy()

    @test RLZoo.nash_conv(p, env) == 11/ 12

    p = get_optimal_kuhn_policy()
    @test RLZoo.nash_conv(p, env) == 0.

    env = OpenSpielEnv("leduc_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 4.747222222222222

    env = OpenSpielEnv("kuhn_poker(players=3)";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 2.0624999999999996

    env = OpenSpielEnv("kuhn_poker(players=4)";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 3.4760416666666663
end


text:
# REAL = double = Cdouble

mutable struct TriangulateIO
  pointlist::Ptr{Cdouble}
  pointattributelist::Ptr{Cdouble}
  pointmarkerlist::Ptr{Cint}
  numberofpoints::Cint
  numberofpointattributes::Cint
  trianglelist::Ptr{Cint}
  triangleattributelist::Ptr{Cdouble}
  trianglearealist::Ptr{Cdouble}
  neighborlist::Ptr{Cint}
  numberoftriangles::Cint
  numberofcorners::Cint
  numberoftriangleattributes::Cint
  segmentlist::Ptr{Cint}
  segmentmarkerlist::Ptr{Cint}
  numberofsegments::Cint
  holelist::Ptr{Cdouble}
  numberofholes::Cint
  regionlist::Ptr{Cdouble}
  numberofregions::Cint
  edgelist::Ptr{Cint}
  edgemarkerlist::Ptr{Cint}
  normlist::Ptr{Cdouble}
  numberofedges::Cint
  TriangulateIO() = new(C_NULL, C_NULL, C_NULL, 0, 0, C_NULL, C_NULL, C_NULL, C_NULL, 0, 0, 0, C_NULL, C_NULL, 0, C_NULL, 0, C_NULL, 0, C_NULL, C_NULL, C_NULL, 0)
end


text:
# Check values of the new multimode polaron codes

@testset verbose = true "MultiplePhonons" begin

# Physical constants
ħ = 1.05457162825e-34;          # kg m2 / s
eV = 1.602176487e-19;           # kg m2 / s2
me = 9.10938188e-31;            # kg
kB =  1.3806504e-23;            # kg m2 / K s2
ε_0 = 8.854E-12                 # Units: C2N−1m−2, permittivity of free space
c = 3e8

MAPI= [
# 96.20813558773261 0.4996300522819191
# 93.13630357703363 1.7139631746083817
# 92.87834578121567 0.60108592692181
# 92.4847918585963 0.0058228799414729
# 92.26701437594754 0.100590086574602
# 89.43972834606603 0.006278895133832249
# 46.89209141511332 0.2460894564364346
# 46.420949316788 0.14174282581124137
# 44.0380222871706 0.1987196948553428
# 42.89702947649343 0.011159939465770681
# 42.67180170168193 0.02557751102757614
# 41.46971205834201 0.012555230726601503
# 37.08982543385215 0.00107488277468418
# 36.53555265689563 0.02126940080871224
# 30.20608114002676 0.009019481779712388
# 27.374810898415028 0.03994453721421388
# 26.363055017011728 0.05011922682554448
# 9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05   # Acoustic modes!
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]

volume = (6.29e-10)^3
ϵ_optic = 4.5
m_eff = 0.12
phonon_freq = MAPI[:, 1]
ir_activity = MAPI[:, 2]
ϵ_static = 24.1

# Ionic dielectric and decomposed alpha

ω = 2π .* phonon_freq
ϵ_ionic = [ϵ_ionic_mode(i, j, volume) for (i, j) in zip(phonon_freq, ir_activity)]

α = [multi_frohlichalpha(ϵ_optic, i, sum(ϵ_ionic), j, m_eff) for (i, j) in zip(ϵ_ionic, phonon_freq)]

@testset "Ionic dielectric and decomposed alphas" begin

    @test ϵ_ionic ≈ [0.2999680470756664, 0.0247387647244569, 0.2543132184018061, 0.16621617310133838, 2.3083204422506296, 3.0707979601813267, 3.2026782087486407, 0.892674135624958, 0.8861579096771846, 0.7209278375756829, 0.40199759805819046, 0.6183279038315278, 0.7666391525823296, 0.1555444994147378, 1.1627710200840813] rtol = 1e-3

    @test α ≈ [0.03401013445306177, 0.002850969846883158, 0.03075081562607006, 0.02275292381052607, 0.33591418423943553, 0.46526818717696034, 0.5046331098089347, 0.14223560721522646, 0.16083871312929882, 0.1622911897190622, 0.09123913334086006, 0.14070972715961402, 0.18159786148348117, 0.039500396977008016, 0.3487735470060312] rtol = 1e-3

    @test sum(α) ≈ 2.663366500992453 rtol = 1e-3

    println('\n', "Ionic dielectric and alphas:")
    println("-------------------------------")
    println("  ω (THz)    ϵ_ionic     αs")
    println("-------------------------------")
    display(hcat(ω, ϵ_ionic, α))
    println("\n-------------------------------")
end

# Variations

v_0, w_0 = var_params(α; v = 0.0, w = 0.0, ω = ω) # Athermal

β = [i .* ħ / (kB * 300) * 1e12 for i in ω]
v, w = var_params(α, β; v = 0.0, w = 0.0, ω = ω) # Thermal

@testset "Multiple mode variations" begin

    @test v_0 ≈ 3.292283619446986 rtol = 1e-3
    @test w_0 ≈ 2.679188425097246 rtol = 1e-3

    @test v ≈ 35.19211042393129 rtol = 1e-3
    @test w ≈ 32.454157668863225 rtol = 1e-3

    println("\nVariational parameters:")
    println("Athermal v = $(v_0[1]) | athermal w = $(w_0[1])")
    println("300K v = $(v[1]) | 300K w = $(w[1])")
end

# Energies

E = multi_F(v_0, w_0, α; ω = ω) * 1000 * ħ / eV * 1e12 # Enthalpy

F = multi_F(v, w, α, β; ω = ω) * 1000 * ħ / eV * 1e12 # Free energy

@testset "Multiple mode energies" begin

    @test E ≈ -19.50612170650821 rtol = 1e-3

    @test F ≈ -42.79764110613318 rtol = 1e-3

    println("\nPolaron energies")
    println("0K energy: $E meV")
    println("300K energy: $F meV")
end

# Mobility

μ = polaron_mobility(β, α, v, w; ω = ω) * eV / (1e12 * me * m_eff) * 100^2

@testset "Multiple mode mobility" begin

    @test μ ≈ 160.50844330430286 rtol = 1e-3

    println("\nMobility at 300K: $μ")
end

# Impedence

Z_dc = polaron_complex_impedence(0, β, α, v, w; ω = ω) / sqrt(ε_0 * ϵ_optic) / c # DC limit

Z = polaron_complex_impedence(50, β, α, v, w; ω = ω) / sqrt(ε_0 * ϵ_optic) / c # AC limit

@testset "Multiple mode impedence" begin

    @test Z_dc ≈ 0.04822198080982958 + 0.0im rtol = 1e-3

    @test Z ≈ 0.015597562999157006 - 0.15169951660903833im rtol = 1e-3

    println("\nComplex Impedence:")
    println("DC limit: $Z_dc")
    println("At 50 THz: $Z")
end

# Conductivity

σ_dc = polaron_complex_conductivity(0, β, α, v, w; ω = ω) / sqrt(ε_0 * ϵ_optic) / c # DC limit

σ = polaron_complex_conductivity(50, β, α, v, w; ω = ω) / sqrt(ε_0 * ϵ_optic) / c # AC limit

@testset "Multiple mode conductivity" begin

    @test σ_dc ≈ 5.783096153975994e-6 - 0.0im rtol = 1e-3

    @test σ ≈ 1.8703663429355958e-7 + 1.8190897521649982e-6im rtol = 1e-3

    println("\nComplex Conductivity:")
    println("DC limit: $σ_dc")
    println("At 50 THz: $σ\n")
end

# Single Mode Tests

Hellwarth_B_freq = HellwarthBScheme(hcat(phonon_freq, ir_activity))

@testset "Hellwarth effective frequency" begin

    @test Hellwarth_B_freq ≈ 2.25 rtol = 1e-3
end

@testset "Single effective mode MAPI" begin

    singlemode_polaron = make_polaron(ϵ_optic, ϵ_static, Hellwarth_B_freq, m_eff, [0.0, 300.0], [0.0, 3.0]; volume = volume, ir_activity = nothing, rtol = 1e-4, verbose = true, threads = true)

    println('\n', singlemode_polaron)

    @test singlemode_polaron.α ≈ 2.393 rtol = 1e-3
    @test singlemode_polaron.v ≈ [3.308644142915268; 19.847591395925644] rtol = 1e-3
    @test singlemode_polaron.w ≈ [2.6633969095604466; 16.948206590039813] rtol = 1e-3
    @test singlemode_polaron.F ≈ [-23.02903831886734, -35.46521250753788] rtol = 1e-3
    @test singlemode_polaron.Z ≈ [0.0 + 0.0im 0.11385135083260453 + 0.0im; 0.018788147630145906 - 0.0320853372512867im 0.11055356479444363 - 0.0023745543208784164im] rtol = 1e-3
    @test singlemode_polaron.σ ≈ [Inf + 0.0im 8.783382829337691 - 0.0im; 13.590340404155473 + 23.208815675232696im 9.04121795743805 + 0.1941942189449604im] rtol = 1e-3
    @test singlemode_polaron.μ ≈ [Inf, 136.42796295778325] rtol = 1e-3
end

# Multimode Tests

@testset "Multiple mode (15) MAPI" begin

    multimode_polaron = make_polaron(ϵ_optic, ϵ_static, phonon_freq, m_eff, [0.0, 300.0], [0.0, 3.0]; volume = volume, ir_activity = ir_activity, rtol = 1e-4, verbose = true, threads = true)

    println('\n', multimode_polaron)

    @test multimode_polaron.α ≈ [0.03401013445306177, 0.002850969846883158, 0.03075081562607006, 0.02275292381052607, 0.33591418423943553, 0.46526818717696034, 0.5046331098089347, 0.14223560721522646, 0.16083871312929882, 0.1622911897190622, 0.09123913334086006, 0.14070972715961402, 0.18159786148348117, 0.039500396977008016, 0.3487735470060312] rtol = 1e-3
    @test multimode_polaron.v ≈ [3.292288128236545; 35.19197149158517] rtol = 1e-3
    @test multimode_polaron.w ≈ [2.679192692280829; 32.45402163343302] rtol = 1e-3
    @test multimode_polaron.F ≈ [-19.506121706507198, -42.79764112300015] rtol = 1e-3
    @test multimode_polaron.Z ≈ [0.0 + 0.0im 0.09684805201838559 + 0.0im; 0.013916089370749993 - 0.02698207327877506im 0.09373671239869472 - 0.00586253350488456im] rtol = 1e-3
    @test multimode_polaron.σ ≈ [Inf + 0.0im 10.325452904413197 - 0.0im; 15.098451823667599 + 29.274570078468706im 10.626611965355764 + 0.6646154649133272im] rtol = 1e-3
    @test multimode_polaron.μ ≈ [Inf, 160.50598091483337] rtol = 1e-3
end

print('\n')

end


text:
@testset "CUCB" begin
  @testset "Algorithm" begin
    n = 3
    ε = .1
    reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]
    instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

    Random.seed!(1)
    n_rounds = 20
    s = simulate(instance, CUCB(), n_rounds)
    @test s.round == n_rounds
    @test sum(values(s.arm_counts)) == n_rounds * n
    @test s.regret ≈ 3.7
    @test s.reward ≈ 36.0
    @test s.regret + s.reward ≈ 39.7

    for arm in keys(s.arm_counts)
      @test s.arm_average_reward[arm] ≈ s.arm_reward[arm] / s.arm_counts[arm]
    end
  end

  @testset "Trace" begin
    n = 3
    ε = .1
    reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]
    instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

    n_rounds = 20
    test_policy_trace(instance, CUCB(), n_rounds, n_rounds * n, 3.7, 36.0, 39.7)
  end
end


text:
function testGlobalNull(X, Y; bootSamples::Int64=300, estimSupport::Union{Array{BitArray},Nothing}=nothing)
    nx, p = size(X)
    ny    = size(Y, 1)

    Sx = Symmetric( cov(X) )
    Sy = Symmetric( cov(Y) )

    eS = estimSupport === nothing ? __initSupport(Sx, Sy, X, Y) : estimSupport

    rp = div((p+1)*p, 2)
    indS = Vector{Vector{Int64}}(undef, rp)

    Δhat = Vector{Float64}(undef, rp)

    # obtain initial estimate
    it = 0
    for col=1:p
        for row=col:p
            it = it + 1
            indS[it] = getLinearSupport(row, col, eS[it])

            tmp = estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, row, col, indS[it])
            Δhat[it] = tmp.p
        end
    end

    # this is not complete yet....
end


text:
module ExpFamily

using Compat

# @compat abstract type EFam end

@compat abstract type EFamily end
@compat abstract type NatParam  <: EFamily end
@compat abstract type MeanParam <: EFamily end

@compat abstract type SuffStats end

const Float = Float64

export
    EFamily,
    natparam,
    meanparam,
    suffstats,
    project,
    loglik,
    uloglik,
    gradloglik

import Base.+, Base.-, Base.*, Base./, Base.-

include("models/tools.jl")
include("models/gaussian.jl")
include("models/gaussian_ops.jl")
include("models/gaussian_loglik.jl")
include("models/gaussian_suffstats.jl")

end # module


text:
"""
    test_constraint_get_ConstraintIndex(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test getting constraints by name.
"""
function test_constraint_get_ConstraintIndex(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, T(1) * x, MOI.GreaterThan(T(1)))
    MOI.set(model, MOI.ConstraintName(), c1, "c1")
    c2 = MOI.add_constraint(model, T(1) * x, MOI.LessThan(T(2)))
    MOI.set(model, MOI.ConstraintName(), c2, "c2")
    F = MOI.ScalarAffineFunction{T}
    @test MOI.get(model, MOI.ConstraintIndex, "c3") === nothing
    @test MOI.get(model, MOI.ConstraintIndex{F,MOI.LessThan{T}}, "c1") ===
          nothing
    @test MOI.get(model, MOI.ConstraintIndex{F,MOI.GreaterThan{T}}, "c2") ===
          nothing
    c1 = MOI.get(model, MOI.ConstraintIndex{F,MOI.GreaterThan{T}}, "c1")
    @test MOI.get(model, MOI.ConstraintIndex, "c1") == c1
    @test MOI.is_valid(model, c1)
    c2 = MOI.get(model, MOI.ConstraintIndex{F,MOI.LessThan{T}}, "c2")
    @test MOI.get(model, MOI.ConstraintIndex, "c2") == c2
    @test MOI.is_valid(model, c2)
    return
end

"""
    test_constraint_ScalarAffineFunction_LessThan(
        model::MOI.ModelLike,
        config::Config,
    )

Add an ScalarAffineFunction-in-LessThan constraint.
"""
function test_constraint_ScalarAffineFunction_LessThan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    )
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0)),
        MOI.LessThan(T(1)),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ScalarAffineFunction_LessThan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1//2],
        ),
    )
    return
end

"""
    test_constraint_ScalarAffineFunction_GreaterThan(
        model::MOI.ModelLike,
        config::Config,
    )

Add an ScalarAffineFunction-in-GreaterThan constraint.
"""
function test_constraint_ScalarAffineFunction_GreaterThan(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.GreaterThan{T},
    )
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0)),
        MOI.GreaterThan(T(1)),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ScalarAffineFunction_GreaterThan),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}) => T[1//2],
        ),
    )
    return
end

"""
    test_constraint_ScalarAffineFunction_EqualTo(
        model::MOI.ModelLike,
        config::Config,
    )

Add an ScalarAffineFunction-in-EqualTo constraint.
"""
function test_constraint_ScalarAffineFunction_EqualTo(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0)),
        MOI.EqualTo(T(1)),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ScalarAffineFunction_EqualTo),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}) => T[1//2],
        ),
    )
    return
end

"""
    test_constraint_ScalarAffineFunction_Interval(
        model::MOI.ModelLike,
        config::Config,
    )

Add an ScalarAffineFunction-in-Interval constraint.
"""
function test_constraint_ScalarAffineFunction_Interval(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.Interval{T},
    )
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0)),
        MOI.Interval(T(1), T(4)),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(3), x)], T(0)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(6),
        variable_primal = [(x, T(2))],
        constraint_primal = [(c, T(4))],
        constraint_dual = [(c, T(-3 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ScalarAffineFunction_Interval),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.Interval{T}) => T[-3//2],
        ),
    )
    return
end

"""
    test_constraint_ScalarAffineFunction_duplicate(
        model::MOI.ModelLike,
        config::Config,
    )

Add a `ScalarAffineFunction`-in-`LessThan` constraint with duplicate terms in
the function.

Taken from https://github.com/JuliaOpt/MathOptInterfaceMosek.jl/issues/41
"""
function test_constraint_ScalarAffineFunction_duplicate(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    objective_function =
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    MOI.set(
        model,
        MOI.ObjectiveFunction{typeof(objective_function)}(),
        objective_function,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(T(1), [x, x]), T(0))
    c = MOI.add_constraint(model, f, MOI.LessThan(T(1)))
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T(1))],
        constraint_dual = [(c, T(-1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ScalarAffineFunction_duplicate),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}) => T[-1//2],
        ),
    )
    return
end

"""
    test_constraint_VectorAffineFunction_duplicate(
        model::MOI.ModelLike,
        config::Config,
    )

Add a `VectorAffineFunction`-in-`Nonpositives` constraint with duplicate terms
in the function.
"""
function test_constraint_VectorAffineFunction_duplicate(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    objective_function =
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    MOI.set(
        model,
        MOI.ObjectiveFunction{typeof(objective_function)}(),
        objective_function,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(T(1), [x, x])),
        T[-1],
    )
    c = MOI.add_constraint(model, f, MOI.Nonpositives(1))
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2),
        variable_primal = [(x, T(1 // 2))],
        constraint_primal = [(c, T[0])],
        constraint_dual = [(c, T[-1//2])],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_VectorAffineFunction_duplicate),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2]),
            MOI.FEASIBLE_POINT,
            (MOI.VectorAffineFunction{T}, MOI.Nonpositives) => [T[-1//2]],
        ),
    )
    return
end

"""
    test_constraint_qcp_duplicate_diagonal(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test a QCP problem with a duplicate diagonal term.

The problem is `max x + 2y | y + x^2 + x^2 <= 1, x >= 0.5, y >= 0.5`.
"""
function test_constraint_qcp_duplicate_diagonal(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{T},
        MOI.LessThan{T},
    )
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(T[1, 2], x), T(0)),
    )
    vc1 = MOI.add_constraint(model, x[1], MOI.GreaterThan(T(1 // 2)))
    # We test this after the creation of every `VariableIndex` constraint
    # to ensure a good coverage of corner cases.
    @test vc1.value == x[1].value
    vc2 = MOI.add_constraint(model, x[2], MOI.GreaterThan(T(1 // 2)))
    @test vc2.value == x[2].value
    MOI.add_constraint(
        model,
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.(T(2), [x[1], x[1]], [x[1], x[1]]),  # quad
            MOI.ScalarAffineTerm.(T(1), [x[2]]),  # affine terms
            T(0),  # constant
        ),
        MOI.LessThan(T(1)),
    )
    _test_model_solution(
        model,
        config;
        objective_value = T(3 // 2),
        variable_primal = [(x[1], T(1 // 2)), (x[2], T(1 // 2))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_qcp_duplicate_diagonal),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2, 1//2]),
        ),
    )
    return
end

"""
    test_constraint_qcp_duplicate_off_diagonal(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test a QCP problem with a duplicate off-diagonal term.

The problem is
`max x + 2y | x^2 + 0.25y*x + 0.25x*y + 0.5x*y + y^2 <= 1, x >= 0.5, y >= 0.5`.

!!! warn
    This problem has an irrational solution! Solvers using rational arithmetic
    should exclude it.
"""
function test_constraint_qcp_duplicate_off_diagonal(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.ScalarQuadraticFunction{T},
        MOI.LessThan{T},
    )
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(T[1, 2], x), T(0)),
    )
    vc1 = MOI.add_constraint(model, x[1], MOI.GreaterThan{T}(T(1 // 2)))
    @test vc1.value == x[1].value
    vc2 = MOI.add_constraint(model, x[2], MOI.GreaterThan{T}(T(1 // 2)))
    @test vc2.value == x[2].value
    MOI.add_constraint(
        model,
        MOI.ScalarQuadraticFunction(
            MOI.ScalarQuadraticTerm.(
                T[2, 1//4, 1//4, 1//2, 2],
                [x[1], x[1], x[2], x[1], x[2]],
                [x[1], x[2], x[1], x[2], x[2]],
            ),  # quad
            MOI.ScalarAffineTerm{T}[],  # affine terms
            T(0),  # constant
        ),
        MOI.LessThan(T(1)),
    )
    x2_solution = (sqrt(T(13)) - T(1)) / T(4)
    _test_model_solution(
        model,
        config;
        objective_value = T(1 // 2) + T(2) * x2_solution,
        variable_primal = [(x[1], T(1 // 2)), (x[2], x2_solution)],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_qcp_duplicate_off_diagonal),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1//2, (sqrt(T(13))-T(1))/T(4)]),
        ),
    )
    return
end

"""
    test_constraint_ZeroOne_bounds(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test a problem with a bounded ZeroOne variable.
"""
function test_constraint_ZeroOne_bounds(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.ZeroOne())
    MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(2) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        variable_primal = [(x, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ZeroOne_bounds),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
    )
    return
end

"""
    test_constraint_ZeroOne_bounds_2(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test a problem with a ZeroOne and binding fractional upper bound.
"""
function test_constraint_ZeroOne_bounds_2(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.ZeroOne())
    MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    MOI.add_constraint(model, x, MOI.LessThan(T(1 // 2)))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(2) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        variable_primal = [(x, T(0))],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_ZeroOne_bounds_2),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[0])),
    )
    return
end

"""
    test_constraint_ZeroOne_bounds_3(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test a problem with a ZeroOne and infeasible fractional bounds.
"""
function test_constraint_ZeroOne_bounds_3(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @requires _supports(config, MOI.optimize!)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.ZeroOne())
    MOI.add_constraint(model, x, MOI.GreaterThan(T(1 // 5)))
    MOI.add_constraint(model, x, MOI.LessThan(T(1 // 2)))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    f = T(2) * x
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == config.infeasible_status
    return
end

function setup_test(
    ::typeof(test_constraint_ZeroOne_bounds_3),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, MOI.INFEASIBLE),
    )
    return
end

"""
    test_constraint_PrimalStart_DualStart_SecondOrderCone(
        model::MOI.ModelLike,
        config::Config{T},
    ) where {T}

Test combining the [`MOI.VariablePrimalStart`](@ref),
[`MOI.ConstraintPrimalStart`](@ref) and [`MOI.ConstraintDualStart`](@ref)
attributes with a `MOI.VectorAffineFunction{T}`-in-`MOI.SecondOrderCone`.
"""
function test_constraint_PrimalStart_DualStart_SecondOrderCone(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        MOI.SecondOrderCone,
    )
    @requires _supports(config, MOI.optimize!)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOIU.operate(vcat, T, x, T(1)),
        MOI.SecondOrderCone(2),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    if MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
        MOI.set(model, MOI.VariablePrimalStart(), x, T(2))
    end
    if MOI.supports(model, MOI.ConstraintPrimalStart(), typeof(c))
        MOI.set(model, MOI.ConstraintPrimalStart(), c, T[2, 2])
    end
    if MOI.supports(model, MOI.ConstraintDualStart(), typeof(c))
        MOI.set(model, MOI.ConstraintDualStart(), c, T[2, -2])
    end
    _test_model_solution(
        model,
        config;
        objective_value = T(1),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c, T[1, 1])],
        constraint_dual = [(c, T[1, -1])],
    )
    return
end

function setup_test(
    ::typeof(test_constraint_PrimalStart_DualStart_SecondOrderCone),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.VectorAffineFunction{T}, MOI.SecondOrderCone) =>
                [T[1, -1]],
        ),
    )
    return
end

"""
    test_constraint_ConstraintPrimalStart(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test the `ConstraintPrimalStart` attribute.
"""
function test_constraint_ConstraintPrimalStart(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        MOI.Nonnegatives,
    )
    @requires MOI.supports(
        model,
        MOI.ConstraintPrimalStart(),
        MOI.ConstraintIndex{MOI.VectorAffineFunction{T},MOI.Nonnegatives},
    )
    x = MOI.add_variable(model)
    f = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x))],
        T[0],
    )
    c = MOI.add_constraint(model, f, MOI.Nonnegatives(1))
    @test MOI.get(model, MOI.ConstraintPrimalStart(), c) === nothing
    MOI.set(model, MOI.ConstraintPrimalStart(), c, T[-1])
    @test MOI.get(model, MOI.ConstraintPrimalStart(), c) == T[-1]
    MOI.set(model, MOI.ConstraintPrimalStart(), c, nothing)
    @test MOI.get(model, MOI.ConstraintPrimalStart(), c) === nothing
    return
end

"""
    test_constraint_ConstraintDualStart(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test the `ConstraintDualStart` attribute.
"""
function test_constraint_ConstraintDualStart(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        MOI.Nonnegatives,
    )
    @requires MOI.supports(
        model,
        MOI.ConstraintDualStart(),
        MOI.ConstraintIndex{MOI.VectorAffineFunction{T},MOI.Nonnegatives},
    )
    x = MOI.add_variable(model)
    f = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x))],
        T[0],
    )
    c = MOI.add_constraint(model, f, MOI.Nonnegatives(1))
    @test MOI.get(model, MOI.ConstraintDualStart(), c) === nothing
    MOI.set(model, MOI.ConstraintDualStart(), c, T[-1])
    @test MOI.get(model, MOI.ConstraintDualStart(), c) == T[-1]
    MOI.set(model, MOI.ConstraintDualStart(), c, nothing)
    @test MOI.get(model, MOI.ConstraintDualStart(), c) === nothing
    return
end

"""
    test_constraint_Indicator_ConstraintName(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test ConstraintName for indicator sets.
"""
function test_constraint_Indicator_ConstraintName(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires(
        MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.GreaterThan{T}},
        ),
    )
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, x[1], MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), x[2])),
        ],
        T[0, 0],
    )
    s = MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(T(1)))
    c = MOI.add_constraint(model, f, s)
    MOI.set(model, MOI.ConstraintName(), c, "my_indicator")
    @test MOI.get(model, MOI.ConstraintName(), c) == "my_indicator"
    return
end

"""
    test_constraint_Indicator_ACTIVATE_ON_ONE(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test ACTIVATE_ON_ONE for indicator sets.
"""
function test_constraint_Indicator_ACTIVATE_ON_ONE(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires(
        MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.GreaterThan{T}},
        ),
    )
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, x[1], MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), x[2])),
        ],
        T[0, 0],
    )
    s = MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(T(1)))
    c = MOI.add_constraint(model, f, s)
    @test MOI.get(model, MOI.ConstraintSet(), c) == s
    @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), f)
    return
end

"""
    test_constraint_Indicator_ACTIVATE_ON_ZERO(
        model::MOI.ModelLike,
        ::Config{T},
    ) where {T}

Test ACTIVATE_ON_ZERO for indicator sets.
"""
function test_constraint_Indicator_ACTIVATE_ON_ZERO(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires(
        MOI.supports_constraint(
            model,
            MOI.VectorAffineFunction{T},
            MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.GreaterThan{T}},
        ),
    )
    x = MOI.add_variables(model, 2)
    MOI.add_constraint(model, x[1], MOI.ZeroOne())
    f = MOI.VectorAffineFunction(
        [
            MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(T(1), x[1])),
            MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(T(2), x[2])),
        ],
        T[0, 0],
    )
    s = MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(MOI.GreaterThan(T(1)))
    c = MOI.add_constraint(model, f, s)
    @test MOI.get(model, MOI.ConstraintSet(), c) == s
    @test isapprox(MOI.get(model, MOI.ConstraintFunction(), c), f)
    return
end


text:
using ForwardNets
using Auto2D

type Tim2DExtractor <: AbstractFeatureExtractor
    feature_means::Vector{Float64} # for standardization
    feature_std::Vector{Float64} # for standardization

    function Tim2DExtractor(filepath::AbstractString,)
        feature_means = vec(h5read(filepath, "initial_obs_mean"))
        feature_std = vec(h5read(filepath, "initial_obs_std"))
        new(feature_means, feature_std)
    end
end

Base.length(::Tim2DExtractor) = N_FEATURES
_standardize(v::Real, μ::Real, σ::Real) = (v - μ)/σ
function AutomotiveDrivingModels.pull_features!{F<:AbstractFloat}(ext::Tim2DExtractor, features::Vector{F}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int)
    #debug = open("debug.log", "a"); println(debug, "AutomotiveDrivingModels.pull_features!{F<:AbstractFloat}(ext::Tim2DExtractor, features::Vector{F}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int)");  close(debug)
    pull_features!(features, rec, roadway, vehicle_index)

    # standardize
    for (i,v) in enumerate(features)
        features[i] = _standardize(v, ext.feature_means[i], ext.feature_std[i])
    end
    features
end

function _pull_W_b_h(filepath::AbstractString, path::AbstractString)
    W_xr = h5read(filepath, joinpath(path, "W_xr:0"))::Matrix{Float32}
    W_xu = h5read(filepath, joinpath(path, "W_xu:0"))::Matrix{Float32}
    W_xc = h5read(filepath, joinpath(path, "W_xc:0"))::Matrix{Float32}
    W_x = vcat(W_xr, W_xu, W_xc)

    W_hr = h5read(filepath, joinpath(path, "W_hr:0"))::Matrix{Float32}
    W_hu = h5read(filepath, joinpath(path, "W_hu:0"))::Matrix{Float32}
    W_hc = h5read(filepath, joinpath(path, "W_hc:0"))::Matrix{Float32}
    W_h = vcat(W_hr, W_hu, W_hc)

    b_r = h5read(filepath, joinpath(path, "b_r:0"))::Vector{Float32}
    b_u = h5read(filepath, joinpath(path, "b_u:0"))::Vector{Float32}
    b_c = h5read(filepath, joinpath(path, "b_c:0"))::Vector{Float32}
    b = vcat(b_r, b_u, b_c)

    (W_x, W_h, b)
end

function _pull_W_b(filepath::AbstractString, path::AbstractString)
    W = h5read(filepath, joinpath(path, "W:0"))::Matrix{Float32}
    b = h5read(filepath, joinpath(path, "b:0"))::Vector{Float32}
    (W,b)
end

function load_gru_driver(
            filepath::AbstractString = POLICY_FILEPATH,
            iteration::Int=-1;
            gru_layer::Bool=true,
            bc_policy::Bool=false,
            action_type::DataType=AccelTurnrate)
    #debug = open("debug.log", "a"); println(debug, "load_gru_driver(filepath=", filepath, " iteration=", iteration, " gru_layer=", gru_layer, ")");  close(debug)
    println("load_gru_driver(filepath=", filepath, " iteration=", iteration, " gru_layer=", gru_layer, ")");

#    if iteration == -1 # load most recent
#        fid = h5open(filepath)
#        for val in keys(read(fid))
#            i = parse(Int, match(r"\d+", val).match)
#            iteration = max(iteration, i)
#        end
#    end

    if gru_layer
         basepath = @sprintf("iter%05d/mlp_policy", iteration)
         layers = sort(collect(keys(h5read(filepath, basepath))))
         layers = layers[1:end-1]
    else
        if bc_policy
            basepath = @sprintf("iter%05d/mlp_policy", iteration)
            layers = sort(collect(keys(h5read(filepath, basepath))))
            layers = layers[1:end-2]
        else
            basepath = @sprintf("iter%05d/mlp_policy/mean_network", iteration)
            layers = sort(collect(keys(h5read(filepath, basepath))))
            layers = layers[1:end-1]
        end
    end

    W, b = _pull_W_b(filepath, joinpath(basepath, layers[1]))

    net = ForwardNet{Float32}()
    push!(net, Variable(:input, Array{Float32}(size(W, 2))))

    # hidden layers
    for layer in layers
        layer_sym = Symbol(layer)
        W, b = _pull_W_b(filepath, joinpath(basepath, layer))
        push!(net, Affine, layer_sym, lastindex(net), size(W, 1))
        copy!(net[layer_sym].W, W)
        copy!(net[layer_sym].b, b)
        push!(net, ELU, Symbol(layer*"ELU"), lastindex(net))
    end

    # Add GRU layer
    if gru_layer
        # output layer
        W, b = _pull_W_b(filepath, joinpath(basepath, "output"))
        push!(net, Affine, :output_layer_mlp, lastindex(net), size(W, 1))
        copy!(net[:output_layer_mlp].W, W)
        copy!(net[:output_layer_mlp].b, b)
        push!(net, Variable(:output_mlp, output(net[:output_layer_mlp])), lastindex(net))
        push!(net, ELU, Symbol("output_ELU"), lastindex(net))

        basepath = @sprintf("iter%05d/gru_policy", iteration)
        W_x, W_h, b = _pull_W_b_h(filepath, joinpath(basepath, "mean_network", "gru"))
        push!(net, GRU, :gru, :output_mlp, size(W_h, 2))
        copy!(net[:gru].W_x, W_x)
        copy!(net[:gru].W_h, W_h)
        copy!(net[:gru].b, b)
        copy!(net[:gru].h_prev, zeros(size(W_h, 2)))
        push!(net, Variable(:output_gru, output(net[:gru])), lastindex(net))

        # Finally, output layer from GRU
        W, b = _pull_W_b(filepath, joinpath(basepath, "mean_network", "output_flat"))
        push!(net, Affine, :output_layer, :gru, size(W, 1))
        copy!(net[:output_layer].W, W)
        copy!(net[:output_layer].b, b)
        push!(net, Variable(:output, output(net[:output_layer])), lastindex(net))
    else
        # Output layer
        if bc_policy
            W, b = _pull_W_b(filepath, joinpath(basepath, "mean_network", "output_flat"))
        else
            W, b = _pull_W_b(filepath, joinpath(basepath, "output"))
        end
        push!(net, Affine, :output_layer, lastindex(net), size(W, 1))
        copy!(net[:output_layer].W, W)
        copy!(net[:output_layer].b, b)
        push!(net, Variable(:output, output(net[:output_layer])), lastindex(net))
        basepath = @sprintf("iter%05d/mlp_policy", iteration)
    end

    if gru_layer
        logstdevs = vec(h5read(filepath, joinpath(basepath, "output_log_std/param:0")))::Vector{Float32}
    else
        if bc_policy
            logstdevs = vec(h5read(filepath, joinpath(basepath, "output_log_std/param:0")))::Vector{Float32}
        else
            logstdevs = vec(h5read(filepath, joinpath(basepath, "output_std_param/param:0")))::Vector{Float32}
        end
    end
    Σ = convert(Vector{Float64}, exp.(logstdevs))
    Σ = Σ.^2

    # extactor
    if gru_layer
        extractor = Tim2DExtractor(joinpath(ROOT_FILEPATH, "julia", "validation", "models", "gail_gru.h5"))
    else
        extractor = Tim2DExtractor(joinpath(ROOT_FILEPATH, "julia", "validation", "models", "gail_mlp.h5"))
    end

    Auto2D.GaussianMLPDriver(action_type, net, extractor, IntegratedContinuous(0.1,1),
                        input = :input, output = :output, Σ = Σ)
end

max_stars_repo_name: yoshihikokuwahara/gail-driver

max_issues_repo_name: yoshihikokuwahara/gail-driver

max_forks_repo_name: yoshihikokuwahara/gail-driver


text:
using ..MaterialParameters: MaterialParamsInfo
import GeoParams.param_info

export  DiffusionCreep,
        SetDiffusionCreep

const AxialCompression, SimpleShear, Invariant = 1,2,3

#=----Diffusion Creep---
Defines diffusion creep law parameters

n is the power-law exponent
r is the water-fugacity exponent
p is the negative defined grain size exponent (value is set negative in the calculation of σ and ε)
A is a   material specific rheological parameter
E is the activation energy
V is the activation volume
R is the universal gas constant
Apparatus defines the appartus type that shall be recreated (Axial Compression, Simple Shear, Invariant)
=#

@with_kw_noshow mutable struct DiffusionCreep{T,N,U1,U2,U3,U4,U5} <: AbstractCreepLaw{T}
    Name::NTuple{N,Char}        = ""
    n::GeoUnit{T,U1}            = 1.0NoUnits         # power-law exponent
    r::GeoUnit{T,U1}            = 0.0NoUnits         # exponent of water-fugacity
    p::GeoUnit{T,U1}            = -3.0NoUnits         # grain size exponent
    A::GeoUnit{T,U2}            = 1.5MPa^(-n-r)*s^(-1)*m^(-p)    # material specific rheological parameter
    E::GeoUnit{T,U3}            = 500kJ/mol          # activation energy
    V::GeoUnit{T,U4}            = 6e-6m^3/mol        # activation volume
    R::GeoUnit{T,U5}            = 8.314J/mol/K       # universal gas constant
    Apparatus::Int32            = AxialCompression # type of experimental apparatus, either AxialCompression, SimpleShear or Invariant
end

DiffusionCreep(args...) = DiffusionCreep(NTuple{length(args[1]), Char}(collect.(args[1])), convert.(GeoUnit,args[2:end-1])..., args[end])

function param_info(s::DiffusionCreep)
    name = String(collect(s.Name))
    eq = L"\tau_{ij} = 2 \eta  \dot{\varepsilon}_{ij}"
    if name == ""
        return MaterialParamsInfo(Equation=eq)
    end
    inf = DiffusionCreep_info[name][2]
    return MaterialParamsInfo(Equation=eq, Comment=inf.Comment, BibTex_Reference=inf.BibTex_Reference)
end

# Calculation routines for linear viscous rheologies
# All inputs must be non-dimensionalized (or converted to consitent units) GeoUnits
function computeCreepLaw_EpsII(TauII, a::DislocationCreep, c::CreepLawVariables)
    @unpack_val n,r,p,A,E,V,R = a
    @unpack_val P,T,f,d       = c

    FT, FE = CorrectionFactor(a);

    return A*(TauII*FT)^n*f^r*d^p*exp(-(E + P*V)/(R*T))/FE;
end

# EpsII .= A.*(TauII.*FT).^n.*f.^r.*d.^p.*exp.(-(E.+P.*V)./(R.*T))./FE; Once we have a
# All inputs must be non-dimensionalized (or converted to consistent units) GeoUnits
function computeCreepLaw_TauII(EpsII, a::DiffusionCreep, c::CreepLawVariables)
    @unpack_val n,r,A,E,V,R = a
    @unpack_val P,T,f,d     = c

    FT, FE = CorrectionFactor(a);

    return A^(-1/n)*(EpsII*FE)^(1/n)*f^(-r/n)*d^(-p/n)*exp((E + P*V)/(n * R*T))/FT;
end

# Print info
function show(io::IO, g::DiffusionCreep)
    print(io, "DiffusionCreep: Name = $(String(collect(g.Name))), n=$(g.n.val), r=$(g.r.val), p=$(g.p.val), A=$(g.A.val), E=$(g.E.val), V=$(g.V.val), Apparatus=$(g.Apparatus)" )
end

# This computes correction factors to go from experimental data to tensor format
# A nice discussion
function CorrectionFactor(a::DiffusionCreep{_T}) where {_T}

    FT = one(_T)
    FE = one(_T)
    if a.Apparatus == AxialCompression
        FT = sqrt(one(_T)*3)               # relation between differential stress recorded by apparatus and TauII
        FE = one(_T)*2/sqrt(one(_T)*3)     # relation between gamma recorded by apparatus and EpsII
    elseif a.Apparatus == SimpleShear
        FT = one(_T)*2                     # it is assumed that the flow law parameters were derived as a function of differential stress, not the shear stress. Must be modidified if it is not the case
        FE = one(_T)*2
    end
    return FT,FE
end

# Add pre-defined creep laws
"""
    SetDiffusionCreep["Name of Diffusion Creep"]
This is a dictionary with pre-defined creep laws
"""
SetDiffusionCreep(name::String) = DiffusionCreep_info[name][1]

# predefined diffusion creep laws are to be added in the dictionary as it is done for dislocation creep laws (see 'DislocationCreep.jl')!

DiffusionCreep_info = Dict([])