text:
module MIPforQP

import JuMP, MathOptInterface
import Mosek, MosekTools
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const VI = MOI.VariableIndex
const SV = MOI.SingleVariable

import LinearAlgebra, SparseArrays
import Gurobi

abstract type AbstractRelaxation end

struct HongboRelaxation <: AbstractRelaxation
    num_layers::Int
end
struct NeuralNetRelaxation <: AbstractRelaxation
    num_layers::Int
end
struct Reformulation{T <: AbstractRelaxation}
    relaxation::T
    soc_lower_bound::Bool
end
Reformulation(relaxation::AbstractRelaxation) = Reformulation(relaxation, true)

function formulate_quadratics!(model::JuMP.Model, reformulation::Reformulation)
    vq_ci, vq_f, vq_s = vectorize_quadratics!(model)
    # The following is a hacky way to determine if the objective is linear or quadratic.
    obj = try
        MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        nothing
    catch InexactError
        MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    end
    if obj === nothing
        MOI.delete(JuMP.backend(model), vq_ci)
        bridge = _relax_quadratic!(JuMP.backend(model), reformulation, vq_f, vq_s)
    else
        bridge = _relax_quadratic!(JuMP.backend(model), reformulation, vq_f, vq_s, obj)
    end
    MOIU.attach_optimizer(model)
    return bridge
end

function vectorize_quadratics!(model::JuMP.Model)
    moi_backend = JuMP.backend(model)
    @assert typeof(moi_backend) <: MOIU.CachingOptimizer
    @assert moi_backend.state in (MOIU.NO_OPTIMIZER, MOIU.EMPTY_OPTIMIZER)
    return vectorize_quadratics!(moi_backend.model_cache.model)
end

function vectorize_quadratics!(cached_model::MOIU.Model{T}) where {T}
    scalar_quad_cons = cached_model.moi_scalarquadraticfunction
    @assert isempty(scalar_quad_cons.moi_equalto)
    @assert isempty(scalar_quad_cons.moi_interval)
    @assert isempty(scalar_quad_cons.moi_semicontinuous)
    @assert isempty(scalar_quad_cons.moi_semiinteger)
    @assert isempty(scalar_quad_cons.moi_zeroone)

    num_gt = length(scalar_quad_cons.moi_greaterthan)
    num_lt = length(scalar_quad_cons.moi_lessthan)

    constants = Array{T}(undef, 0)
    affine_terms = Array{MOI.VectorAffineTerm{T}}(undef, 0)
    quad_terms = Array{MOI.VectorQuadraticTerm{T}}(undef, 0)

    func_count = 0
    for (q_ci, q_f, q_s) in scalar_quad_cons.moi_greaterthan
        func_count += 1
        push!(constants, -(q_f.constant - q_s.lower))
        for aff in q_f.affine_terms
            push!(affine_terms, MOI.VectorAffineTerm{T}(func_count, MOI.ScalarAffineTerm{T}(-aff.coefficient, aff.variable_index)))
        end
        for quad in q_f.quadratic_terms
            push!(quad_terms, MOI.VectorQuadraticTerm{T}(func_count, MOI.ScalarQuadraticTerm{T}(-quad.coefficient, quad.variable_index_1, quad.variable_index_2)))
        end
    end

    @assert func_count == num_gt

    for (q_ci, q_f, q_s) in scalar_quad_cons.moi_lessthan
        func_count += 1
        push!(constants, q_f.constant - q_s.upper)
        for aff in q_f.affine_terms
            push!(affine_terms, MOI.VectorAffineTerm{T}(func_count, MOI.ScalarAffineTerm{T}(aff.coefficient, aff.variable_index)))
        end
        for quad in q_f.quadratic_terms
            push!(quad_terms, MOI.VectorQuadraticTerm{T}(func_count, MOI.ScalarQuadraticTerm{T}(quad.coefficient, quad.variable_index_1, quad.variable_index_2)))
        end
    end

    @assert func_count == num_gt + num_lt

    vector_quad_cons = cached_model.moi_vectorquadraticfunction.moi_nonpositives
    @assert isempty(vector_quad_cons)
    constraint_index = MOI.ConstraintIndex{MOI.VectorQuadraticFunction{T},MOI.Nonpositives}(1)
    vector_func = MOI.VectorQuadraticFunction{T}(affine_terms, quad_terms, constants)
    vector_set = MOI.Nonpositives(num_gt + num_lt)
    retval = (constraint_index, vector_func, vector_set)
    push!(vector_quad_cons, retval)
    empty!(scalar_quad_cons.moi_greaterthan)
    empty!(scalar_quad_cons.moi_lessthan)
    return retval
end

struct MIPforQPBridge
    reformulation::Reformulation
    squared_vars::Dict{VI,VI}
    binaries::Dict{VI,Union{Nothing,Vector{VI}}}
end
MIPforQPBridge(reformulation::Reformulation) = MIPforQPBridge(reformulation, Dict{VI,VI}(), Dict{VI,Vector{VI}}())

function _indices_and_coefficients(
    I::AbstractVector{Int}, J::AbstractVector{Int}, V::AbstractVector{Float64},
    indices::AbstractVector{Int}, coefficients::AbstractVector{Float64},
    f::MOI.ScalarQuadraticFunction, canonical_index::Dict{VI,Int}
)
    variables_seen = 0
    for (i, term) in enumerate(f.quadratic_terms)
        vi_1 = term.variable_index_1
        vi_2 = term.variable_index_2
        if !haskey(canonical_index, vi_1)
            variables_seen += 1
            canonical_index[vi_1] = variables_seen
        end
        if !haskey(canonical_index, vi_2)
            variables_seen += 1
            canonical_index[vi_2] = variables_seen
        end

        I[i] = canonical_index[term.variable_index_1]
        J[i] = canonical_index[term.variable_index_2]
        V[i] =  term.coefficient
        # Gurobi returns a list of terms. MOI requires 0.5 x' Q x. So, to get
        # from
        #   Gurobi -> MOI => multiply diagonals by 2.0
        #   MOI -> Gurobi => multiply diagonals by 0.5
        # Example: 2x^2 + x*y + y^2
        #   |x y| * |a b| * |x| = |ax+by bx+cy| * |x| = 0.5ax^2 + bxy + 0.5cy^2
        #           |b c|   |y|                   |y|
        #   Gurobi needs: (I, J, V) = ([0, 0, 1], [0, 1, 1], [2, 1, 1])
        #   MOI needs:
        #     [SQT(4.0, x, x), SQT(1.0, x, y), SQT(2.0, y, y)]
        if I[i] == J[i]
            V[i] *= 0.5
        end
    end
    for (i, term) in enumerate(f.affine_terms)
        indices[i] = term.variable_index.value
        coefficients[i] = term.coefficient
    end
    return
end

function _indices_and_coefficients(f::MOI.ScalarQuadraticFunction, canonical_index::Dict{VI,Int})
    f_canon = MOI.Utilities.canonical(f)
    nnz_quadratic = length(f_canon.quadratic_terms)
    nnz_affine = length(f_canon.affine_terms)
    I = Vector{Int}(undef, nnz_quadratic)
    J = Vector{Int}(undef, nnz_quadratic)
    V = Vector{Float64}(undef, nnz_quadratic)
    indices = Vector{Int}(undef, nnz_affine)
    coefficients = Vector{Float64}(undef, nnz_affine)
    _indices_and_coefficients(I, J, V, indices, coefficients, f_canon, canonical_index)
    return indices, coefficients, I, J, V
end

function _get_bounds(model, x::VI)
    c_lt = MOI.ConstraintIndex{MOI.SingleVariable,MOI.LessThan{Float64}}(x.value)
    c_gt = MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{Float64}}(x.value)
    c_int = MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{Float64}}(x.value)
    if MOI.is_valid(model, c_int)
        @assert !MOI.is_valid(c_lt) && !MOI.is_valid(c_gt)
        int::MOI.Interval{Float64} = MOI.get(model, MOI.ConstraintSet(), c_int)
        return int.lower, int.upper
    elseif MOI.is_valid(model, c_lt)
        lt::MOI.LessThan{Float64} = MOI.get(model, MOI.ConstraintSet(), c_lt)
        if MOI.is_valid(model, c_gt)
            gt_1::MOI.GreaterThan{Float64} = MOI.get(model, MOI.ConstraintSet(), c_gt)
            return gt_1.lower, lt.upper
        else
            return -Inf, lt.upper
        end
    elseif MOI.is_valid(model, c_gt)
        gt_2::MOI.GreaterThan{Float64} = MOI.get(model, MOI.ConstraintSet(), c_gt)
        return gt_2.lower, Inf
    else
        return -Inf, Inf
    end
end

function _map_to_interval_if_needed(model, x::VI, y::VI, l::Float64, u::Float64)
    if l == 0 && u == 1
        return x, y
    else
        Δ = u - l
        @assert Δ > 0
        x̃_v, x̃_c = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0))
        ỹ_v, ỹ_c = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0))
        MOIU.normalize_and_add_constraint(model, SV(x) - (l + Δ * SV(x̃_v)), MOI.EqualTo(0.0))
        MOIU.normalize_and_add_constraint(model, SV(y) - (l^2 + 2l * Δ * SV(x̃_v) + Δ^2 * SV(ỹ_v)), MOI.EqualTo(0.0))
        return x̃_v, ỹ_v
    end
end

function reformulate_unary_quadratic_term!(model::MOI.ModelLike, reformulation::Reformulation{HongboRelaxation}, _x::VI, _y::VI, l::Float64, u::Float64)
    num_layers = reformulation.relaxation.num_layers
    ξ_v = MOI.add_variables(model, num_layers)
    η_v = MOI.add_variables(model, num_layers)
    λ_1_v, λ_1_c = MOI.add_constrained_variables(model, [MOI.GreaterThan(0.0) for i in 1:num_layers])
    λ_2_v, λ_2_c = MOI.add_constrained_variables(model, [MOI.GreaterThan(0.0) for i in 1:num_layers])
    z_v, z_c = MOI.add_constrained_variables(model, [MOI.ZeroOne() for i in 1:num_layers])
    x = SV(_x)
    y = SV(_y)
    ξ = [SV(ξ_v[i]) for i in 1:num_layers]
    η = [SV(η_v[i]) for i in 1:num_layers]
    λ_1 = [SV(λ_1_v[i]) for i in 1:num_layers]
    λ_2 = [SV(λ_2_v[i]) for i in 1:num_layers]
    z = [SV(z_v[i]) for i in 1:num_layers]

    θ_min = (if l > 0
        atan(l^2 - 1) / (2l)
    elseif l == 0
        -π / 2
    else
        atan(l^2 - 1) / (2l) - π
    end)
    θ_max = (if u > 0
        atan(u^2 - 1) / (2u)
    elseif u == 0
        -π / 2
    else
        atan(u^2 - 1) / (2u) - π
    end)
    @assert -3π / 2 < θ_min < π / 2
    @assert -3π / 2 < θ_max < π / 2
    θ_d = θ_max - θ_min
    θ_mid = (θ_max + θ_min) / 2

    ν = num_layers
    # RLP inequality (1)
    MOIU.normalize_and_add_constraint(model, y - (l + u) * x + l * u, MOI.LessThan(0.0))
    # Equation (11)
    MOIU.normalize_and_add_constraint(model, ξ[ν] * cos(θ_d / 2^(ν + 1)) + η[ν] * sin(θ_d / 2^(ν + 1)) - (y + 1.0) / 2.0 * cos(θ_d / 2^(ν + 1)), MOI.GreaterThan(0.0))
    if !reformulation.soc_lower_bound
        # Equation 12
        MOIU.normalize_and_add_constraint(model, ξ[ν] * cos(θ_d / 2^ν) + η[ν] * sin(θ_d / 2^ν) - (y + 1.0) / 2.0, MOI.LessThan(0.0))
        # Equation 13
        MOIU.normalize_and_add_constraint(model, ξ[ν] - (y + 1.0) / 2.0, MOI.LessThan(0.0))
    end
    # Equation block (20)
    MOIU.normalize_and_add_constraint(model, ξ[1] - x * cos(θ_mid) - (y / 2.0 - 0.5) * sin(θ_mid), MOI.EqualTo(0.0))
    C = max(l^2, u^2) / 2 + 0.5
    MOIU.normalize_and_add_constraint(model, (1.0λ_2[1] - 1.0λ_1[1]) * C + x * sin(θ_mid) - (y / 2.0 - 0.5) * cos(θ_mid), MOI.EqualTo(0.0))
    MOI.add_constraint(model, η[1] - (1.0λ_1[1] + 1.0λ_2[1]) * C, MOI.EqualTo(0.0))
    MOIU.normalize_and_add_constraint(model, 1.0λ_1[1] - (1.0 - z[1]), MOI.LessThan(0.0))
    MOI.add_constraint(model, 1.0λ_2[1] - 1.0z[1], MOI.LessThan(0.0))
    for j in 1:(ν - 1)
        # Equation block 21
        MOI.add_constraint(model, ξ[j + 1] - ξ[j] * cos(θ_d / 2^(j + 1)) - η[j] * sin(θ_d / 2^(j + 1)), MOI.EqualTo(0.0))
        C_j = C * sin(θ_d / 2^(j + 1))
        MOI.add_constraint(model, (1.0λ_2[j + 1] - 1.0λ_1[j + 1]) * C_j + ξ[j] * sin(θ_d / 2^(j + 1)) - η[j] * cos(θ_d / 2^(j + 1)), MOI.EqualTo(0.0))
        MOI.add_constraint(model, η[j + 1] - (1.0λ_1[j + 1] + 1.0λ_2[j + 1]) * C_j, MOI.EqualTo(0.0))
        MOIU.normalize_and_add_constraint(model, 1.0λ_1[j] - (1.0 - z[j + 1]), MOI.LessThan(0.0))
        MOI.add_constraint(model, 1.0λ_2[j] - 1.0z[j + 1], MOI.LessThan(0.0))
    end
