text:
type class_homogeneous_residuals <: abstract_residuals
	r_G::Float64
	r_D::Array{Float64,1}
	r_P::Array{Float64,1}

	r_norm::Float64
	r_D_norm::Float64
	r_G_norm::Float64
	r_P_norm::Float64

	r_D_norm_scaled::Float64
	r_G_norm_scaled::Float64
	r_P_norm_scaled::Float64

	mu::Float64
	scaled_mu::Float64

  val_c::Float64

	primal_norm::Float64
	dual_norm::Float64

	primal_infeas_sign::Int64
	dual_infeas_sign::Int64
	primal_infeas_norm::Float64
	dual_infeas_norm::Float64

  function class_homogeneous_residuals()
      return new();
  end
end

function finite_difference_merit_function(res::class_homogeneous_residuals, theta::class_theta, nlp, newton_solver, intial_point::class_variables, direction::class_variables)
		new_vars = deepcopy(intial_point)

		alpha = 1e-4;
		move!(new_vars, alpha, direction);
		orginal_mf = current_merit_function(res)
		val_mf = proximal_merit_function!(deepcopy(res), nlp, new_vars, newton_solver, intial_point);

		finite_dif = (val_mf - orginal_mf)/alpha

		computed_dif = merit_function_derivative(res, theta)

		@show finite_dif, computed_dif
end

function update_residuals!(res::class_homogeneous_residuals, nlp_eval::internal_AbstractNLPEvaluator, vars::class_variables, newt::abstract_newton_direction)
		update_residuals!(res, nlp_eval, vars, newt, vars)
end