end

function reformulate_unary_quadratic_term!(model::MOI.ModelLike, reformulation::Reformulation{NeuralNetRelaxation}, _x::VI, _y::VI, l::Float64, u::Float64)
    x, y = _map_to_interval_if_needed(model, _x, _y, l, u)
    MOI.add_constraint(model, 1.0SV(x) - SV(y), MOI.GreaterThan(0.0))
    return _formulate_zero_one_quadratic!(model, x, y, reformulation.relaxation.num_layers, !reformulation.soc_lower_bound)
end

function _formulate_zero_one_quadratic!(model::MOI.ModelLike, x::VI, y::VI, num_layers::Int, impose_lower_bound::Bool)
    z_v, z_c = MOI.add_constrained_variables(model, [MOI.ZeroOne() for i in 1:num_layers])
    g_v, g_c = MOI.add_constrained_variables(model, [MOI.Interval(0.0, 1.0) for i in 1:num_layers])

    gv = Dict(i => SV(g_v[i]) for i in 1:num_layers)
    gv[0] = SV(x)

    agg_aff = 1.0 * SV(y) - SV(x) + sum(gv[s] / 2.0^(2s) for s in 1:num_layers)

    if impose_lower_bound
        agg_c = MOI.add_constraint(model, agg_aff, MOI.EqualTo(0.0))
    else
        agg_c = MOI.add_constraint(model, agg_aff, MOI.LessThan(0.0))
    end

    for i in 1:num_layers
        α_i = SV(z_v[i])
        MOI.add_constraint(model, gv[i] - 2.0gv[i - 1], MOI.LessThan(0.0))
        MOIU.normalize_and_add_constraint(model, gv[i] + 2.0gv[i - 1] - 2.0, MOI.LessThan(0.0))
        MOI.add_constraint(model, gv[i] + 2.0gv[i - 1] - 2.0 * α_i, MOI.GreaterThan(0.0))
        MOIU.normalize_and_add_constraint(model, gv[i] - 2.0gv[i - 1] + 2.0 * α_i, MOI.GreaterThan(0.0))
    end
    return z_v
end

function _get_Q_matrix(q_func::MOI.ScalarQuadraticFunction, canonical_index::Dict{VI,Int}, n::Int)
    indices, coefficients, I, J, V = _indices_and_coefficients(q_func, canonical_index)
    Q = Matrix{Float64}(SparseArrays.sparse(I, J, V, n, n))
    Q = 1 / 2 * (Q + Q')

    @assert n == size(Q, 1) == size(Q, 2)
    @assert LinearAlgebra.issymmetric(Q)
    return Q
end

function _compute_diagonal_shift(Q::Matrix{Float64})
    n = size(Q, 1)
    @assert n == size(Q, 2)
    sdp_model = JuMP.Model(JuMP.with_optimizer(Mosek.Optimizer))
    δ = JuMP.@variable(sdp_model, [i = 1:n])
    JuMP.@objective(sdp_model, Min, sum(δ))
    JuMP.@SDconstraint(sdp_model, Q + LinearAlgebra.diagm(δ .- 1e-4) >= 0)
    JuMP.optimize!(sdp_model)
    @assert JuMP.primal_status(sdp_model) == MOI.FEASIBLE_POINT
    return JuMP.value.(δ)
end

# Assumption: Constraints we wish to reformulate have been bridged to the form
# `VectorQuadraticFunction`-in-`Nonpositives`.
function _relax_quadratic!(
    model,
    reformulation::Reformulation,
    f::MOI.VectorQuadraticFunction{Float64},
    s::MOI.Nonpositives,
    obj::Union{Nothing,MOI.ScalarQuadraticFunction{Float64}}=nothing
)
    squared_vars = Dict{VI,VI}()
    binaries = Dict{VI,Union{Nothing,Vector{VI}}}()
    # Walk all constraints in vector and populate squared_vars dict
    sqt_terms = vcat([t -> t.scalar_term for t in f.quadratic_terms], obj.quadratic_terms)
    for st in sqt_terms
        for x in (st.variable_index_1, st.variable_index_2)
            if !haskey(squared_vars, x)
                l, u = _get_bounds(model, x)
                y_lb = l < 0 < u ? 0.0 : min(l^2, u^2)
                y_ub = max(l^2, u^2)
                y_v, y_c = MOI.add_constrained_variable(model, MOI.Interval{Float64}(y_lb, y_ub))
                squared_vars[x] = y_v
                if reformulation.soc_lower_bound
                    # Alternatively, could pass SOC constraint with
                    # MOI.add_constraint(model, MOI.VectorOfVariables([y_v, x]), MOI.SecondOrderCone(2))
                    MOI.add_constraint(model, 1.0SV(y_v) - 1.0SV(x) * SV(x), MOI.GreaterThan(0.0))
                end
                z_v = reformulate_unary_quadratic_term!(model, reformulation, x, y_v, l, u)
                binaries[x] = z_v
            end
        end
    end

    n = length(keys(squared_vars))
    q_funcs = MOIU.scalarize(f)
    canonical_index = Dict{VI,Int}()
    for q_func in q_funcs
        Q = _get_Q_matrix(q_func, canonical_index, n)
        δ_value = _compute_diagonal_shift(Q)
        quad_shift = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], [MOI.ScalarQuadraticTerm(2δ_value[canonical_index[xi]], xi, xi) for xi in keys(canonical_index)], 0.0)
        aff_shift = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-δ_value[canonical_index[xi]], squared_vars[xi]) for xi in keys(canonical_index)], 0.0)
        shifted_q_func = MOIU.canonical(q_func + quad_shift + aff_shift)
        MOIU.normalize_and_add_constraint(model, shifted_q_func, MOI.LessThan(0.0))
    end
    # Now handle objective
    if obj !== nothing
        Q = _get_Q_matrix(obj, canonical_index, n)
        δ_value = _compute_diagonal_shift(Q)
        quad_shift = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], [MOI.ScalarQuadraticTerm(2δ_value[canonical_index[xi]], xi, xi) for xi in keys(canonical_index)], 0.0)
        aff_shift = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-δ_value[canonical_index[xi]], squared_vars[xi]) for xi in keys(canonical_index)], 0.0)
        shifted_q_func = MOIU.canonical(obj + quad_shift + aff_shift)
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), shifted_q_func)
    end

    return MIPforQPBridge(reformulation, squared_vars, binaries)
end

include("box_qp_util.jl")

end # module

max_stars_repo_name: joehuchette/compact-mips-for-qp

max_issues_repo_name: joehuchette/compact-mips-for-qp

max_forks_repo_name: joehuchette/compact-mips-for-qp


text:
abstract type Objective end

# reduce Base.show() to Base.string()
Base.show(io::IO, obj::Objective) = print(io, string(obj))
Base.show(io::IO, ::MIME"application/prs.juno.inline", obj::Objective) = print(io, string(obj))
Base.string(obj::Objective) = @sprintf "ToDo: implement Base.string for %s" typeof(obj)

mutable struct ExpectedScoreObjective <: Objective
    score::Score
    orientation::Symbol
    function ExpectedScoreObjective(score, orientation)
        @assert (orientation .== (:minimise, :maximise)) |> any "orientation must be ':minimise' or ':maximise'"
        new(score, orientation)
    end
end
minimise(score::Score) = ExpectedScoreObjective(score, :minimise)
maximise(score::Score) = ExpectedScoreObjective(score, :maximise)

function Base.string(obj::ExpectedScoreObjective)
    if obj.orientation == :minimise
        return @sprintf "argmin: expected %s" string(obj.score)
    end
    if obj.orientation == :maximise
        return @sprintf "argmax: expected %s" string(obj.score)
    end
end

evaluate(obj::ExpectedScoreObjective, args...) = evaluate(obj.score, args...)
(obj::ExpectedScoreObjective)(args...) = evaluate(obj.score, args...)
integrand_x1(obj::ExpectedScoreObjective, args...; kwargs...) = integrand_x1(obj.score, args...; kwargs...)

update!(obj::ExpectedScoreObjective, x::TI, n::TI) where {TS<:Score,TI<:Integer} = update!(obj.score, x, n)

function update(obj::ExpectedScoreObjective, prior::Prior)
    obj       = deepcopy(obj)
    obj.score = update(obj.score, prior)
    return obj
end

function add!(
        JuMP_model_and_indicator_variables::Tuple,
        objective::ExpectedScoreObjective,
        problem::Problem;
        partial::Tuple{TI,TI} = (0, 0)
    ) where {TI<:Integer}
    m, ind = JuMP_model_and_indicator_variables
    sense  = (objective.orientation == :minimise) ? MOI.OptimizationSense(0) : MOI.OptimizationSense(1)
    # score integrand must include a factor pmf(x1, n1, prior) but in many
    # cases it is more effective to integrate that in the score calculation
    @objective(m, sense,
        sum(
            integrand_x1(
                objective, x1, n1, n2, c2;
                partial_stage_one = partial
            ) *
            ind[(n1, x1, n2, c2)]
            for (n1, x1, n2, c2) in grid(problem)
        )
    )
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ b9f84364-6b75-11eb-183f-9358dcc0e4d0
begin
	using Test
	using LinearAlgebra
end

# ╔═╡ 9c44adc6-6b7c-11eb-20b5-11ad83fee84b
ArrayFloatOrIntDim2 = Array{<:Union{Float64, Int64}, 2}

# ╔═╡ 95a22bc2-6b77-11eb-2cd5-c50658a942b9
md"""
## Gram - Schmidt
"""

# ╔═╡ 3d8a1668-6b70-11eb-293c-ab5ce07a3ce5
begin
	function GramSchmidt(A::ArrayFloatOrIntDim2)
		m, n = size(A)
		Q = zeros(m, n)
		R = zeros(n, n)
		for j = 1:n
			v = A[:,j]
			for i = 1:j-1
				R[i, j] = Q[:, i]' * v
				v = v - R[i, j] * Q[:, i]
			end
			R[j, j] = norm(v)
			Q[:,j] = v / R[j,j]
		end
		return (Q, R)
	end
end

# ╔═╡ 7356cb64-6b71-11eb-01f7-91f369a9fb74
begin
	A = [
		1 0
		1 1
		1 3
		1 4
	]
	Q, R = GramSchmidt(A)

	@test A == round.(Int, Q * R)
end

# ╔═╡ 78c7863c-6b7c-11eb-0be7-cbfb3c368065
md"""
## Householder (probably wrong)
"""

# ╔═╡ 83f357d4-6b7c-11eb-11b6-5ddf77041fe5
begin
	function householder(A::ArrayFloatOrIntDim2)
		A = convert(Array{Float64}, A)
		m, n = size(A)
		U = zeros(m, n)
		for k = 1:n
			w = A[k:m, k]
			wₙ = norm(w)
			w[1] = (w[1] - wₙ)
			u = w / wₙ
			U[k:m, k] = u
			A[k:m, k:n] = A[k:m, k:n] - 2 * u * (u' * A[k:m, k:n])
		end
		R = triu(A[:, 1:n])
		return U, R
	end
end

# ╔═╡ 43133fa8-6b7d-11eb-331e-73d6eabd7d83
begin
	U, Rₕ = householder(A)
end

# ╔═╡ Cell order:
# ╠═b9f84364-6b75-11eb-183f-9358dcc0e4d0
# ╠═9c44adc6-6b7c-11eb-20b5-11ad83fee84b
# ╟─95a22bc2-6b77-11eb-2cd5-c50658a942b9
# ╠═3d8a1668-6b70-11eb-293c-ab5ce07a3ce5
# ╠═7356cb64-6b71-11eb-01f7-91f369a9fb74
# ╟─78c7863c-6b7c-11eb-0be7-cbfb3c368065
# ╠═83f357d4-6b7c-11eb-11b6-5ddf77041fe5
# ╠═43133fa8-6b7d-11eb-331e-73d6eabd7d83


text:
using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
import ModelingToolkit: Interval, infimum, supremum

@parameters x y t
@variables u(..)
Dx=Differential(x)^1
Dy=Differential(y)^1
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# 2D PDE
eq  = Dxx(u(x,y)) + Dyy(u(x,y))+ u(x,y) ~ 0

# Boundary conditions

bcs = [Dxx(u(-1,y)) ~ -u(-1,y), Dxx(u(1,y)) ~ -u(1,y),
        Dyy(u(x,-1)) ~ -u(x,-1)+2,Dy(u(range(0, step=0.05, stop=1),0)) ~ 0,Dyy(u(x,1))~ -u(x,1)]
    #
# Space and time domains
domains = [x ∈ Interval(-1.0,1.0),
           y ∈ Interval(-1.0,1.0)]

# Neural network
dim = 2 # number of dimensions
chain = FastChain(FastDense(dim,8,Flux.σ),FastDense(8,8,Flux.σ),FastDense(8,1))
# Initial parameters of Neural network
initθ = Float64.(DiffEqFlux.initial_params(chain))

# Discretization
dx = 0.05
discretization = PhysicsInformedNN(chain,GridTraining(dx),init_params =initθ)

@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
prob = discretize(pde_system,discretization)

#Optimizer
opt = Optim.BFGS()

#Callback function
cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob, opt, cb = cb, maxiters=1500)
phi = discretization.phi

using Plots

xs,ys = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains]
analytic_sol_func(x,y) = y^2*sin(pi*x)

u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u = abs.(u_predict .- u_real)

p1 = plot(xs, ys, u_real, linetype=:surface,title = "analytic");
p2 = plot(xs, ys, u_predict, linetype=:surface,title = "predict");
p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
plot(p1,p2,p3)