function update_residuals!(res::class_homogeneous_residuals, nlp_eval::internal_AbstractNLPEvaluator, vars::class_variables, newt::abstract_newton_direction, intial_point::class_variables)
			try
				start_advanced_timer("residuals");

        nlp_vals = newt.nlp_vals;
        update_nlp_cache!(nlp_vals, nlp_eval, vars)

				res.mu = mu(newt,vars);
        res.val_c = nlp_vals.val_c

				val_x_scaled = x_scaled(vars);
        val_s = s(vars);
        val_kappa = kappa(vars);
        val_tau = tau(vars);

				prox_mod_x = -newt.delta_mod[1:n(vars)] .* (x(vars) - x(intial_point))
				res.r_D = s(vars) - tau(vars) * nlp_vals.val_gradlag + prox_mod_x;

				prox_mod_tau = -newt.delta_mod[n(vars)+1] * (tau(vars) - tau(intial_point))
				res.r_G = kappa(vars) + dot(nlp_vals.val_gradlag, x(vars)) + dot(nlp_vals.val_a, y(vars)) + prox_mod_tau;

				prox_mod_y = - newt.delta_mod[(n(vars)+2):(n(vars) + 1 + m(vars))] .* (y(vars) - y(intial_point))
				res.r_P = -tau(vars) * nlp_vals.val_a + prox_mod_y;

				res.r_D_norm = norm(res.r_D,1);
				res.r_G_norm = abs(res.r_G);
				res.r_P_norm = norm(res.r_P,1);
				res.r_norm = res.r_D_norm + res.r_G_norm + res.r_P_norm

				# other stuff
				scale = kappa(vars) + tau(vars);
				res.r_D_norm_scaled = res.r_D_norm/scale;
				res.r_G_norm_scaled = res.r_G_norm/scale;
				res.r_P_norm_scaled = res.r_P_norm/scale;

				res.scaled_mu = res.mu/(kappa(vars) + tau(vars));

				res.primal_norm = norm(res.r_P,1)/tau(vars);
				res.dual_norm = norm(res.r_D,1)/tau(vars);

				primal_infeas_obj = dot(nlp_vals.val_b, y(vars));
				dual_infeas_obj = -dot(x_scaled(vars), nlp_vals.val_gradc);
				res.primal_infeas_sign = sign(primal_infeas_obj);
				res.dual_infeas_sign = sign(dual_infeas_obj);

				res.primal_infeas_norm = norm(s(vars) + nlp_vals.val_jac_a' * y(vars),1)/abs(primal_infeas_obj);
 				res.dual_infeas_norm = norm(nlp_vals.val_a,1)/abs(dual_infeas_obj);

				pause_advanced_timer("residuals");
        #println(res.dual_infeas_norm, norm(nlp_vals.val_a,1)/nlp_vals.val_c)
			catch e
				  println("ERROR in class_residuals.update")
				  throw(e)
			end
end

max_forks_repo_name: ohinder/Everest


text:
# convergence/runtests.jl
# this file runs all the convergence tests

include("calc_line.jl")
include("calc_line2.jl")

include("./p1/conservative/runtests.jl")
include("./p1/conservative_3dg/runtests.jl")
include("./p1/conservative_dg/runtests.jl")
include("./p1/entropy/runtests.jl")
include("./p1/source_term/runtests.jl")

include("./p2/conservative/runtests.jl")
include("./p2/conservative_3dg/runtests.jl")
include("./p2/conservative_dg/runtests.jl")
include("./p2/entropy/runtests.jl")

include("./p3/conservative/runtests.jl")
include("./p3/conservative_dg/runtests.jl")
include("./p3/entropy/runtests.jl")

include("./p4/conservative/runtests.jl")
include("./p4/conservative_dg/runtests.jl")
include("./p4/entropy/runtests.jl")

function run_convergence_tests()
  @testset "----- testing convergence rates -----" begin

    original_dir = pwd()
    cd(dirname(@__FILE__))

    start_dir = pwd()
    println("start_dir = ", start_dir)
  #=
    # test p1 elements
    cd("./p1/conservative")
    run_convergence_p1_conservative()
  =#
    cd(start_dir)
    cd("./p1/conservative_dg")
    test_convergence_p1_dg()

    cd(start_dir)
    cd("./p1/source_term")
    test_convergence_p1_source()

    cd(start_dir)
    cd("./p1/conservative_3dg")
    test_convergence_p1_3dg()

  #=
    cd(start_dir)
    cd("./p1/entropy")
    test_convergence_p1_entropy()
  =#

  #----------------------------------------------------------------------------
  # p=2 elements
  #=
    cd(start_dir)
    cd("./p2/conservative")a
    test_convergence_p2_conservative()
  =#
    cd(start_dir)
    cd("./p2/conservative_dg")
    test_convergence_p2_dg()

    cd(start_dir)
    cd("./p2/conservative_3dg")
    test_convergence_p2_3dg()

  #=
    cd(start_dir)
    cd("./p2/entropy")
    test_convergence_p2_entropy()
  =#

  #----------------------------------------------------------------------------
  # p=3 elements
  #=
    cd(start_dir)
    cd("./p3/conservative")
    test_convergence_p3_conservative()
  =#
    cd(start_dir)
    cd("./p3/conservative_dg")
    test_convergence_p3_dg()

  #=
    cd(start_dir)
    cd("./p3/entropy")
    test_convergence_p3_entropy()
  =#

  #----------------------------------------------------------------------------
  # p=4 elements
  #=
    # test p4 elements
    cd(start_dir)
    cd("./p4/conservative")
    test_convergence_p4_conservative()
  =#
    cd(start_dir)
    cd("./p4/conservative_dg")
    test_convergence_p4_dg()

  #=
    cd(start_dir)
    cd("./p4/entropy")
    test_convergence_p4_entropy()
  =#
    # return to original directory
    cd(start_dir)

    cd(original_dir)

  end  # end facts block

  return nothing

end

#run_convergence_tests()
add_func1!(EulerTests, run_convergence_tests, [TAG_CONVERGENCE, TAG_NLSOLVERS, TAG_SHORTTEST])

ext: jl


text:
abstract type AbstractPartition{N} <: OnlineStat{N} end

nobs(o::AbstractPartition) = isempty(o.parts) ? 0 : sum(nobs, o.parts)

#-----------------------------------------------------------------------# Part
"""
    Part(stat, a, b)

`stat` summarizes a Y variable over an X variable's range `a` to `b`.
"""
mutable struct Part{T, O <: OnlineStat} <: OnlineStat{VectorOb}
    stat::O
    a::T
    b::T
end
# Part(o::O, ab::T) where {O<:OnlineStat, T} = Part{T,O}(o, ab, ab)
nobs(o::Part) = nobs(o.stat)
Base.show(io::IO, o::Part) = print(io, "Part $(o.a) to $(o.b) | $(o.stat)")
function _merge!(o::Part, o2::Part)
    _merge!(o.stat, o2.stat)
    o.a = min(o.a, o2.a)
    o.b = max(o.b, o2.b)
    o
end
Base.in(x, o::Part) = (o.a ≤ x ≤ o.b)
Base.isless(o::Part, o2::Part) = o.b < o2.a
value(o::Part) = value(o.stat)

midpoint(o::Part{<:Number}) = middle(o.a, o.b)

function midpoint(o::Part)
    if o.a == o.b
        return o.a
    else
        return (o.a:o.b)[round(Int, length(o.a:o.b) / 2)]
    end
end

midpoint(o::Part{<:TimeType}) = o.a + fld(o.b - o.a, 2)

width(o::Part) = o.b - o.a

isfull(o::Part{Int}) = (nobs(o) == o.b - o.a + 1)

function _fit!(p::Part, xy)
    x, y = xy
    x in p || error("$x ∉ [$(o.a), $(o.b)]")
    _fit!(p.stat, y)
end

#-----------------------------------------------------------------------# Partition
"""
    Partition(stat, nparts=100)

Split a data stream into `nparts` where each part is summarized by `stat`.

# Example

    o = Partition(Extrema())
    fit!(o, cumsum(randn(10^5)))

    using Plots
    plot(o)
"""
struct Partition{T, O <: OnlineStat{T}} <: AbstractPartition{T}
    parts::Vector{Part{Int, O}}
    b::Int  # max partition size
    init::O
end
function Partition(o::O, b::Int=100) where {T, O<:OnlineStat{T}}
    Partition{T, O}(Part{Int, O}[], b, o)
end
function _fit!(o::Partition, y)
    isempty(o.parts) && push!(o.parts, Part(copy(o.init), 1, 1))
    lastpart = last(o.parts)
    n = nobs(o)
    if (n + 1) ∈ lastpart
        _fit!(lastpart, (n+1, y))
    else
        stat = fit!(copy(o.init), y)
        push!(o.parts, Part(stat, n + 1, n + nobs(lastpart)))
    end
    (length(o.parts) > o.b) && isfull(o.parts[end]) && merge_next!(o.parts)
end
function _merge!(o::Partition, o2::Partition)
    n = nobs(o)
    for p in o2.parts
        push!(o.parts, Part(copy(p.stat), p.a + n, p.b + n))
    end
    while length(o.parts) > o.b
        merge_next!(o.parts)
    end
    o
end
function merge_next!(parts::Vector{<:Part})
    n = nobs(parts[1])
    ind = 1
    for (i, p) in enumerate(parts)
        if nobs(p) < n
            ind = i
            break
        end
    end
    if ind == length(parts)
        ind -= 1
    end
    merge!(parts[ind], parts[ind+1])
    deleteat!(parts, ind + 1)
end

#-----------------------------------------------------------------------# IndexedPartition
"""
    IndexedPartition(T, stat, b=100)

Summarize data with `stat` over a partition of size `b` where the data is indexed by a
variable of type `T`.

# Example

    o = IndexedPartition(Float64, Hist(10))
    fit!(o, eachrow(randn(10^4, 2)))

    using Plots
    plot(o)
"""
struct IndexedPartition{IN, O<:OnlineStat{IN}, T} <: AbstractPartition{VectorOb}
    parts::Vector{Part{T, O}}
    b::Int
    init::O
end
function IndexedPartition(T::Type, o::O, b::Int=100) where {IN, O<:OnlineStat{IN}}
    IndexedPartition{IN, O, T}(Part{T, O}[], b, o)
end
function _fit!(o::IndexedPartition{I,O,T}, xy) where {I,O,T}
    x, y = xy
    isempty(o.parts) && push!(o.parts, Part(copy(o.init), T(x), T(x)))
    addpart = true
    for p in o.parts
        if x in p
            addpart = false
            _fit!(p, xy)
            break
        end
    end
    addpart && push!(o.parts, Part(fit!(copy(o.init), y), x, x))
    length(o.parts) > o.b && merge_nearest!(sort!(o.parts))
end

function merge_nearest!(parts::Vector{<:Part})
    diff = get_diff(parts[1], parts[2])
    ind = 1
    for i in 2:(length(parts) - 1)
        newdiff = get_diff(parts[i], parts[i+1])
        if newdiff < diff
            diff = newdiff
            ind = i
        end
    end
    merge!(parts[ind], parts[ind + 1])
    deleteat!(parts, ind + 1)
end

get_diff(a::Part{T}, b::Part{T}) where {T<:Dates.TimeType} =
    (Dates.value(b.a) - Dates.value(a.b)) * (nobs(a) + nobs(b)) / 2

get_diff(a::Part{T}, b::Part{T}) where {T<:Number} = (b.a - a.b) * (nobs(a) + nobs(b)) / 2

function _merge!(o::IndexedPartition, o2::IndexedPartition)
    # If there's any overlap, merge
    for p2 in o2.parts
        pushpart = true
        for p in o.parts
            if (p2.a ∈ p) || (p2.b ∈ p)
                pushpart = false
                merge!(p, p2)
                break
            end
        end
        pushpart && push!(o.parts, p2)
    end
    # merge parts that overlap
    for i in reverse(2:length(o.parts))
        p1, p2 = o.parts[i-1], o.parts[i]
        if p1.a > p2.b
            # info("hey I deleted something at $i")
            merge!(p1, p2)
            deleteat!(o.parts, i)
        end
    end
    # merge until there's b left
    sort!(o.parts)
    while length(o.parts) > o.b
        merge_nearest!(o.parts)
    end
    o
end


text:
using SymbolicManipulations
using SymbolicUtils
using Test

# Taken from SymbolicUtils
macro eqtest(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end
macro eqtest_broken(expr)
    @assert expr.head == :call && expr.args[1] in [:(==), :(!=)]
    if expr.args[1] == :(==)
        :(@test_broken isequal($(expr.args[2]), $(expr.args[3])))
    else
        :(@test_broken !isequal($(expr.args[2]), $(expr.args[3])))
    end |> esc
end
SymbolicUtils.show_simplified[] = false

@testset "SymbolicManipulations.jl" begin
    @testset "together" begin
        @syms x y z

        @eqtest together(x/z + z/y) == (x*y + z^2) / (y*z)
        @eqtest together(x/z + (1 + z)/y) == (x*y + z*(1 + z)) /(y*z)
        @eqtest together((x + y)/z + (1 + z)/y) == (y*(x + y) + z*(1 + z)) /(y*z)
    end

    @testset "count_ops" begin
        @syms x y z

        @test SymbolicManipulations.count_ops(x) == 0
        @test SymbolicManipulations.count_ops(1) == 0
        @test SymbolicManipulations.count_ops(x + y) == 1
        @test SymbolicManipulations.count_ops(x + y + z) == 2

        @test SymbolicManipulations.count_ops(2x) == 1
        @test SymbolicManipulations.count_ops(2x + 0) == 1
        @test SymbolicManipulations.count_ops(2x + 1) == 2
        @test SymbolicManipulations.count_ops(2x + y + 1) == 3

        @test SymbolicManipulations.count_ops(x^2) == 1
        @test SymbolicManipulations.count_ops(x^2 * y^3) == 3

        @test SymbolicManipulations.count_ops(x^(2y + z)) == 3
    end

    @testset "fullsimplify" begin
        @syms x

        @eqtest fullsimplify(cos(2x) + 2sin(x)^2) == 1
        @eqtest fullsimplify(exp(x) * (exp(x) - 1)^-2) == 1//4*(sinh(1//2 * x))^-2
    end

    @testset "integrate" begin
        @syms x y

        @eqtest integrate(x^2, x) == x^3 // 3
        @eqtest integrate(x^y, x) == x^(y + 1) / (y + 1)
        @eqtest integrate(y^x, y) == y^(x + 1) / (x + 1)
        @eqtest integrate(5x, x) == 5//2 * x^2
        @eqtest integrate(x + x^2, x) == 1//2 * x^2 + 1//3 * x^3
        @eqtest integrate((x + x^2) / x, x) == x + 1//2 * x^2
        @eqtest integrate((x + y)^2, x) == y^2 * x + y * x^2 + 1//3 * x^3

        @eqtest_broken integrate(sin(x)^2 * cos(x), x) == sin(x)^3 // 3
        @eqtest_broken integrate((x^2 + x) / sqrt(x), x) == 2//15 * x^(3//2) * (3x + 5)
        @eqtest_broken integrate(x^4 / (1 - x^2)^(5//2), x) == asin(x) + x//3 * (4x^2 - 3) / (1 - x^2)^(3/2)
        @eqtest_broken integrate(sec(x)^2 / (1 + sec(x)^2 - 3tan(x)), x) == log(2cos(x) - sin(x)) - log(cos(x) - sin(x))
        @eqtest_broken integrate((sin(x) + cos(x))^2, x) == x - 1//2 * cos(2x)
        @eqtest_broken integrate(tan(x) * sec(x)^2, x) == sec(x)^2 // 2
        @eqtest_broken integrate(sin(x) * cos(x), x) == -1//2 * cos(x)^2

        # SAINT could not solve these:
        @eqtest_broken integrate(x * sqrt(1 + x), x) == 2//15 * (1 + x)^(3//2) * (3x - 2)
        @eqtest_broken integrate(cos(sqrt(x)), x) == 2cos(sqrt(x)) + 2sqrt(x)*sin(sqrt(x))
    end
end


text: # import MbedTLS


text:
# TODO Move the constants to CLIMAParameters or LandParameters
Base.@kwdef struct phys{T<:Number}
    visc_0::T = 13.3e-06;            # Kinematic viscosity at 0C and 1013.25 hPa (m2/s)
    Dh_0::T = 18.9e-06;              # Molecular diffusivity (heat) at 0C and 1013.25 hPa (m2/s)
    Dv_0::T = 21.8e-06;              # Molecular diffusivity (H2O) at 0C and 1013.25 hPa (m2/s)
    Dc_0::T = 13.8e-06;              # Molecular diffusivity (CO2) at 0C and 1013.25 hPa (m2/s)
    Wtoμmole_s::T = 4.57;            # converts Watts/m2 to μmole/s
    Cdryleaf::T = 1396;              # Specific heat of dry leaf at constant pressure (J/kg/K)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
"""
**ErrorAnalysis** is a collection of tools to calculate errors for Monte Carlo data.
"""
module ErrorAnalysis

include("binning.jl")
include("Jackknife.jl")
include("autocorrtimes.jl") # old

"""
    iswithinerrorbars(a, b, δ[, print=false])

Checks whether numbers `a` and `b` are equal up to given error `δ`.
Will print `x ≈ y + k·δ` for `print=true`.

Is equivalent to `isapprox(a,b,atol=δ,rtol=zero(b))`.
"""
function iswithinerrorbars(a::T, b::S, δ::Real, print::Bool=false) where T<:Number where S<:Number
  equal = isapprox(a,b,atol=δ,rtol=zero(δ))
  if print && !equal
    out = a>b ? abs(a-(b+δ))/δ : -abs(a-(b-δ))/δ
    println("x ≈ y + ",round(out,4),"·δ")
  end
  return equal
end

"""
    iswithinerrorbars(A::AbstractArray{T<:Number}, B::AbstractArray{T<:Number}, Δ::AbstractArray{<:Real}[, print=false])

Elementwise check whether `A` and `B` are equal up to given real error matrix `Δ`.
Will print `A ≈ B + K.*Δ` for `print=true`.
"""
function iswithinerrorbars(A::AbstractArray{T}, B::AbstractArray{S},
                           Δ::AbstractArray{<:Real}, print::Bool=false) where T<:Number where S<:Number
  size(A) == size(B) == size(Δ) || error("A, B and Δ must have same size.")

  R = iswithinerrorbars.(A,B,Δ,false)
  allequal = all(R)

  if print && !all(R)
    if T<:Real && S<:Real
      O = similar(A, promote_type(T,S))
      for i in eachindex(O)
        a = A[i]; b = B[i]; δ = Δ[i]
        O[i] = R[i] ? 0.0 : round(a>b ? abs(a-(b+δ))/δ : -abs(a-(b-δ))/δ,4)
      end
      println("A ≈ B + K.*Δ, where K is:")
      display(O)
    else
      warn("Unfortunately print=true is only supported for real input.")
    end
  end

  return allequal
end

# Exports
export iswithinerrorbars

export binning_error
export jackknife_error

export plot_binning_R
export plot_error

# export tau_binning
# export tau_integrated
# export tau_fitting

# export Neff_binning
# export Neff_integrated
# export Neff_fitting

# export tau_binning_error
# export tau_integrated_error
# export tau_fitting_error

end # module


text:
function model_evaluation(parsable_trees=parsable_trees; random_transp=false)
    if random_transp
        for (i,t) in enumerate(parsable_trees)
            transp = sample_chord_transposition()
            parsable_trees[i] = map(transp, t)
        end
    end

    vals = [
        single_fold_cross_validation(parsable_trees, make_jazz_grammar, true),
        single_fold_cross_validation(parsable_trees, make_jazz_grammar, false),
        single_fold_cross_validation(parsable_trees, make_plain_PCFG_jazz_grammar, true),
        single_fold_cross_validation(parsable_trees, make_plain_PCFG_jazz_grammar, false),
        single_fold_cross_validation_baseline(parsable_trees)
    ]
end

function sample_chord_transposition()
    i = rand(1:12)
    c -> begin
        m = match(r"(\d+)(.*)", c)
        string(mod(parse(m[1]) + i, 12), m[2])
    end
end


text:
module ManifoldLearning

    import Base: show, summary
    import SparseArrays: AbstractSparseArray, SparseMatrixCSC, spzeros, spdiagm, findnz
    import Statistics: mean, std
    import StatsBase: StatsBase, fit, standardize
    import MultivariateStats: fit, outdim, projection, transform, KernelPCA, principalvars
    import LinearAlgebra: eigvals, mul!, svd, qr, Symmetric, eigen, eigen!, diagm, tr,
                          rmul!, I, schur, opnorm
    import LightGraphs: neighbors, nv, add_edge!, connected_components, vertices,
                        dijkstra_shortest_paths, induced_subgraph, weights
    import SimpleWeightedGraphs: SimpleWeightedGraph

    export

    # Transformation types
    AbstractDimensionalityReduction,
    Isomap,             # Type: Isomap model
    HLLE,               # Type: Hessian Eigenmaps model
    LLE,                # Type: Locally Linear Embedding model
    LTSA,               # Type: Local Tangent Space Alignment model
    LEM,                # Type: Laplacian Eigenmaps model
    DiffMap,            # Type: Diffusion maps model

    ## common interface
    outdim,             # the output dimension of the transformation
    fit,                # perform the manifold learning
    transform,          # the projection matrix
    eigvals,            # eigenvalues from the spectral analysis
    neighbors,          # the number of nearest neighbors used for building local
    vertices,           # vertices of largest connected component
    projection          # the projection matrix (deprecated)

    abstract type AbstractDimensionalityReduction end
    vertices(R::AbstractDimensionalityReduction) = Int[]

    const Projection{T <: Real} = AbstractMatrix{T}

    include("utils.jl")
    include("isomap.jl")
    include("hlle.jl")
    include("lle.jl")
    include("ltsa.jl")
    include("lem.jl")
    include("diffmaps.jl")

    show(io::IO, ::MIME"text/plain", R::T) where {T<:AbstractDimensionalityReduction} = summary(io, R)
    function show(io::IO, R::T) where {T<:AbstractDimensionalityReduction}
        summary(io, R)
        io = IOContext(io, :limit=>true)
        println(io)
        println(io, "connected component: ")
        Base.show_vector(io, vertices(R))
        println(io)
        println(io, "eigenvalues: ")
        Base.show_vector(io, eigvals(R))
        println(io)
        println(io, "projection:")
        Base.print_matrix(io, transform(R), "[", ",","]")
    end

    # deprecated functions
    @deprecate transform(DimensionalityReduction, X; k=k, d=d) fit(DimensionalityReduction, X; k=k, maxoutdim=d)
    @deprecate projection(DimensionalityReductionModel) transform(DimensionalityReductionModel)
end


text:
using JuliaForDataScience
using Test

@testset "JuliaForDataScience.jl" begin
    # Write your own tests here.
end


text:
@component struct EditorMoveGuides
    es::NTuple{4,Entity} #for now only x,y,z axis and sphere
end
@component struct EditorRotateGuides
    es::NTuple{3,Entity}
end

@component @with_kw struct Movable
    axis_thickness::Float32 = 0.2f0
    axis_length::Float32    = 5.0f0
    font_size::Float32      = 0.3f0
end

@component @with_kw struct Rotatable
    ring_radii::Float32 = 5.0f0
    rotation_speed::Float32 = 0.02f0 / 200
    font_size::Float32 = 0.3f0
end

struct Editor <: System end

function Overseer.requested_components(::Editor)
    return (Movable, Rotatable, InstancedDefaultVao, EditorMoveGuides, EditorRotateGuides)
end

function Overseer.update(::Editor, m::AbstractLedger)
    poly_vao      = m[InstancedDefaultVao]
    movable       = m[Movable]
    rotatable     = m[Rotatable]
    spatial       = m[Spatial]
    rotation_comp = m[Rotation]
    move_guides   = m[EditorMoveGuides]
    rotate_guides = m[EditorRotateGuides]
    text          = m[Text]
    text_vao      = m[TextVao]
    line_vao      = m[LineVao]

    poly = m[PolygonGeometry] #arrows are here
    selectable = m[Selectable]
    camera = singleton(m, Camera3D) #this we need to know how much a dragged mouse correspond to a movement in the world
    camera_entity = Entity(m[Camera3D].indices.packed[1])
    eye = spatial[camera_entity].position
    keyboard = singleton(m, Keyboard)
    mouse = singleton(m, Mouse)
    rotated = false
    moved = false

    set_visibility = b -> begin
        for g in (move_guides, rotate_guides)
            for e in @entities_in(g)
                for guide_entity in g[e].es
                    if guide_entity in poly_vao
                        poly_vao[guide_entity].visible = b
                    end
                    if guide_entity in text_vao
                        text_vao[guide_entity].visible = b
                    end
                    if guide_entity in line_vao
                        line_vao[guide_entity].visible = b
                    end
                end
            end
        end
    end

    if pressed(keyboard) && keyboard.button ∈ CTRL_KEYS
        set_visibility(true)
        #setup move guides
        for e in @entities_in(movable && spatial && !move_guides)
            mov            = movable[e]
            editor_guides  = assemble_axis_arrows(spatial[e].position; axis_length = mov.axis_length, thickness = mov.axis_thickness)
            guide_entities = map(x -> Entity(m, x..., Selectable()), editor_guides)
            move_guides[e] = EditorMoveGuides(guide_entities)
        end

        #setup rotation_comp guides
        for e in @entities_in(rotatable && spatial && !rotate_guides)
            if !in(e, rotation_comp)
                rotation_comp[e] = Rotation(Z_AXIS, 0.0f0)
            end

            editor_guides = assemble_orientation_sphere(spatial[e].position;
                                                        radius = rotatable[e].ring_radii)
            guide_entities = map(x -> Entity(m, x..., Selectable(; color_modifier = 1.5f0)),
                                 editor_guides)
            rotate_guides[e] = EditorRotateGuides(guide_entities)
        end

        #sanitize positions of move_guides
        for e in @entities_in(move_guides)
            guide_entities = move_guides[e].es
            t_spat = spatial[e]
            for g in guide_entities
                spatial[g] = t_spat
            end
            pos = t_spat.position
            text[guide_entities[1]] = Text(;
                                           str = "Position: (x=$(pos[1]), y=$(pos[2]), z=$(pos[3]))",
                                           font_size = movable[e].font_size,
                                           offset = Vec3f0(0, 0, -2.0f0))
        end

        #sanitize positions and orientiations of rotation_comp guides
        for e in @entities_in(rotate_guides)
            rot = rotation_comp[e]
            guide_entities = rotate_guides[e].es
            for guide in guide_entities
                spatial[guide] = spatial[e]
            end
            ax = direction(rot)
            rotation_comp[guide_entities[1]] = rot #this is the in plane one, so should rotate along the direction of the object
            # rotation_comp[guide_entities[2]] = Rotation(rot.q * rotation(Y_AXIS, X_AXIS))#this is the out of plane one, so should rotate along the perpendicular direction of the object
            rotation_comp[guide_entities[2]] = Rotation(rot.q * rotation(X_AXIS, Z_AXIS))
            rotation_comp[guide_entities[3]] = Rotation(rot.q * rotation(Y_AXIS, Z_AXIS))
            #this is the out of plane one, so should rotate along the perpendicular direction of the object
            text[guide_entities[1]] = Text(;
                                           str = "Orientation: (x=$(ax[1]), y=$(ax[2]), z=$(ax[3]))",
                                           font_size = rotatable[e].font_size,
                                           offset = Vec3f0(0, 0, -2.0f0))
        end

        #handle moving
        for parent_entity in @entities_in(move_guides)
            for e in move_guides[parent_entity].es[2:end]
                if selectable[e].selected
                    if pressed(mouse)
                        mouse_drag = (mouse.dx * camera.right - mouse.dy * camera.up) *
                                     camera.translation_speed / 2
                        guide_direction = direction(rotation_comp[e])
                        move = guide_direction * dot(guide_direction, mouse_drag)

                        t_spat = spatial[parent_entity]
                        spatial[parent_entity] = Spatial(t_spat;
                                                         position = t_spat.position + move)
                        moved = true
                    end
                end
            end
        end

        #handle rotating
        for parent_entity in @entities_in(rotate_guides)
            rotate_entities = rotate_guides[parent_entity].es
            for e in rotate_entities
                if selectable[e].selected
                    if pressed(mouse)
                        rotation_axis = direction(rotation_comp[e])
                        ray = Ray(mouse, singleton(m, Canvas), camera, eye)
                        plane = Plane(spatial[e].position, rotation_axis)
                        closest_point = intersect(ray, plane)
                        tangent = cross(rotation_axis, closest_point - spatial[e].position)
                        mouse_drag = (mouse.dx * camera.right - mouse.dy * camera.up) *
                                     rotatable[parent_entity].rotation_speed

                        guide_rotation = Quaternions.qrotation(rotation_axis,
                                                               dot(tangent, mouse_drag))

                        t_rot = rotation_comp[parent_entity]
                        rotation_comp[parent_entity] = Rotation(guide_rotation * t_rot.q)
                        rotated = true
                    end
                end
            end
        end
    else
        set_visibility(false)
    end
    rotated && push!(singleton(m, UpdatedComponents), Rotation)
    return moved && push!(singleton(m, UpdatedComponents), Spatial)
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ de098278-8e74-11eb-13e3-49c1b6e06e7d
using Pkg; Pkg.activate("MLJ_env", shared=true)

# ╔═╡ 10d6d90e-8d43-11eb-216c-8f94a3946fab
begin
	push!(LOAD_PATH, "./src")
	using YaDeepLearning
end

# ╔═╡ ac463e7a-8b59-11eb-229e-db560e17c5f5
begin
	using Test
	using Random
	using Distributions
	using PlutoUI
	using Plots
	using Flux.Data.MNIST
	using Images
end

# ╔═╡ 8c80e072-8b59-11eb-3c21-a18fe43c4536
md"""
## Deep Learning (DL) (cont'ed)

ref. from book **"Data Science from Scratch"**, Chap 19

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ 3972ef40-8d4f-11eb-3369-6dd94e7cfc13
Layers.Linear <: AbstractLayers.AbstractLayer, Layers.Linear <: AbstractLayers.AL

# ╔═╡ 6b3c5460-8dab-11eb-3046-331119a0dc9b
# on using: using vs import and provide extension
# cf. https://towardsdatascience.com/julias-big-problem-with-namespace-996d2e9ed71e

# ╔═╡ e7373726-8b59-11eb-2a2b-b5138e4f5268
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ f5ee64b2-8b59-11eb-2751-0778efd589cd
html"""

"""

# ╔═╡ 81290d1c-8ce2-11eb-3340-337957fd81b7
md"""
#### Dropout

Droput is a technique for regularizing neural network (to prevent overfitting). It works as follows:
  - at training time we randomly turn off a set of neurons (by replacing their output with 0) given a fixed probability. This seems to prevent overfitting by forcing the neural network to not rely  on any particular neuron.

  - at evaluation time, we do not drop the neurons, but we scale down the outputs (uniformly) by using the same fraction (used at training time).
"""

# ╔═╡ 8ff1bb20-8ce2-11eb-1de6-fd84daec8930
md"""
Let us define a helper function on initialize the binary mask we need to apply dorpout at training time.
"""

# ╔═╡ d3ee2138-8ce2-11eb-0b29-659a3be01512
function prand(dt::DataType, shape::NTuple{N, Integer};
               keep_proba::Float64=.5)::BitArray where {N}
	rand(dt, shape) .≥ (1. - keep_proba)
end

# ╔═╡ 13a0940e-8ce4-11eb-231c-f331a607203c
begin
	using .TensorDT: Tensor
	import .AbstractLayers: AbstractLayer, forward, backward, parms, ∇parms

	mutable struct Dropout <: AbstractLayer
		keep_prob::Float64
		train::Bool
		mask::BitArray
		_type::Symbol
		#
		Dropout(;keep_prob=.5, train=true) =
			new(keep_prob, train, BitArray(undef), :regularization)
	end

	function forward(self::Dropout, input::Tensor)::Tensor
		if self.train
			self.mask = prand(eltype(input), size(input);
				keep_proba=self.keep_prob)
			input .* self.mask
		else
			input * self.keep_prob # (1. - )
		end
	end

	function backward(self::Dropout, ∇p::Tensor)::Tensor
		!self.train &&
			throws(RuntimeError("Do NOT all backward when NOT in training mode"))
		∇p .* self.mask
	end

	parms(_self::Dropout) = []
	∇parms(_self::Dropout) = []
end

# ╔═╡ d3d074e4-8ce2-11eb-0307-d932886e4b10
md"""
Then let us define our dropout  layer.
"""

# ╔═╡ 80ddf0f4-8ce2-11eb-3046-331119a0dc9b
html"""

"""

# ╔═╡ 80a71408-8ce2-11eb-3978-75a4a2df9116
md"""
#### MNIST
"""

# ╔═╡ d1ed1e36-8ce5-11eb-34ee-4db60fb0db8a
const NROWS, NCOLS, NLABELS = 28, 28, 10

# ╔═╡ d2b198b6-8ce2-11eb-170e-0f17904c9f2c
begin
  ## Load labels and images from Flux.Data.MNIST
  train_labels, train_imgs = MNIST.labels(:train), MNIST.images(:train);
  test_labels, test_imgs = MNIST.labels(:test), MNIST.images(:test)
  size(train_labels), size(train_imgs), size(test_labels), size(test_imgs)
end

# ╔═╡ d37f37b8-8ce8-11eb-2c00-3f98ca407f41
size(train_imgs[1]), size(test_imgs[1])

# ╔═╡ 5bf7729e-8dfd-11eb-070f-9b7ec0746bd7
typeof(train_imgs[1])

# ╔═╡ 394093fa-8cea-11eb-0071-eff3045a012b
plot(train_imgs[1]), train_labels[1]

# ╔═╡ 435c4444-8dfd-11eb-24e8-5543f905f199
# sum((train_imgs[10])),
minimum(Float32.(train_imgs[1]))

# ╔═╡ e95983ba-8ceb-11eb-38fd-ed92cdcf754c
mosaicview(train_imgs[1:100]; nrow=5, ncol=20, npad=5)

# ╔═╡ d3a749a2-8cec-11eb-1f06-b568f244b576
md"""
Each image is  28x28 pixels, but our linear layer can only deal with one-dimensional input, therefore we need to flatten them.

We will also recenter both train and test sets by using the mean calculated over the training set.
"""

# ╔═╡ 31d3226a-8cf4-11eb-0897-39989ba76b58
begin
	## reshape (flatten) and cast from Gray{Normed{UInt8,8}} to Float32
	_train_imgs = [Float32.(img) for img ∈ (reshape.(train_imgs, NROWS * NCOLS))]
	_test_imgs = [Float32.(img) for img ∈ (reshape.(test_imgs, NROWS * NCOLS))]

	## calc. avg
	img_avg = Float32.(sum(_train_imgs) / length(_train_imgs) / (NROWS * NCOLS))
	@assert size(_train_imgs[1]) == size(img_avg) == (NROWS * NCOLS,)

	## Re-center, re-scale
	λ = img -> Float32.(img .- img_avg)
	ntrain_imgs = λ.(_train_imgs); _train_imgs = nothing
	ntest_imgs = λ.(_test_imgs); _test_imgs = nothing

  	## ntrain_imgs ≡ Array{Array{Float32,1},1}
	@assert size(ntrain_imgs) == size(train_imgs)
	@assert size(ntest_imgs) == size(test_imgs)
end

# ╔═╡ 6e7e7896-8cf8-11eb-062e-99492ec8cff8
md"""
Let us write a one hot encoder for our 10 outputs (digits from 0 to 9)
"""

# ╔═╡ 70f707c0-8cf6-11eb-11c2-73d7b28f7a0c
function one_hot_encoder(label::Integer; nlabels=NLABELS, DT=Float64)::Tensor
	# @assert label < nlabels
	a = zeros(DT, nlabels)
	a[label + 1] = one(DT)
	a
end

# ╔═╡ 70d84470-8cf6-11eb-23b6-c7ba506d8552
begin
  @test one_hot_encoder(3) == [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.]
  @test one_hot_encoder(2; nlabels=5) == [0., 0., 1., 0., 0.]
end

# ╔═╡ c451451a-8cf7-11eb-183b-e358c9d618e0
begin
	ntrain_labels = one_hot_encoder.(train_labels)
	ntest_labels = one_hot_encoder.(test_labels);
	ntrain_labels[1:3]
end

# ╔═╡ ebd531e6-8cf7-11eb-228f-d530e7a48151
md"""
Let adjust our (previoulsy defined) train_loop function.
"""

# ╔═╡ 30bdc788-8d01-11eb-2f8f-9fd79593ddb8
begin

const AL, ALoss, NAOpt = AbstractLayers.AL, Loss.ALoss, AbstractOptimizers.NAOpt

function train_loop(model::AL, imgs::AA, labels::AA, loss_fn::ALoss;
    epochs=1, opt_fn::NAOpt=nothing, verbose=false)
	#
	rec_avg_loss, rec_acc = [], []
	ncorr, totloss = 0, 0.0
	#
	for e ∈ 1:epochs, t in 1:length(imgs)
		ŷ = Layers.forward(model, imgs[t])
		y = labels[t]
		## NOTE argmax => CartesianIndex
		argmax(ŷ)[1] == argmax(y) && (ncorr += 1)
		totloss += Loss.loss(loss_fn, ŷ, y)
		#
		if !isnothing(opt_fn)
			∇p = Loss.∇loss(loss_fn, ŷ, y)
			Layers.backward(model, ∇p)
			Optimizers.astep(opt_fn, model)
		end
		avg_loss = totloss / (e * t)
		acc = ncorr / (e * t)
		verbose && t % 5000 == 0 && (@show "$(e)/$(t)", avg_loss, acc)
		push!(rec_avg_loss, avg_loss)
		push!(rec_acc, acc)
	end
	(rec_avg_loss, rec_acc)
end

end

# ╔═╡ 84946000-8d02-11eb-3160-571efee8fb0b
md"""
Let us create a baseline to train a mluticlass logistic regression (which is a single linear layer followed by a softmax). One pass through the 60000 training examples should be enough to learn the model.
"""

# ╔═╡ d05f1c40-8d17-11eb-2724-4d989c4c6b92
begin
	lr_model = Layers.Linear(NROWS * NCOLS, NLABELS)
	lr_loss = Loss.SoftmaxXEntropy()
	lr_opt = Optimizers.MomentumGD(;η=.01, α=.99)
end

# ╔═╡ 848fd812-8d02-11eb-01d6-75965b08bcc5
with_terminal() do
	# Train
	global (tr1_avg_loss, tr1_acc) = train_loop(lr_model, ntrain_imgs, ntrain_labels,
		lr_loss; opt_fn=lr_opt, verbose=true)

	# Test model
	global (te1_avg_loss, te1_acc) = train_loop(lr_model, ntrain_imgs, ntrain_labels,		lr_loss)
	# rouglhy 90% of accuracy
end

# ╔═╡ 7d35f85e-8e07-11eb-228c-81d28a444b52
begin
	ntr, nte = 60_000, 10_000
	lg1 = grid(2, 1, heights=[0.5, 0.5])
	a = plot(tr1_avg_loss[1:ntr], label="train avg loss", leg=:bottomright)
	b = plot(tr1_acc[1:ntr], label="train acc.")
	plot(a, b, layout=lg1)
end

# ╔═╡ 6e7c18b4-8e09-11eb-0238-2515d39a89cd
begin
	lg2 = grid(2, 1, heights=[0.5, 0.5])
	a₂ = plot(te1_avg_loss[1:nte], label="test avg loss", leg=:bottomright)
	b₂ = plot(te1_acc[1:nte], label="test acc.")
	plot(a₂, b₂, layout=lg2)
end

# ╔═╡ 40d46802-8e0f-11eb-2e47-53096e24dbd8
md"""
And now let us define and train our neural network
"""

# ╔═╡ 80dff352-8d02-11eb-06b1-5f5e325046f5
begin
	dropout₁ = Dropout(;keep_prob=0.9)
	dropout₂ = Dropout(;keep_prob=0.8)

	m_model = Layers.Sequential([
		Layers.Linear(NROWS * NCOLS, 3 * NLABELS),
		dropout₁,
		Activations.Tanh(),
		Layers.Linear(3 * NLABELS, NLABELS),
		dropout₂,
		Activations.Tanh(),
		Layers.Linear(NLABELS, NLABELS)
	])
	m_loss = Loss.SoftmaxXEntropy()
	m_opt = Optimizers.MomentumGD(;η=.01, α=.99)
	with_terminal() do
		## Train
		global (tr_avg_loss, tr_acc) = train_loop(m_model, ntrain_imgs,
			ntrain_labels, m_loss; epochs=5, opt_fn=m_opt)

		## Test
		dropout₁.train = dropout₂.train = false
		global (te_avg_loss, te_acc) = train_loop(m_model, ntrain_imgs,
			ntrain_labels, m_loss, verbose=true)
	end
end

# ╔═╡ f1af4256-8e02-11eb-0238-2515d39a89cd
begin
	## 2 epochs for training...
	x₁ = plot(tr_avg_loss[1:ntr], label="train avg loss", leg=:bottomright)
	x₂ = plot(tr_avg_loss[ntr+1:end], label="train avg loss", leg=:bottomright)

	y₁ = plot(tr_acc[1:ntr], label="train acc.", leg=:bottomright)
	y₂ = plot(tr_acc[ntr+1:end], label="train acc.")
	plot(x₁, x₂, y₁, y₂, layout=(2, 2))
end

# ╔═╡ 43a52692-8e10-11eb-2043-8f0be195f58a
begin
	x₃ = plot(te_avg_loss[1:ntr], label="test1 avg loss", leg=:bottomright)
	y₃ = plot(te_acc[1:ntr], label="test2 acc.", leg=:bottomright)
	plot(x₃, y₃, layout=(2, 1))
end

# ╔═╡ Cell order:
# ╠═8c80e072-8b59-11eb-3c21-a18fe43c4536
# ╠═de098278-8e74-11eb-13e3-49c1b6e06e7d
# ╠═ac463e7a-8b59-11eb-229e-db560e17c5f5
# ╠═10d6d90e-8d43-11eb-216c-8f94a3946fab
# ╠═3972ef40-8d4f-11eb-3369-6dd94e7cfc13
# ╠═6b3c5460-8dab-11eb-3046-331119a0dc9b
# ╟─e7373726-8b59-11eb-2a2b-b5138e4f5268
# ╟─f5ee64b2-8b59-11eb-2751-0778efd589cd
# ╟─81290d1c-8ce2-11eb-3340-337957fd81b7
# ╟─8ff1bb20-8ce2-11eb-1de6-fd84daec8930
# ╠═d3ee2138-8ce2-11eb-0b29-659a3be01512
# ╟─d3d074e4-8ce2-11eb-0307-d932886e4b10
# ╠═13a0940e-8ce4-11eb-231c-f331a607203c
# ╟─80ddf0f4-8ce2-11eb-3046-331119a0dc9b
# ╟─80a71408-8ce2-11eb-3978-75a4a2df9116
# ╠═d1ed1e36-8ce5-11eb-34ee-4db60fb0db8a
# ╠═d2b198b6-8ce2-11eb-170e-0f17904c9f2c
# ╠═d37f37b8-8ce8-11eb-2c00-3f98ca407f41
# ╠═5bf7729e-8dfd-11eb-070f-9b7ec0746bd7
# ╠═394093fa-8cea-11eb-0071-eff3045a012b
# ╠═435c4444-8dfd-11eb-24e8-5543f905f199
# ╠═e95983ba-8ceb-11eb-38fd-ed92cdcf754c
# ╟─d3a749a2-8cec-11eb-1f06-b568f244b576
# ╠═31d3226a-8cf4-11eb-0897-39989ba76b58
# ╟─6e7e7896-8cf8-11eb-062e-99492ec8cff8
# ╠═70f707c0-8cf6-11eb-11c2-73d7b28f7a0c
# ╠═70d84470-8cf6-11eb-23b6-c7ba506d8552
# ╠═c451451a-8cf7-11eb-183b-e358c9d618e0
# ╟─ebd531e6-8cf7-11eb-228f-d530e7a48151
# ╠═30bdc788-8d01-11eb-2f8f-9fd79593ddb8
# ╟─84946000-8d02-11eb-3160-571efee8fb0b
# ╠═d05f1c40-8d17-11eb-2724-4d989c4c6b92
# ╠═848fd812-8d02-11eb-01d6-75965b08bcc5
# ╠═7d35f85e-8e07-11eb-228c-81d28a444b52
# ╠═6e7c18b4-8e09-11eb-0238-2515d39a89cd
# ╟─40d46802-8e0f-11eb-2e47-53096e24dbd8
# ╠═80dff352-8d02-11eb-06b1-5f5e325046f5
# ╠═f1af4256-8e02-11eb-0238-2515d39a89cd
# ╠═43a52692-8e10-11eb-2043-8f0be195f58a

max_stars_repo_path: Data Science From Scratch/19-DL-p2.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/19-DL-p2.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/19-DL-p2.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
function SW07()
    # Parameters (with values)
    # ε_w = 10
    # ρ_ga = 0.51
    # ε_p = 10
    # l_bar = 0
    # Π_bar = 0.7
    # B = 0.742
    # μ_w = 0
    # μ_p = 0
    # α = 0.24
    # ψ = 0.2696
    # φ = 6.0144
    # δ = 0.025
    # σ_c = 1.5
    # λ = 0.6361
    # ϕ_p = 1.5
    # ι_w = 0.3243
    # ξ_w = 0.8087
    # ι_p = 0.47
    # ξ_p = 0.6
    # σ_l = 1.9423
    # ϕ_w = 1.5
    # r_π = 1.488
    # r_Δy = 0.2347
    # r_y = 0.0593
    # ρ = 0.8762
    # ρ_a = 0.9977
    # ρ_b = 0.5799
    # ρ_g = 0.9957
    # ρ_i = 0.7165
    # ρ_r = 0
    # ρ_p = 0
    # ρ_w = 0
    # γ_bar = 0.3982
    # gy_ss = 0.18

    # Define symbolic ones
    @parameters ε_w ρ_ga ε_p l_bar Π_bar B μ_w μ_p α ψ φ δ σ_c λ ϕ_p ι_w ξ_w ι_p ξ_p σ_l ϕ_w r_π r_Δy r_y ρ ρ_a ρ_b ρ_g ρ_i ρ_r ρ_p ρ_w γ_bar gy_ss # 34 parameters
    # states
    @variables y_f_m y_m k_f k c_f_m c_m i_f_m i_m π_m w_m r_m ε_a b ε_g ε_i ε_r ε_pm η_p_aux ε_wm η_w_aux # 20 variables
    # controls
    @variables k_s_f k_s r_f r r_k_f r_k z_f z w_f w l_f l q_f q y_f y i_f i c_f c π μ_pm # 22 variables
    # Markovize
    @make_markov y_f_m y_m k_f k c_f_m c_m i_f_m i_m π_m w_m r_m ε_a b ε_g ε_i ε_r ε_pm η_p_aux ε_wm η_w_aux k_s_f k_s r_f r r_k_f r_k z_f z w_f w l_f l q_f q y_f y i_f i c_f c π μ_pm # 42 variables

    x_sym = [y_f_m, y_m, k_f, k, c_f_m, c_m, i_f_m, i_m, π_m, w_m, r_m, ε_a, b, ε_g, ε_i,
             ε_r, ε_pm, η_p_aux, ε_wm, η_w_aux]
    y_sym = [k_s_f, k_s, r_f, r, r_k_f, r_k, z_f, z, w_f, w, l_f, l, q_f, q, y_f, y, i_f, i,
             c_f, c, π, μ_pm]
    p = [ε_w, ρ_ga, ε_p, l_bar, Π_bar, B, μ_w, μ_p, α, ψ, φ, δ, σ_c, λ, ϕ_p, ι_w, ξ_w, ι_p,
         ξ_p, σ_l, ϕ_w, r_π, r_Δy, r_y, ρ, ρ_a, ρ_b, ρ_g, ρ_i, ρ_r, ρ_p, ρ_w, γ_bar, gy_ss]
    p_f = []

    # Pre-compute variables
    π_bar = 1 + Π_bar / 100
    γ = 1 + γ_bar / 100
    β = 1 / (1 + B / 100)
    β_bar = β * γ^-σ_c
    R_star = π_bar / (β * γ^-σ_c)
    R_k_star = 1 / (β * γ^-σ_c) - (1 - δ)
    w_star = (α^α * (1 - α)^(1 - α) / (ϕ_p * R_k_star^α))^(1 / (1 - α))
    k_1 = 1 - (1 - δ) / γ
    i_k = (1 - (1 - δ) / γ) * γ
    l_k = (1 - α) / α * R_k_star / w_star
    k_y = ϕ_p * l_k^(α - 1)
    i_y = i_k * k_y
    c_y = 1 - gy_ss - i_k * k_y
    z_y = R_k_star * k_y
    w_hlc = 1 / ϕ_w * (1 - α) / α * R_k_star * k_y / c_y

    # Defining H
    H = [α * r_k_f + (1 - α) * w_f - ε_a, (1 - ψ) / ψ * r_k_f - z_f,
         w_f + l_f - k_s_f - r_k_f, k_f + z_f - k_s_f,
         1 / (1 + β * γ) * (i_f_m + β * γ * i_f_p + 1 / (γ^2 * φ) * q_f) + ε_i - i_f,
         -r_f +
         b / ((1 - λ / γ) / (σ_c * (1 + λ / γ))) +
         R_k_star / (R_k_star + 1 - δ) * r_k_f_p +
         (1 - δ) / (R_k_star + 1 - δ) * q_f_p - q_f,
         λ / γ / (1 + λ / γ) * c_f_m +
         1 / (1 + λ / γ) * c_f_p +
         (σ_c - 1) * w_hlc / (σ_c * (1 + λ / γ)) * (l_f - l_f_p) -
         (1 - λ / γ) / (σ_c * (1 + λ / γ)) * r_f + b - c_f,
         c_y * c_f + i_y * i_f + ε_g + z_y * z_f - y_f,
         ϕ_p * (α * k_s_f + (1 - α) * l_f + ε_a) - y_f,
         σ_l * l_f + 1 / (1 - λ / γ) * c_f - λ / γ / (1 - λ / γ) * c_f_m - w_f,
         (1 - k_1) * k_f + k_1 * i_f + k_1 * γ^2 * φ * ε_i - k_f_p,
         # Sticky Part
         α * r_k + (1 - α) * w - ε_a - μ_pm, (1 - ψ) / ψ * r_k - z, w + l - k_s - r_k,
         k + z - k_s, 1 / (1 + β * γ) * (i_m + β * γ * i_p + 1 / (γ^2 * φ) * q) + ε_i - i,
         -r +
         π_p +
         b / ((1 - λ / γ) / (σ_c * (1 + λ / γ))) +
         R_k_star / (R_k_star + 1 - δ) * r_k_p +
         (1 - δ) / (R_k_star + 1 - δ) * q_p - q,
         λ / γ / (1 + λ / γ) * c_m +
         1 / (1 + λ / γ) * c_p +
         (σ_c - 1) * w_hlc / (σ_c * (1 + λ / γ)) * (l - l_p) -
         (1 - λ / γ) / (σ_c * (1 + λ / γ)) * (r - π_p) + b - c,
         c_y * c + i_y * i + ε_g + z_y * z - y, ϕ_p * (α * k_s + (1 - α) * l + ε_a) - y,
         1 / (1 + β * γ * ι_p) * (β * γ * π_p +
                                  ι_p * π_m +
                                  (1 - ξ_p) * (1 - β * γ * ι_p) / (ξ_p * ((ϕ_p - 1) * ε_p + 1)) * μ_pm) +
         ε_pm - π,
         1 / (1 + β * γ) * w_m + β * γ / (1 + β * γ) * w_p + ι_w / (1 + β * γ) * π_m -
         (1 + β * γ * ι_w) / (1 + β * γ) * π +
         β * γ / (1 + β * γ) * π_p +
         (1 - ξ_w) * (1 - β * γ * ξ_w) / ((1 + β * γ) * ξ_w) / ((ϕ_w - 1) * ε_w + 1) *
         (σ_l * l + 1 / (1 - λ / γ) * c - λ / γ / (1 - λ / γ) * c_m - w) +
         ε_wm - w,
         r_π * (1 - ρ) * π +
         r_y * (1 - ρ) * (y - y_f) +
         r_Δy * (y - y_f - y_m + y_f_m) +
         ρ * r_m +
         ε_r, (1 - k_1) * k + k_1 * i + k_1 * γ^2 * φ * ε_i - k_p,
         # Shock Transition
         ρ_a * ε_a - ε_a_p, ρ_b * b - b_p, ρ_g * ε_g - ε_g_p, ρ_i * ε_i - ε_i_p,
         ρ_r * ε_r - ε_r_p, ρ_p * ε_pm + η_p_aux_p - μ_p * η_p_aux - ε_pm_p, 0 - η_p_aux_p,
         ρ_w * ε_wm + η_w_aux_p - μ_w * η_w_aux - ε_wm_p, 0 - η_w_aux_p,
         # Variable Mapping
         y_f - y_f_m_p, y - y_m_p, c_f - c_f_m_p, c - c_m_p, i_f - i_f_m_p, i - i_m_p,
         π - π_m_p, w - w_m_p, r - r_m_p]

    x̄ = [y_f_m_ss ~ 0, y_m_ss ~ 0, k_f_ss ~ 0, k_ss ~ 0, c_f_m_ss ~ 0, c_m_ss ~ 0,
          i_f_m_ss ~ 0, i_m_ss ~ 0, π_m_ss ~ 0, w_m_ss ~ 0, r_m_ss ~ 0, ε_a_ss ~ 0,
          b_ss ~ 0, ε_g_ss ~ 0, ε_i_ss ~ 0, ε_r_ss ~ 0, ε_pm_ss ~ 0, η_p_aux_ss ~ 0,
          ε_wm_ss ~ 0, η_w_aux_ss ~ 0]
    ȳ = [k_s_f_ss ~ 0, k_s_ss ~ 0, r_f_ss ~ 0, r_ss ~ 0, r_k_f_ss ~ 0, r_k_ss ~ 0,
          z_f_ss ~ 0, z_ss ~ 0, w_f_ss ~ 0, w_ss ~ 0, l_f_ss ~ 0, l_ss ~ 0, q_f_ss ~ 0,
          q_ss ~ 0, y_f_ss ~ 0, y_ss ~ 0, i_f_ss ~ 0, i_ss ~ 0, c_f_ss ~ 0, c_ss ~ 0,
          π_ss ~ 0, μ_pm_ss ~ 0]

    x̄_iv = deepcopy(x̄)
    ȳ_iv = deepcopy(ȳ)

    n_ϵ = 2 # TODO: change
    n_x = length(x_sym)
    n_y = length(y_sym)
    n_p = length(p)
    n_p = length(p_f)
    Γ = zeros(Num, n_ϵ, n_ϵ) # TODO: change, also make sure it is not a float64 matrix
    Γ[1, 1] = 1
    Γ[2, 2] = 1
    η = zeros(n_x, n_ϵ) # TODO: change

    return H, (; x = x_sym, y = y_sym, x̄, ȳ, Γ, η, p_f, p, x̄_iv, ȳ_iv), "SW07"
end


text:
using PyPlot
using JLD
using Statistics, StatsBase
using ColorSchemes
using Printf

@load "/Users/milan/cams/entropy/information_all_unstructured.jld"
ICgridded = load("/Users/milan/cams/entropy/information_all.jld","IC")

# normalize
IC = mean(IC,dims=2)[:,1,:]
ICgridded = mean(ICgridded,dims=2)[:,1,:]

IC[varnames.=="t",:] .= ICgridded[varnames.=="t",:]

varnames[varnames .== "c"] .= "ch4_c"

nvars = length(varnames)

## sort into groups
aero = Array(1:15)
ozone = [35,44,45]
methane = [25,26,41]
dynamics = [33,34,54,55,53]
clouds = [22,27,28,31,32,51]
hydro = vcat(Array(36:40),46)
nitro = [42,43,52]
oopp = Array(47:50)
ces = vcat(Array(17:21),[23,24,16])
co12 = [29,30]

grouped = vcat(aero,co12,clouds,methane,ces,hydro,dynamics,nitro,ozone,oopp)
groups = [aero,co12,clouds,methane,ces,hydro,dynamics,nitro,ozone,oopp]

# sort
varnames = varnames[grouped]
IC = IC[grouped,:]

## sum over 3 dimensions
# normalise such that max entropy = 1
ICsum = copy(IC)
ICsumnan = copy(ICsum)
ICsumnan[iszero.(ICsumnan)] .= NaN

## bit start
mininf = [argmin(ICsum[i,11:22])+9 for i in 1:nvars]

IC_noiseremoved = copy(ICsum)
for i in 1:nvars
    IC_noiseremoved[i,mininf[i]:end] .= 0
    # IC_noiseremoved[i,IC_noiseremoved[i,:] .< 0.01] .= 0
end

ICcsum = cumsum(IC_noiseremoved,dims=2)
ICcsum_norm = copy(ICcsum)
for i in 1:nvars
    ICcsum_norm[i,:] ./= ICcsum_norm[i,end]
end

inflevel = 0.98
infbits = [argmax(ICcsum_norm[i,:] .> inflevel) for i in 1:nvars]

# manual adjustment
for varname in ["aermr03","aermr08","aermr10","aermr16","aermr18",
                "c2h4","c5h8","no","so2","ole","par"]
    infbits[varnames.==varname] .= 1+8+3
end

# for cloud
infbits[varnames.=="cc"] .= 1+8+5
for varname in ["ciwc","clwc","crwc","cswc"]
    infbits[varnames.==varname] .= 1+8+3
end

# for aerosols
for varname in ["aermr04","aermr05","aermr06"]
    infbits[varnames.==varname] .= 1+8+5
end

# don't go below 0
infbits[infbits .< 9] .= 9

@save "/Users/milan/cams/entropy/keepbits_98.jld" infbits varnames

infbitsx = copy(vec(hcat(infbits,infbits)'))
infbitsy = copy(vec(hcat(Array(0:nvars-1),Array(1:nvars))'))
## plotting
ioff()
fig,ax1 = subplots(1,1,figsize=(8,10),sharey=true)
ax1.invert_yaxis()
tight_layout(rect=[0.06,0.08,0.93,0.98])
pos = ax1.get_position()
cax = fig.add_axes([pos.x0,0.06,pos.x1-pos.x0,0.02])

ax1right = ax1.twinx()
ax1right.invert_yaxis()

# information
cmap = ColorMap(ColorSchemes.turku.colors).reversed()
cmap_array = cmap(0:1/100000:1)
cmap_array[1,:] = [1,1,1,1]

pcm = ax1.pcolormesh(ICsumnan,vmin=0,vmax=1;cmap)
cbar = colorbar(pcm,cax=cax,orientation="horizontal")
cbar.set_label("information content [bit]")

# useful bit ranges
ax1.plot(vcat(infbits,infbits[end]),Array(0:nvars),"C1",ds="steps-pre",zorder=10)

# for legend only
ax1.fill_betweenx([-1,-1],[-1,-1],[-1,-1],color="w",edgecolor="k",label="unused bits")

# grey shading
ax1.fill_betweenx(infbitsy,infbitsx,fill(32,length(infbitsx)),alpha=0.5,color="grey",label="noise / false information")

ax1.axvline(1,color="k",lw=1,zorder=3)
ax1.axvline(9,color="k",lw=1,zorder=3)

for (ig,group) in enumerate(groups)
    y = sum([length(g) for g in groups[1:ig]])
    ax1.axhline(y,color="w",lw=1.5,zorder=2)
end

ax1.set_title("Bitwise information content",loc="left",fontweight="bold")

ax1.set_xlim(0,32)
ax1.set_ylim(nvars,0)
ax1right.set_ylim(nvars,0)

ax1.set_yticks(Array(1:nvars).-0.5)
ax1right.set_yticks(Array(1:nvars).-0.5)
ax1.set_yticklabels(varnames)
ax1right.set_yticklabels([@sprintf "%4.1f" i for i in ICcsum[:,end]])
ax1right.set_ylabel("information per value [bit]")
# ax1.set_ylabel("variable")

ax1.text(infbits[1]+0.1,0.8,"$(infbits[1]-9) significant bits",fontsize=8,color="k")
for i in 2:nvars
    ax1.text(infbits[i]+0.1,(i-1)+0.8,"$(infbits[i]-9)",fontsize=8,color="k")
end

ax1.set_xticks([1,9])
ax1.set_xticklabels([])
ax1.text(0,nvars+2,"sign",rotation=90)
ax1.text(2,nvars+1.2,"exponent")
ax1.text(10,nvars+1.2,"significant bits")

ax1.legend(loc=(0.64,0.21))

savefig("/Users/milan/git/Elefridge.jl/plots/bitinformation_all.png",dpi=100)
close(fig)


text:
@testset "JSON" begin
    io = IOBuffer()
    val, e = readjson(io)
    @test val === nothing
    @test e isa Empty

    io = IOBuffer("{bad]")
    val, e = readjson(io)
    @test val === nothing
    @test e !== nothing

    io = IOBuffer("[1,2,3]")
    val, e = readjson(io)
    @test val == [1, 2, 3]
    @test e === nothing

    # writejson can throw.
    io = IOBuffer()
    @test_throws ArgumentError writejson(io, @__MODULE__)

    io = IOBuffer()
    writejson(io, [1, 2, 3])
    @test String(take!(io)) == "[1,2,3]"

    io = IOBuffer()
    d = Dict("foo" => 0, "bar" => [1, 2, 3])
    writejson(io, d)
    @test JSON.parse(String(take!(io))) == d

    # The try version just returns the exception and backtrace.
    io = IOBuffer()
    e, bt = trywritejson(io, @__MODULE__)
    @test e isa ArgumentError

    io = IOBuffer()
    e, bt = trywritejson(io, [1, 2, 3])
    @test e === nothing && bt === nothing
end


text:
using Schrodinger
using Test, LinearAlgebra, SparseArrays
using Schrodinger: gaussianpulse

println("Testing Optimal Control Functions...")

function opt_2Q_QFT(n)
    # Hamiltonian & controls
    Hd = 0.5 * (σx⊗σx + σy⊗σy + σz⊗σz)
    Hc = [0.5*σx⊗σ0, 0.5*σy⊗σ0, 0.5*σ0⊗σx, 0.5*σ0⊗σy]
    t = 6.0
    # Initial controls
    tvec = range(-t/2,t/2,length=n)
    ui = reduce(hcat, @. [2*gaussianpulse(tvec,[[t/4,t,2π/t,0,π]]),
                          2*gaussianpulse(tvec,[[t/4,t,2π/t,-π/2,π]]),
                          3*gaussianpulse(tvec,[[t/5,t,4π/t,0,π]]),
                          3*gaussianpulse(tvec,[[t/5,t,4π/t,-π/2,π]])])
    # Create objective function
    Ut = Operator([0.5  0.5    0.5  0.5
                   0.5  0.5im -0.5 -0.5im
                   0.5 -0.5    0.5 -0.5
                   0.5 -0.5im -0.5  0.5im],
                  (2,2))
    O = NormPSU(Ut,Hd,Hc,t,n)
    return grape(O,ui)
end

grape_res = opt_2Q_QFT(200)

@testset "Optimal Control Tests" begin
@test 1 - gate_fidelity(grape_res.Uf, grape_res.Ut) < 1E-8
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

given a random X variable taking values in the finite set {α₁,...,αₙ}, compute
the distribution minimizing a given convex spectral function over all distributions
satisfying some prior information (expressed using convex constraints)

adapted from Boyd and Vandenberghe, "Convex Optimization", section 7.2

p ∈ ℝᵈ is the probability variable
minimize    f(p)            (note: enforces p ≥ 0)
subject to  Σᵢ pᵢ = 1       (probability distribution)
            gⱼ(D p) ≤ kⱼ ∀j (prior info as convex constraints)
            B p = b         (prior info as equalities)
            C p ≤ c         (prior info as inequalities)
where f and gⱼ are different convex spectral functions
=#

struct NonparametricDistrJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    d::Int
    num_spec::Int # number of spectral cones
    use_EFs::Bool # use standard cone extended formulations for spectral cones
end

function build(inst::NonparametricDistrJuMP{T}) where {T <: Float64}
    d = inst.d
    @assert d >= 2
    p0 = rand(T, d)
    p0 ./= sum(p0)

    # pick random spectral cones (or EFs)
    if inst.use_EFs
        exts = [VecNegGeomEFExp(), VecNegGeomEFPow(), VecInvEF(), VecNegLogEF(),
            VecNegEntropyEF(), VecPower12EF(1.5), VecNeg2SqrtEF()]
    else
        exts = [VecNegGeom(), VecInv(), VecNegLog(), VecNegEntropy(),
            VecPower12(1.5), VecNeg2Sqrt()]
    end
    @assert all(is_domain_pos, exts) # all domains must be positive
    @assert 1 <= inst.num_spec <= length(exts)
    Random.seed!(inst.num_spec)
    exts = Random.shuffle!(exts)[1:inst.num_spec]

    model = JuMP.Model()
    JuMP.@variable(model, p[1:d])
    JuMP.@constraint(model, sum(p) == 1)

    # linear prior constraints
    B = randn(T, round(Int, sqrt(d - 2)), d)
    b = B * p0
    JuMP.@constraint(model, B * p .== b)
    C = randn(T, round(Int, log(d - 1)), d)
    c = C * p0
    JuMP.@constraint(model, C * p .<= c)

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(exts[1], d, vcat(1.0 * epi, p), model)

    # convex constraints
    con_aff = Vector{Tuple{T, Matrix{T}}}()
    for ext in exts[2:end]
        D = rand(T, d, d)
        val_p0 = get_val(D * p0, ext)
        push!(con_aff, (val_p0, D))
        add_homog_spectral(ext, d, vcat(val_p0, D * p), model)
    end

    # save for use in tests
    model.ext[:exts] = exts
    model.ext[:con_aff] = con_aff
    model.ext[:p_var] = p

    return model
end

function test_extra(inst::NonparametricDistrJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective and constraints
    tol = eps(T)^0.2
    exts = model.ext[:exts]
    con_aff = model.ext[:con_aff]
    p_opt = JuMP.value.(model.ext[:p_var])
    d = length(p_opt)
    @test sum(p_opt) ≈ 1 atol=tol rtol=tol
    @test minimum(p_opt) >= -tol
    p_opt = pos_only(p_opt)
    obj_result = get_val(p_opt, exts[1])
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    for (i, ext) in enumerate(exts[2:end])
        (val_p0, D) = con_aff[i]
        @test val_p0 >= get_val(D * p_opt, ext) - tol
    end
    return
end


text:
using JuMP, Cbc

simMod = Model(Cbc.Optimizer)

@variable(simMod, x1>=0)
@variable(simMod, x2>=0)

@constraints(simMod, begin
    2*x1+x2<=4
    x1+2*x2<=4
    x1+x2>=1
end)

@objective(simMod, Max, -10x1-4x2)

res = JuMP.optimize!(simMod)

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
# Command: Accept Alert
"""
    accept(session::Session)::Nothing

Accepts the currently displayed alert dialog.

!!! note

    This command uses the JSON Wire Protocol instead of the current W3c WebDriver API.
"""
function accept(session::Session)::Nothing
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/accept_alert",
        [("Content-Type" => "application/json; charset=utf-8")],
    )
    @assert response.status == 200
end


text:
using NameResolution
using Test
using PrettyPrint

@testset "share freevar" begin
    println("""test case:
    function f(x)
        y = 1 + x
        g -> begin
            y = 2
            y + g
        end
    end
    """)
    ana = top_analyzer()
    enter!(ana, :f)
    is_local!(ana, :x)
    enter!(ana, :x)

    enter!(ana, :y)
    require!(ana, :x)
    lambda = child_analyzer!(ana)

    is_local!(lambda, :g)
    enter!(lambda, :g)

    enter!(lambda, :y)
    require!(lambda, :y)
    require!(lambda, :g)

    run_analyzer(ana)
    print("f ")
    pprint(ana.solved)
    println()
    print("lambda ")
    pprint(lambda.solved)
    println()

    @test lambda.solved.freevars[:y] === ana.solved.bounds[:y]
    @test lambda.solved.freevars[:y].is_mutable.x === true
    @test lambda.solved[:kk] === :kk
end

@testset "local" begin
    println("""test case:
    function f(x)
        y = 1 + x
        g -> begin
            local y = 2
            y + g
        end
    end
    """)
    ana = top_analyzer()
    enter!(ana, :f)
    is_local!(ana, :x)
    enter!(ana, :x)

    enter!(ana, :y)
    require!(ana, :x)
    lambda = child_analyzer!(ana)

    is_local!(lambda, :g)
    enter!(lambda, :g)

    is_local!(lambda, :y)
    enter!(lambda, :y)
    require!(lambda, :y)
    require!(lambda, :g)

    run_analyzer(ana)
    print("f ")
    pprint(ana.solved)
    println()
    print("lambda ")
    pprint(lambda.solved)
    println()

    @test haskey(lambda.solved.bounds, :y)
    @test lambda.solved.bounds[:y].is_mutable.x === false
end

@testset "local + mutable" begin
    println("""test case:
    function f(x)
        y = 1 + x
        g -> begin
            local y = 2
            y += 1
            y + g
        end
    end
    """)
    ana = top_analyzer()
    enter!(ana, :f)
    is_local!(ana, :x)
    enter!(ana, :x)

    enter!(ana, :y)
    require!(ana, :x)
    lambda = child_analyzer!(ana)

    is_local!(lambda, :g)
    enter!(lambda, :g)

    is_local!(lambda, :y)
    enter!(lambda, :y)
    enter!(lambda, :y)
    require!(lambda, :y)
    require!(lambda, :g)

    run_analyzer(ana)
    print("f ")
    pprint(ana.solved)
    println()
    print("lambda ")
    pprint(lambda.solved)
    println()

    @test haskey(lambda.solved.bounds, :y)
    @test lambda.solved.bounds[:y].is_mutable.x === true
    @test lambda.solved[:y].is_mutable.x === true
end

@testset "scope extrusion" begin
    @test_throws Any begin
        ana = top_analyzer()
        !is_local(ana, :x)
        !is_global(ana, :x)
        enter!(ana, :x)
        run_analyzer(ana)
    end
end

@testset "deeper scope" begin
    ana = top_analyzer()
    is_local!(ana, :x)
    enter!(ana, :x)
    shallow = child_analyzer!(ana, true)
    deep = child_analyzer!(shallow, true)
    require!(deep, :x)
    run_analyzer(ana)

    @test haskey(shallow.solved.freevars, :x)
end

@testset "uninitialized" begin
    ana = top_analyzer()
    is_local!(ana, :x)
    require!(ana, :x)
    run_analyzer(ana)
    @test haskey(ana.solved.bounds, :x)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

export @var

struct Binding
    mod::Module
    var::Symbol

    function Binding(m::Module, v::Symbol)
        # Normalise the binding module for module symbols so that:
        #   Binding(Base, :Base) === Binding(Main, :Base)
        m = nameof(m) === v ? parentmodule(m) : m
        new(Base.binding_module(m, v), v)
    end
end

bindingexpr(x) = Expr(:call, Binding, splitexpr(x)...)

defined(b::Binding) = isdefined(b.mod, b.var)
resolve(b::Binding) = getfield(b.mod, b.var)

function splitexpr(x::Expr)
    isexpr(x, :macrocall) ? splitexpr(x.args[1]) :
    isexpr(x, :.)         ? (x.args[1], x.args[2]) :
    error("Invalid @var syntax `$x`.")
end
splitexpr(s::Symbol) = Expr(:macrocall, getfield(Base, Symbol("@__MODULE__")), nothing), quot(s)
splitexpr(r::GlobalRef) = r.mod, quot(r.name)
splitexpr(other)     = error("Invalid @var syntax `$other`.")

macro var(x)
    esc(bindingexpr(x))
end

function Base.show(io::IO, b::Binding)
    if b.mod === Main
        print(io, b.var)
    else
        print(io, b.mod, '.', Base.isoperator(b.var) ? ":" : "", b.var)
    end
end

aliasof(b::Binding)     = defined(b) ? (a = aliasof(resolve(b), b); defined(a) ? a : b) : b
aliasof(d::DataType, b) = Binding(d.name.module, d.name.name)
aliasof(λ::Function, b) = (m = typeof(λ).name.mt; Binding(m.module, m.name))
aliasof(m::Module,   b) = Binding(m, nameof(m))
aliasof(other,       b) = b

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
@testset "Functions and Calls" begin
    tmpdir = mktempdir()
    atexit(() -> rm(tmpdir; force = true, recursive = true))

    model_filename = joinpath(tmpdir, "model.jld2")
    saving_filename = joinpath(tmpdir, "saving.jl")
    loading_filename = joinpath(tmpdir, "loading.jl")

    saving_contents = """
    module TestModule
        struct MyHardStruct{T,T2,T3}
            f::T
            u::T2
            p::T3
        end
        function (m::MyHardStruct)(t)
            m.u[2]*t + (1-t)*m.u[1]
        end
        export MyHardStruct
    end
    using .TestModule
    function f(u,p,t)
        [u[1],u[2]] .* p
    end
    u = [[1.0,3.0],[5.0,5.0]]
    p = (1.0,2.0)
    m = MyHardStruct(f,u,p)

    using JLD2
    @save "out.jld2" m
    """

    loading_contents = """
    module TestModule
        struct MyHardStruct{T,T2,T3}
            f::T
            u::T2
            p::T3
        end
        function (m::MyHardStruct)(t)
            m.u[2]*t + (1-t)*m.u[1]
        end
        export MyHardStruct
    end
    using .TestModule, JLD2, Test
    function f(u,p,t)
        [u[1],u[2]] .* p
    end
    JLD2.@load "out.jld2" m
    @test m(0.5) == [3.0,4.0]
    @test m.f(u[1],m.p,0.0) == [1.0,6.0]
    """

    rm(mycustomload_filename; force = true, recursive = true)
    rm(model_filename; force = true, recursive = true)
    rm(saving_filename; force = true, recursive = true)
    rm(loading_filename; force = true, recursive = true)

    open(saving_filename, "w") do io
        println(io, saving_contents)
    end
    open(loading_filename, "w") do io
        println(io, loading_contents)
    end

    saving_cmd = `$(Base.julia_cmd()) $(saving_filename)`
    loading_cmd = `$(Base.julia_cmd()) $(loading_filename)`

    rm(model_filename; force = true, recursive = true)

    @test success(saving_cmd)
    @test success(loading_cmd)

    rm(tmpdir; force = true, recursive = true)
end


text:
export logsoftmax, logsumexp, softmax, logp
using AutoGrad: AutoGrad, @primitive1, value

"""
    softmax(x; dims=:)
    logsoftmax(x; dims=:)

Treat entries in `x` as as unnormalized log probabilities and return normalized (log)
probabilities, i.e.

    softmax(x; dims) = exp.(x) ./ sum(exp.(x); dims=dims)
    logsoftmax(x; dims) = x .- log.(sum(exp.(x); dims=dims))

For numerical stability `x = x .- maximum(x,dims=dims)` is performed before exponentiation.

`dims` is an optional argument, if not specified the normalization is over the whole `x`,
otherwise the normalization is performed over the given dimensions.  In particular, if `x` is
a matrix, `dims=1` normalizes columns of `x` and `dims=2` normalizes rows of `x`.

"""
softmax, logsoftmax, logp

"""
    logsumexp(x;dims=:)

Compute `log(sum(exp(x);dims))` in a numerically stable manner.

`dims` is an optional argument, if not specified the summation is over the whole `x`,
otherwise the summation is performed over the given dimensions.  In particular if `x` is a
matrix, `dims=1` sums columns of `x` and `dims=2` sums rows of `x`.

"""
logsumexp

# Separate user facing methods from underlying implementations
logsoftmax(x; dims=:) = _logsoftmax(x, dims=dims)
softmax(x; dims=:)    = _softmax(x, dims=dims)
logsumexp(x; dims=:)  = _logsumexp(x, dims=dims)
const logp = logsoftmax

function _logsoftmax(x; dims)
    isa(x, Number) && return zero(x)
    isempty(x) && return similar(x)
    x = x .- maximum(x,dims=dims)
    return (x .- log.(sum(exp.(x); dims=dims)))
end

function ∇logsoftmax(x,y,dy; dims)
    return (dy - exp.(y).*sum(dy; dims=dims))
end

function _softmax(x; dims)
    isa(x, Number) && return one(x)
    isempty(x) && return similar(x)
    x = x .- maximum(x; dims=dims)
    x = exp.(x)
    return (x ./ sum(x; dims=dims))
end

function ∇softmax(x,y,dy; dims)
    return (y .* dy .- y .* sum(y .* dy; dims=dims))
end

function _logsumexp(x; dims)
    isa(x, Number) && return x
    isempty(x) && throw(ArgumentError("reducing over an empty collection is not allowed"))
    xmax = maximum(x,dims=dims)
    return (xmax + log.(sum(exp.(x .- xmax),dims=dims)))
end

function ∇logsumexp(x,y,dy; dims)
    return (dy .* exp.(x .- y))
end

@primitive1 _logsoftmax(x;dims),dy,y  ∇logsoftmax(x,y,dy,dims=dims)
@primitive1 _softmax(x;dims),dy,y     ∇softmax(x,y,dy,dims=dims)
@primitive1 _logsumexp(x;dims),dy,y   ∇logsumexp(x,y,dy;dims=dims)

# Math for the softmax loss: x is unnormalized input, p is target probabilities, q is
# estimated probabilities. Read left column down, right column (loss gradients) back up.
#
# x			dx = -p + qz/z = -p + exp(logq)
# xmax  = max(x,1)	-sum(db)=0
# logqz = x .- xmax	-p + qz/z
# qz    = exp(logqz)	rep(1/z)
# z     = sum(qz,1)	1/z
# logz  = log(z)	sum(p)=1
# logq  = logqz.-logz	-p
# plogq = p .* logq	-1
# loss  = -sum(plogq)	1


text:
function generate_training_sets(domains,dx,eqs,eltypeθ)
    if dx isa Array
        dxs = dx
    else
        dxs = fill(dx,length(domains))
    end
    spans = [infimum(d.domain):dx:supremum(d.domain) for (d,dx) in zip(domains,dxs)]
    train_set = adapt(eltypeθ, hcat(vec(map(points -> collect(points), Iterators.product(spans...)))...))
end

function get_loss_function_(loss,initθ,pde_system,strategy::GridTraining)
    eqs = pde_system.eqs
    if !(eqs isa Array)
        eqs = [eqs]
    end
    domains = pde_system.domain
    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars, pde_system.depvars)
    eltypeθ = eltype(initθ)
    parameterless_type_θ =  DiffEqBase.parameterless_type(initθ)
    dx = strategy.dx
    train_set = generate_training_sets(domains,dx,eqs,eltypeθ)
    get_loss_function(loss,train_set,eltypeθ,parameterless_type_θ,strategy)
end

function get_bounds_(domains,eqs,eltypeθ,dict_indvars,dict_depvars,strategy)
    dict_span = Dict([Symbol(d.variables) => [infimum(d.domain), supremum(d.domain)] for d in domains])
    args = get_argument(eqs,dict_indvars,dict_depvars)

    bounds= map(args) do pd
        span = map(p -> get(dict_span, p, p), pd)
        map(s -> adapt(eltypeθ,s), span)
    end
    bounds
end

function get_loss_function_(loss,initθ,pde_system,strategy::StochasticTraining)
    eqs = pde_system.eqs
    if !(eqs isa Array)
        eqs = [eqs]
    end
    domains = pde_system.domain

    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars, pde_system.depvars)

    eltypeθ = eltype(initθ)
    parameterless_type_θ =  DiffEqBase.parameterless_type(initθ)

    bound = get_bounds_(domains,eqs,eltypeθ,dict_indvars,dict_depvars,strategy)[1]

    get_loss_function(loss,bound,eltypeθ,parameterless_type_θ,strategy)
end

function get_loss_function_(loss,initθ,pde_system,strategy::QuasiRandomTraining)
    eqs = pde_system.eqs
    if !(eqs isa Array)
        eqs = [eqs]
    end
    domains = pde_system.domain

    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars, pde_system.depvars)

    eltypeθ = eltype(initθ)
    parameterless_type_θ =  DiffEqBase.parameterless_type(initθ)

    bound = get_bounds_(domains,eqs,eltypeθ,dict_indvars,dict_depvars,strategy)[1]

    get_loss_function(loss,bound,eltypeθ,parameterless_type_θ,strategy)
end

function get_bounds_(domains,eqs,eltypeθ,dict_indvars,dict_depvars,strategy::QuadratureTraining)
    dict_lower_bound = Dict([Symbol(d.variables) => infimum(d.domain) for d in domains])
    dict_upper_bound = Dict([Symbol(d.variables) => supremum(d.domain) for d in domains])

    args = get_argument(eqs,dict_indvars,dict_depvars)

    lower_bounds= map(args) do pd
        span = map(p -> get(dict_lower_bound, p, p), pd)
        map(s -> adapt(eltypeθ,s), span)
    end
    upper_bounds= map(args) do pd
        span = map(p -> get(dict_upper_bound, p, p), pd)
        map(s -> adapt(eltypeθ,s), span)
    end
    bound = lower_bounds,upper_bounds
end

function get_loss_function_(loss,initθ,pde_system,strategy::QuadratureTraining)
    eqs = pde_system.eqs
    if !(eqs isa Array)
        eqs = [eqs]
    end
    domains = pde_system.domain

    depvars,indvars,dict_indvars,dict_depvars = get_vars(pde_system.indvars, pde_system.depvars)

    eltypeθ = eltype(initθ)
    parameterless_type_θ =  DiffEqBase.parameterless_type(initθ)

    bound = get_bounds_(domains,eqs,eltypeθ,dict_indvars,dict_depvars,strategy)
    lb,ub = bound
    get_loss_function(loss,lb[1],ub[1],eltypeθ,parameterless_type_θ,strategy)
end

"""
the method that trains a neural network using the results from one already obtained prediction.

Arguments:
* `loss`: the body of loss function,
* `initθ`: the initial parameter of the neural network,,
* `pde_system`: PDEs are defined using the ModelingToolkit.jl,
* `strategy`: determines which training strategy will be used.

"""
function neural_adapter(loss,initθ,pde_system,strategy)
    loss_function__ = get_loss_function_(loss,initθ,pde_system,strategy)

    function loss_function_(θ,p)
        loss_function__(θ)
    end
    f_ = OptimizationFunction(loss_function_, Optimization.AutoZygote())
    prob = Optimization.OptimizationProblem(f_, initθ)
end

"""
the method that trains a neural network using the results from many already obtained predictions.

Arguments:
* `loss`: the body of loss functions,
* `initθ`: the initial parameter of the neural network,
* `pde_system`: PDEs are defined using the ModelingToolkit.jl,
* `strategy`: determines which training strategy will be used.

"""
function neural_adapter(losses::Array,initθ,pde_systems::Array,strategy)
    loss_functions_ = map(zip(losses,pde_systems)) do (l,p)
        get_loss_function_(l,initθ,p,strategy)
    end
    loss_function__ = θ -> sum(map(l->l(θ) ,loss_functions_))
    function loss_function_(θ,p)
        loss_function__(θ)
    end

    f_ = OptimizationFunction(loss_function_, Optimization.AutoZygote())
    prob = Optimization.OptimizationProblem(f_, initθ)
end


text:
using Pkg.Artifacts
using polymake_jll

# adapted from https://github.com/JuliaLang/julia/pull/38797
function full_artifact_dir(m::Module)
   artifacts_toml = joinpath(dirname(dirname(Base.pathof(m))), "StdlibArtifacts.toml")

   # If this file exists, it's a stdlib JLL and we must download the artifact ourselves
   if isfile(artifacts_toml)
       # we need to remove the _jll for the artifact name
       meta = artifact_meta(string(m)[1:end-4], artifacts_toml)
       hash = Base.SHA1(meta["git-tree-sha1"])
       if !artifact_exists(hash)
           dl_info = first(meta["download"])
           download_artifact(hash, dl_info["url"], dl_info["sha256"])
       end
       return artifact_path(hash)
   else
      # Otherwise, we can just use the artifact directory given to us by the module
      return m.artifact_dir
   end
end

function prepare_deps_tree()
   mutable_artifacts_toml = joinpath(dirname(pathof(polymake_jll)), "..", "MutableArtifacts.toml")
   polymake_tree = "polymake_tree"
   polymake_tree_hash = artifact_hash(polymake_tree, mutable_artifacts_toml)

   # create a directory tree for polymake with links to dependencies
   # looking similiar to the tree in the build environment
   # for compiling wrappers at run-time
   polymake_tree_hash = create_artifact() do art_dir
      mkpath(joinpath(art_dir,"deps"))
      deps = [ polymake_jll.FLINT_jll,
               polymake_jll.GMP_jll,
               polymake_jll.MPFR_jll,
               polymake_jll.PPL_jll,
               polymake_jll.Perl_jll,
               polymake_jll.bliss_jll,
               polymake_jll.boost_jll,
               polymake_jll.cddlib_jll,
               polymake_jll.lrslib_jll,
               polymake_jll.normaliz_jll ]
      for dep in deps
         symlink(full_artifact_dir(dep), joinpath(art_dir,"deps","$dep"))
      end
      for dir in readdir(polymake_jll.artifact_dir)
         symlink(joinpath(polymake_jll.artifact_dir,dir), joinpath(art_dir,dir))
      end
   end

   bind_artifact!(mutable_artifacts_toml,
      polymake_tree,
      polymake_tree_hash;
      force=true
   )

   # Point polymake to our custom tree
   return artifact_path(polymake_tree_hash)
end

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
function computer_test()
    @testset "Computer" begin
        comp = current_computer();
        @test (comp.name == :LH)  &&  !is_remote(comp)  &&  run_local()

        compList = deepcopy(FilesLH.defaultCompList);
        compNew = Computer(name = :newComp, homeDir = "/Users/lutz/Documents",
            isRemote = false);

        # Add new computer that does not exist
        doReplace = add_computer(compNew, compList = compList);
        @test !doReplace
        # Add again, now it exists
        doReplace = add_computer(compNew, compList = compList);
        @test doReplace

        # Remove computer
        compExists = remove_computer(compNew, compList = compList);
        @test compExists
        # Remove again; it no longer exists
        compExists = remove_computer(compNew, compList = compList);
        @test !compExists
    end
end

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


text:
module MeasureToolbox # TODO make this an independent package

# Import the necessary packages.
import FastGaussQuadrature
import JuMP
import Distributions
const JuMPC = JuMP.Containers
using ..InfiniteOpt

# include jl files here
include("integrals.jl")
include("expectations.jl")
include("support_sums.jl")

end # end module


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra, Random
using Distributions, PDMats, StatsBase

struct test_dist <: Distribution{Univariate, Continuous} end
Distributions.sampler(d::test_dist) = Distributions.sampler(Distributions.Normal(0,1))

struct test_batsampler{T} <: BAT.BATSampler{T, Continuous} end

function Random.rand!(rng::AbstractRNG, s::test_batsampler, x::Integer)
    return 0.5
end

Base.length(s::test_batsampler{Multivariate}) = 2
function Random.rand!(rng::AbstractRNG, s::test_batsampler, x::AbstractVector{<:Real})
    for i in eachindex(x)
        x[i] = 1
    end
    return x
end

@testset "bat_sampler" begin
    @test BAT.bat_sampler(@inferred test_dist()) == Distributions.Normal(0,1)

    @test rand!(test_batsampler{Univariate}(), 1) == 0.5
    x = zeros(2,3)
    rand!(test_batsampler{Multivariate}(), x)
    @test x == ones(2,3)
end


text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
######################################################################################################
using JLD2;
using Distributed;
using Printf;
@everywhere include("runit.jl");
@everywhere include("../utilities/thresholds.jl");
include("../utilities/experiment_helpers.jl");

# setup of the 2nd experiment from Optimal Best Arm Identification with Fixed Confidence (Garivier and Kaufmann 2016)
dist = Bernoulli();
μ = [.3, .21, .2, .19, .18];
δs = (0.1, 0.01, 0.001, 0.0001);
βs = GK16.(δs); # Recommended in Section 6 of Wouter Koolen's paper
N = 3000;
seed = 1234;
pep = BestArm(dist);
# sampling rules to be compared
srs = [
    FWSampling(),
    TrackAndStop(DTracking),
    DaBomb(CTracking, nanswers(pep, μ)),
    Menard(CTracking, 1/oracle(pep, μ)[1]),
    OptimisticTrackAndStop(CTracking),
    RoundRobin()
];

println("μ=$μ, N=$N");
# compute
@time data = pmap(
    ((sr,i),) -> runit(seed+i, sr, μ, pep, βs),
    Iterators.product(srs, 1:N)
);

dump_stats(pep, μ, δs, βs, srs, data, N);
# save
@save isempty(ARGS) ? "BAI1.dat" : ARGS[1]  dist μ pep srs data δs βs N seed

# visualise by loading viz_bai.jl


text:
function associate(state::SlamState, z, R, gate1, gate2)
    #
    # Simple gated nearest-neighbour data-association. No clever feature
    # caching tricks to speed up association, so computation is O(N), where
    # N is the number of features in the state.

    x = state.x
    P = state.cov
    n_observables = size(z, 1)

    zf = Array{AbstractFloat}(n_observables, 0)
    zn = Array{AbstractFloat}(n_observables, 0)
    idf = Array{Int}(1, 0)

    Nxv = 3 # number of vehicle pose components [x, y, phi]
    Nf = round(Int, (length(x) - Nxv)/2) # number of features already in map

    # linear search for nearest-neighbour, no clever tricks (like a quick
    # bounding-box threshold to remove distant features; or, better yet,
    # a balanced k-d tree lookup). TODO: implement clever tricks.
    for i = 1:size(z,2)
        jbest = 0
        nbest = Inf
        outer = Inf

        # search for neighbours
        for j = 1:Nf
            nis, nd = compute_association(x, P, z[:,i], R, j)
            ingate = 0
            if nis < gate1
                if nd < nbest # if within gate, store nearest-neighbour
                    ingate = 1
                end
            end
            if ingate == 1 # if within gate, store nearest-neighbour
                nbest = nd
                jbest = j
            elseif nis < outer # else store best nis value
                outer = nis
            end
        end
        # add nearest-neighbour to association list
        if jbest != 0
            zf = [zf  z[:,i]]
            idf = [idf jbest]
        elseif outer > gate2 # z too far to associate, but far enough to be a new feature
            zn = [zn z[:,i]]
        end
    end
    zf, idf, zn
end

function compute_association(x, P, z, R, idf)
    # Return normalised innovation squared (Mahalanobis distance) and
    # normalised distance
    zp, H = predict_observation(x, idf)
    v = z - zp
    v[2] = mpi_to_pi(v[2])
    S = H*P*H' + R
    nis = dot(v, inv(S)*v)
    nd = nis + log(det(S))
    nis, nd
end


text:
### estimators for π₀

"""
    estimate_pi0(PValues, Pi0Estimator)

Estimate π₀, the fraction of tests under the null hypothesis

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, StoreyBootstrap())
0.0
julia> estimate_pi0(pvals, FlatGrenander())
0.42553191489361697
```

```jldoctest
julia> subtypes(Pi0Estimator)
10-element Array{Union{DataType, UnionAll},1}:
 MultipleTesting.BUM
 MultipleTesting.CensoredBUM
 MultipleTesting.ConvexDecreasing
 MultipleTesting.FlatGrenander
 MultipleTesting.LeastSlope
 MultipleTesting.Oracle
 MultipleTesting.RightBoundary
 MultipleTesting.Storey
 MultipleTesting.StoreyBootstrap
 MultipleTesting.TwoStep
```

# See also

`Pi0Estimator`s:

[`Storey`](@ref)
[`StoreyBootstrap`](@ref)
[`LeastSlope`](@ref)
[`Oracle`](@ref)
[`TwoStep`](@ref)
[`RightBoundary`](@ref)
[`CensoredBUM`](@ref)
[`BUM`](@ref)
[`FlatGrenander`](@ref)
[`ConvexDecreasing`](@ref)

"""
function estimate_pi0 end

function estimate_pi0(pValues::AbstractVector{T}, method::M) where {T<:AbstractFloat, M<:Pi0Estimator}
    estimate_pi0(PValues(pValues), method)
end

## Storey estimator

"""
Storey's π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, Storey())
0.22222222222222224

julia> estimate_pi0(pvals, Storey(0.4))
0.33333333333333337
```

# References

Storey, J.D., Taylor, J.E., and Siegmund, D. (2004). Strong control,
conservative point estimation and simultaneous conservative consistency of false
discovery rates: a unified approach. Journal of the Royal Statistical Society:
Series B (Statistical Methodology) 66, 187–205.

"""
struct Storey <: Pi0Estimator
    λ::Float64

    function Storey(λ)
        isin(λ, 0, 1) || throw(DomainError("λ must be in [0, 1]"))
        return new(λ)
    end
end

Storey() = Storey(0.1)

function estimate_pi0(pValues::PValues{T}, pi0estimator::Storey) where T<:AbstractFloat
    lambda = pi0estimator.λ
    pi0 = (sum(pValues .>= lambda) / length(pValues)) / (1 - lambda)
    pi0 = clamp(pi0, 0, 1)
    return pi0
end

## Storey bootstrap estimator

"""
Storey's closed-form bootstrap π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, StoreyBootstrap())
0.0

julia> estimate_pi0(pvals, StoreyBootstrap(0.1:0.1:0.9, 0.2))
0.0
```

# References

Robinson, D. (2016). Original Procedure for Choosing λ.
http://varianceexplained.org/files/pi0boot.pdf

Storey, J.D., Taylor, J.E., and Siegmund, D. (2004). Strong control,
conservative point estimation and simultaneous conservative consistency of false
discovery rates: a unified approach. Journal of the Royal Statistical Society:
Series B (Statistical Methodology) 66, 187–205.

"""
struct StoreyBootstrap <: Pi0Estimator
    λseq::Vector{Float64}
    q   ::Float64

    function StoreyBootstrap(λseq, q)
        isin(λseq, 0, 1) || throw(DomainError("λseq must be in [0, 1]"))
        isin(q, 0, 1) || throw(DomainError("q must be in [0, 1]"))
        return new(λseq, q)
    end
end

StoreyBootstrap() = StoreyBootstrap(0.05:0.05:0.95, 0.1)

function estimate_pi0(pValues::PValues{T}, pi0estimator::StoreyBootstrap) where T<:AbstractFloat
    lambdas = pi0estimator.λseq
    q = pi0estimator.q
    n = length(pValues)
    w = [sum(pValues .>= l) for l in lambdas]  # TODO: check if >= or >
    pi0 = w ./ n ./ (1 .- lambdas)
    min_pi0 = quantile(pi0, q)
    mse = (w ./ (n.^2 .* (1 .- lambdas).^2 )) .* (1 .- w/n) + (pi0 .- min_pi0).^2
    pi0 = clamp(pi0[argmin(mse)], 0, 1)
    return pi0
end

## Least SLope (LSL) estimator

"""
Least Slope (LSL) π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, LeastSlope())
1.0
```

# References

Benjamini, Y., and Hochberg, Y. (2000). On the Adaptive Control of the False
Discovery Rate in Multiple Testing With Independent Statistics. Journal of
Educational and Behavioral Statistics 25, 60–83.

"""
struct LeastSlope <: Pi0Estimator
end

function estimate_pi0(pValues::PValues{T}, pi0estimator::LeastSlope) where T<:AbstractFloat
    n = length(pValues)
    pValues = sort_if_needed(pValues)
    s0 = lsl_slope(1, n, pValues)
    sx = 0
    for i in 2:n
        s1 = lsl_slope(i, n, pValues)
        if (s1 - s0) < 0
            sx = s1
            break
        end
        s0 = s1
    end
    pi0 = min( 1/sx + 1, n ) / n
    return pi0
end

function lsl_slope(i::Integer, n::Integer, pValue::AbstractVector{T}) where T<:AbstractFloat
    s = (1 - pValue[i]) / (n - i + 1)
    return s
end

# alternative, vectorized version
# used for comparison and compactness
function lsl_pi0_vec(pValues::AbstractVector{T}) where T<:AbstractFloat
    n = length(pValues)
    pValues = sort_if_needed(pValues)
    s = (1 .- pValues) ./ (n:-1:1)
    d = diff(s) .< 0
    idx = something(findfirst(d), 0) + 1
    pi0 = min( 1/s[idx] + 1, n ) / n
    return pi0
end

## Oracle

"""
Oracle π₀

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, Oracle(0.5)) # a bit boring...
0.5
```
"""
struct Oracle <: Pi0Estimator
    π0::Float64

    function Oracle(π0)
        isin(π0, 0, 1) || throw(DomainError("π0 must be in [0, 1]"))
        return new(π0)
    end
end

Oracle() = Oracle(1.0)

function estimate_pi0(pValues::PValues{T}, pi0estimator::Oracle) where T<:AbstractFloat
    pi0estimator.π0
end

## Two-Step estimator: Benjamini, Krieger and Yekutieli (2006)

"""
Two-step π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, TwoStep())
0.2

julia> estimate_pi0(pvals, TwoStep(0.05, BenjaminiLiu()))
0.2

```

# References

Benjamini, Y., Krieger, A.M., and Yekutieli, D. (2006). Adaptive linear step-up
procedures that control the false discovery rate. Biometrika 93, 491–507.

"""
struct TwoStep <: Pi0Estimator
    α::Float64
    adjustment::PValueAdjustment

    function TwoStep(α, method)
        isin(α, 0, 1) || throw(DomainError("α must be in [0, 1]"))
        return new(α, method)
    end
end

TwoStep() = TwoStep(0.05)

TwoStep(α) = TwoStep(α, BenjaminiHochberg())

function estimate_pi0(pValues::PValues{T}, pi0estimator::TwoStep) where T<:AbstractFloat
    alpha = pi0estimator.α
    padj = adjust(pValues, pi0estimator.adjustment)
    pi0 = sum(padj .>= (alpha/(1+alpha))) / length(padj)
    return pi0
end

# RightBoundary procedure as defined in Definition 2 of Liang and Nettleton, 2012

"""
Right boundary π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, RightBoundary())
0.2127659574468085

julia> estimate_pi0(pvals, RightBoundary(0.1:0.1:0.9))
0.25
```

# References

Liang, K., and Nettleton, D. (2012). Adaptive and dynamic adaptive procedures
for false discovery rate control and estimation. Journal of the Royal
Statistical Society: Series B (Statistical Methodology) 74, 163–182.

"""
struct RightBoundary <: Pi0Estimator
    λseq::Vector{Float64}

    function RightBoundary(λseq)
        isin(λseq, 0, 1) || throw(DomainError("λseq must be in [0, 1]"))
        return new(λseq)
    end
end

# λseq used in Liang, Nettleton 2012
RightBoundary() = RightBoundary([0.02:0.02:0.1; 0.15:0.05:0.95])

function estimate_pi0(pValues::PValues{T}, pi0estimator::RightBoundary) where T<:AbstractFloat
    n = length(pValues)
    λseq = sort_if_needed(pi0estimator.λseq)
    # make sure we catch p-values equal to 1 despite left closure
    # use closed=:left because we have been using >= convention in this package
    # note that original paper uses > convention.
    h = fit(Histogram, pValues, [λseq; Inf], closed=:left)
    pi0_estimates = reverse(cumsum(reverse(h.weights)))./(1 .- λseq)./n
    pi0_decrease = diff(pi0_estimates) .>= 0
    pi0_decrease[end] = true
    pi0 = pi0_estimates[something(findfirst(pi0_decrease), 0) + 1]
    pi0 = clamp(pi0, 0, 1)
    return pi0
end

## Censored BUM

"""
Censored Beta-Uniform Mixture (censored BUM) π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, CensoredBUM())
0.21052495526400936

```

# References

Markitsis, A., and Lai, Y. (2010). A censored beta mixture model for the
estimation of the proportion of non-differentially expressed genes.
Bioinformatics 26, 640–646.

"""
struct CensoredBUM <: Pi0Estimator
    γ0::Float64
    λ::Float64
    xtol::Float64
    maxiter::Int64

    function CensoredBUM(γ0, λ, xtol, maxiter)
        isin(γ0, 0, 1) || throw(DomainError("γ0 must be in [0, 1]"))
        isin(λ, 0, 1) || throw(DomainError("λ must be in [0, 1]"))
        isin(xtol, 0, 1) || throw(DomainError("xtol must be in [0, 1]"))
        maxiter > 0 || throw(DomainError("maxiter must be a positive number"))
        return new(γ0, λ, xtol, maxiter)
    end
end

CensoredBUM() = CensoredBUM(0.5, 0.05, 1e-6, 10000)
CensoredBUM(γ0, λ) = CensoredBUM(γ0, λ, 1e-6, 10000)

struct CensoredBUMFit <: Pi0Fit
    π0::Float64
    param::Vector{Float64}
    is_converged::Bool
end

function fit(pi0estimator::CensoredBUM, pValues::AbstractVector{T}; kw...) where T<:AbstractFloat
    π0, param, is_converged = cbum_pi0(pValues, pi0estimator.γ0, pi0estimator.λ,
                                       pi0estimator.xtol, pi0estimator.maxiter;
                                       kw...)
    return CensoredBUMFit(π0, param, is_converged)
end

function estimate_pi0(pValues::PValues{T}, pi0estimator::CensoredBUM) where T<:AbstractFloat
    estimate_pi0(fit(pi0estimator, pValues))
end

function estimate_pi0(pi0fit::CensoredBUMFit)
    pi0 = pi0fit.is_converged ? pi0fit.π0 : NaN
    return pi0
end

function cbum_pi0(pValues::AbstractVector{T},
                  γ0::AbstractFloat = 0.5, λ::AbstractFloat = 0.05,
                  xtol::AbstractFloat = 1e-6, maxiter::Integer = 10000) where T<:AbstractFloat
    n = length(pValues)
    idx_right = pValues .>= λ
    n2 = sum(idx_right)
    n1 = n - n2
    sz = (1-γ0)*n
    szr = (1-γ0)*n2
    szl = sz - szr
    pr = pValues[idx_right]
    lpr = log.(pr)
    ll = log(λ)
    zr = fill(1-γ0, n2)
    pi0_old = γ0 = α = γ = Inf
    for i in 1:maxiter
        γ = 1 - sz/n
        α = -szr / ( ll * szl + sum(zr .* lpr) )
        γ = clamp(γ, 0, 1)
        α = clamp(α, 0, 1)
        xl = (1-γ) * (λ^α)
        szl = (xl ./ (γ*λ + xl)) * n1
        xr = (1-γ) * α * pr.^(α-1)
        zr = xr ./ (γ .+ xr)
        szr = sum(zr)
        sz = szl + szr
        pi0_new = γ + (1-γ)*α
        if abs(pi0_new - pi0_old) <= xtol
            return pi0_new, [γ, α], true
        end
        γ0 = γ
        pi0_old = pi0_new
    end
    return NaN, [γ, α], false
end

function cbum_pi0_naive(pValues::AbstractVector{T},
                        γ0::AbstractFloat = 0.5, λ::AbstractFloat = 0.05,
                        xtol::AbstractFloat = 1e-6, maxiter::Integer = 10000) where T<:AbstractFloat
    n = length(pValues)
    z = fill(1-γ0, n)
    idx_left = pValues .< λ
    idx_right = .!idx_left
    pi0_old = γ0 = α = γ = Inf
    # compute constant values only once
    lpr = log.(pValues[idx_right])
    ll = log(λ)
    for i in 1:maxiter
        γ = sum(1 .- z) / n
        α = -sum(z[idx_right])
        α = α / ( ll * sum(z[idx_left]) + sum(z[idx_right] .* lpr) )
        xl = (1-γ) * (λ^α)
        z[idx_left] .= xl / (γ*λ + xl)
        xr = (1-γ) * α * pValues[idx_right].^(α-1)
        z[idx_right] = xr ./ (γ .+ xr)
        pi0_new = γ + (1-γ)*α
        if abs(pi0_new - pi0_old) <= xtol
            return pi0_new, [γ, α], true
        end
        γ0 = γ
        pi0_old = pi0_new
    end
    return NaN, [γ, α], false
end

## BUM

"""
Beta-Uniform Mixture (BUM) π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, BUM())
0.22802795505154264
```

# References

Pounds, S., and Morris, S.W. (2003). Estimating the occurrence of false
positives and false negatives in microarray studies by approximating and
partitioning the empirical distribution of p-values. Bioinformatics 19,
1236–1242.

"""
struct BUM <: Pi0Estimator
    γ0::Float64
    xtol::Float64
    maxiter::Int64

    function BUM(γ0, xtol, maxiter)
        isin(γ0, 0, 1) || throw(DomainError("γ0 must be in [0, 1]"))
        isin(xtol, 0, 1) || throw(DomainError("xtol must be in [0, 1]"))
        maxiter > 0 || throw(DomainError("maxiter must be a positive number"))
        return new(γ0, xtol, maxiter)
    end
end

BUM() = BUM(0.5, 1e-6, 10000)

BUM(y0::Float64) = BUM(y0, 1e-6, 10000)

struct BUMFit <: Pi0Fit
    π0::Float64
    param::Vector{Float64}
    is_converged::Bool
end

function fit(pi0estimator::BUM, pValues::AbstractVector{T}; kw...) where T<:AbstractFloat
    π0, param, is_converged = cbum_pi0(pValues, pi0estimator.γ0, eps(),
                                       pi0estimator.xtol, pi0estimator.maxiter;
                                       kw...)
    return BUMFit(π0, param, is_converged)
end

function estimate_pi0(pValues::PValues{T}, pi0estimator::BUM) where T<:AbstractFloat
    estimate_pi0(fit(pi0estimator, pValues))
end

function estimate_pi0(pi0fit::BUMFit)
    pi0 = pi0fit.is_converged ? pi0fit.π0 : NaN
    return pi0
end

## Longest constant interval in the Grenander estimator

"""
Flat Grenander π₀ estimator

Estimates π₀ by finding the longest constant interval in the Grenander estimator.

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, FlatGrenander())
0.42553191489361697
```