max_stars_repo_name: mohsenhos/helmholtzpde

max_issues_repo_name: mohsenhos/helmholtzpde

max_forks_repo_name: mohsenhos/helmholtzpde


text:
module Report

import Base.write
import Base.close

type Markdown
  filename::String
  writeorappend::ASCIIString
  figurefolder::String
  iostream::IO
end
Markdown(filename) = Markdown(filename, "w", "", open(filename,"w"))
Markdown(filename, writeorappend) = Markdown(filename, writeorappend,"", open(filename,writeorappend))
Markdown(filename, writeorappend, figurefolder) = Markdown(filename, writeorappend,figurefolder, open(filename,writeorappend))
close(doc::Markdown) = close(doc.iostream)

type Header
  level::Integer
  text::String
end

type Blockquote
  text::String
end

type Author
  name::String
  email::String
  address::String
  affiliation::String
  website::String
end

type Figure
  url::String
  caption::String
end

type Code
  language::String
  code::String
  identifier::String
  numberLines::Bool
  startFrom::Integer
end
Code(language, code) = Code(language, code, "", false, 1)

type Table
  nrows::Integer
  ncolumns::Integer
  header::Array{String}
  data::Array{Any}
  caption::String
end

type Link
  text::String
  url::String
end

type Paragraph
  text::String
end

include("reportwriter.jl")

end


text:
module ClusteringDifferences

using DifferencesBase
using Distances

import Clustering: FuzzyCMeansResult, KmeansResult, KmedoidsResult

export
    # clustering.jl
    AbstractClustering,
    HierarchicalClustering,
    PartitionalClustering,
    assignments,
    constraints,
    features,
    parameters,
    weights,
    # difference.jl
    AbstractClusteringDifference,
    PartitionalClusteringDifference,
    backwarddiff,
    forwarddiff,
    # kmeans.jl
    kmeans
    # pckmeans.jl
    #pckmeans

include("clustering.jl")
include("difference.jl")
include("kmeans.jl")
#include("pckmeans.jl")

end # module


text:
module HelioseismicKernels

using BipolarSphericalHarmonics
using BipolarSphericalHarmonics: kronindex
using DSP: hilbert
using DelimitedFiles
using Distributed
using DualNumbers
using EllipsisNotation
using FastGaussQuadrature
using FFTW
using FileIO
using FITSIO
using HybridArrays
using JLD2
using LegendrePolynomials
using LinearAlgebra
BLAS.set_num_threads(1)
using LoopVectorization
using MPI
using NamedArrays
using NumericalIntegration
using OffsetArrays
using OffsetArrays: no_offset_view
using ParallelUtilities
using ParallelUtilities: ProductSplit, extremaelement
using Polynomials
using Printf
using ProgressMeter
using Rotations
using SparseArrays
using SphericalHarmonicModes
using SphericalHarmonicArrays
using SphericalHarmonicArrays: shmodes
using SphericalHarmonics
using StaticArrays
using StructArrays
using SuiteSparse
using TimerOutputs
using UnPack
using DistributedArrays
using VectorSphericalHarmonics
using WignerD
using WignerSymbols

export TravelTimes, Amplitudes, SoundSpeed, Flow
export Point2D, Point3D
export los_radial, los_earth

firstshmodes(x) = first(shmodes(x))

# remove wrappers to get the underlying array
reinterpret_as_float(arr::OffsetArray) = reinterpret_as_float(no_offset_view(arr))
# convert complex arrays to real
reinterpret_as_float(arr::Array{<:Complex}) = reinterpret(real(eltype(arr)), arr)
# real arrays may be returned as is
reinterpret_as_float(arr::Array{<:Real}) = arr

reinterpret_as_complex(arr::OffsetArray) = reinterpret_as_complex(no_offset_view(arr))
reinterpret_as_complex(arr::Array{<:Complex}) = arr
reinterpret_as_complex(arr::Array{<:Real}) = reinterpret(complex(eltype(arr)), arr)

function save_to_fits_and_return(filename, arr::Array; kwargs...)
	filepath = joinpath(SCRATCH_kerneldir[], filename)
	FITS(filepath, "w") do f
		outarr = reinterpret_as_float(arr)
		FITSIO.write(f, outarr; kwargs...)
	end
	arr
end

function save_to_fits_and_return(filename, val::Number; kwargs...)
	save_to_fits_and_return(filename, [val]; kwargs...)
end

function Powspec(ω)
	σ = 2π*0.4e-3
	ω0 = 2π*3e-3
	exp(-(ω-ω0)^2/(2σ^2))
end

#################################################################
# Cross covariances and changes in cross covariances
#################################################################
abstract type SeismicMeasurement end
struct TravelTimes <: SeismicMeasurement end
struct Amplitudes <: SeismicMeasurement end

abstract type PerturbationParameter end
struct SoundSpeed <: PerturbationParameter end
struct Flow <: PerturbationParameter end

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

const SCRATCH = Ref("")
const SCRATCH_kerneldir = Ref("")

function __init__()
    SCRATCH[] = get(ENV, "SCRATCH",
                isdir(joinpath("/scratch", ENV["USER"])) ?
                joinpath("/scratch", ENV["USER"]) : pwd())

    SCRATCH[] = expanduser(SCRATCH[])
    SCRATCH_kerneldir[] = joinpath(SCRATCH[], "kernels")

    if !isdir(SCRATCH_kerneldir[])
        @info "Creating $(SCRATCH_kerneldir[])"
        mkdir(SCRATCH_kerneldir[])
    end
end

function rank_size(comm)
	rank = MPI.Comm_rank(comm)
	np = MPI.Comm_size(comm)
	rank, np
end
function productsplit(iters, comm)
	rank, np = rank_size(comm)
	ProductSplit(iters, np, rank + 1)
end

include("$(@__DIR__)/pointsonasphere.jl")
include("$(@__DIR__)/finite_difference.jl")
include("$(@__DIR__)/continuous_FFT.jl")
include("$(@__DIR__)/directions.jl")
include("$(@__DIR__)/timer_utils.jl")
include("$(@__DIR__)/greenfn.jl")
include("$(@__DIR__)/crosscov.jl")
include("$(@__DIR__)/kernel.jl")
include("$(@__DIR__)/kernel3D.jl")
include("$(@__DIR__)/measurementsvalidation.jl")

end # module


text:
# Written in Julia by André Ferrari
# January 2015
#
# Translation of matlab files bvn.m and bvnu.m authored by:
#
# Alan Genz
# Department of Mathematics
# Washington State University
# Pullman, WA 99164-3113
# Email : alangenz@wsu.edu
#
# Original source available from
# http://www.math.wsu.edu/faculty/genz/software/matlab/bvn.m
#
# Copyright (C) 2013, Alan Genz,  All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided the following conditions are met:
#   1. Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#   2. Redistributions in binary form must reproduce the above copyright
#      notice, this list of conditions and the following disclaimer in
#      the documentation and/or other materials provided with the
#      distribution.
#   3. The contributor name(s) may not be used to endorse or promote
#      products derived from this software without specific prior
#      written permission.
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#

function bvn(xl::Number,xu::Number,yl::Number,yu::Number,r::Number)
    p = bvnu(xl,yl,r) - bvnu(xu,yl,r) - bvnu(xl,yu,r) + bvnu(xu,yu,r)
    p = max(0.0, min(p,1.0))
end

phid(z::Number) = erfc(-z/sqrt(2.0))/2.0

function bvnu(dh::Number,dk::Number,r::Number)
    if dh ==  Inf || dk ==  Inf
        p = 0.0
    elseif dh == -Inf
        if dk == -Inf
            p = 1.0
        else
            p = phid(-dk)
        end
    elseif dk == -Inf
        p = phid(-dh)
    elseif r == 0.0
        p = phid(-dh)*phid(-dk)
    else
        tp = 2.0*pi
        h = dh
        k = dk
        hk = h*k
        bvn = 0.0
        if abs(r) < 0.3
            w = [0.1713244923791705, 0.3607615730481384, 0.4679139345726904]
            x = [0.9324695142031522, 0.6612093864662647, 0.2386191860831970]
        elseif abs(r) < 0.75
            w = [.04717533638651177, 0.1069393259953183, 0.1600783285433464]
            w = [w, 0.2031674267230659, 0.2334925365383547, 0.2491470458134029]
            x = [0.9815606342467191, 0.9041172563704750, 0.7699026741943050]
            x = [x, 0.5873179542866171, 0.3678314989981802, 0.1252334085114692]
        else
            w = [.01761400713915212, .04060142980038694, .06267204833410906]
            w = [w, .08327674157670475, 0.1019301198172404, 0.1181945319615184]
            w = [w, 0.1316886384491766, 0.1420961093183821, 0.1491729864726037]
            w = [w, 0.1527533871307259]
            x = [0.9931285991850949, 0.9639719272779138, 0.9122344282513259]
            x = [x, 0.8391169718222188, 0.7463319064601508, 0.6360536807265150]
            x = [x, 0.5108670019508271, 0.3737060887154196, 0.2277858511416451]
            x = [x, 0.07652652113349733]
        end
        w = [w, w]
        x = [1.0-x, 1.0+x]

        if abs(r) < 0.925
            hs = ( h*h + k*k )/2
            asr = asin(r)/2.0
            sn = sin(asr*x)
            bvn = dot(exp((sn*hk-hs)./(1-sn.^2)),w)
            bvn = bvn*asr/tp + phid(-h)*phid(-k)
        else
            if r < 0.0
                k = -k
                hk = -hk
            end
            if abs(r) < 1.0
                as = 1.0-r^2
                a = sqrt(as)
                bs = (h-k)^2
                asr = -( bs/as + hk )/2.0
                c = (4.0-hk)/8.0
                d = (12.0-hk)/80.0
                if asr > -100.0
                    bvn = a*exp(asr)*(1.0-c*(bs-as)*(1-d*bs)/3.0+c*d*as^2)
                end
                if hk  > -100.0
                    b = sqrt(bs)
                    sp = sqrt(tp)*phid(-b/a)
                    bvn = bvn - exp(-hk/2.0)*sp*b*( 1.0 - c*bs*(1-d*bs)/3.0 )
                end
                a = a/2.0
                xs = (a*x).^2
                asr = -( bs./xs + hk )/2.0
                ix = find( asr .> -100.0 )
                xs = xs[ix]
                sp = ( 1.0 + c*xs.*(1.0+5.0*d*xs))
                rs = sqrt(1.0-xs)
                ep = exp( -(hk/2.0)*xs./(1+rs).^2 )./rs
                bvn = (a*dot( exp(asr[ix]).*(sp-ep),w[ix]) - bvn)/tp
            end

            if r > 0.0
                bvn =  bvn + phid(-max(h,k))
            elseif h >= k
                bvn = -bvn
            else
                if h < 0.0
                    L = phid(k)-phid(h)
                else
                    L = phid(-h)-phid(-k)
                end
                bvn =  L - bvn
            end
        end
        p = max(0.,min(1.,bvn))
    end

    return p

end


text:
@json struct Tag
    links::NamedTuple
    name::String
    target::Commit
end
@json_struct Tag

endpoint(::BitbucketAPI, ::typeof(get_tags), workspace::String, repo::String) = Endpoint(:GET, "/repositories/$workspace/$repo/refs/tags")
@not_implemented(::BitbucketAPI, ::typeof(get_tags), ::UUID)
into(::BitbucketAPI, ::typeof(get_tags)) = Page{get_tags, Tag}


text:
function helper_fn(z, attempt)
	single_data_name, test_fold = z
	proc_id = myid();

	helper_start_time = time()

	transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/synthetic_data_assortments/", single_data_name, "_transaction_counts.csv") ) )
	assortments = convert(Array{Int64}, readcsv( string("../tcm_data/synthetic_data_assortments/", single_data_name, "_assortments.csv") ) )

	grand_M = Int64(50)
	folds = Int64[ 1*ones(grand_M/5); 2*ones(grand_M/5); 3*ones(grand_M/5); 4*ones(grand_M/5); 5*ones(grand_M/5) ]

	train_assortments = assortments[ find(folds .!= test_fold), :]
	test_assortments = assortments[ find(folds .== test_fold), :]

	train_transaction_counts = transaction_counts[ find(folds .!= test_fold), :]
	test_transaction_counts = transaction_counts[ find(folds .== test_fold), :]

	num_train_transactions = sum(train_transaction_counts)
	num_test_transactions = sum(test_transaction_counts)

	num_train_assortments = size(train_assortments, 1)
	num_test_assortments = size(test_assortments, 1)

	N = 10

	time_limit = 1200

	perTol = 1e-6

	## RANKING-BASED MODEL ##

	ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking(N, train_assortments, train_transaction_counts, time_limit, perTol)
	ranking_K = length(ranking_lambda)

	test_ranking_predict, test_ranking_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, test_assortments)
	ranking_KL, ranking_noinf_count, ranking_KL_avg = tcm_evaluateKL(N, test_ranking_predict, test_transaction_counts)

	avg_CS_size, weighted_avg_CS_size, max_CS_size = tcm_rankingStats(N, ranking_lambda, orderings)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"ranking", ",",
							ranking_loglik, ",",
							ranking_K, ",",
							perTol, ",",
							ranking_KL, ",",
							ranking_noinf_count, ",",
							ranking_KL_avg, ",",
							ranking_elapsed_time, ",",
							avg_CS_size, ",",
							weighted_avg_CS_size, ",",
							max_CS_size)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(ranking_outcsvhandle, row_string);
		flush(ranking_outcsvhandle);
	end

	## LC-MNL MODEL ##

	K_list = [5]
	numInitializations = 5

	for K in K_list
		LCMNL_u_val, LCMNL_class_probs, LCMNL_loglik, LCMNL_grand_elapsed_time = tcm_estimateMMNL(N, K, train_assortments, train_transaction_counts, numInitializations, time_limit, perTol)

		test_LCMNL_predict, test_LCMNL_predict_by_class = tcm_predictMMNL(N, K, LCMNL_class_probs, LCMNL_u_val, test_assortments)
		LCMNL_KL, LCMNL_noinf_count, LCMNL_KL_avg = tcm_evaluateKL(N, test_LCMNL_predict, test_transaction_counts)

		row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"LCMNL", ",",
							LCMNL_loglik, ",",
							K, ",",
							perTol, ",",
							numInitializations, ",",
							LCMNL_KL, ",",
							LCMNL_noinf_count, ",",
							LCMNL_KL_avg, ",",
							LCMNL_grand_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(MNL_outcsvhandle, row_string);
			flush(MNL_outcsvhandle);
		end
	end

	## MNL MODEL ##

	MNL_u_val, MNL_loglik, MNL_elapsed_time = tcm_estimateMNL(N, train_assortments, train_transaction_counts)
	test_MNL_predict = tcm_predictMNL(N, MNL_u_val, test_assortments)
	MNL_KL, MNL_noinf_count, MNL_KL_avg = tcm_evaluateKL(N, test_MNL_predict, test_transaction_counts)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"MNL", ",",
							MNL_loglik, ",",
							1, ",",
							0, ",",
							1, ",",
							MNL_KL, ",",
							MNL_noinf_count, ",",
							MNL_KL_avg, ",",
							MNL_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(MNL_outcsvhandle, row_string);
		flush(MNL_outcsvhandle);
	end

	## HALO-MNL MODEL ##

	HALOMNL_u_val, HALOMNL_alpha_val, HALOMNL_loglik, HALOMNL_elapsed_time = tcm_estimateHALOMNL(N, train_assortments, train_transaction_counts)

	test_HALOMNL_predict = tcm_predictHALOMNL(N, HALOMNL_u_val, HALOMNL_alpha_val, test_assortments)

	HALOMNL_KL, HALOMNL_noinf_count, HALOMNL_KL_avg = tcm_evaluateKL(N, test_HALOMNL_predict, test_transaction_counts)

	row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"HALOMNL", ",",
							HALOMNL_loglik, ",",
							1, ",",
							0, ",",
							1, ",",
							HALOMNL_KL, ",",
							HALOMNL_noinf_count, ",",
							HALOMNL_KL_avg, ",",
							HALOMNL_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(MNL_outcsvhandle, row_string);
		flush(MNL_outcsvhandle);
	end

	## FOREST, COLD-START ##

	depth_list = [4]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_createBasicForest(N)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_loglik, forest_elapsed_time = tcm_estimateForest(N, train_assortments, train_transaction_counts, time_limit, perTol, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		test_forest_predict, test_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_KL, forest_noinf_count, forest_KL_avg = tcm_evaluateKL(N, test_forest_predict, test_transaction_counts)

		simple_avg_leaves, weighted_avg_leaves, max_leaves = tcm_leafStats(lambda, forest_isLeaf)
		simple_avg_depth, weighted_avg_depth, max_depth = tcm_depthStats(lambda, forest_left, forest_right, forest_isLeaf)

		row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cold", ",",
							depth_limit, ",",
							forest_loglik, ",",
							K, ",",
							perTol, ",",
							forest_KL, ",",
							forest_noinf_count, ",",
							forest_KL_avg, ",",
							forest_elapsed_time, ",",
							simple_avg_leaves, ",",
							weighted_avg_leaves, ",",
							max_leaves, ",",
							simple_avg_depth, ",",
							weighted_avg_depth, ",",
							max_depth)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end

		@show forest_KL
	end

	## FOREST, WARM-START ##

	depth_list = [4]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_convertRankingToForest(N, orderings)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_loglik, forest_elapsed_time = tcm_estimateForest(N, train_assortments, train_transaction_counts, time_limit, perTol, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		test_forest_predict, test_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_KL, forest_noinf_count, forest_KL_avg = tcm_evaluateKL(N, test_forest_predict, test_transaction_counts)

		simple_avg_leaves, weighted_avg_leaves, max_leaves = tcm_leafStats(lambda, forest_isLeaf)
		simple_avg_depth, weighted_avg_depth, max_depth = tcm_depthStats(lambda, forest_left, forest_right, forest_isLeaf)

		row_string = string(single_data_name,  ",",
							test_fold, ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-ranking-warm", ",",
							depth_limit, ",",
							forest_loglik, ",",
							K, ",",
							perTol, ",",
							forest_KL, ",",
							forest_noinf_count, ",",
							forest_KL_avg, ",",
							forest_elapsed_time, ",",
							simple_avg_leaves, ",",
							weighted_avg_leaves, ",",
							max_leaves, ",",
							simple_avg_depth, ",",
							weighted_avg_depth, ",",
							max_depth)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	helper_end_time = time()

	return helper_end_time - helper_start_time;

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
"""
Types used to build the structure for the 1d1d scrap bulk modeling
"""

"""
    M1d1dMeshSettings

Settings for the mesh of the 1d1d modelling. Mesh is generated in the simulation step, from bulk, scrap and shaft properties itself.

# Arguments
- `no_cells_bulk::Integer`: number of cells over bulk length
- `no_cells_scrap_piece::Union{Nothing, Integer}`: number of cells over scrap piece
"""
struct M1d1dMeshSettings
    no_cells_bulk::Integer
    no_cells_scrap_piece::Union{Nothing, Integer}
end

#= #############################################################################
Simulation Setting Types
############################################################################# =#

"""
    SimSettingsSubModels

Settings for activation and special options for the submodels

"""
struct SimSettingsSubModels{T <: AbstractFloat}
    scrap_discretization_model_enabled::Bool    # defaults to false
    min_no_of_scrap_cells::Integer              # defaults to 3
    cross_radiation_model_enabled::Bool         # defaults to false
    cross_radiation_constant_c::T               # defaults to ?
    furnace_side_rad_model::Union{Nothing,HotSideRadiationModel}

    SimSettingsSubModels(;scrap_discretization_model_enabled = false,
            min_no_of_scrap_cells = 5,
            cross_radiation_model_enabled = false,
            cross_radiation_constant_c = 1.0,
            furnace_side_rad_model = nothing) =
            new{typeof(cross_radiation_constant_c)}(scrap_discretization_model_enabled,
                    min_no_of_scrap_cells, cross_radiation_model_enabled,
                    cross_radiation_constant_c, furnace_side_rad_model)
end

"""
SimSettings:
Struct representing the simulation settings for a phase in the
simulation of the 1d1d model.

All uinits musst be given in SI units, temperatures in K

dt_save:
interval time for results beeing written to Vector{ResultsTimeStep**ScrapDisc}
"""
struct SimSettings{T <: AbstractFloat}
    duration::T
    submodels::SimSettingsSubModels{T}
    dt_save::T # defaults  to duration?

    SimSettings(
                duration::T2,
                submodels::SimSettingsSubModels{T2}
                ) where {T2 <: AbstractFloat}  = new{T2}(duration, submodels, duration)

    SimSettings(
        duration::T2, submodels::SimSettingsSubModels{T2}, dt_save::T2
        ) where {T2 <: AbstractFloat} = new{T2}(duration, submodels, dt_save)
end

#= #############################################################################
Scrap Bulk Types
############################################################################# =#

# """
# ScrapMixture:
# Desribes the mixture of different ScrapPiece(s). Here mass_frac[i] is the mass
# fraction of a ScrapPiece (pieces[i]).
# (sum(mass_frac[:] musst be equal to 1)
#TODO: Make this work in the future
# """
# mutable struct ScrapMixture
#     pieces::Vector{ScrapPiece}
#     mass_frac::Vector{Float64}
# end

"""
CharScrapPiece:
A mutable struct CharScrapPiece (charcterisitic piece of scrap) is asumed
to be a flat shaped object which has a charcterisitic
thickness that can represent the heating throughout it's body.
Also it has a surf_area which ist mostly relevant for convection and
radiation transport.
"""
struct CharScrapPiece{T <: AbstractFloat}
    char_thickness::T
    char_cross_section_area::T
    thermal_props::SolidThermalMaterial
    surf_area::T
    vol::T

    # initialize from dimensions and mat
    CharScrapPiece( T::T2,
                    W::T2,
                    L::T2,
                    mat::SolidThermalMaterial) where {T2<:AbstractFloat}=
                new{T2}(
                        char_thickness_quader([T, W, L]),
                        char_cross_section_area_quader([T, W, L]),
                        mat,
                        surface_quader(T,W,L),
                        volume_quader(T,W,L)
                    )
end

"""
Bulk:
Struct representing the properties and geometry of the scrap bulk

scrap_pieces_contact_factor: factor which defines the average amount of scrap
piece surface, which has direct contact with another scrap piece.
All units musst be given in SI units
"""
struct Bulk{T <: AbstractFloat}
    mass_flow::T
    T_charged::T
    length::T
    cross_sect_surf_area::T
    density::T
    porosity::T
    scrap_pieces_contact_factor::T
   # mat::ScrapMixture # Future see ScrapMixture
    scrap_props::CharScrapPiece
end

#= #############################################################################
Off-gas Types
############################################################################# =#

"""
OffGas:
Struct representing the off-gas

All uinits musst be given in SI units, temperatures in K
"""
struct OffGas{T <: AbstractFloat}
    mass_flow::T
    T_in::T
    thermal_props::IncompressibleFluidThermalMaterial
end

#= #############################################################################
Final Model Type
############################################################################# =#

"""
M1d1d:
Model parent type
"""
struct M1d1d
    sim_phase::Vector{SimSettings}
    heat_transfer::Vector{HeatTransfer}
    bulk::Vector{Bulk}
    offgas::Vector{OffGas}
    mesh_settings::Vector{M1d1dMeshSettings}
end


text:
using Mimi
using Distributions
using Test

include("test-model/test-model.jl")
using .TestModel

m = create_model()
mm1 = create_marginal_model(create_model())
mm2 = create_marginal_model(create_model())

simdef = @defsim begin
    share = Uniform(0, 1)
    save(emissions.E_Global)
end

# Test running one MarginalModel

function post_trial1(si, trialnum, ntimesteps, tup)
    @test length(si.models) == 1
    @test si.models[1] isa MarginalModel
end

si = run(simdef, mm1, 2, post_trial_func = post_trial1, results_in_memory = true)
@test all(iszero, si.results[1][:emissions, :E_Global][!, :E_Global])   # Test that the marginal emission saved from the MarginalModel are zeros (because there's no difference between mm.base and mm.marginal)

# Test running a vector of MarginalModels

function post_trial2(si, trialnum, ntimesteps, tup)
    @test length(si.models) == 2
    @test si.models[1] isa MarginalModel
    @test si.models[2] isa MarginalModel
end

si = run(simdef, [mm1, mm2], 2, post_trial_func = post_trial2, results_in_memory = true)
@test all(iszero, si.results[1][:emissions, :E_Global][!, :E_Global])  # Test that the regular model has non-zero emissions
@test all(iszero, si.results[2][:emissions, :E_Global][!, :E_Global])   # Test that the marginal emission saved from the MarginalModel are zeros (because there's no difference between mm.base and mm.marginal)

# Test running a vector of a Model and a MarginalModel

function post_trial3(si, trialnum, ntimesteps, tup)
    @test length(si.models) == 2
    @test si.models[1] isa Model
    @test si.models[2] isa MarginalModel
end

si = run(simdef, [m, mm1], 2, post_trial_func = post_trial3, results_in_memory = true)
@test all(!iszero, si.results[1][:emissions, :E_Global][!, :E_Global])  # Test that the regular model has non-zero emissions
@test all(iszero, si.results[2][:emissions, :E_Global][!, :E_Global])   # Test that the marginal emission saved from the MarginalModel are zeros (because there's no difference between mm.base and mm.marginal)


text:
using LightGraphs

import MathProgBase

mutable struct NodeData{S}
    nlds::NLDS{S}
    npath::Dict{Int, Int}

    # Feasibility cuts
    fcuts::CutStore{S}
    # Optimality cuts
    ocuts::CutStore{S}

    function NodeData{S}(nlds::NLDS{S}, nvars_a) where S
        new{S}(nlds, Dict{Int, Int}(), CutStore{S}(nvars_a), CutStore{S}(nvars_a))
    end
end

NodeData(nlds::NLDS{S}, parent) where {S} = NodeData{S}(nlds, parent)

function Base.show(io::IO, data::NodeData)
    println(io, "Node of $(data.nlds.nx) variables")
end

# Mutable for setprobability!
mutable struct Transition{S} <: SOI.AbstractTransition
    source::Int
    target::Int
    σ::Int # FIXME NLDS is the only one who needs to map out_transitions to 1:n when it uses AveragedCut, it should have an internal dictionary
    proba::S
    childT::Union{Nothing, AbstractMatrix{S}}
    function Transition(source::Int, target::Int, σ::Int, proba::S, childT) where S
        new{S}(source, target, σ, proba, childT)
    end
end
SOI.get(::SOI.AbstractStochasticProgram, ::SOI.Source, tr::Transition) = tr.source
SOI.get(::SOI.AbstractStochasticProgram, ::SOI.Target, tr::Transition) = tr.target
SOI.get(::SOI.AbstractStochasticProgram, ::SOI.Probability, tr::Transition) = tr.proba

Base.:(==)(t1::Transition, t2::Transition) = t1.source == t2.source && t1.target == t2.target && t1.σ == t2.σ
Base.hash(t::Transition, u::UInt) = hash(t.source, hash(t.target, hash(t.σ, u)))

"""
    StochasticProgram{S, TT}

StochasticProgram of coefficient type `S` and transition type `TT`.
"""
mutable struct StochasticProgram{S, TT} <: SOI.AbstractStochasticProgram
    out_transitions::Vector{Vector{TT}} # out_transitions[i] : outgoing transitions for node i
    data::Vector{NodeData{S}}           # data[i] : data for node i
    num_stages::Int
    function StochasticProgram{S, TT}(num_stages) where {S, TT}
        new{S, TT}(Vector{TT}[], NodeData{S}[], num_stages)
    end
end
StochasticProgram{S}(num_stages) where S = StochasticProgram{S, Transition{S}}(num_stages)

nodedata(sp::StochasticProgram, node::Int) = sp.data[node]

SOI.get(sp::StochasticProgram, ::SOI.NumberOfStages) = sp.num_stages

SOI.get(sp::StochasticProgram, ::SOI.NodeObjectiveValueBound, node) = getobjectivebound(nodedata(sp, node).nlds)
SOI.set!(sp::StochasticProgram, ::SOI.TransitionObjectiveValueBound, tr::Transition, θlb) = setθbound!(nodedata(sp, SOI.get(sp, SOI.Source(), tr)).nlds, edgeid(sp, tr), θlb)
SOI.get(sp::StochasticProgram, ::SOI.Dimension, node) = nodedata(sp, node).nlds.nx

SOI.get(sp::StochasticProgram, ::SOI.OutTransitions, node::Int) = sp.out_transitions[node]
SOI.get(::StochasticProgram{S, TT}, ::SOI.TransitionType) where {S, TT} = TT

SOI.get(::StochasticProgram, ::SOI.MasterNode) = 1

# If the graph is not a tree, this will loop if I don't use a num_stages limit
function SOI.get(sp::StochasticProgram, nop::SOI.NumberOfPathsFrom, node)
    len = nop.length
    @assert len >= 0
    if iszero(len) || iszero(outdegree(sp, node))
        1
    else
        npath = nodedata(sp, node).npath
        if !(len in keys(npath))
            npath[len] = sum(map(tr -> SOI.get(sp, SOI.NumberOfPathsFrom(len-1), SOI.get(sp, SOI.Target(), tr)), SOI.get(sp, SOI.OutTransitions(), node)))
        end
        npath[len]
    end
end

function SOI.add_scenario_node!(sp::StochasticProgram{S}, data::NodeData) where S
    push!(sp.out_transitions, Transition{S}[])
    push!(sp.data, data)
    @assert length(sp.out_transitions) == length(sp.data)
    length(sp.data)
end

function SOI.add_scenario_transition!(sp::StochasticProgram, parent, child, proba, childT=nothing)
    tr = Transition(parent, child, outdegree(sp, parent)+1, proba, childT)
    push!(sp.out_transitions[parent], tr)
    data = nodedata(sp, parent)
    empty!(data.npath)
    childdata = nodedata(sp, child)
    add_scenario_transition!(data.nlds, childdata.fcuts, childdata.ocuts, proba, childT)
    @assert length(data.nlds.childFC) == length(data.nlds.proba) == outdegree(sp, parent)
    tr
end

function SOI.set!(sp::StochasticProgram, ::SOI.Probability, tr, proba)
    tr.proba = proba
    data = nodedata(sp, SOI.get(sp, SOI.Source(), tr))
    setprobability!(data.nlds, edgeid(sp, tr), proba)
end

edgeid(sp::StochasticProgram, tr::Transition) = tr.σ

function SOI.get(sp::StochasticProgram, ::SOI.Solution, node)
    getsolution(nodedata(sp, node).nlds)
end

function SOI.set!(sp::StochasticProgram, ::SOI.SourceSolution, tr, sol::Solution)
    data = nodedata(sp, SOI.get(sp, SOI.Source(), tr))
    if tr.childT !== nothing
        T = tr.childT
        x = T * sol.x
        xuray = sol.xuray
        if xuray !== nothing
            xuray = T * sol.xuray
        end
    else
        x = sol.x
        xuray = sol.xuray
    end
    setparentx(nodedata(sp, SOI.get(sp, SOI.Target(), tr)).nlds, x, xuray, sol.objvalxuray)
end

function SOI.getbellmanvalue(sp::StochasticProgram, tr::SOI.AbstractTransition, sol::Solution)
    @assert length(sol.θ) == outdegree(sp, SOI.get(sp, SOI.Source(), tr))
    SOI.getbellmanvalue(sol, edgeid(sp, tr))
end

function SOI.getbellmanvalue(sp::StochasticProgram, node, sol::Solution)
    @assert length(sol.θ) == 1
    SOI.getbellmanvalue(sol, 1)
end

SOI.get(sp::StochasticProgram, ::SOI.NeedAllSolutions, node) = needallsolutions(SOI.get(sp, CutGenerator(), node))

"""
    CutGenerator <: AbstractNodeAttribute

The cut generator of the node.
"""
struct CutGenerator <: SOI.AbstractNodeAttribute end

SOI.get(sp::StochasticProgram, ::CutGenerator, node) = nodedata(sp, node).nlds.cutgen
function SOI.set!(sp::StochasticProgram, ::CutGenerator, node, cutgen::AbstractOptimalityCutGenerator)
    nodedata(sp, node).nlds.cutgen = cutgen
end


text:
function plot( obj::tikzObj, x::Vector{Float64}, y::Vector{Float64};
               lw = 0.0, smooth = false, color = "", linestyle="", style=tikzStyle() )

    @assert length(x) == length(y)

    n = length(x)

    # print start
    @printf( obj.fID, "\\draw[" )
    if lw > 0.0
        @printf( obj.fID, "line width=%3.2fpt, ", lw )
    end
    if length(color) > 0
        @printf( obj.fID, "%s, ", color )
    end
    if linestyle == "--" || linestyle == "dashed"
        @printf( obj.fID, "dashed, " )
    elseif length(linestyle) > 0
        @printf( obj.fID, "%s, ", linestyle )
    end
    if style.set
        @printf( obj.fID, "%s, ", style.name )
    end

    @printf( obj.fID, "] plot[" )
    if smooth || obj.smooth
        @printf( obj.fID, "smooth, " )
    end
    @printf( obj.fID, "] coordinates { %%\n" )

    # print coordinates
    for jj in 1:n-1
        @printf( obj.fID, "  (%6.4f,%6.4f) %%\n", x[jj], y[jj] )
    end

    # print last point and close
    @printf( obj.fID, "  (%6.4f,%6.4f) };\n\n", x[n], y[n] )

end

function plot( obj::tikzObj, objAx::axisObj, x::Vector{Float64}, y::Vector{Float64};
               lw = 0.0, smooth = false, color = "", linestyle="", style=tikzStyle() )

    @assert length(x) == length(y)

    @printf( obj.fID, "  \\addplot [ " )

    # add options
    if lw > 0.0
        @printf( obj.fID, "line width=%3.2fpt, ", lw )
    end
    if length(color) > 0
        @printf( obj.fID, "%s, ", color )
    end
    if linestyle == "--" || linestyle == "dashed"
        @printf( obj.fID, "dashed, " )
    elseif length(linestyle) > 0
        @printf( obj.fID, "%s, ", linestyle )
    end
    if style.set
        @printf( obj.fID, "%s, ", style.name )
    end
    # end options
    @printf( obj.fID, "  ]\n " )
    @printf( obj.fID, "    table[row sep=crcr]{%%\n " )

    # coordinates
    for jj in 1:length(x)
        @printf( obj.fID, "    %6.4f %6.4f\\\\ \n ", x[jj], y[jj] )
    end

    # close plot
    @printf( obj.fID, "  };\n\n" )

end

max_stars_repo_name: mopg/julia-tikz

max_issues_repo_name: mopg/julia-tikz

max_forks_repo_name: mopg/julia-tikz


text:
# using Pkg
# Pkg.activate(".")
using Weave

weave("plots_ising.jmd", doctype = "md2html", out_path="build/")

max_stars_repo_name: edwinb-ai/IsingRBM

max_issues_repo_name: edwinb-ai/IsingRBM

max_forks_repo_name: edwinb-ai/IsingRBM


text:
# necessary packages for the framework
    using CSV
    using DelimitedFiles
    using DataFrames
    using JuMP
    using Random
    using Shapefile
    using Plots
    using StatsPlots
    using Dates
    using Distributions
    using Gurobi

# necessary functions for the framework
    include("part1_adjacency_matrix.jl")
    include("part1_workload_calculation.jl")
    include("part1_sets_m_n.jl")
    include("part1_districting_model.jl")
    include("part1_normalize.jl")
    include("part1_equalize.jl")
    include("part1_plot_districts.jl")
    include("part2_small_functions.jl")
    include("part2_capacity_heuristic.jl")
    include("part2_ressource_flow_matrix.jl")
    include("part2_fill_queues.jl")
    include("part2_incident_dispatch.jl")
    include("part2_determine_patrol_location.jl")
    include("part2_fill_exchange_queue.jl")
    include("part2_fulfill_exchange_queue.jl")
    include("part2_prepare_next_minute.jl")
    include("part2_simulation.jl")
    include("evaluation_ressource_status.jl")
    include("evaluation_prepare_results.jl")
    include("evaluation_plot_results.jl")
    include("part2_prepare_stage.jl")


text:
export embedded_number_field

################################################################################
#
#  Types
#
################################################################################

mutable struct EmbeddedNumField{S, E} <: Field
  field::S
  embedding::E
end

mutable struct EmbeddedNumFieldElem{T} <: FieldElem
  parent
  element::T

  function EmbeddedNumFieldElem{T}(parent, element::T) where {T}
    @assert Hecke.parent(element) === number_field(parent)
    return new{T}(parent, element)
  end
end

function (E::EmbeddedNumField{S})(x::T) where {S, T <: NumFieldElem}
  if elem_type(S) === T
    @assert parent(x) === number_field(E)
    return EmbeddedNumFieldElem{T}(E, x)
  else
    return E(number_field(E)(x))
  end
end

function (E::EmbeddedNumField)(x::EmbeddedNumFieldElem)
  @assert parent(x) === E
  return x
end

function (E::EmbeddedNumField)(x::RingElement)
  return E(number_field(E)(x))
end

function hash(x::EmbeddedNumFieldElem, h::UInt)
  return xor(hash(data(x), h), 0x5cccaf5d1346dc53%UInt)
end

function Base.deepcopy_internal(x::EmbeddedNumFieldElem, id::IdDict)
  return parent(x)(Base.deepcopy_internal(data(x), id))
end

number_field(K::EmbeddedNumField) = K.field

embedding(K::EmbeddedNumField) = K.embedding

parent(x::EmbeddedNumFieldElem{T}) where {T} = x.parent::parent_type(x)

elem_type(::Type{EmbeddedNumField{S, E}}) where {S, E} = EmbeddedNumFieldElem{elem_type(S)}

elem_type(K::EmbeddedNumField{S, E}) where {S, E} = elem_type(EmbeddedNumField{S, E})

parent_type(::Type{EmbeddedNumFieldElem{T}}) where {T} = EmbeddedNumField{parent_type(T), embedding_type(parent_type(T))}

parent_type(x::EmbeddedNumFieldElem{T}) where {T} = parent_type(EmbeddedNumFieldElem{T})

data(x::EmbeddedNumFieldElem) = x.element

function embedded_field(K::SimpleNumField, i::NumFieldEmb)
  @assert number_field(i) === K
  E = EmbeddedNumField(K, i)
  return E, E(gen(K))
end

function embedded_field(K::NonSimpleNumField, i::NumFieldEmb)
  @assert number_field(i) === K
  E = EmbeddedNumField(K, i)
  return E, E.(gens(K))
end

(E::EmbeddedNumField)() = zero(E)

unary_ops = [:(-)]

binary_ops = [:(+), :(*), :(-), :(div), :(divexact), :(//)]

for b in unary_ops
  @eval begin
    function ($b)(x::EmbeddedNumFieldElem)
      return parent(x)($(b)(data(x)))
    end
  end
end

for b in binary_ops
  @eval begin
    function ($b)(x::EmbeddedNumFieldElem, y::EmbeddedNumFieldElem)
      return parent(x)($(b)(data(x), data(y)))
    end
  end
end

function ==(x::EmbeddedNumFieldElem, y::EmbeddedNumFieldElem)
  return ==(data(x), data(y))
end

iszero(x::EmbeddedNumFieldElem) = iszero(data(x))

isone(x::EmbeddedNumFieldElem) = isone(data(x))

isunit(x::EmbeddedNumFieldElem) = isunit(data(x))

zero(E::EmbeddedNumField) = E(zero(number_field(E)))

one(E::EmbeddedNumField) = E(one(number_field(E)))

# Now the ordering

function isless(x::EmbeddedNumFieldElem, y::EmbeddedNumFieldElem)
  i = embedding(parent(x))
  # Need to exclude equality
  if x == y
    return false
  end
  p = 32
  xe, ye = i(data(x), p), i(data(y), p)
  while overlaps(xe, ye)
    p = Int(floor(p * 1.42))
    xe, ye = i(data(x), p), i(data(y), p)
  end
  if real(xe) < real(ye)
    return true
  else
    @assert real(xe) > real(ye)
    return false
  end
end

isless(x::EmbeddedNumFieldElem, y) = isless(x, parent(x)(y))

isless(x, y::EmbeddedNumFieldElem) = isless(parent(y)(x), y)

# Support comparing with floats

function isless(x::EmbeddedNumFieldElem, y::AbstractFloat)
  i = embedding(parent(x))
  # Need to exclude equality
  p = 32
  xe = i(data(x), p)
  # check if y is "equal" to x as a rational
  if isrational(data(x))
    xq = QQ(data(x))
    d = denominator(xq)
    if isone(d)
      yd = y
    else
      yd = d * y
    end
    if isinteger(yd) && numerator(xq) == yd
      return false
    end
  end

  while contains(real(xe), BigFloat(y))
    p = Int(floor(p * 1.42))
    xe = i(data(x), p)
  end
  if real(xe) < y
    return true
  else
    @assert real(xe) > y
    return false
  end
end

function isless(y::AbstractFloat, x::EmbeddedNumFieldElem)
  i = embedding(parent(x))
  # Need to exclude equality
  p = 32
  xe = i(data(x), p)
  # check if y is "equal" to x as a rational
  if isrational(data(x))
    xq = QQ(data(x))
    d = denominator(xq)
    if isone(d)
      yd = y
    else
      yd = d * y
    end
    if isinteger(yd) && numerator(xq) == yd
      return false
    end
  end

  while contains(real(xe), BigFloat(y))
    p = Int(floor(p * 1.42))
    xe = i(data(x), p)
  end
  if real(xe) > y
    return true
  else
    @assert real(xe) < y
    return false
  end
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, E::EmbeddedNumField)
  print(io, "Embedded field\n$(number_field(E))\nat\n$(embedding(E))")
end

# I overload this to get the value at the embedding
# (But I don't want this in the general printing routines, e.g., over polynomial
# rings.)
function Base.show(io::IO, ::MIME"text/plain", x::EmbeddedNumFieldElem)
  print(io, "$(data(x))")
  a = Float64(real(embedding(parent(x))(data(x), 32)))
  print(io, @sprintf " (%.2f)" a)
end

function Base.show(io::IO, x::EmbeddedNumFieldElem)
  print(io, "$(data(x))")
end

function AbstractAlgebra.expressify(x::EmbeddedNumFieldElem; context = nothing)
  AbstractAlgebra.expressify(data(x), context = context)
end

################################################################################
#
#  Constructors
#
################################################################################

function embedded_number_field(f::Union{fmpq_poly, fmpz_poly}, r::Union{AbstractFloat, Tuple}, var = "a")
  K, a = number_field(f, var)
  r = real_embedding(K, r)
  return embedded_field(K, r)
end

function embedded_number_field(f::Vector{<:Union{fmpq_poly, fmpz_poly}}, r::Vector{<:Union{AbstractFloat, Tuple}}, var = "a")
  K, a = number_field(f, var)
  r = real_embedding(K, r)
  return embedded_field(K, r)
end

################################################################################
#
#  Creation of elements
#
################################################################################

(E::EmbeddedNumField)(v::Vector) = E(number_field(E)(v))

################################################################################
#
#  Promote rule
#
################################################################################

function AbstractAlgebra.promote_rule(::Type{EmbeddedNumFieldElem{T}}, ::Type{S}) where {T <: NumFieldElem, S <: RingElement}
  if T === S
    return Union{}
  else
    if AbstractAlgebra.promote_rule(T, S) === T
      return EmbeddedNumFieldElem{T}
    else
      return Union{}
    end
  end
end

function AbstractAlgebra.promote_rule(::Type{EmbeddedNumFieldElem{T}}, ::Type{EmbeddedNumFieldElem{T}}) where {T <: NumFieldElem}
  return EmbeddedNumFieldElem{T}
end


text:
module TestProperties
using Test
using CovarianceFunctions
using CovarianceFunctions: input_trait, DotProductInput, IsotropicInput, StationaryLinearFunctionalInput, GenericInput
using CovarianceFunctions: EQ, RQ, Exp, Dot, ExponentialDot, Poly, Line, Cos

using LinearAlgebra

@testset "properties" begin
    dot_kernels = [Dot(), Dot()^3, ExponentialDot(), Line(1.), Poly(5, 1.)]
    for k in dot_kernels
        @test input_trait(k) isa DotProductInput
    end

    iso_kernels = [EQ(), RQ(1.), Exp()]
    for k in iso_kernels
        @test input_trait(k) isa IsotropicInput
    end

    k = CovarianceFunctions.NeuralNetwork()
    @test input_trait(k) isa GenericInput

    # testing that constant kernels don't confuse the input_trait inference
    @test input_trait(1*EQ() + 1) isa IsotropicInput
    @test input_trait(1*EQ() + 2 + RQ(1.)*1) isa IsotropicInput

    @test input_trait(1*Dot() + 1) isa DotProductInput
    @test input_trait(1*Dot() + 2 + Dot()^2*1) isa DotProductInput

    w = randn()
    @test input_trait(1*Cos(w) + 1) isa StationaryLinearFunctionalInput
    @test input_trait(1*Cos(w) + 2 + Cos(w)^2*1) isa StationaryLinearFunctionalInput
end

end


text:
##############################################################################
##
## NA's via the NAtype
##
## Inspirations:
##  * R's NA's
##  * Panda's discussion of NA's:
##    http://pandas.pydata.org/pandas-docs/stable/missing_data.html
##  * NumPy's analysis of the issue:
##    https://github.com/numpy/numpy/blob/master/doc/neps/missing-data.rst
##
## NAtype is a composite type representing missingness:
## * An object of NAtype can be generated by writing NA
##
##############################################################################

type NAtype
end

const NA = NAtype()

Base.show(io::IO, x::NAtype) = print(io, "NA")

type NAException <: Exception
    msg::String
end
NAException() = NAException("NA found")

Base.length(x::NAtype) = 1
Base.size(x::NAtype) = ()
Base.ndims(x::NAtype) = 0

isna(x::NAtype) = true
isna(x::Any) = false

# TODO: Rethink this rule
Base.promote_rule{T}(::Type{T}, ::Type{NAtype} ) = T

# TODO: Try to remove this line
baseval(::Type{NAtype}) = NA


text:
module Trivia

import HTTP
import JSON
using Base64: base64decode
using Random: shuffle!

export Token, request_questions, present, prompt, play_round

# This is the main source file for the Trivia.jl package. It defines the `Trivia` Julia module,
# imports a few dependencies (HTTP, JSON, Base64, and Random), and lists a few helpful types
# and methods that will be exported.

# All of the actual work is split out into the following files:
include("requests.jl")
include("tokens.jl")
include("questions.jl")
include("interactive.jl")

end # module


text:
mutable struct NfRelOrdToFqMor{T, S, U} <: Map{NfRelOrd{T, S, U}, FqFiniteField, HeckeMap, NfRelOrdToFqMor}
  header::MapHeader{NfRelOrd{T, S, U}, FqFiniteField}
  poly_of_the_field::fq_poly
  P::NfRelOrdIdl{T, S, U}

  function NfRelOrdToFqMor{T, S, U}(O::NfRelOrd{T, S, U}, P::NfRelOrdIdl{T, S, U}) where {T, S, U}
    z = new{T, S, U}()
    z.P = P
    p = minimum(P, copy = false)
    F, mF = ResidueField(order(p), p)
    mmF = extend(mF, nf(order(p)))
    Fx, = PolynomialRing(F, "x", cached = false)
    if is_index_divisor(O, p)
      A, OtoA = AlgAss(O, P, p)
      AtoO = pseudo_inv(OtoA)
      x = rand(A)
      h = minpoly(x)
      while degree(h) != dim(A)
        x = rand(A)
        h = minpoly(x)
      end
      # F and base_ring(h) are the same as in "==" but not as in "==="
      hh = Fx()
      ccall((:fq_poly_set, libflint), Nothing, (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}), hh, h, F)
      z.poly_of_the_field = hh
      FF, mFF = field_extension(hh)

      M = zero_matrix(F, dim(A), dim(A))
      t = one(A)
      for i = 1:dim(A)
        for j = 1:dim(A)
          M[j, i] = t.coeffs[j]
        end
        t = t*x
      end
      Minv = inv(M)

      function _image_index_div(a::NfRelOrdElem)
        b = OtoA(a)
        bb = zero_matrix(F, dim(A), 1)
        for i = 1:dim(A)
          bb[i, 1] = b.coeffs[i]
        end
        @assert mod(AtoO(A([ bb[i, 1] for i = 1:dim(A) ])), P) == mod(a, P)
        bb = Minv*bb
        g = Fx([ bb[i, 1] for i = 1:dim(A) ])
        return mFF(g)
      end

      function _preimage_index_div(a::fq)
        g = pseudo_inv(mFF)(a)
        c = zero_matrix(F, dim(A), 1)
        for i = 1:dim(A)
          c[i, 1] = coeff(g, i - 1)
        end
        c = M*c
        b = A([ c[i, 1] for i = 1:dim(A) ])
        return AtoO(b)
      end
      z.header = MapHeader(O, FF, _image_index_div, _preimage_index_div)
    else
      h = P.non_index_div_poly
      # F and base_ring(h) are the same as in "==" but not as in "==="
      hh = Fx()
      ccall((:fq_poly_set, libflint), Nothing, (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}), hh, h, F)
      z.poly_of_the_field = hh
      FF, mFF = field_extension(hh)

      function _image(x::NfRelOrdElem)
        f = parent(nf(O).pol)(elem_in_nf(x))
        if iszero(f)
          ff = Fx()
        else
          ff = Fx([ mmF(coeff(f, i)) for i = 0:degree(f) ])
        end
        return image(mFF, ff)
      end

      function _preimage(x::fq)
        f = preimage(mFF, x)
        immF = pseudo_inv(mmF)
        y = nf(O)([ immF(coeff(f, i)) for i = 0:degree(f) ])
        return O(y)
      end
      z.header = MapHeader(O, FF, _image, _preimage)
    end
    return z
  end