# References

Langaas, M., Lindqvist, B.H., and Ferkingstad, E. (2005). Estimating the
proportion of true null hypotheses, with application to DNA microarray data.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67,
555–572.

"""
struct FlatGrenander <: Pi0Estimator
end

function estimate_pi0(pValues::PValues{T}, pi0estimator::FlatGrenander) where T<:AbstractFloat
    p, f, F = grenander(pValues)
    pi0 = longest_constant_interval(p, f)
    return pi0
end

function longest_constant_interval(p::AbstractVector{T}, f::AbstractVector{T}) where T<:AbstractFloat
    p = [0.0; p]
    f = [Inf; f]
    i2 = length(f)
    Δp_max = Δp = 0.0
    pi0 = 1
    for i1 in (length(f)-1):-1:1
        if f[i2] ≈ f[i1] # within constant interval
            Δp = p[i2] - p[i1]
        else
            if Δp >= Δp_max
                Δp_max = Δp
                pi0 = f[i2]
            end
            i2 = i1
            Δp = 0.0
        end
        if f[i1] > 1
            break
        end
    end
    return pi0
end

## Convex Decreasing π₀ estimator

"""
Convex Decreasing π₀ estimator

# Examples

```jldoctest
julia> pvals = PValues([0.001, 0.002, 0.01, 0.03, 0.5]);

julia> estimate_pi0(pvals, ConvexDecreasing())
0.013007051336745304
```