end

function extend(f::NfRelOrdToFqMor{T, S}, K::NfRel{T}) where {T, S}
  nf(domain(f)) != K && error("Number field is not the number field of the order")

  g = NfRelToFqMor{T}()
  g.header.domain = K
  g.header.codomain = f.header.codomain

  O = domain(f)
  P = f.P
  u = anti_uniformizer(P)

  function _image(x::NfRelElem{T})
    m = denominator(x, O)
    a = O(m*x)
    b = O(m)
    l = valuation(m, P)
    if l != 0
      a = O(elem_in_nf(a)*u^l)
      b = O(elem_in_nf(b)*u^l)
    end
    return f(a)//f(b)
  end

  function _preimage(x::fq)
    return elem_in_nf(preimage(f, x))
  end

  g.header.image = _image
  g.header.preimage = _preimage

  return g
end

mutable struct RelOrdToAlgAssMor{S, T} <: Map{S, AlgAss{T}, HeckeMap, RelOrdToAlgAssMor}
  header::MapHeader{S, AlgAss{T}}

  function RelOrdToAlgAssMor{S, T}(O::S, A::AlgAss{T}, _image::Function, _preimage::Function) where { S <: Union{ NfRelOrd, AlgAssRelOrd }, T }
    z = new{S, T}()
    z.header = MapHeader(O, A, _image, _preimage)
    return z
  end
end

function RelOrdToAlgAssMor(O::Union{ NfRelOrd, AlgAssRelOrd }, A::AlgAss{T}, _image, _preimage) where {T}
  return RelOrdToAlgAssMor{typeof(O), T}(O, A, _image, _preimage)
end

mutable struct RelOrdQuoMap{T1, T2, T3, S} <: Map{T1, RelOrdQuoRing{T1, T2, T3}, HeckeMap, RelOrdQuoMap}
  header::MapHeader{T1, RelOrdQuoRing{T1, T2, T3}}

  function RelOrdQuoMap{T1, T2, T3, S}(O::T1, Q::RelOrdQuoRing{T1, T2, T3}) where { T1, T2, T3, S }
    z = new{T1, T2, T3, S}()

    _image = function (x::S)
      return Q(x)
    end

    _preimage = function (x::RelOrdQuoRingElem{T1, T2, T3, S})
      return x.elem
    end

    z.header = MapHeader(O, Q, _image, _preimage)
    return z
  end
end

function RelOrdQuoMap(O::T1, Q::RelOrdQuoRing{T1, T2, T3}) where { T1, T2, T3 }
  S = elem_type(O)
  return RelOrdQuoMap{T1, T2, T3, S}(O, Q)
end

mutable struct NfRelOrdToRelFinFieldMor{S, T} <: Map{S, RelFinField{T}, HeckeMap, NfRelOrdToRelFinFieldMor}
  header::MapHeader{S, RelFinField{T}}
  poly_of_the_field
  P
  map_subfield::Union{NfOrdToFqMor, NfRelOrdToRelFinFieldMor}

    function NfRelOrdToRelFinFieldMor{S, T}(O::S, P, mapsub::NfOrdToFqMor) where {S, T <: fq}
    z = new{S, T}()
    z.P = P
    z.map_subfield = mapsub
    p = minimum(P, copy = false)
    FK, mK = codomain(mapsub), mapsub
    mmK = extend(mK, nf(order(p)))
    FKx, = PolynomialRing(FK, "x", cached = false)
    if is_index_divisor(O, p)
      A, OtoA = AlgAss(O, P, p)
      AtoO = pseudo_inv(OtoA)
      x = rand(A)
      h = minpoly(x)
      while degree(h) != dim(A)
        x = rand(A)
        h = minpoly(x)
      end
      hh = FKx()
      ccall((:fq_poly_set, libflint), Nothing, (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}), hh, h, FK)
      z.poly_of_the_field = hh

      FE = RelFinField(hh, :v)
      FEabs, FEabstoFE = Hecke.absolute_field(FE, cached = false)
      FE2, mE2 = field_extension(hh)
      FE2toFEabs = hom(FE2, FEabs, gen(FEabs))
      mE = compose(mE2, compose(FE2toFEabs, FEabstoFE))

      M = zero_matrix(FK, dim(A), dim(A))
      t = one(A)
      for i = 1:dim(A)
        for j = 1:dim(A)
          M[j, i] = t.coeffs[j]
        end
        t = t*x
      end
      Minv = inv(M)

      function _image_index_div(a::NfRelOrdElem)
        b = OtoA(a)
        bb = zero_matrix(FK, dim(A), 1)
        for i = 1:dim(A)
          bb[i, 1] = b.coeffs[i]
        end
        @assert mod(AtoO(A([ bb[i, 1] for i = 1:dim(A) ])), P) == mod(a, P)
        bb = Minv*bb
        g= FKx([ bb[i, 1] for i = 1:dim(A) ])
        return mE(g)
      end

      function _preimage_index_div(a::RelFinFieldElem)
        g = pseudo_inv(mE)(a)
        c = zero_matrix(FK, dim(A), 1)
        for i = 1:dim(A)
          c[i,1] = coeff(g, i-1)
        end
        c = M*c
        b = A([c[i,1] for i=1:dim(A)])
        return AtoO(b)
      end
      z.header = MapHeader(O, FE, _image_index_div, _preimage_index_div)
    else
      h = P.non_index_div_poly
      hh = FKx()
      ccall((:fq_poly_set, libflint), Nothing, (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}), hh, h, FK)
      z.poly_of_the_field = hh

      FE = RelFinField(hh, :v)
      FEabs, FEabstoFE = Hecke.absolute_field(FE, cached = false)
      FE2, mE2 = field_extension(hh)
      FE2toFEabs = hom(FE2, FEabs, gen(FEabs))
      mE = compose(mE2, compose(FE2toFEabs, FEabstoFE))

      function _image(x::NfRelOrdElem)
        f = parent(nf(O).pol)(elem_in_nf(x))
        if iszero(f)
          ff = FKx()
        else
          ff = FKx([ mmK(coeff(f, i)) for i =0:degree(f)])
        end
        return image(mE, ff)
      end

      function _preimage(x::RelFinFieldElem)
        f = preimage(mE, x)
        immK = pseudo_inv(mmK)
        y = nf(O)([ immK(coeff(f,i)) for i=0:degree(f)])
        return O(y)
      end
      z.header = MapHeader(O, FE, _image, _preimage)
    end
    return z
  end

  function NfRelOrdToRelFinFieldMor{S, T}(O::S, P, mapsub::NfRelOrdToRelFinFieldMor) where {S, T <: RelFinFieldElem}
    z = new{S, T}()
    z.P = P
    z.map_subfield = mapsub
    p = minimum(P, copy = false)
    FK, mK = codomain(mapsub), mapsub
    mmK = extend(mK, nf(order(p)))
    FKx, = PolynomialRing(FK, "x", cached = false)
    FKabs, FKabstoFK = Hecke.absolute_field(FK, cached = false)
    FKtoFKabs = pseudo_inv(FKabstoFK)
    FKabsz, _ = PolynomialRing(FKabs, "z", cached = false)
    if is_index_divisor(O, p)
      A, OtoA = AlgAss(O, P, p)
      AtoO = pseudo_inv(OtoA)
      x = rand(A)
      h = minpoly(x)
      while degree(h) != dim(A)
        x = rand(A)
        h = minpoly(x)
      end
      hh = FKabsz()
      ccall((:fq_poly_set, libflint), Nothing, (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}), hh, h, FKabs)
      h = map_coefficients(FKabstoFK, hh)
      h = FKx(collect(coefficients(h)))
      z.poly_of_the_field = h

      FE = RelFinField(h, :v)
      FEabs, FEabstoFE = Hecke.absolute_field(FE, cached = false)
      FKxtoFKabsz = MapFromFunc(f -> FKabsz(FKtoFKabs.(collect(coefficients(f)))), FKx, FKabsz)
      FE2, mE2 = field_extension(hh)
      FE2toFEabs = hom(FE2, FEabs, gen(FEabs))
      mE = compose(FKxtoFKabsz,compose(mE2, compose(FE2toFEabs, FEabstoFE)))

      M = zero_matrix(FK, dim(A), dim(A))
      t = one(A)
      for i = 1:dim(A)
        for j = 1:dim(A)
          M[j, i] = t.coeffs[j]
        end
        t = t*x
      end
      Minv = inv(M)

      function _image_index_div(a::NfRelOrdElem)
        b = OtoA(a)
        bb = zero_matrix(FK, dim(A), 1)
        for i = 1:dim(A)
          bb[i, 1] = b.coeffs[i]
        end
        @assert mod(AtoO(A([ bb[i, 1] for i = 1:dim(A) ])), P) == mod(a, P)
        bb = Minv*bb
        g= FKx([ bb[i, 1] for i = 1:dim(A) ])
        return mE(g)
      end

      function _preimage_index_div(a::RelFinFieldElem)
        g = pseudo_inv(mE)(a)
        c = zero_matrix(FK, dim(A), 1)
        for i = 1:dim(A)
          c[i,1] = coeff(g, i-1)
        end
        c = M*c
        b = A([c[i,1] for i=1:dim(A)])
        return AtoO(b)
      end
      z.header = MapHeader(O, FE, _image_index_div, _preimage_index_div)
    else
      h = P.non_index_div_poly
      hh = FKabsz()
      ccall((:fq_poly_set, libflint), Nothing, (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}), hh, h, FKabs)
      h = map_coefficients(FKabstoFK, hh)
      h = FKx(collect(coefficients(h)))
      z.poly_of_the_field = h

      FE = RelFinField(h, :v)
      FEabs, FEabstoFE = Hecke.absolute_field(FE, cached = false)
      FKxtoFKabsz = MapFromFunc(f -> FKabsz(FKtoFKabs.(collect(coefficients(f)))), FKx, FKabsz)
      FE2, mE2 = field_extension(hh)
      FE2toFEabs = hom(FE2, FEabs, gen(FEabs))
      mE = compose(FKxtoFKabsz,compose(mE2, compose(FE2toFEabs, FEabstoFE)))

      function _image(x::NfRelOrdElem)
        f = parent(nf(O).pol)(elem_in_nf(x))
        if iszero(f)
          ff = FKx()
        else
          ff = FKx([ mmK(coeff(f, i)) for i =0:degree(f)])
        end
        return image(mE, ff)
      end

      function _preimage(x::RelFinFieldElem)
        f = preimage(mE, x)
        immK = pseudo_inv(mmK)
        y = nf(O)([ immK(coeff(f,i)) for i=0:degree(f)])
        return O(y)
      end
      z.header = MapHeader(O, FE, _image, _preimage)
    end
    return z
  end