# References

Langaas, M., Lindqvist, B.H., and Ferkingstad, E. (2005). Estimating the
proportion of true null hypotheses, with application to DNA microarray data.
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67,
555–572.

"""
struct ConvexDecreasing <: Pi0Estimator
    gridsize::Int64
    xtol::Float64
    maxiter::Int64

    function ConvexDecreasing(gridsize, xtol, maxiter)
        gridsize > 0 || throw(DomainError("gridsize must be a positive number"))
        isin(xtol, 0, 1) || throw(DomainError("xtol must be in [0, 1]"))
        maxiter > 0 || throw(DomainError("maxiter must be a positive number"))
        return new(gridsize, xtol, maxiter)
    end
end

ConvexDecreasing() = ConvexDecreasing(100, 1e-5, 10000)

struct ConvexDecreasingFit <: Pi0Fit
    π0::Float64
    param::Vector{Float64}
    is_converged::Bool
end

function fit(pi0estimator::ConvexDecreasing, pValues::AbstractVector{T}) where T<:AbstractFloat
    π0, param, is_converged = convex_decreasing(pValues,
                                                pi0estimator.gridsize,
                                                pi0estimator.xtol,
                                                pi0estimator.maxiter)
    return ConvexDecreasingFit(π0, param, is_converged)
end

function estimate_pi0(pValues::PValues{T}, pi0estimator::ConvexDecreasing) where T<:AbstractFloat
    estimate_pi0(fit(pi0estimator, pValues))
end

function estimate_pi0(pi0fit::ConvexDecreasingFit)
    pi0 = pi0fit.is_converged ? pi0fit.π0 : NaN
    return pi0
end

function convex_decreasing(pValues::AbstractVector{T},
                           gridsize::Integer = 100,
                           xtol::AbstractFloat = 1e-5,
                           maxiter::Integer = 10000) where T<:AbstractFloat

    n = length(pValues)
    p = sort(pValues)
    dx = 1 / gridsize
    t = collect(dx:dx:1)
    x = collect(0:dx:1)
    f = fill(one(T), gridsize+1)
    f_p = fill(one(T), n)
    theta = dx * find_theta(t, p)
    f_theta = triangular_weighting(x, theta)
    f_theta_p = triangular_weighting(p, theta)
    idx_lower = round.(Int, round.(gridsize .* p, RoundDown) .+ 1)
    p_upper = round.(gridsize .* p, RoundUp)/gridsize
    idx_upper = round.(Int, gridsize .* p_upper) .+ 1
    px = p_upper .- p
    thetas = T[]
    pi0_new = pi0_old = Inf
    for j in 1:maxiter
        if sum((f_p.-f_theta_p)./f_p) > 0
            ε = 0.0
        else
            l = 0.0
            u = 1
            while abs(u-l) > xtol
                ε = (l+u)/2.0
                if decide(f_p, f_theta_p, ε) < 0
                    l = ε
                else
                    u = ε
                end
            end
        end
        @. f = (1-ε)*f + ε*f_theta
        @. f_p = (1-ε)*f_p + ε*f_theta_p
        theta = dx * find_theta(t, p, f_p)
        f_theta .= triangular_weighting(x, theta)
        pi0_new = f[end]
        if abs(pi0_new - pi0_old) <= xtol
            return pi0_new, thetas, true
        end
        pi0_old = pi0_new
        f_theta_p .= triangular_weighting(p, theta)
        if sum(f_theta_p ./ f_p) < sum(1 ./ f_p)
            theta = 0.0
            f_theta .= fill(1, size(f_theta))
            f_theta_p .= fill(1, size(f_theta_p))
        end
        if !(theta in thetas)
            append!(thetas, theta)
            sort!(thetas)
        end
    end
    return NaN, thetas, false
end

function find_theta(t::Vector{Float64}, p::Vector{Float64})
    return argmax( [theta.^-2 * sum(theta .- p[p .< theta]) for theta in t] )
end

function find_theta(t::Vector{Float64}, p::Vector{Float64}, f_p::Vector{Float64})
    return argmax( [theta.^-2 * sum( (theta .- p) .* (p .< theta) ./ f_p ) for theta in t] )
end

function decide(f_p::Vector{Float64}, f_theta_p::Vector{Float64}, ε::Float64)
    idx = f_p .> 0
    return sum( @. ( f_p[idx] - f_theta_p[idx] ) / ( (1-ε) * f_p[idx] + ε * f_theta_p[idx] ) )
end

function triangular_weighting(x::Vector{T}, mid::T) where T<:AbstractFloat
    return @. 2 / mid^2 * (mid-x) * (x<mid)
end