end

mutable struct NfRelToRelFinFieldMor{S, T} <: Map{S, Hecke.RelFinField{T}, HeckeMap, NfRelToRelFinFieldMor}
  header::MapHeader{S, Hecke.RelFinField{T}}

  function NfRelToRelFinFieldMor{S, T}() where {S <: NfRel, T <: FinFieldElem}
    z = new{S, T}()
    z.header = MapHeader{S, Hecke.RelFinField{T}}()
    return z
  end
end

function extend(f::NfRelOrdToRelFinFieldMor{S, T}, E::NfRel) where {S, T}
  nf(domain(f)) != E && error("Number field is not the number field of the order")

  g = NfRelToRelFinFieldMor{typeof(E), T}()
  g.header.domain = E
  g.header.codomain = f.header.codomain

  OE = domain(f)
  P = f.P
  u = anti_uniformizer(P)

  function _image(x::NfRelElem)
    m = denominator(x, OE)
    a = OE(m*x)
    b = OE(m)
    l = valuation(m, P)
    if l != 0
      a = OE(elem_in_nf(a)*u^l)
      b = OE(elem_in_nf(b)*u^l)
    end
    return f(a)//f(b)
  end

  function _preimage(x::RelFinFieldElem)
    return elem_in_nf(preimage(f, x))
  end

  g.header.image = _image
  g.header.preimage = _preimage

  return g
end


text:
## The birthday paradox in theory and in practice.

## Let's simulate the classic birthday paradox.
# In this simulation, we will draw birthdays at random for a group
# of N people, and see if there are any two birthdays that
# are the same.

"""
N is the number of people

This function returns true if two people had the same birthday.
"""
function simple_birthday_simulation(N::Int)
    birthdays = rand(1:365, N)
    # so each birthday is an integer from 1 to 365.
    # we need to check if all of the birthdays are unique
    return !allunique(birthdays) # ! reverses true -> false, false -> true
end

simple_birthday_simulation(20)

## Let's look at what happens as we vary the number of people.
ntrials = 100
prob_dup_birthday = zeros(0) # this is an empty array
for people=5:40
    numdup = 0
    for t=1:ntrials
        numdup += simple_birthday_simulation(people)
    end
    push!(prob_dup_birthday, numdup/ntrials)
end
# make a quick people - by - prob table
display([5:40 prob_dup_birthday])
# but also make a plot
using Plots
plot(5:40, prob_dup_birthday)
xlabel!("People")
ylabel!("Prob. of shared birthday")

## Some ideas to improve
# IDEA -- add confidence intervals
# IDEA -- add "exact" curve based on prob.

## Let's use real data.
using DelimitedFiles
birthdays = readdlm("1-julia-intro/birthday.txt")
plot(birthdays[:,2])

## First step, let's simulate a random birthday from data.
using StatsBase
bday_weights = fweights(Int.(birthdays[:,2]))
sample(bday_weights)

## Let's do the same demo
sample_bdays = [sample(fweights([2,4,3,1])) for i=1:1000]
sum(sample_bdays .== 4)

## Let's implement the sampling ourselves too...
"""
N is the number of people
data is an array of length number of days with the number of people
that were born on that day.

Internally, what this does is convert data to a cum. prob. distribution
and then we use a random [0,1] sample to find where in the CDF we fall.

The idea of how this works is:
give probs [0.2 0.4 0.3 0.1] (i.e. we have 4 birthdays)
then we are going to make an array
[0.2, 0.6, 0.9, 1.0]
and then we are going to guess a random [0,1] variable
and see which bin it falls into.
if the rv. is less than 0.2, we will call it birthday 1
if the rv. is less between 0.2, 0.6, we will call it birthday 2
if the rv. is less between 0.6, 0.9, we will call it birthday 3
if the rv. is less between 0.9, 1, we will call it birthday 4
but this will work for even larger lists too.
"""
function sample_birthdays(N::Int, data)
    probs = data/sum(data)
    cumdist = cumsum(probs)
    bdays = zeros(Int, N)
    for i=1:N
        rv = rand()
        for j=1:length(cumdist)
            if rv <= cumdist[j]
                bdays[i] = j
                break # need to stop when we find the birthday
            end
        end
    end
    return bdays
end
bdays = sample_birthdays(1000, [2,4,3,1])
# Let's check their frequency
sum(bdays .== 3) # count the number of "2" birthdays in our test
## IDEA - use binary search to find the right birthday faster.

## Let's simulate the data-based shared birthday probability
function data_birthday_simulation(N::Int, data)
    birthdays = sample_birthdays(N, data)
    # so each birthday is an integer from 1 to 365.
    # we need to check if all of the birthdays are unique
    return !allunique(birthdays) # ! reverses true -> false, false -> true
end

data_birthday_simulation(20, Int.(birthdays[:,2]))

ntrials = 100
prob_dup_birthday = zeros(0) # this is an empty array
for people=5:40
    numdup = 0
    for t=1:ntrials
        numdup += data_birthday_simulation(people, Int.(birthdays[:,2]))
    end
    push!(prob_dup_birthday, numdup/ntrials)
end
# make a quick people - by - prob table
display([5:40 prob_dup_birthday])
# but also make a plot
using Plots
plot(5:40, prob_dup_birthday)
xlabel!("People")
ylabel!("Prob. of shared birthday based on data")

## TODO Do a test to show that this curve is higher than the uniform
# case.

## Use Julia's sampler
function data_birthday_simulation_julia(N::Int, data)
    w = fweights(data)
    birthdays = [sample(w) for i=1:N]
    # so each birthday is an integer from 1 to 365.
    # we need to check if all of the birthdays are unique
    return !allunique(birthdays) # ! reverses true -> false, false -> true
end
ntrials = 100
prob_dup_birthday = zeros(0) # this is an empty array
for people=5:40
    numdup = 0
    for t=1:ntrials
        numdup += data_birthday_simulation_julia(people, Int.(birthdays[:,2]))
    end
    push!(prob_dup_birthday, numdup/ntrials)
end
# make a quick people - by - prob table
display([5:40 prob_dup_birthday])
# but also make a plot
using Plots
plot(5:40, prob_dup_birthday)
xlabel!("People")
ylabel!("Prob. of shared birthday based on data")


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

abstract type SpatialVolume{T<:Real} end

export SpatialVolume

Base.eltype(b::SpatialVolume{T}) where T = T

Base.rand(rng::AbstractRNG, vol::SpatialVolume) =
    rand!(rng, vol, Vector{float(eltype(vol))}(ndims(vol)))

Base.rand(rng::AbstractRNG, vol::SpatialVolume, n::Integer) =
    rand!(rng, vol, Matrix{float(eltype(vol))}(ndims(vol), n))

doc"""
    log_volume(vol::SpatialVolume)

Get the logarithm of the volume of the space in `vol`.
"""
function log_volume end
export log_volume

doc"""
    fromuhc!(Y::VecOrMat, X::VecOrMat, vol::SpatialVolume)

Bijective transformation of coordinates `X` within the unit hypercube to
coordinates `Y` in `vol`. If `X` and `Y` are matrices, the transformation is
applied to the column vectors. Use `Y === X` to transform in-place.

Use `inv(fromuhc!)` to get the the inverse transformation.
"""
function fromuhc! end
export fromuhc!

function inv_fromuhc! end

Base.inv(::typeof(fromuhc!)) = inv_fromuhc!
Base.inv(::typeof(inv_fromuhc!)) = fromuhc!

doc"""
    fromuhc(X::VecOrMat, vol::SpatialVolume)

Bijective transformation from unit hypercube to `vol`. See `fromuhc!`.

Use `inv(fromuhc)` to get the the inverse transformation.
"""
function fromuhc(X::VecOrMat, vol::SpatialVolume)
    fromuhc!(similar(X), X, vol)
end
export fromuhc

function inv_fromuhc(X::VecOrMat, vol::SpatialVolume)
    inv_fromuhc!(similar(X), X, vol)
end

Base.inv(::typeof(fromuhc)) = inv_fromuhc
Base.inv(::typeof(inv_fromuhc)) = fromuhc

struct HyperRectVolume{T<:Real} <: SpatialVolume{T}
    # ToDo: Use origin and widths instead of lo and hi, as in GeometryTypes.HyperRectangle?
    lo::Vector{T}
    hi::Vector{T}

    function HyperRectVolume{T}(lo::Vector{T}, hi::Vector{T}) where {T<:Real}
        (indices(lo) != indices(hi)) && throw(ArgumentError("lo and hi must have the same indices"))
        new{T}(lo, hi)
    end
end

export HyperRectVolume

HyperRectVolume(lo::Vector{T}, hi::Vector{T}) where {T<:Real} = HyperRectVolume{T}(lo, hi)

Base.ndims(vol::HyperRectVolume) = size(vol.lo, 1)

function Base.isempty(vol::HyperRectVolume)
        @inbounds for i in eachindex(vol.lo, vol.hi)
             (vol.lo[i] > vol.hi[i]) && return true
        end
        isempty(vol.lo)
end

Base.similar(vol::HyperRectVolume) = HyperRectVolume(similar(vol.lo), similar(vol.hi))

Base.in(x::AbstractVector, vol::HyperRectVolume) =
    _all_lteq(vol.lo, x, vol.hi)

# TODO: Remove?
# function Base.in(X::AbstractMatrix, vol::HyperRectVolume, j::Integer)
#     lo = vol.lo
#     hi = vol.hi
#     for i in indices(X, 1)
#         (lo[i] <= X[i, j] <= hi[i]) || return false
#     end
#     return true
# end

# ToDo:
# Base.in(x::Matrix, vol::HyperRectVolume) =

function Base.intersect(a::HyperRectVolume, b::HyperRectVolume)
    c = similar(a)
    c.lo .= max.(a.lo, b.lo)
    c.hi .= min.(a.hi, b.hi)
    c
end

function Base.rand!(rng::AbstractRNG, vol::HyperRectVolume, x::StridedVecOrMat{<:Real})
    rand!(rng, x)
    x .= x .* (vol.hi - vol.lo) .+ vol.lo # TODO: Avoid memory allocation
end

function log_volume(vol::HyperRectVolume{T}) where {T}
    R = promote_type(Float64, T)
    S = promote_type(Double{Float64}, T)
    s = zero(S)
    hi = vol.hi
    lo = vol.lo
    @assert indices(hi) == indices(lo)
    @inbounds @simd for i in eachindex(hi)
        d = max(zero(R), R(hi[i]) - R(lo[i]))
        s += JuliaLibm.log(d)
    end
    R(s)
end

# ToDo:
# log_intersect_volume(a::HyperRectVolume{T}, b::HyperRectVolume{T}) where {T} = ...

doc"""
    fromuhc!(Y::VecOrMat, X::VecOrMat, vol::HyperRectVolume)

Linear bijective transformation from unit hypercube to hyper-rectangle `vol`.
"""
function fromuhc!(Y::VecOrMat, X::VecOrMat, vol::HyperRectVolume)
    _all_in_ui(X) || throw(ArgumentError("X not in unit hypercube"))
    Y .= unsafe_fromui.(X, vol.lo, vol.hi)
end

function inv_fromuhc!(Y::VecOrMat, X::VecOrMat, vol::HyperRectVolume)
    X in vol || throw(ArgumentError("X not in vol"))
    Y .= unsafe_inv_fromui.(X, vol.lo, vol.hi)
end


text:
export plotModel, cutBoundaryLayer

"""
function plotModel

	plots model on a regular mesh

	Requirement: PyPlot must be installed.

	Parameters:
	m 				:: model. 2D or 3D Array.
	includeMeshInfo	:: A boolean indicating whether to include distances info (in KMs) in the image. Relevant only in 2D.
	cutPad			:: A flag for cutting the outmost layer. If cutPad == 0, there's no cutting.
	limits			:: Maximum and minimum values in the colorbar.
	filename        :: A filename for automatically saving the image as a file.
"""
function plotModel(m::Union{Array{Float64},Array{Float32}};includeMeshInfo::Bool=false,M_regular = [],cutPad::Int64 = 0,limits = [],filename="",figTitle="")

if limits!=[]
	vmin = limits[1];
	vmax = limits[2];
else
	vmin = minimum(m);
	vmax = maximum(m);
end

limits = tuple([vmax,vmin]...);
if cutPad > 0
	m,M_regular = cutBoundaryLayer(m,M_regular,cutPad);
end
	if length(size(m))==2
		T = m';
		imshow(T, clim = limits,cmap = "jet"); colorbar();
		if figTitle != ""
			title(figTitle)
		end
		if includeMeshInfo
			Omega = M_regular.domain;
			tics = 0:2:floor(Int64,Omega[2])
			# xticks(linspace(0,size(m,1)*tics[end]/Omega[2],length(tics)),tics);
			xlabel("Lateral distance (km)");
			tics = 0:floor(Int64,Omega[4])
			# println(tics)
			# println(linspace(0,size(m,2)*tics[end]/Omega[4],length(tics)))
			# yticks(linspace(0,size(m,2)*tics[end]/Omega[4],length(tics)),tics);
			ylabel("Depth (km)")
		end
		if filename != ""
			savefig(string(filename[1:end-4],".png"));
		end
	elseif length(size(m))==3
		lin = zeros(Int64,16);
		v = m;
		for k=1:16
			lin[k] = Int(round(k*(size(m,2)/16)));
		end
		for k=1:16
#			subplot(string("44",k))
			subplot(4,4,k)
			pic = reshape(v[:,lin[k],:],size(m,1),size(m,3))';
			imshow(pic,clim = limits); title(string(figTitle,": frame",lin[k]));colorbar()
		end
		if filename != ""
			savefig(string(filename[1:end-4],".png"));
		end
	end
return;
end

"""
function cutBoundaryLayer

	Removes the boundary layer from the model and returns the cut model and a corresponding new mesh.
	Used mostly for plotting - sometimes we prefer not to show the boundary of the model.

	The "boundary" here is all the boundaries of the regular domain except the top surface.

	Parameters:
	m 				:: model. 2D or 3D Array.
	M_regular		:: A corresponding regular mesh.
	pad			    :: The thickness of the layer in cells.

"""
function cutBoundaryLayer(m::Array{Float64},M_regular,pad::Int64)
if pad<=0
	return m,M_regular;
end
Omega = M_regular.domain;
if length(size(m))==2
	mnew = m[pad+1:end-pad,1:end-pad];
	nnew = collect(size(mnew)).-1;
	OmegaNew = [Omega[1],Omega[2] - 2*pad*M_regular.h[1],Omega[3],Omega[4]-pad*M_regular.h[2]];
	MshNew = getRegularMesh(OmegaNew,nnew);
elseif length(size(m))==3
	mnew = m[pad+1:end-pad,pad+1:end-pad,1:end-pad];
	nnew = collect(size(mnew)).-1;
	OmegaNew = [Omega[1],Omega[2] - 2*pad*M_regular.h[1],Omega[3],Omega[4]-2*pad*M_regular.h[2],Omega[4],Omega[5]-pad*M_regular.h[3]];
	MshNew = getRegularMesh(OmegaNew,nnew);
end
return mnew,MshNew
end

max_stars_repo_name: JuliaInv/jInvVisPyPlot

max_issues_repo_name: JuliaInv/jInvVisPyPlot

max_forks_repo_name: JuliaInv/jInvVisPyPlot


text:
# Producing the Julia Microbenchmarks plot

# Load the required Julia packages
using Base.MathConstants
using CSV
using DataFrames
using Gadfly
using StatsBase

# Load benchmark data from file
benchmarks =
    CSV.read("benchmarks.csv", DataFrame; header = ["language", "benchmark", "time"])

# Capitalize and decorate language names from datafile
dict = Dict(
    "c" => "C",
    "julia" => "Julia",
    "lua" => "LuaJIT",
    "fortran" => "Fortran",
    "java" => "Java",
    "javascript" => "JavaScript",
    "matlab" => "Matlab",
    "mathematica" => "Mathematica",
    "python" => "Python",
    "octave" => "Octave",
    "r" => "R",
    "rust" => "Rust",
    "go" => "Go",
);
benchmarks[!, :language] = [dict[lang] for lang in benchmarks[!, :language]]

# Normalize benchmark times by C times
ctime = benchmarks[benchmarks[!, :language] .== "C", :]
benchmarks = innerjoin(benchmarks, ctime, on = :benchmark, makeunique = true)
select!(benchmarks, Not(:language_1))
rename!(benchmarks, :time_1 => :ctime)
benchmarks[!, :normtime] = benchmarks[!, :time] ./ benchmarks[!, :ctime];

# Compute the geometric mean for each language
langs = [];
means = [];
priorities = [];
for lang in benchmarks[!, :language]
    data = benchmarks[benchmarks[!, :language] .== lang, :]
    gmean = geomean(data[!, :normtime])
    push!(langs, lang)
    push!(means, gmean)
    if (lang == "C")
        push!(priorities, 1)
    elseif (lang == "Julia")
        push!(priorities, 2)
    else
        push!(priorities, 3)
    end
end

# Add the geometric means back into the benchmarks dataframe
langmean = DataFrame(language = langs, geomean = means, priority = priorities)
benchmarks = innerjoin(benchmarks, langmean, on = :language)

# Put C first, Julia second, and sort the rest by geometric mean
sort!(benchmarks, [:priority, :geomean]);
sort!(langmean, [:priority, :geomean]);

p = plot(
    benchmarks,
    x = :language,
    y = :normtime,
    color = :benchmark,
    Scale.y_log10,
    Guide.ylabel(nothing),
    Guide.xlabel(nothing),
    Theme(
        guide_title_position = :left,
        colorkey_swatch_shape = :circle,
        minor_label_font = "Georgia",
        major_label_font = "Georgia",
    ),
)

draw(SVG("benchmarks.svg", 9inch, 9inch / golden), p)


text:
abstract type AbstractCoronaModel{T} end

sample_position(m::AbstractMetricParams{T}, model::AbstractCoronaModel{T}, N) where {T} =
    error("Not implemented for $(typeof(model)).")
sample_velocity(
    m::AbstractMetricParams{T},
    model::AbstractCoronaModel{T},
    sampler::AbstractDirectionSampler,
    us,
    N,
) where {T} = error("Not implemented for $(typeof(model)).")

@with_kw struct LampPostModel{T} <: AbstractCoronaModel{T}
    @deftype T
    h = 5.0
    θ = 0.01
    ϕ = 0.0
end

function sample_position(::AbstractMetricParams{T}, model::LampPostModel{T}, N) where {T}
    u = @SVector [T(0.0), model.h, model.θ, model.ϕ]
    fill(u, N)
end

function sample_velocity(
    m::AbstractMetricParams{T},
    ::AbstractCoronaModel{T},
    sampler::AbstractDirectionSampler,
    us,
    N,
) where {T}
    map(1:N) do i
        θ, ϕ = sample_angles(sampler, i, N)
        r, t, p = vector_to_local_sky(m, us[i], θ, ϕ)
        @SVector [T(0.0), r, t, p]
    end
end

export LampPostModel


text:
struct DiscretePOMDP
    # TODO: Use sparse matrices?
    T::Array{Float64, 3} # T(s,a,s′)
    R::Array{Float64, 2} # R(s,a) = ∑_s' R(s,a,s')*T(s,a,s′)
    O::Array{Float64, 3} # O(a,s′,o)
    γ::Float64
end

n_states(pomdp::DiscretePOMDP) = size(pomdp.T, 1)
n_actions(pomdp::DiscretePOMDP) = size(pomdp.T, 2)
n_observations(pomdp::DiscretePOMDP) = size(pomdp.O, 3)
discount(pomdp::DiscretePOMDP) = pomdp.γ

ordered_states(pomdp::DiscretePOMDP) = collect(1:n_states(pomdp))
ordered_actions(pomdp::DiscretePOMDP) = collect(1:n_actions(pomdp))
ordered_observations(pomdp::DiscretePOMDP) = collect(1:n_observations(pomdp))

transition(pomdp::DiscretePOMDP, s::Int, a::Int) = Categorical(pomdp.T[s,a,:])
observation(pomdp::DiscretePOMDP, a::Int, s′::Int) = Categorical(pomdp.O[a,s′,:])
reward(pomdp::DiscretePOMDP, s::Int, a::Int) = pomdp.R[s,a]

reward(pomdp::DiscretePOMDP, b::Vector{Float64}, a::Int) = sum(reward(pomdp,s,a)*b[s] for s in ordered_states(pomdp))

function POMDP(pomdp::DiscretePOMDP; γ::Float64=discount(pomdp))
    return POMDP(
        γ,
        ordered_states(pomdp),
        ordered_actions(pomdp),
        ordered_observations(pomdp),
        (s,a, s′=nothing) -> begin
            S′ = transition(pomdp, s, a)
            if s′ == nothing
                return S′
            end
            return pdf(S′, s′)
        end,
        (s,a) -> reward(pomdp, s, a),
        (a,s′,o=nothing) -> begin
            O′ = observation(pomdp, a, s′)
            if o == nothing
                return O′
            end
            return pdf(O′, o)
        end,
        (s, a)->begin
            s′ = rand(transition(pomdp,s,a))
            r = reward(pomdp, s, a)
            o = rand(observation(pomdp,a,s′))
            return (s′, r, o)
        end
    )
end