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

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

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

const SUPPORTED_OBJECTIVES = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64},
    MOI.ScalarQuadraticFunction{Float64}
}

const SUPPORTED_SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64}
}

const SUPPORTED_SCALAR_FUNCTIONS = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64}
}

const SUPPORTED_VECTOR_FUNCTIONS = Union{
    MOI.VectorOfVariables,
    MOI.VectorAffineFunction{Float64},
}

const SUPPORTED_VECTOR_SETS = Union{
    MOI.Zeros,
    MOI.Nonpositives,
    MOI.Nonnegatives,
    MOI.SecondOrderCone,
    MOI.PositiveSemidefiniteConeTriangle,
}

"""
    diff_optimizer(optimizer_constructor)::Optimizer

Creates a `DiffOpt.Optimizer`, which is an MOI layer with an internal optimizer
and other utility methods. Results (primal, dual and slack values) are obtained
by querying the internal optimizer instantiated using the
`optimizer_constructor`. These values are required for find jacobians with respect to problem data.

One define a differentiable model by using any solver of choice. Example:

```julia
julia> using DiffOpt, GLPK

julia> model = diff_optimizer(GLPK.Optimizer)
julia> model.add_variable(x)
julia> model.add_constraint(...)

julia> backward_quad(model)  # for convex quadratic models

julia> backward_quad(model)  # for convex conic models
```
"""
function diff_optimizer(optimizer_constructor)::Optimizer
    return Optimizer(MOI.instantiate(optimizer_constructor, with_bridge_type=Float64))
end

Base.@kwdef struct QPCache
    problem_data::Tuple{
        Matrix{Float64}, Vector{Float64}, # Q, q
        Matrix{Float64}, Vector{Float64}, # G, h
        Matrix{Float64}, Vector{Float64}, # A, b
        Int, Vector{VI}, # nz, var_list
        Int, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}}, # nineq_le, le_con_idx
        Int, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}}, # nineq_ge, ge_con_idx
        Int, Vector{MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}}, # nineq_sv_le, le_con_sv_idx
        Int, Vector{MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}}, # nineq_sv_ge, ge_con_sv_idx
        Int, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}}, # neq, eq_con_idx
        Int, Vector{MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}}, # neq_sv, eq_con_sv_idx
    }
    inequality_duals::Vector{Float64}
    equality_duals::Vector{Float64}
    lhs::SparseMatrixCSC{Float64, Int}
end

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

const CACHE_TYPE = Union{
    Nothing,
    QPCache,
    ConicCache,
}

mutable struct Optimizer{OT <: MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT
    primal_optimal::Vector{Float64}  # solution
    var_idx::Vector{VI}
    gradient_cache::CACHE_TYPE

    function Optimizer(optimizer_constructor::OT) where {OT <: MOI.ModelLike}
        new{OT}(
            optimizer_constructor,
            zeros(0),
            Vector{VI}(),
            nothing,
        )
    end
end

function MOI.add_variable(model::Optimizer)
    model.gradient_cache = nothing
    vi = MOI.add_variable(model.optimizer)
    push!(model.var_idx, vi)
    return vi
end

function MOI.add_variables(model::Optimizer, N::Int)
    model.gradient_cache = nothing
    return VI[MOI.add_variable(model) for i in 1:N]
end

function MOI.add_constraint(model::Optimizer, f::SUPPORTED_SCALAR_FUNCTIONS, s::SUPPORTED_SCALAR_SETS)
    model.gradient_cache = nothing
    return MOI.add_constraint(model.optimizer, f, s)
end

function MOI.add_constraint(model::Optimizer, vf::SUPPORTED_VECTOR_FUNCTIONS, s::SUPPORTED_VECTOR_SETS)
    model.gradient_cache = nothing
    return MOI.add_constraint(model.optimizer, vf, s)
end

function MOI.add_constraints(model::Optimizer, f::AbstractVector{F}, s::AbstractVector{S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    return CI{F, S}[MOI.add_constraint(model, f[i], s[i]) for i in eachindex(f)]
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES}, f::SUPPORTED_OBJECTIVES)
    model.gradient_cache = nothing
    MOI.set(model.optimizer, attr, f)
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, sense)
end

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

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

function MOI.optimize!(model::Optimizer)
    model.gradient_cache = nothing
    MOI.optimize!(model.optimizer)

    # do not fail. interferes with MOI.Tests.linear12test
    if !in(MOI.get(model.optimizer, MOI.TerminationStatus()),  (MOI.LOCALLY_SOLVED, MOI.OPTIMAL))
        @warn "problem status: $(MOI.get(model.optimizer, MOI.TerminationStatus()))"
        return
    end
    # save the solution
    model.primal_optimal = MOI.get(model.optimizer, MOI.VariablePrimal(), model.var_idx)
    return
end

const QP_SET_TYPES = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
}

const QP_FUNCTION_TYPES = Union{
    MOI.SingleVariable,
    MOI.ScalarAffineFunction{Float64},
}

const QP_OBJECTIVE_TYPES = Union{
    MOI.ScalarAffineFunction{Float64},
    MOI.ScalarQuadraticFunction{Float64},
    MOI.SingleVariable,
}

"""
    backward(model::Optimizer, params...)

Wrapper method for the backward pass. Checks the model/problem type and dispatches a call to
[`backward_quad`](@ref) or [`backward_conic`](@ref) accordingly.
"""
function backward(model::Optimizer, params...)
    # check if the model is a QP
    isQP = true

    # check constraints
    con_types = MOI.get(model.optimizer, MOI.ListOfConstraints())
    for (func, set) in con_types

        if !(func <: QP_FUNCTION_TYPES) || !(set <: QP_SET_TYPES)
            isQP = false
        end
    end

    # check objective
    if !(MOI.get(model.optimizer, MOI.ObjectiveFunctionType()) <: QP_OBJECTIVE_TYPES)
        isQP = false
    end

    if isQP
        return backward_quad(model, params...)
    else
        return backward_conic(model, params...)
    end
end

"""
    backward_quad(model::Optimizer, params::Array{String}, dl_dz::Array{Float64})

Method to differentiate optimal solution `z` and return
product of jacobian matrices (`dz / dQ`, `dz / dq`, etc) with
the backward pass vector `dl / dz`

The method computes the product of
1. jacobian of problem solution `z*` with respect to
    problem parameters `params` recieved as method arguments
2. a backward pass vector `dl / dz`, where `l` can be a loss function

Note that this method *does not returns* the actual jacobians.

For more info refer eqn(7) and eqn(8) of https://arxiv.org/pdf/1703.00443.pdf
"""
function backward_quad(model::Optimizer, params::Vector{String}, dl_dz::Vector{Float64})
    z = model.primal_optimal
    if model.gradient_cache !== nothing
        (
            Q, q, G, h, A, b, nz, var_list,
            nineq_le, le_con_idx,
            nineq_ge, ge_con_idx,
            nineq_sv_le, le_con_sv_idx,
            nineq_sv_ge, ge_con_sv_idx,
            neq, eq_con_idx,
            neq_sv, eq_con_sv_idx,
        ) = model.gradient_cache.problem_data
        λ = model.gradient_cache.inequality_duals
        ν = model.gradient_cache.equality_duals
        LHS = model.gradient_cache.lhs

    else # full computation
        problem_data = get_problem_data(model.optimizer)
        (
            Q, q, G, h, A, b, nz, var_list,
            nineq_le, le_con_idx,
            nineq_ge, ge_con_idx,
            nineq_sv_le, le_con_sv_idx,
            nineq_sv_ge, ge_con_sv_idx,
            neq, eq_con_idx,
            neq_sv, eq_con_sv_idx,
        ) = problem_data

        # separate λ, ν

        λ = MOI.get.(model.optimizer, MOI.ConstraintDual(), le_con_idx)
        append!(
            λ,
            MOI.get.(model.optimizer, MOI.ConstraintDual(), ge_con_idx),
        )
        append!(
            λ,
            MOI.get.(model.optimizer, MOI.ConstraintDual(), le_con_sv_idx),
        )
        append!(
            λ,
            MOI.get.(model.optimizer, MOI.ConstraintDual(), ge_con_sv_idx),
        )
        ν = MOI.get.(model.optimizer, MOI.ConstraintDual(), eq_con_idx)
        append!(
            ν,
            MOI.get.(model.optimizer, MOI.ConstraintDual(), eq_con_sv_idx)
        )

        LHS = create_LHS_matrix(z, λ, Q, G, h, A)
        model.gradient_cache = QPCache(
            problem_data,
            λ,
            ν,
            LHS,
        )
    end

    nineq_total = nineq_le + nineq_ge + nineq_sv_le + nineq_sv_ge
    RHS = [dl_dz; zeros(neq + neq_sv + nineq_total)]
    partial_grads = -LHS \ RHS
    dz = partial_grads[1:nz]

    if nineq_total > 0
        dλ = partial_grads[nz+1:nz+nineq_total]
    end
    @debug begin
        @show neq+neq_sv
    end
    @assert length(eq_con_idx) + length(eq_con_sv_idx) == neq+neq_sv
    if neq+neq_sv > 0
        dν = partial_grads[nz+nineq_total+1:nz+nineq_total+neq+neq_sv]
    end

    grads = []
    for param in params
        if param == "Q"
            push!(grads, 0.5 * (dz * z' + z * dz'))
        elseif param == "q"
            push!(grads, dz)
        elseif param == "G"
            push!(grads, Diagonal(λ) * dλ * z' - λ * dz')
        elseif param == "h"
            push!(grads, -Diagonal(λ) * dλ)
        elseif param == "A"
            push!(grads, dν * z' - ν * dz')
        elseif param == "b"
            push!(grads, -dν)
        else
            push!(grads, [])
        end
    end
    return grads
end

# `MOI.supports` methods

function MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}})
    return true
end

function MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}})
    return true
end

function MOI.supports(::Optimizer, ::MOI.AbstractModelAttribute)
    return true
end

function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction)
    return false
end

function MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<: SUPPORTED_OBJECTIVES})
    return true
end

MOI.supports_constraint(::Optimizer, ::Type{<: SUPPORTED_SCALAR_FUNCTIONS}, ::Type{<: SUPPORTED_SCALAR_SETS}) = true
MOI.supports_constraint(::Optimizer, ::Type{<: SUPPORTED_VECTOR_FUNCTIONS}, ::Type{<: SUPPORTED_VECTOR_SETS}) = true
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS} = true
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_VECTOR_FUNCTIONS, S<:SUPPORTED_VECTOR_SETS} = true

MOI.get(model::Optimizer, attr::MOI.SolveTime) = MOI.get(model.optimizer, attr)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.optimizer)
    empty!(model.primal_optimal)
    empty!(model.var_idx)
    model.gradient_cache = nothing
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           isempty(model.primal_optimal) &&
           model.gradient_cache === nothing
           isempty(model.var_idx)
end

# now supports name too
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = true #!copy_names

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    model.gradient_cache = nothing
    return MOIU.default_copy_to(model.optimizer, src, copy_names)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return MOI.get(model.optimizer, MOI.TerminationStatus())
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::VI, value::Float64)
    model.gradient_cache = nothing
    MOI.set(model.optimizer, MOI.VariablePrimalStart(), vi, value)
end

function MOI.supports(model::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return MOI.supports(model.optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(model.optimizer, attr)
    return MOI.get(model.optimizer, attr, vi)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.is_valid(model::Optimizer, v::VI)
    return v in model.var_idx
end

MOI.is_valid(model::Optimizer, con::CI) = MOI.is_valid(model.optimizer, con)

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F <: SUPPORTED_VECTOR_FUNCTIONS, S <: SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintBasisStatus, ci::CI{F,S}) where {F <: SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintBasisStatus(), ci)
end

# helper methods to check if a constraint contains a Variable
function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.SingleVariable, S}) where {S <: SUPPORTED_SCALAR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v == func.variable
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.ScalarAffineFunction{Float64}, S}) where {S <: SUPPORTED_SCALAR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.variable_index, func.terms)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorOfVariables, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v in func.variables
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorAffineFunction{Float64}, S}) where {S <: SUPPORTED_VECTOR_SETS}
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.scalar_term.variable_index, func.terms)
end

function MOI.delete(model::Optimizer, v::VI)
    model.gradient_cache = nothing
    # delete in inner solver
    MOI.delete(model.optimizer, v)

    # delete from var_idx
    filter!(x -> x ≠ v, model.var_idx)
end

# for array deletion
function MOI.delete(model::Optimizer, indices::Vector{VI})
    model.gradient_cache = nothing
    for i in indices
        MOI.delete(model, i)
    end
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.AbstractFunctionModification
)
    model.gradient_cache = nothing
    MOI.modify(
        model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        chg
    )
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S <: SUPPORTED_SCALAR_SETS}
    model.gradient_cache = nothing
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.modify(model::Optimizer, ci::CI{F, S}, chg::MOI.AbstractFunctionModification) where {F<:MOI.VectorAffineFunction{Float64}, S <: SUPPORTED_VECTOR_SETS}
    model.gradient_cache = nothing
    MOI.modify(model.optimizer, ci, chg)
end

"""
    π(v::Vector{Float64}, model::MOI.ModelLike, conic_form::MatOI.GeometricConicForm, index_map::MOIU.IndexMap)

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

"""
    Dπ(v::Vector{Float64}, model, conic_form::MatOI.GeometricConicForm, index_map::MOIU.IndexMap)

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

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

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

# Map the rows corresponding to `F`-in-`S` constraints and store it in `x`.
function _map_rows!(f::Function, x::Vector, model, conic_form::MatOI.GeometricConicForm, index_map::MOIU.DoubleDicts.IndexWithType{F, S}, map_mode, k) where {F, S}
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        r = MatOI.rows(conic_form, index_map[ci])
        k += 1
        _assign_mapped!(x, f(ci, r), r, k, map_mode)
    end
    return k
end

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

"""
    map_rows(f::Function, model, conic_form::MatOI.GeometricConicForm, index_map::MOIU.IndexMap, map_mode::Union{Nested{T}, Flattened{T}})

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

"""
    backward_conic(model::Optimizer, dA::Matrix{Float64}, db::Vector{Float64}, dc::Vector{Float64})

Method to differentiate optimal solution `x`, `y`, `s` given perturbations related to
conic program parameters `A`, `b`, `c`. This is similar to [`backward_quad`](@ref) method
but it this *does returns* the actual jacobians.

For theoretical background, refer Section 3 of Differentiating Through a Cone Program, https://arxiv.org/abs/1904.09043
"""
function backward_conic(model::Optimizer, dA::AbstractMatrix{<:Real}, db::AbstractVector{<:Real}, dc::AbstractVector{<:Real})
    if !in(
            MOI.get(model, MOI.TerminationStatus()), (MOI.LOCALLY_SOLVED, MOI.OPTIMAL)
        )
        error("problem status: ", MOI.get(model.optimizer, MOI.TerminationStatus()))
    end

    if model.gradient_cache !== nothing
        M = model.gradient_cache.M
        vp = model.gradient_cache.vp
        Dπv = model.gradient_cache.Dπv
        (x, y, s) = model.gradient_cache.xys
        A = model.gradient_cache.A
        b = model.gradient_cache.b
        c = model.gradient_cache.c

        m = A.m
        n = A.n
        N = m + n + 1
        # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
        (u, v, w) = (x, y - s, 1.0)
    else

        # fetch matrices from MatrixOptInterface
        cone_types = unique!([S for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraints())])
        # We use `MatOI.OneBasedIndexing` as it is the same indexing as used by `SparseMatrixCSC`
        # so we can do an allocation-free conversion to `SparseMatrixCSC`.
        conic_form = MatOI.GeometricConicForm{Float64, MatOI.SparseMatrixCSRtoCSC{Float64, Int, MatOI.OneBasedIndexing}, Vector{Float64}}(cone_types)
        index_map = MOI.copy_to(conic_form, model)

        # fix optimization sense
        if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
            conic_form.sense = MOI.MIN_SENSE
            conic_form.c = -conic_form.c
        end

        A = convert(SparseMatrixCSC{Float64, Int}, conic_form.A)
        b = conic_form.b
        c = conic_form.c

        # programs in tests were cross-checked against `diffcp`, which follows SCS format
        # hence, some arrays saved during `MOI.optimize!` are not same across all optimizers
        # specifically there's an extra preprocessing step for `PositiveSemidefiniteConeTriangle` constraint for SCS/Mosek

        # get x,y,s
        x = model.primal_optimal
        s = map_rows((ci, r) -> MOI.get(model, MOI.ConstraintPrimal(), ci), model.optimizer, conic_form, index_map, Flattened{Float64}())
        y = map_rows((ci, r) -> MOI.get(model, MOI.ConstraintDual(), ci), model.optimizer, conic_form, index_map, Flattened{Float64}())

        # pre-compute quantities for the derivative
        m = A.m
        n = A.n
        N = m + n + 1
        # NOTE: w = 1 systematically since we asserted the primal-dual pair is optimal
        (u, v, w) = (x, y - s, 1.0)

        # find gradient of projections on dual of the cones
        Dπv = Dπ(v, model.optimizer, conic_form, index_map)

        # Q = [
        #      0   A'   c;
        #     -A   0    b;
        #     -c' -b'   0;
        # ]
        # M = (Q- I) * B + I
        # with B =
        # [
        #  I    .   .
        #  .  Dπv   .
        # .    .   I
        # ]

        M = [
            spzeros(n,n)  (A' * Dπv)    c
            -A               -Dπv + I    b
            -c'            -b' * Dπv    0
        ]
        # find projections on dual of the cones
        vp = π(v, model.optimizer, conic_form, index_map)

        model.gradient_cache = ConicCache(
            M = M,
            vp = vp,
            Dπv = Dπv,
            xys = (x, y, s),
            A = A,
            b = b,
            c = c,
        )
    end # cache

    # dQ * [u, vp, max(0, w)]
    RHS = [dA' * vp + dc; -dA * u + db; -dc ⋅ u - db ⋅ vp]

    dz = if norm(RHS) <= 1e-4
        RHS .= 0
    else
        lsqr(M, RHS)
    end

    du, dv, dw = dz[1:n], dz[n+1:n+m], dz[n+m+1]
    dx = du - x * dw
    dy = Dπv * dv - y * dw
    ds = Dπv * dv - dv - s * dw
    return -dx, -dy, -ds
end

MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{MOI.VariableIndex}) = true

function MOI.get(model::Optimizer, ::MOI.VariableName, v::VI)
    return MOI.get(model.optimizer, MOI.VariableName(), v)
end

function MOI.set(model::Optimizer, ::MOI.VariableName, v::VI, name::String)
    MOI.set(model.optimizer, MOI.VariableName(), v, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    return MOI.get(model.optimizer, MOI.VariableIndex, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, con::CI, name::String)
    MOI.set(model.optimizer, MOI.ConstraintName(), con, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    return MOI.get(model.optimizer, MOI.ConstraintIndex, name)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, con::CI)
    return MOI.get(model.optimizer, MOI.ConstraintName(), con)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{F,S})
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{VF, VS}}) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{VF,VS})
end

function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    MOI.set(model.optimizer, MOI.Name(), name)
end

function MOI.get(model::Optimizer, ::Type{CI{F, S}}, name::String) where {F<:SUPPORTED_SCALAR_FUNCTIONS, S<:SUPPORTED_SCALAR_SETS}
    return MOI.get(model.optimizer, CI{F,S}, name)
end

function MOI.get(model::Optimizer, ::Type{CI{VF, VS}}, name::String) where {VF<:SUPPORTED_VECTOR_FUNCTIONS, VS<:SUPPORTED_VECTOR_SETS}
    return MOI.get(model.optimizer, CI{VF,VS}, name)
end

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Real, Nothing})
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
end
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

MOI.supports(model::Optimizer, ::MOI.Silent) = MOI.supports(model.optimizer, MOI.Silent())

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    MOI.set(model.optimizer, MOI.Silent(), value)
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return MOI.get(model.optimizer, MOI.Silent())
end

function MOI.get(model::Optimizer, ::MOI.SolverName)
    return MOI.get(model.optimizer, MOI.SolverName())
end


text:
module WebDAV

import HTTP
using EzXML
using Base.Filesystem: readdir, mkdir, isdir, isfile
using Base: open
using Base64
import Base: download

namespace = Dict("d" => "DAV:")

struct Server
    url
    headers
end

escape_no_slash(url) = join(HTTP.escapeuri.(split(url,'/')),'/')

"""
    server = WebDAV.Server(url,username,password)

`url` is for example https://server/remote.php/webdav for an NextCloud instance
"""
function Server(url::AbstractString,username::AbstractString,password::AbstractString)
    userinfo = username * ":" * password
    headers = [("Authorization", "Basic $(base64encode(userinfo))")]
    return Server(url,headers)
end

function upload(s::Server,stream::IOStream,remotepath::AbstractString)
    remotepath_escaped = escape_no_slash(remotepath)
    r = HTTP.request("PUT", s.url * "/" * remotepath_escaped,s.headers, stream);
end

function upload(s::Server,localpath::AbstractString,remotepath::AbstractString)
    open(localpath) do stream
        upload(s,stream,remotepath)
    end
end

function download(s::Server,remotepath::AbstractString,localpath::AbstractString)
    open(localpath,"w") do file
        open(s,remotepath,"r") do stream
            while !eof(stream)
                write(file,readavailable(stream))
            end
        end
    end
end

function download(s::Server,remotepath::AbstractString)
    localpath = tempname()
    download(s,remotepath,localpath)
    return localpath
end

# Filesystem-like API

function Base.open(s::Server,remotepath::AbstractString,
                   mode::AbstractString = "r")

    remotepath_escaped = escape_no_slash(remotepath)
    if mode == "r"
        io = Base.BufferStream()
        r = HTTP.request("GET", s.url * "/" * remotepath_escaped,s.headers,response_stream = io)
        return io
    else
        error("unsupported mode $(mode)")
    end
end

function Base.open(f::Function,s::Server,remotepath::AbstractString,
                   mode::AbstractString = "r")

    remotepath_escaped = s.url * "/" * escape_no_slash(remotepath)
    @debug "remotepath_escaped: $remotepath_escaped"
    if mode == "r"
        r = HTTP.open("GET", remotepath_escaped,s.headers) do io
            f(io)
        end
    elseif mode == "w"
        r = HTTP.open("PUT", remotepath_escaped,s.headers) do io
            f(io)
        end
    else
        error("unsupported mode $(mode)")
    end

    return nothing
end

function properties(s::Server,dir)
    r = HTTP.request("PROPFIND", s.url * "/" * escape_no_slash(dir), s.headers; status_exception = false);

    body = String(r.body)
    @debug "properties status: $(r.status)"
    @debug "properties body: $body"
    if r.status == 404
        return nothing,r.status
    else
        doc = EzXML.parsexml(body)
        return doc,r.status
    end
end

function Base.Filesystem.readdir(s::Server,dir::AbstractString=".")
    doc,status = properties(s,dir)
    found = true

    if status == 404
        found = false
    elseif status == 207
        # 207 Multi-Status (WebDAV; RFC 4918)
        status = nodecontent.(findall("/d:multistatus/d:response/d:propstat/d:status",root(doc),namespace))

        if any(contains("404"),status)
            found = false
        end
    end

    if !found
        error("directory $(dir) not found on server $(s.url)")
    end

    path = HTTP.URI(s.url).path
    response = findall("d:response",root(doc),namespace)
    list = Vector{String}(undef,length(response))
    for i = 1:length(response)
        url = nodecontent(findfirst("d:href",response[i],namespace))

        if startswith(url,path)
            list[i] = url[length(path)+1:end]
        else
            list[i] = url
        end
    end

    return HTTP.unescapeuri.(list)
end

function Base.Filesystem.mkdir(s::Server,dir::AbstractString)
    escaped_dir = escape_no_slash(dir)
    if !endswith(dir,"/")
        escaped_dir = escaped_dir * "/"
    end

    r = HTTP.request("MKCOL", s.url * "/" * escaped_dir, s.headers);
    return nothing
end

"""
    rm(server::WebDAV.Server,path)

Removes the `path` on the WebDAV `server`.
"""
function Base.Filesystem.rm(s::Server, dir::AbstractString)
    escaped_dir = escape_no_slash(dir)
    if isdir(s,dir) && !endswith(dir,"/")
        escaped_dir = escaped_dir * "/"
    end

    r = HTTP.request("DELETE", s.url * "/" * escaped_dir, s.headers);
    return nothing
end

function Base.Filesystem.isdir(s::Server, dir::AbstractString)
    doc,status = properties(s,dir)

    # not found
    if status == 404
        return false
    end

    resourcetype = findall("d:response[1]/d:propstat/d:prop/d:resourcetype/*",root(doc),namespace)

    return "collection" in nodename.(resourcetype)
end

function Base.Filesystem.isfile(s::Server, dir::AbstractString)
    doc,status = properties(s,dir)

    # not found
    if status == 404
        return false
    elseif status == 207
        # 207 Multi-Status (WebDAV; RFC 4918)
        status = nodecontent.(findall("/d:multistatus/d:response/d:propstat/d:status",root(doc),namespace))

        if any(contains("404"),status)
            return false
        end
    end

    resourcetype = findall("d:response[1]/d:propstat/d:prop/d:resourcetype/*",root(doc),namespace)

    return !("collection" in nodename.(resourcetype))
end

export download, upload
end # module

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the moon position.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorne, CA.
#
#   [2] Meeus, J (1998). Astronomical algorithms. Willmann-Bell, Inc, Richmond,
#       VA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export moon_position_i

"""
    moon_position_i(JD_TDB::Number[, model])

Compute the Moon position represented in the IAU-76/FK5 MOD (mean-equator,
mean-equinox of date) at the Julian Day `JD_TDB` (Barycentric Dynamical Time).

The `model` must be `Val(:Meeus)` or `Val(:Vallado)`. `Val(:Meeus)` uses the
algorithm in [2, p. 337] that provides an accuracy of 10" in the longitude and
4" in the latitude (the reference does not mention the timespan). `Val(:fast)`
uses the algorithm in [1, p. 288] that is 10x faster than `Val(:Meeus)` but can
lead to errors of 0.3° in longitude and 0.2° in latitude.

"""
moon_position_i(JD_TDB::Number) = moon_position_i(JD_TDB, Val(:Meeus))

function moon_position_i(JD_TDB::Number, ::Val{:Meeus})
    # Number of Julian centuries from J2000 epoch.
    T_TDB = (JD_TDB - JD_J2000)/36525.0

    # Moon's mean latitude reffered to the mean equinox of data [deg].
    L´ = @evalpoly(T_TDB, 218.3164477, 481_267.88123421, -0.0015786, 1/538841, -1/65_194_000)

    # Mean elongation of the Moon [deg].
    D = @evalpoly(T_TDB, 297.8501921, 445_267.1114034, -0.0018819, 1/545868, -1/113_065_000)

    # Sun's mean anomaly [deg].
    M = @evalpoly(T_TDB, 357.5291092, 35_999.0502909, -0.0001536, 1/24_490_000)

    # Moon's mean anomaly [deg].
    M´ = @evalpoly(T_TDB, 134.9633964, 477_198.8675055, 0.0087414, 1/69_699, -1/14_712_000)

    # Moon's argument of latitude (mean distance of the Moon from its ascending
    # node) [deg].
    F = @evalpoly(T_TDB, 93.2720950, 483_202.0175233, -0.0036539, -1/3_526_000, 1/863_310_000)

    # Obliquity of the ecliptic [deg].
    ϵ = @evalpoly(T_TDB, 23.439_291, -0.013_004_2, -1.64e-7, +5.04e-7)

    # Additional arguments required for the algorithm [deg].
    A₁ = @evalpoly(T_TDB, 119.75, 131.849)
    A₂ = @evalpoly(T_TDB,  53.09, 479_264.290)
    A₃ = @evalpoly(T_TDB, 313.45, 481_266.484)

    # Convert everything to [rad] and limit the angles between [0, 2π].
    L´ = mod2pi(deg2rad(L´))
    D  = mod2pi(deg2rad(D))
    M  = mod2pi(deg2rad(M))
    M´ = mod2pi(deg2rad(M´))
    F  = mod2pi(deg2rad(F))
    ϵ  = mod2pi(deg2rad(ϵ))
    A₁ = mod2pi(deg2rad(A₁))
    A₂ = mod2pi(deg2rad(A₂))
    A₃ = mod2pi(deg2rad(A₃))

    # Term used to corret the arguments of the angle M that depends on the
    # Earth's orbit eccentricity around the Sun.
    E  = @evalpoly(T_TDB, 1, -0.002516, 0.0000074)
    E² = E * E

    # Compute the sum of the terms in the tables 47.A and 47.B [2].

    # Auxiliary variables to simplify the expressions.
    tab = _tab_47a

    # Longitude and distance of the Moon.
    num_terms = size(tab)[2]

    Σl = 0.0
    Σr = 0.0

    @inbounds for k = 1:num_terms
        aD  = tab[1, k]
        aM  = tab[2, k]
        aM´ = tab[3, k]
        aF  = tab[4, k]

        arg = mod2pi(aD * D + aM * M + aM´ * M´ + aF * F)

        # Check if we need to apply the correction `E`.
        E_corr = if ((aM == 1) || (aM == -1))
            E
        elseif ((aM == 2) || (aM == -2))
            E²
        else
            one(E)
        end

        # Compute the quantities.
        sin_arg, cos_arg = sincos(arg)

        Σl += tab[5, k] * E_corr * sin_arg
        Σr += tab[6, k] * E_corr * cos_arg
    end

    # Auxiliary variables to simplify the expressions.
    tab = _tab_47b

    # Latitude of the Moon.
    num_terms = size(tab)[2]

    Σb = 0.0

    @inbounds for k = 1:num_terms
        aD  = tab[1, k]
        aM  = tab[2, k]
        aM´ = tab[3, k]
        aF  = tab[4, k]

        arg = mod2pi(aD * D + aM * M + aM´ * M´ + aF * F)

        # Check if we need to apply the correction `E`.
        E_corr = if ((aM == 1) || (aM == -1))
            E
        elseif ((aM == 2) || (aM == -2))
            E²
        else
            one(E)
        end

        Σb += tab[5, k] * E_corr * sin(arg)
    end

    # Apply the corrections to the terms.
    Σl +=  3958sin(A₁) + 1962sin(L´ - F) + 318sin(A₂)
    Σb += -2235sin(L´) + 382sin(A₃) + 175sin(A₁ - F) + 175sin(A₁ + F) +
           127sin(L´ - M´) - 115sin(L´ + M´)

    # Convert to [rad].
    Σl = mod2pi(deg2rad(Σl / 1_000_000))
    Σb = mod2pi(deg2rad(Σb / 1_000_000))

    # Compute the Moon coordinates [rad] and [m].
    λ = mod2pi(L´ + Σl)
    β = Σb
    Δ = 385_000.56e3 + Σr

    # Compute the Moon vector in MOD [m].
    #
    # Notice that λ and β provide us the geocentric latitude and longitude of
    # the Moon w.r.t. the mean equator of date in the ecliptic plane. Hence, we
    # need also to rotate the mean ecliptic to obtain the vector in the MOD.
    sin_λ, cos_λ = sincos(λ)
    sin_β, cos_β = sincos(β)
    sin_ϵ, cos_ϵ = sincos(ϵ)

    r_moon_i = SVector(
        Δ * cos_β * cos_λ,
        Δ * (cos_ϵ * cos_β * sin_λ - sin_ϵ * sin_β),
        Δ * (sin_ϵ * cos_β * sin_λ + cos_ϵ * sin_β)
    )

    return r_moon_i
end

function moon_position_i(JD_TDB::Number, ::Val{:Vallado})
    # Number of Julian centuries from J2000 epoch.
    T_TDB = (JD_TDB - JD_J2000)/36525.0

    # Auxiliary computation to improve performance.
    sin1, cos1 = sincos(deg2rad(134.9 + 477_198.85T_TDB))
    sin2, cos2 = sincos(deg2rad(259.2 - 413_335.38T_TDB))
    sin3, cos3 = sincos(deg2rad(235.7 + 890_534.23T_TDB))
    sin4, cos4 = sincos(deg2rad(269.9 + 954_397.70T_TDB))
    sin5       = sin(deg2rad(357.5 +  35_999.05T_TDB))
    sin6       = sin(deg2rad(186.6 + 966_404.05T_TDB))

    # Ecliptic latitude of the Moon [deg].
    λₑ = 218.32 + 481_267.8813T_TDB +
         6.29sin1 - 1.27sin2 + 0.66sin3 + 0.21sin4 - 0.19sin5 - 0.11sin6

    # Ecliptic longitude of the Moon [deg].
    ϕₑ = 5.13sin(deg2rad( 93.3 + 483_202.03T_TDB)) +
         0.28sin(deg2rad(228.2 + 960_400.87T_TDB)) -
         0.28sin(deg2rad(318.3 +   6_003.18T_TDB)) -
         0.17sin(deg2rad(217.6 - 407_332.20T_TDB))

    # Parallax [deg].
    P = 0.9508 + 0.0518cos1 + 0.0095cos2 + 0.0078cos3 + 0.0028cos4

    # Obliquity of the ecliptic [deg].
    ϵ = @evalpoly(T_TDB, 23.439_291, -0.013_004_2, -1.64e-7, +5.04e-7)

    # Convert to radians and limit to the interval [0,2π].
    λₑ = mod2pi(deg2rad(λₑ))
    ϕₑ = mod2pi(deg2rad(ϕₑ))
    P  = mod2pi(deg2rad( P))
    ϵ  = mod2pi(deg2rad( ϵ))

    # Compute the distance from Earth to the Moon [m].
    r = R0 / sin(P)

    # Auxiliary variables.
    sin_λ, cos_λ = sincos(λₑ)
    sin_ϕ, cos_ϕ = sincos(ϕₑ)
    sin_ϵ, cos_ϵ = sincos(ϵ)

    # Compute the Moon vector represented in MOD (IAU-76/KF5 mean-equator,
    # mean-equinox of date).
    r_mod = SVector(
        r * (cos_ϕ * cos_λ),
        r * (cos_ϵ * cos_ϕ * sin_λ - sin_ϵ * sin_ϕ),
        r * (sin_ϵ * cos_ϕ * sin_λ + cos_ϵ * sin_ϕ)
    )

    return r_mod
end


text:
#===============================================================================

The radical of n, rad(n), is the product of distinct prime factors of n. For
example, 504 = 2^3 × 3^2 × 7, so rad(504) = 2 × 3 × 7 = 42.

We shall define the triplet of positive integers (a, b, c) to be an abc-hit if:

    gcd(a, b) = gcd(a, c) = gcd(b, c) = 1
    a < b
    a + b = c
    rad(abc) < c

For example, (5, 27, 32) is an abc-hit, because:

    gcd(5, 27) = gcd(5, 32) = gcd(27, 32) = 1
    5 < 27
    5 + 27 = 32
    rad(4320) = 30 < 32

It turns out that abc-hits are quite rare and there are only thirty-one abc-hits
for c < 1000, with ∑c = 12523.

Find ∑c for c < 120000.

--------------------------------------------------------------------------------

Analyse
-------

1) a < b; a + b = c < bound
    a)  ⇒ 2a < a + b = c < bound - 1 ⇒ 1 ≤ a < bound ÷ 2
    b)  ⇒ b < bound - a

2) if a and b are coprime, then c = a+b is coprime with a and with b
    => gcd(a, b) = 1 ⇒ gcd(a, c) = gcd(b, c) = 1

3) gcd(a, b) == 1 ⇔ gcd(rad(a), rad(b)) == 1

4) if a and b are coprime, the rad(ab) = rad(a) rad(b)

5) rad(n) = 1 if n = 1
          ≥ 2 if n > 1

6) So the conditions can be simplifed to
        gcd(a, b) = 1
        1 ≤ a < bound ÷ 2
        a < b < bound - a
        c = a + b
        rad(a) rad(b) rad(c) < c

7) Precomputing rad(n) for all n in 2..bound-2

case a=1
--------
    gcd(1,b) = 1 is granted
    2 ≤ b ≤ bound-2
    c = 1+b
    rad(b) rad(c) < c ⇒ rad(b) < (1+b) ÷ 2
                      ⇒ rad(b) rad(1+b) < 1+b

    Therefore the constraints are:
        1 < b < bound-1
        c = b+1
        rad(b) rad(b+1) < (b+1)

case a>1
--------
    a < b => a+1 ≤ b,
    2a+1 ≤ a + b = c < bound => a < (bound-1)÷2
                             => a ≤ (bound-1)÷2 -1
    rad(a) ≥ 2, rad(b) ≥ 2, rad(c) ≥ 2
    rad(a) rad(b) rad(a+b) < a+b ⇒ rad(a) < c ÷ 4
                                 ⇒ rad(b) < c ÷ (2 rad(a))

    Therefore the constraints are:
        2   ≤ a ≤ (bound-1)÷2 -1
        a+1 ≤ b ≤ bound - 1 - a
        gcd(a,b) = 1
        c = a + b
        rad(a) < c ÷ 4                  < (bound-1) ÷ 4
        rad(b) < c ÷ (2 rad(a))         < (bound-1) ÷ 4
        rad(c) < c ÷ (rad(a) rad(b))    < (bound-1) ÷ 4

Conclusion
==========

The principle behind the algorithm is to precompute rad(n) for 2 ≤ n ≤ bound-, and to proceed to all possible tests before computing gcd(rad(a), rad(b)).

Note that precomputing gcd(rad(a), rad(b)) for all possible a and b is not efficient.

===============================================================================#

import Primes

function euler(bound::Int)
    #precimpute the radicals
    rad = Array{Int}(bound-1)
    for n in 1:bound-1
        rad[n]=prod(keys(Primes.factor(n)))
    end

    res::Int = 0

    # case a=1
    for b in 2:bound-2
        c=b+1
        if (rad[b]*rad[c] < c) res+=c end
    end

    # case a>1
    for a in 2:(bound-1) ÷ 2 -1
        rad_a = rad[a]
        if (rad_a ≥ (bound-1) ÷ 4) continue end

        for b in a+1:bound-1-a
            rad_b = rad[b]
            c=a+b
            if (rad_a * rad_b * rad[c] ≥ c) continue end
            if (gcd(rad_a, rad_b) == 1) res += c end
        end
    end

    res
end

println("  1_000: $(euler(  1_000))")
@time println("120_000: $(euler(120_000))")

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
using Test
using GraphIO.EdgeList
using GraphIO.EdgeList: IntEdgeListFormat

@testset "EdgeList" begin
    for g in values(digraphs)
        readback_test(EdgeListFormat(), g)
        readback_test(IntEdgeListFormat(), g)
    end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
module O

using Plots

struct Wireframe{F, X, Y}
    f::F
    x::X
    y::Y
end

vecwithnan(X) = [
    vec([X ; fill(NaN, 1, size(X,  2))])
    vec([X'; fill(NaN, 1, size(X', 2))])
]

@recipe function f(wf::Wireframe)
    f, x, y = wf.f, wf.x, wf.y
    X, Y = reim(complex.(x', y))
    XX, YY = vecwithnan(X), vecwithnan(Y)
    ZZ = f.(XX, YY)
    XX, YY, ZZ
end

@recipe function f(wf::Wireframe{<:AbstractMatrix})
    Z, x, y = wf.f, wf.x, wf.y
    X, Y = reim(complex.(x', y))
    XX, YY, ZZ = vecwithnan(X), vecwithnan(Y), vecwithnan(Z)
    XX, YY, ZZ
end

end

# %%
using Plots
x = y = range(-1, 1, 51)
f(x, y) = x == y == 0 ? zero(x) : x^2*y/(x^2 + y^2)
plot(O.Wireframe(f, x, y); label="")

# %%
Z = f.(x', y)
plot(O.Wireframe(Z, x, y); label="")

# %%
wireframe(x, y, f)

# %%

max_stars_repo_path: 0009/my Wireframe recipe WIP.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0009/my Wireframe recipe WIP.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0009/my Wireframe recipe WIP.jl

max_forks_repo_name: genkuroki/public


text:
"""
The Yarn ClusterManager for Julia.

ugi: UserGroupInformation representing the user to represent (default is current logged in user)
yarnhost: yarn host name
rmport: yarn resource manager port
schedport: yarn scheduler port
launch_timeout: number of seconds to wait until all workers are launched (default is 60, provide larger value for more workers)
keep_connected: if false, YarnManager will disconnect from the cluster once all workers are removed (default true, use false if you do not plan to add more workers after removing all existing workers)
"""
struct YarnManager <: ClusterManager
    ugi::UserGroupInformation
    clnt::Union{YarnClient,Nothing}
    am::YarnAppMaster
    app::Union{YarnApp,ContainerIdProto}
    launch_timeout::Integer
    keep_connected::Bool

    function YarnManager(; kwargs...)
        params = Dict(kwargs)
        paramkeys = keys(params)
        unmanaged = (:unmanaged in paramkeys) ? params[:unmanaged] : true
        launch_timeout  = (:launch_timeout in paramkeys) ? params[:launch_timeout] : 60
        keep_connected  = (:keep_connected in paramkeys) ? params[:keep_connected] : true
        @debug("YarnManager constructor", params)

        if unmanaged
            ugi             = (:ugi             in paramkeys) ? params[:ugi]            : UserGroupInformation()
            rmport          = (:rmport          in paramkeys) ? params[:rmport]         : 8032
            yarnhost        = (:yarnhost        in paramkeys) ? params[:yarnhost]       : "localhost"
            schedport       = (:schedport       in paramkeys) ? params[:schedport]      : 8030

            clnt    = YarnClient(yarnhost, rmport, ugi)
            am      = YarnAppMaster(yarnhost, schedport, ugi)
            app     = submit(clnt, am)
        else
            ugi     = UserGroupInformation()
            am      = YarnAppMaster(ugi)
            clnt = nothing
            register(am)
            app = parse_container_id(ENV["CONTAINER_ID"])
        end

        new(ugi, clnt, am, app, launch_timeout, keep_connected)
    end
end

function YarnManager(fn::Function; kwargs...)
    yam = YarnManager(; kwargs...)
    try
        fn(yam)
    finally
        disconnect(yam)
    end
end

function disconnect(yarncm::YarnManager)
    ret = unregister(yarncm.am, true)
end

function show(io::IO, yarncm::YarnManager)
    print(io, "YarnManager for ")
    show(io, yarncm.am)
end

function setup_worker(host, port)
    @debug("YarnManager setup_worker", host, port)

    c = connect(IPv4(host), port)
    if :wait_connected in names(Base; all=true) # < Julia 1.3
        Base.wait_connected(c)
    else
        Sockets.wait_connected(c) # >= Julia 1.3
    end
    # identify container id so that rmprocs can clean things up nicely if required
    serialize(c, ENV["CONTAINER_ID"])

    redirect_stdout(c)
    redirect_stderr(c)

    start_worker(c, ENV["JULIA_YARN_KUKI"])
end

_envdict(envdict::Dict) = envdict
function _envdict(envhash::Base.EnvDict)
    envdict = Dict{String,String}()
    for (n,v) in envhash
        envdict[n] = v
    end
    envdict
end

function _currprocname()
    # Julia throws ENOBUFS sometimes while fetching process title
    pt = try
        Sys.get_process_title()
    catch
        "julia"
    end
    p = joinpath(Sys.BINDIR, pt)
    isfile(p) && (return p)

    haskey(ENV, "_") && contains(ENV["_"], "julia") && (return ENV["_"])

    "julia"
end

function container_start(manager::YarnManager, cmd::String, env::Dict{String,String}, ipaddr::IPv4, port::UInt16, cid::ContainerIdProto)
    try
        env["JULIA_YARN_KUKI"] = cluster_cookie()

        initargs = "using Elly; Elly.setup_worker($(ipaddr.host), $(port))"

        clc = launchcontext(cmd="$cmd -e '$initargs'", env=env)
        # TODO: do we need to add tokens into clc?

        @debug("YarnManager container_start", initargs, clc)
        container_start(manager.am, cid, clc)
    catch ex
        @error("error starting yarn container", exception=ex)
        rethrow(ex)
    end
end

container_stop(manager::YarnManager, cid::ContainerIdProto) = container_stop(manager.am, cid)
container_release(manager::YarnManager, cid::ContainerIdProto) = container_release(manager.am, cid)

function launch(manager::YarnManager, params::Dict, instances_arr::Array, c::Condition)
    @debug("YarnManager launch", params)

    paramkeys   = keys(params)
    np          = (:np       in paramkeys)  ? params[:np]               : 1
    cmd         = (:exename  in paramkeys)  ? params[:exename]          : _currprocname()
    appenv      = (:env      in paramkeys)  ? _envdict(params[:env])    : Dict{String,String}()
    mem         = (:mem      in paramkeys)  ? params[:mem]              : YARN_CONTAINER_MEM_DEFAULT
    cpu         = (:cpu      in paramkeys)  ? params[:cpu]              : YARN_CONTAINER_CPU_DEFAULT
    loc         = (:loc      in paramkeys)  ? params[:loc]              : YARN_CONTAINER_LOCATION_DEFAULT
    priority    = (:priority in paramkeys)  ? params[:priority]         : YARN_CONTAINER_PRIORITY_DEFAULT

    acceptors = Task[]
    ipaddr = getipaddr()
    (port, server) = listenany(ipaddr, 11000)
    accept_task = @async begin
        while length(acceptors) < np
            sock = accept(server)
            push!(acceptors, @async begin
                config = WorkerConfig()
                config.io = sock
                # keep the container id in userdata to be used with rmprocs if required
                cidstr = deserialize(sock)
                @debug("got container id string", cidstr)
                cid = parse_container_id(cidstr)
                @debug("got container id", cid)
                config.userdata = Dict(:container_id => cid)
                push!(instances_arr, config)
            end)
        end
    end

    on_alloc = (cid) -> container_start(manager, cmd, appenv, ipaddr, port, cid)
    callback(manager.am, on_alloc, nothing)

    container_allocate(manager.am, convert(Int, np); mem=mem, cpu=cpu, loc=loc, priority=priority)

    timedwait(()->istaskdone(accept_task), Float64(manager.launch_timeout))
    nconnected = 0
    nsuccess = 0
    while !isempty(acceptors)
        acceptor = pop!(acceptors)
        nconnected += 1
        try
            wait(acceptor)
            nsuccess += 1
        catch ex
            @error("Could not decipher worker connection", exception=ex)
        end
    end

    (nconnected == np) || @warn("Not all workers connected back", nconnected, nsuccess, nrequested=np)
    notify(c)
    nothing
end

function manage(manager::YarnManager, id::Integer, config::WorkerConfig, op::Symbol)
    @debug("YarnManager manage", id, op, container=config.userdata, nprocs=nprocs())
    # This function needs to exist, but so far we don't do anything
    if op in [:deregister, :finalize, :interrupt]
        if (config.userdata !== nothing) && haskey(config.userdata, :container_id)
            cid = config.userdata[:container_id]
            cidstr = container_id_string(cid)
            if cidstr in manager.am.containers.active
                container_stop(manager, cid)
                container_release(manager, cid)
            end
        end
        !manager.keep_connected && (1 == nprocs()) && disconnect(manager)
    end
    nothing
end


text:
mutable struct LazyNode{T}
	data::Union{T, Nothing}
	realised::Bool
	f::Union{Function, Nothing}
	LazyNode{T}(f::Function) where T = new{T}(nothing, false, f)
	LazyNode(v::T) where T = new{T}(v, true, nothing)
end

export LazyArray
struct LazyArray{T, N}
	data::Array{LazyNode{T}, N}
	LazyArray{T, N}(dims::Int...) where T where N =
		new{T, N}(Array{LazyNode{T}, N}(undef, dims))
end

function LazyArray{T}(dims::Int...) where T
	LazyArray{T, length(dims)}(dims...)
end

function LazyArray(l::Array{T}) where T
	r = LazyArray{T}(size(l)...)
	r.data .= map(x -> LazyNode(x), l)
	r
end

LazyArray{T, N}(l::Array) where T where N = LazyArray(l)

export LazyVector
LazyVector{T} = LazyArray{T, 1}

LazyVector(l::Array{T}) where T = LazyArray(l)

import Base.get
export get

get(n::LazyNode) = n.realised ? n.data :
	(n.realised = true; n.data = n.f(); n.f = nothing; n.data)

get(l::LazyArray, indices::Int...) = get(l.data[indices...])

lazyget(l::LazyArray{T}, indices::Int...) where T =
	LazyNode{T}(() -> (l[indices...]))

import Base.getindex
export getindex

getindex(l::LazyArray, indices::Int...) = get(l, indices...)

function getindex(l::LazyArray{T}, range::AbstractRange) where T
	r = LazyVector{T}(length(range))
	for (i, j) in zip(1 : length(range), range)
		r[i] = lazyget(l, j)
	end
	r
end

getindex(l::LazyArray, ::Colon) = l[1:end]

set(n::LazyNode{T}, v::T) where T =
	(n.data = v; n.realised = true; n.f = nothing)

set(n::LazyNode, f::Function) = (n.realised = false; n.f = f)

set(l::LazyArray{T}, v::Union{T, Function}, indices::Int...) where T =
	set(l.data[indices...], v)

set(l::LazyArray, v::LazyNode, indices::Int...) = (l.data[indices...] = v)

import Base.setindex!
export setindex!

setindex!(l::LazyArray, v, indices::Int...) = set(l, v, indices...)

import Base.length
export length

length(l::LazyArray) = length(l.data)

import Base.size
export size
size(l::LazyArray) = size(l.data)

import Base.lastindex
export lastindex
lastindex(l::LazyArray) = length(l)

import Base.iterate
export iterate

iterate(iter::LazyArray, state::Int=1) = state > length(iter) ?
	nothing : (iter[state], state + 1)

import Base.show
export show

show(io::IO, n::LazyNode) = print(io, n.realised ? string(n.data) : "pending")

show(io::IO, l::LazyArray) = print(io, string(l.data))

import Base.map
export map

function map(f::Function, l::LazyArray{T}) where T
	r = LazyArray{T}(size(l)...)
	for i in 1 : length(r)
		r.data[i] = LazyNode{T}(() -> f(l[i]))
	end
	r
end

import Base.map!
export map!

function map!(f::Function, l::LazyArray{T}) where T
	ghost = copy(l.data)
	for i in 1 : length(l)
		l.data[i] = LazyNode{T}(() -> f(get(ghost[i])))
	end
end

export realize!
realize!(l::LazyArray) = for i in l end


text:
using PyCall
anki_vector = pyimport("anki_vector")
util = pyimport("anki_vector.util")

args = anki_vector.util.parse_command_args()
@pywith anki_vector.Robot(args.serial) as robot begin
  robot.behavior.drive_off_charger()

  import Pkg
  Pkg.activate(".")

  using DistancePlan
  # include("../src/GraphDistancePlan.jl")
  # using GraphDistancePlan
  import Makie
  using StaticArrays

  scene = Makie.Scene()

  robot_theta = pi / 4
  robot_x = -250.
  robot_y = -250.

  B = 300.

  obstacles = [
      Box(SVector{2}([-150. -150.]), SVector{2}([150. 150.])),

      Box(SVector{2}([-B -B]), SVector{2}([-B B])), # Left
      Box(SVector{2}([B -B]), SVector{2}([B B])), # Right
      Box(SVector{2}([-B B]), SVector{2}([B B])), # Top
      Box(SVector{2}([-B -B]), SVector{2}([B -B])), # Bottom
     ]
  bounds = Box(SVector{2}([-B -B]), SVector{2}([B B]))
  goal = SVector{2}([250. 250.])
  origin = SVector{2}([robot_x robot_y])

  println("Planning...")
  tree, path = evt_goal(origin, bounds, obstacles, goal, min_radius=30.)
  push!(path, goal)
  println("Planning complete!")
  println("path = ", path)

  DistancePlan.draw!(origin)
  DistancePlan.draw!(bounds)
  DistancePlan.draw!(tree)
  DistancePlan.draw!.(obstacles)
  for point in path
    DistancePlan.draw!(point, color="purple")
  end
  Makie.plot!([x for (x, y) in path], [y for (x, y) in path], color="black")

  Makie.save("drive_vector.png", scene, resolution=(1000, 1000))

  for (x, y) in path
    println("going to (", x, ", ", y, ")")
    dx = x - robot_x
    dy = y - robot_y
    ego_theta = mod2pi(atan(dy, dx) - robot_theta)
    if ego_theta > pi
      ego_theta = ego_theta - 2pi
    end
    println("turning ", ego_theta, "rads")
    robot.behavior.turn_in_place(util.radians(ego_theta))
    println("driving ", sqrt(dx ^ 2 + dy ^ 2), "mm")
    robot.behavior.drive_straight(util.distance_mm(sqrt(dx ^ 2 + dy ^ 2)),
                                  util.speed_mmps(50))
    robot_x = x
    robot_y = y
    robot_theta = mod2pi(robot_theta + ego_theta)
    println("now at (", robot_x, ", ", robot_y, ")")
  end
end

max_stars_repo_name: krzentner/DistancePlan

max_issues_repo_name: krzentner/DistancePlan

max_forks_repo_name: krzentner/DistancePlan


text:
using URIParser
using FilePaths

@testset "URI" begin
    @test string(URI(p"/foo/bar")) == "file:///foo/bar"
    @test string(URI(p"/foo foo/bar")) == "file:///foo%20foo/bar"
    @test_throws ArgumentError URI(p"foo/bar")
    @test string(URI(WindowsPath("C:\\foo\\bar"))) == "file:///C:/foo/bar"
    @test string(URI(p"/foo/bar", query="querypart", fragment="fragmentpart")) == "file:///foo/bar?querypart#fragmentpart"
end


text:
module Rotations

using StaticArrays
using LinearAlgebra: norm, rank, svd, Diagonal
using SHIPs: _mrange, IntS
using SHIPs.SphericalHarmonics: index_y

export ClebschGordan, CoeffArray, single_B

"""
`ClebschGordan: ` storing precomputed Clebsch-Gordan coefficients; see
`?clebschgordan` for the convention that is use.
"""
struct ClebschGordan{T}
	# maxL::Int
	# cg::Array{T, 3}   # rewrite as Dict!
	vals::Dict{Tuple{IntS, IntS, IntS}, T}
end

"""
`CoeffArray: ` storing recursively precomputed coefficients for a
rotation-invariant basis.
"""
struct CoeffArray{T}
   vals::Vector{Dict}
   cg::ClebschGordan{T}
end

# ----------------------------------------------------------------------
#     ClebschGordan code
# ----------------------------------------------------------------------

cg_conditions(j1,m1,j2,m2,j3,m3) =
	cg_l_condition(j1, j2, j3)   &&
	cg_m_condition(m1, m2, m3)   &&
	(abs(m1) <= j1) && (abs(m2) <= j2) && (abs(m3) <= j3)

cg_l_condition(j1, j2, j3) = (abs(j1-j2) <= j3 <= j1 + j2)

cg_m_condition(m1, m2, m3) = (m3 == m1 + m2)

"""
`clebschgordan(j1, m1, j2, m2, j3, m3, T=Float64)` :

A reference implementation of Clebsch-Gordon coefficients based on

https://hal.inria.fr/hal-01851097/document
Equation (4-6)

This heavily uses BigInt and BigFloat and should therefore not be employed
for performance critical tasks.

The ordering of parameters corresponds to the following convention:
```
clebschgordan(j1, m1, j2, m2, j3, m3) = C_{j1m1j2m2}^{j3m3}
```
where
```
   D_{m1k1}^{l1} D_{m2k2}^{l2}}
	=
	∑_j  C_{l1m1l2m2}^{j(m1+m2)} C_{l1k1l2k2}^{j2(k1+k2)} D_{(m1+m2)(k1+k2)}^{j}
```
"""
function clebschgordan(j1, m1, j2, m2, j3, m3, T=Float64)
	if !cg_conditions(j1, m1, j2, m2, j3, m3)
		return zero(T)
	end

   N = (2*j3+1) *
       factorial(big(j1+m1)) * factorial(big(j1-m1)) *
       factorial(big(j2+m2)) * factorial(big(j2-m2)) *
       factorial(big(j3+m3)) * factorial(big(j3-m3)) /
       factorial(big( j1+j2-j3)) /
       factorial(big( j1-j2+j3)) /
       factorial(big(-j1+j2+j3)) /
       factorial(big(j1+j2+j3+1))

   G = big(0)
   # 0 ≦ k ≦ j1+j2-j3
   # 0 ≤ j1-m1-k ≤ j1-j2+j3   <=>   j2-j3-m1 ≤ k ≤ j1-m1
   # 0 ≤ j2+m2-k ≤ -j1+j2+j3  <=>   j1-j3+m2 ≤ k ≤ j2+m2
   lb = (0, j2-j3-m1, j1-j3+m2)
   ub = (j1+j2-j3, j1-m1, j2+m2)
   for k in maximum(lb):minimum(ub)
      bk = big(k)
      G += (-1)^k *
           binomial(big( j1+j2-j3), big(k)) *
           binomial(big( j1-j2+j3), big(j1-m1-k)) *
           binomial(big(-j1+j2+j3), big(j2+m2-k))
   end

   return T(sqrt(N) * G)
end

ClebschGordan(T=Float64) =
	ClebschGordan(Dict{Tuple{IntS,IntS,IntS}, T}())

_cg_key(j1, m1, j2, m2, j3, m3) =
	IntS.((index_y(j1,m1), index_y(j2,m2), index_y(j3,m3)))

function (cg::ClebschGordan{T})(j1, m1, j2, m2, j3, m3) where {T}
	if !cg_conditions(j1,m1,j2,m2,j3,m3)
		return zero(T)
	end
	key = _cg_key(j1, m1, j2, m2, j3, m3)
	if haskey(cg.vals, key)
		return cg.vals[key]
	end
	val = clebschgordan(j1, m1, j2, m2, j3, m3, T)
	cg.vals[key] = val
	return val
end

# ----------------------------------------------------------------------
#     CoeffArray code
# ----------------------------------------------------------------------

dicttype(N::Integer) = dicttype(Val(N))

dicttype(::Val{N}) where {N} =
   Dict{Tuple{SVector{N,Int8}, SVector{N,Int8}, SVector{N,Int8}}, Float64}

CoeffArray(T=Float64) = CoeffArray(Dict[], ClebschGordan(T))

function get_vals(A::CoeffArray, valN::Val{N}) where {N}
	if length(A.vals) < N
		for n = length(A.vals)+1:N
			push!(A.vals, dicttype(n)())
		end
	end
   return A.vals[N]::dicttype(valN)
end

_key(ll::StaticVector{N}, mm::StaticVector{N}, kk::StaticVector{N}) where {N} =
      (SVector{N, Int8}(ll), SVector{N, Int8}(mm), SVector{N, Int8}(kk))

function (A::CoeffArray)(ll::StaticVector{N},
                         mm::StaticVector{N},
                         kk::StaticVector{N}) where {N}
   if       sum(mm) != 0 ||
            sum(kk) != 0 ||
            !all(abs.(mm) .<= ll) ||
            !all(abs.(kk) .<= ll)
      return 0.0
   end
   vals = get_vals(A, Val(N))  # this should infer the type!
   key = _key(ll, mm, kk)
   if haskey(vals, key)
      val  = vals[key]
   else
      val = _compute_val(A, key...)
      vals[key] = val
   end
   return val
end

function (A::CoeffArray)(ll::StaticVector{1},
                         mm::StaticVector{1},
                         kk::StaticVector{1})
   if ll[1] == mm[1] == kk[1] == 0
      return 1.0
   else
      return 0.0
   end
end

function (A::CoeffArray)(ll::StaticVector{2},
                         mm::StaticVector{2},
                         kk::StaticVector{2})
   if ll[1] != ll[2] || sum(mm) != 0 || sum(kk) != 0
      return 0.0
   else
      return 8 * pi^2 / (2*ll[1]+1) * (-1)^(mm[1]-kk[1])
   end
end

function _compute_val(A::CoeffArray, ll::StaticVector{N},
                                     mm::StaticVector{N},
                                     kk::StaticVector{N}) where {N}
	val = 0.0
   llp = ll[1:N-2]
   mmp = mm[1:N-2]
   kkp = kk[1:N-2]
   for j = abs(ll[N-1]-ll[N]):(ll[N-1]+ll[N])
      if abs(kk[N-1]+kk[N]) > j || abs(mm[N-1]+mm[N]) > j
         continue
      end
		cgk = try
			A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
		catch
			@show (ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
			0.0
		end
		cgm = A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N])
		if cgk * cgm  != 0
			val += cgk * cgm * A( SVector(llp..., j),
								       SVector(mmp..., mm[N-1]+mm[N]),
								       SVector(kkp..., kk[N-1]+kk[N]) )
		end
   end
   return val
end

_len_mrange(ll) = sum(_ -> 1, _mrange(ll))

function basis(A::CoeffArray{T}, ll; ordered=false) where {T}
	CC = compute_Al(A, ll, Val(ordered))
	svdC = svd(CC)
	rk = rank(Diagonal(svdC.S))
	return svdC.U[:, 1:rk]
end

compute_Al(ll::SVector{N}; ordered = false) where {N} =
		compute_Al(CoeffArray(N, sum(ll)), ll; ordered=ordered)

compute_Al(A::CoeffArray, ll::SVector{N}; ordered = false) where {N} =
		compute_Al(A, ll, Val(ordered))

# unordered
function compute_Al(A::CoeffArray{T}, ll::SVector, ::Val{true}) where {T}
	num_mm_sorted = sum(mm -> issorted(mm), _mrange(ll))
	@show num_mm_sorted
	num_mm = _len_mrange(ll)
   CC = zeros(T, num_mm, num_mm_sorted)
	im = 0
   for mm in _mrange(ll)
		if issorted(mm)
			im += 1
			for (ik, kk) in enumerate(_mrange(ll))
		      CC[ik, im] = A(ll, mm, kk)
			end
		end
	end
   return CC
end

# ordered
function compute_Al(A::CoeffArray{T}, ll::SVector, ::Val{false}) where {T}
	len = _len_mrange(ll)
   CC = zeros(T, len, len)
   for (im, mm) in enumerate(_mrange(ll)), (ik, kk) in enumerate(_mrange(ll))
      CC[ik, im] = A(ll, mm, kk)
   end
   return CC
end

function single_B(A::CoeffArray{T}, ll::SVector) where {T}
	MM = collect(_mrange(ll))
	Is = sortperm([ (sum(mm .!= 0), norm(mm)) for mm in MM ])
	MM = MM[Is]
	CC = zeros(T, length(MM))
	for mm in MM
		for (ik, kk) in enumerate(_mrange(ll))
			CC[ik] = A(ll, mm, kk)
		end
		if norm(CC) > 0
			CC ./= norm(CC)
			break
		end
	end
	return CC
end

end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using NiLang, NiLang.AD, Test
using NiLang.AD: hessian_numeric

@testset "hessian" begin
    h1 = hessian_backback(PlusEq(*), (0.0, 2.0, 3.0); iloss=1)
    h2 = hessian_numeric(PlusEq(*), (0.0, 2.0, 3.0); iloss=1)
    @test h1 ≈ h2

    @i function test(a,b,c,d)
        a += b*c
        a += b^d
        c += b/d
        ROT(a, c, d)
        b += d ^ 2
        a += c * d
    end
    h1 = hessian_backback(test, (0.0, 2.0, 1.0, 3.0); iloss=1)
    h2 = hessian_numeric(test, (0.0, 2.0, 1.0, 3.0); iloss=1)
    @show h2
    @test isapprox(h1, h2, atol=1e-8)
end


text:
@info("ExcelWriter: running build.jl")
using Pkg

try
    println("")
    println("Building PyCall ...")
    ENV["PYTHON"] = ""
    Pkg.build("Conda")
    Pkg.build("PyCall")
    using PyCall

    global const pyModPandas = PyNULL()
    global const pyModxlsxwriter = PyNULL()
    global const pyModnumpy = PyNULL()

    @info("ExcelWriter build: adding pandas...")
    copy!(pyModPandas, pyimport_conda("pandas", "pandas"))
    @info("ExcelWriter build: adding numpy...")
    copy!(pyModnumpy, pyimport_conda("numpy","numpy"))
    @info("ExcelWriter build: adding xlsxwriter...")
    copy!(pyModxlsxwriter, pyimport_conda("xlsxwriter", "xlsxwriter"))
    @info("ExcelWriter build: testing pyModPandas.ExcelWriter...")
    writer = pyModPandas.ExcelWriter("blabla_deletme.xlsx", engine="xlsxwriter")
    @show typeof(writer)
    println("If this worked, then pandas and xlsxwriter should be installed propoerly");
    println("Done building PyCall.")

catch e

    @warn("ExcelWriter: Build script failed.")
    println("ExcelWriter: Build script failed.")
    @show e

end
# exit()

@info("ExcelWriter: build script finished.")


text:
const BoolTypes = Union{Bool}
const IntTypes = Union{Int8, Int16, Int32, Int64, Int128}
const UIntTypes = Union{UInt8, UInt16, UInt32, UInt64, UInt128}
const IntegerTypes = Union{BoolTypes, IntTypes, UIntTypes, Ptr}
const FloatingTypes = Union{Float16, Float32, Float64}
const ScalarTypes = Union{IntegerTypes, FloatingTypes}

type_length(::Type{Vec{N,T}}) where {N,T} = N
type_size(::Type{Vec{N,T}}) where {N,T} = (N,)
type_size(::Type{Vec{N,T}}, n::Integer) where {N,T} = (N,)[n]

# Element-wise access

@generated function vsetindex(v::AbstractSIMDVector{N,T}, x::Number, ::Val{I}) where {N,T,I}
    @assert isa(I, Integer)
    1 <= I <= N || throw(BoundsError())
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    vtyp = "<$N x $typ>"
    decls = String[]
    instrs = String[]
    push!(instrs, "%res = insertelement $vtyp %0, $typ %1, $ityp $(I-1)")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}, T}, extract_data(v), T(x)))
    end
end

@generated function vsetindex(v::AbstractSIMDVector{N,T}, x::Number, i::Int) where {N,T}
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    vtyp = "<$N x $typ>"
    decls = String[]
    instrs = String[]
    push!(instrs, "%res = insertelement $vtyp %0, $typ %2, $ityp %1")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        @boundscheck 1 <= i <= N || throw(BoundsError())
        Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,VE{T}}, Tuple{NTuple{N,VE{T}}, Int, T},
            extract_data(v), i-1, T(x)))
    end
end

@inline Base.setindex(v::AbstractStructVec{N,T}, x::Number, i) where {N,T} = setindex(v, x, i)
@inline Base.setindex(v::AbstractStructVec{N,T}, x::Number, i::Integer) where {N,T} = setindex(v, x, Int(i))

@inline getvalindex(v::AbstractSIMDVector{N,T}, ::Val{I}) where {N,T,I} = extract_data(v)[I].value

@inline function vbroadcast(::Type{Vec{N,T}}, s::S) where {N,T,S<:ScalarTypes}
    @inbounds ntuple(i -> VE{T}(s), Val(N))
end
@inline function vbroadcast(::Type{Vec{N,T}}, s::Core.VecElement{T}) where {N,T}
    @inbounds ntuple(i -> s, Val(N))
end
@inline vbroadcast(::Type{Vec{N,T}}, s::Vec{N,T}) where {N,T} = s
@inline function vbroadcast(::Type{Vec{N,T}}, ptr::Ptr{T}) where {N,T}
    s = Core.VecElement{T}(VectorizationBase.load(ptr))
    @inbounds ntuple(_ -> s, Val(N))
end
@inline function vbroadcast(::Type{Vec{N,T}}, ptr::Ptr) where {N,T}
    s = Core.VecElement{T}(VectorizationBase.load(Base.unsafe_convert(Ptr{T},ptr)))
    @inbounds ntuple(_ -> s, Val(N))
end
@inline function svbroadcast(::Type{SVec{N,T}}, s::S) where {N,T,S<:ScalarTypes}
    @inbounds SVec(ntuple(i -> VE{T}(s), Val(N)))
end
@inline svbroadcast(::Type{SVec{N,T}}, s::SVec{N,T}) where {N,T} = s
@inline function vbroadcast(::Type{SVec{N,T}}, s::S) where {N,T,S<:ScalarTypes}
    @inbounds SVec(ntuple(i -> VE{T}(s), Val(N)))
end
@inline function vbroadcast(::Type{SVec{N,T}}, s::Core.VecElement{T}) where {N,T}
    @inbounds SVec(ntuple(i -> s, Val(N)))
end
@inline vbroadcast(::Type{SVec{N,T}}, s::Vec{N,T}) where {N,T} = SVec(s)
@inline vbroadcast(::Type{SVec{N,T}}, s::SVec{N,T}) where {N,T} = s

@inline vone(::Type{Vec{N,T}}) where {N,T} = ntuple(i -> Core.VecElement(one(T)), Val(N))
@inline vzero(::Type{Vec{N,T}}) where {N,T} = ntuple(i -> Core.VecElement(zero(T)), Val(N))
@inline vone(::Type{SVec{N,T}}) where {N,T} = SVec(ntuple(i -> Core.VecElement(one(T)), Val(N)))
@inline vzero(::Type{SVec{N,T}}) where {N,T} = SVec(ntuple(i -> Core.VecElement(zero(T)), Val(N)))
@inline vone(::Type{T}) where {T} = one(T)
@inline vzero(::Type{T}) where {T} = zero(T)

@inline Vec{N,T}(v::Vararg{T,N}) where {T,N} = ntuple(n -> VE(v[n]), Val(N))

@inline function pirate_convert(::Type{Vec{N,T}}, xs::NTuple{N,T}) where {N,T<:ScalarTypes}
    @inbounds ntuple(i -> VE(xs[i]), Val(N))
end
@inline function pirate_convert(::Type{Vec{N,T1}}, xs::Vec{N,T2}) where {N,T1<:ScalarTypes,T2<:ScalarTypes}
    @inbounds ntuple(i -> VE(T1(xs[i].value)), Val(N))
end
@generated function Base.convert(::Type{SVec{N,T}}, xs::NTuple{N,T}) where {N,T}
    quote
        $(Expr(:meta,:inline))
        @inbounds SVec((Base.Cartesian.@ntuple $N n -> VE(xs[n])))
    end
end
@generated function Base.convert(::Type{SVec{N,T1}}, xs::SVec{N,T2}) where {N,T1,T2}
    quote
        $(Expr(:meta,:inline))
        @inbounds SVec((Base.Cartesian.@ntuple $N n -> VE{T1}(xs[n])))
    end
end
@inline Base.convert(::Type{SVec{N,T}}, x::SVec{N,T}) where {N,T} = x

# @inline similar(s::T, ::Vec{N,T}, ::Vec{N,T}) = ntuple(i -> VE{T}(s), Val(N))
# @inline similar(s::T, ::Vec{N,T}, ::Vec{N,T}) = ntuple(i -> VE{T}(s), Val(N))

# Floating point formats

int_type(::Type{Float16}) = Int16
int_type(::Type{Float32}) = Int32
int_type(::Type{Float64}) = Int64
# int_type(::Type{Float128}) = Int128
# int_type(::Type{Float256}) = Int256

uint_type(::Type{Float16}) = UInt16
uint_type(::Type{Float32}) = UInt32
uint_type(::Type{Float64}) = UInt64
# uint_type(::Type{Float128}) = UInt128
# uint_type(::Type{Float256}) = UInt256

significand_bits(::Type{Float16}) = 10
significand_bits(::Type{Float32}) = 23
significand_bits(::Type{Float64}) = 52
# significand_bits(::Type{Float128}) = 112
# significand_bits(::Type{Float256}) = 136

exponent_bits(::Type{T}) where {T<:FloatingTypes} =
    8*sizeof(T) - 1 - significand_bits(T)
sign_bits(::Type{T}) where {T<:FloatingTypes} = 1

significand_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(uint_type(T)(1) << significand_bits(T) - 1)
exponent_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(uint_type(T)(1) << exponent_bits(T) - 1) << significand_bits(T)
sign_mask(::Type{T}) where {T<:FloatingTypes} =
    uint_type(T)(1) << (significand_bits(T) + exponent_bits(T))

@generated function vecbool_to_unsigned(vb::Vec{N,Bool}) where {N}
    Nout = VectorizationBase.nextpow2(max(N, 8))
    Utype = VectorizationBase.mask_type(Nout)
    btype = llvmtype(Bool)
    vbtyp = "<$N x $btype>"
    decls = String[]
    instrs = String["%maskb = trunc $vbtyp %0 to <$N x i1>"]
    u1name = Nout == N ? "%uout" : "%uabridged"
    push!(instrs, "$u1name = bitcast <$N x i1> %maskb to i$N")
    if Nout > N
        push!(instrs, "%uout = zext i$N %uabridged to i$Nout")
    end
    push!(instrs, "ret i$Nout %uout")
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall($(join(instrs, "\n")), $Utype, Tuple{Vec{$N,Bool}}, vb)
    end
end


text:
type OIM<:FiniteHorizonAgent
    policy::Matrix{UInt64}
    ns::Array{UInt64, 4}
    phat::Array{Float64, 4}
    n::Array{UInt64, 3}
    rewards::Array{Float64, 3}
    optR::Vector{Float64}
    Ve::Matrix{Float64}
    Vr::Matrix{Float64}
    δ::Float64
    ϵ::Float64
    rewards_known::Bool
end

function OIM(S, A, H, rewards, δ, ϵ, keep_updating=false)
    policy = zeros(UInt64, S, H)
    ns = zeros(UInt64, S+1, S, A, H)
    phat = zeros(S+1, S, A, H)
    # pretend we have seen one transition to optimistic state S+1 from each state
    n = ones(UInt64, S, A, H)
    ns[end, :, :, :] = 1
    max_rewards = vec(maximum(reshape(rewards, S*A, H), 1))
    # rough guess of required optimistic rewards, exact value will require transferring entire theory to our setting
    #optR = 8 * cumsum(max_rewards).^2 * H^2 / ϵ^2 * log(S * A * H / ϵ / δ)
    optR = cumsum(max_rewards).^2 * H  / ϵ * log(S * A * H / δ)

    m = OIM(policy, ns, phat, n, rewards, optR, zeros(S, H+1),zeros(S, H+1), δ, ϵ, true)
    rand!(m.policy, 1:A) # Initialize with uniformly random policy
    m
end

function update_policy!(m::OIM)
    S = nS(m)
    H = horizon(m)
    A = nA(m)

    Qe = zeros(A)
    Qr = zeros(A)

    m.phat .= m.ns ./ reshape(m.n, 1, S, A, H)
    for t=H:-1:1
        for s ∈ 1:S
            for a ∈ 1:A
                Qe[a] = m.optR[t] * m.phat[S+1, s, a, t] + sum(m.phat[1:S, s, a, t] .* m.Ve[1:S,t+1])
                Qr[a] = m.rewards[s,a,t] + sum(m.phat[1:S, s, a, t] .* m.Vr[1:S,t+1])
            end
            bestA = indmax(Qe+Qr)
            m.policy[s, t] = bestA
            m.Ve[s, t] = Qe[bestA]
            m.Vr[s, t] = Qr[bestA]
        end
    end
end


text:
# Test functionality with continuous features.
@testset "TuRF - Continuous Features" begin
    data = rand(1000, 10)
    for idx1 = 1:size(data, 2) - 1
        for idx2 = idx1+1:size(data, 2)
            target = Int64.(data[:, idx1] .> data[:, idx2])
            weights = Relief.turf(data, target, 5, rba=Relief.relieff)
            @test all(weights[idx1] .>= weights[(1:end .!= idx1) .& (1:end .!= idx2)])
            @test all(weights[idx2] .>= weights[(1:end .!= idx1) .& (1:end .!= idx2)])
        end
    end
end

# Test functionality with discrete features.
@testset "TuRF - Discrete Features" begin
    data = rand([0, 1, 2, 3], 1000, 10)
    for idx1 = 1:size(data, 2) - 1
        for idx2 = idx1+1:size(data, 2)
            target = Int64.(data[:, idx1] .> data[:, idx2])
            weights = Relief.turf(data, target, 5, rba=Relief.relieff)
            @test all(weights[idx1] .>= weights[(1:end .!= idx1) .& (1:end .!= idx2)])
            @test all(weights[idx2] .>= weights[(1:end .!= idx1) .& (1:end .!= idx2)])
        end
    end
end

# Test functionality with default RBA (ReliefF).
@testset "TuRF - Default RBA" begin
    data = rand([0, 1, 2, 3], 1000, 10)
    idx1, idx2 = 1, 2
    target = Int64.(data[:, idx1] .> data[:, idx2])
    weights = Relief.turf(data, target, 5)
    @test all(weights[idx1] .>= weights[(1:end .!= idx1) .& (1:end .!= idx2)])
    @test all(weights[idx2] .>= weights[(1:end .!= idx1) .& (1:end .!= idx2)])

end


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

using Markdown
using InteractiveUtils

# ╔═╡ e3725bf4-a3dc-11eb-2ebb-bf94a858f39a
using SparkSQL, DataFrames, Decimals, Dates

# ╔═╡ cccca8da-6ee8-4ce1-9f46-fd13dd6c209b
using StatsPlots

# ╔═╡ c031a91e-ef69-4b98-9b7a-66347f65f908
# Stock Analysis Example with Julia & Spark using the SparkSQL.jl package.

# ╔═╡ b0c17e9b-21b9-44cc-924c-35f656a10793
# SparkSQL and DataFrames packages are required; Decimals, and Dates are needed if those types are present in your data.

# ╔═╡ 863e724b-44a8-4dc6-b91f-4b596a167568
# Use initJVM to initialize the Java Virtual Machine.

# ╔═╡ b2203415-7169-41dc-8f80-bc843f5ed06c
SparkSQL.initJVM()

# ╔═╡ 337c4b04-6ae5-495a-aa47-22bc7bdb66b4
# SparkSession is the entry point to using the Apache Spark Dataset & Dataframe API. Pass the the url for Spark, and your app name as arguments.  (Additional arguments are supported, see SparkSQL.jl package help for more details)

# ╔═╡ 58199e82-0651-414c-863d-06cee4ce61ca
sprk = SparkSQL.SparkSession("spark://localhost:7077", "Julia Spark Stock App");

# ╔═╡ 689913c0-af7e-4616-aafa-560c3537d43d
# From Julia, instruct Spark to query data using a CSV file. The example file does not have column headers, so "_c0, _c1, ... _cN" is used to refer the columns positionally to alias column names.

# ╔═╡ d5ce8110-6afa-4a31-8135-50e6d9fc938f
loadData = sql(sprk, "SELECT _c0 TICKER, _c1 DATE,_c2 OPEN,_c3 HIGH,_c4 LOW,_c4 CLOSE,_c6 VOL FROM CSV.`sampleData.csv`;");

# ╔═╡ 8d1dd0b1-f422-4ed6-b540-f164d574a559
# Create a temporary database view in Spark with Julia to provide

# ╔═╡ 8a78b16e-316e-4639-ba72-a4fd6c326822
createOrReplaceTempView(loadData, "rawdata");

# ╔═╡ 112e6688-220e-4a88-8fa5-711c0b45f1d7
# Prep the data in Spark from Julia by casting the string values into the correct data types.  Notice we use the database view name we just created in the FROM clause.

# ╔═╡ 6b228ad4-60f8-4104-8c71-ac88d6444f27
 prepData= sql(sprk,"SELECT TICKER, to_date(DATE, 'yyyymmdd' ) AS DATE, cast(OPEN AS double) OPEN, cast(LOW as double) LOW, cast(HIGH AS double) AS HIGH, cast(CLOSE as double) as CLOSE, cast(VOL as integer) as VOL FROM rawdata")

# ╔═╡ 719cd252-573e-408c-8694-974e5a455860
#Create a database view for the prepared data, "stocks".

# ╔═╡ c489dd64-6abc-417f-a640-2fe235c8a94b
createOrReplaceTempView(prepData,"stocks")

# ╔═╡ aeb0d1df-70c9-47e0-a815-fffca74bda73
#Create an aggregate query in Spark from Julia to get the average, min and max close price for IBM stock by year.

# ╔═╡ d99f6039-6bc7-4d38-aa12-9cf0660a59ee
queryTickerIBM = sql(sprk, "SELECT YEAR(DATE) AS TRADE_YEAR, AVG(CLOSE) AVG_CLOSE, MIN(CLOSE) MIN_CLOSE, MAX(CLOSE) MAX_CLOSE FROM stocks WHERE TICKER = 'IBM' AND YEAR(DATE) > 2007 GROUP BY YEAR(DATE) ORDER BY YEAR(DATE) DESC")

# ╔═╡ 9e2a359f-18e0-43e5-aa02-0f3884e93f24
# Move the query results from Spark to a Julia DataFrame with the "toJuliaDF" function.

# ╔═╡ ddfc9503-e132-4832-869e-9db1ee79b273
ibmTicker = toJuliaDF(queryTickerIBM)

# ╔═╡ 80e4f0ea-ded8-49f4-ac0a-181e23de530c
# Create a new Spark query to get the daily trade prices after 2007.

# ╔═╡ f0886653-6f2a-4e0a-98fa-afa2b5e26417
queryIBMdailyPrice = sql(sprk, "SELECT TICKER, DATE, OPEN, HIGH, LOW, CLOSE, VOL FROM stocks WHERE TICKER = 'IBM' AND YEAR(DATE) > 2007")

# ╔═╡ e8d16d90-e0cb-49d3-a62f-34091ed7a9c0
# Move the daily stock price query data from Spark to a Julia DataFrame.

# ╔═╡ e01e724f-198b-400a-832d-38186d587dc2
ibmTickerDaily = toJuliaDF(queryIBMdailyPrice)

# ╔═╡ 76ec754a-f397-43c9-a34a-780aa0f6fa82
# The data is now in a Julia DataFrame. Let's create a chart.

# ╔═╡ 08d03e24-2e08-4d99-9d0a-c5bad13d3153
@df ibmTickerDaily density(:CLOSE, label="IBM Stock Price Distribution 2008-2020", legend=:topleft)

# ╔═╡ 97e46d43-ac63-44da-90ed-dc959a422147
# Add data to the Julia DataFrame:

# ╔═╡ eb35660f-a0c8-4f76-9529-5a23cc627fb5
push!(ibmTickerDaily, ["IBM", Dates.Date("9999-1-30", dateformat"yyyy-mm-dd"), 100, 100, 100, 100, 1000000])

# ╔═╡ f0f0bb02-7c25-4f23-98b7-5708f16015e7
# Send the Julia DataFrame data back to Apache Spark.

# ╔═╡ e3aee473-c606-4977-b616-bdaeda545c3d
moveToSpark = toSparkDS(sprk, ibmTickerDaily)

# ╔═╡ 3ca39b89-016a-4126-a378-7b2678049eba
# The data is in Spark, create database view to use it with SQL.

# ╔═╡ 92459661-2d84-44ff-a2fa-8a9cce682716
createOrReplaceTempView(moveToSpark,"ticker_refresh")

# ╔═╡ 2e250ef2-6f53-4519-96a9-b37ed8f099e7
# Split the data back into columns and assign column names. (Note: "," is the default delimiter. To use a different delimiter, see SparkSQL.jl help for instructions.)

# ╔═╡ ae425c7d-69c5-4a2a-82d0-a98ec213d4a5
wrangleData = sql(sprk, "SELECT
	split(value, ',' )[0] AS TICKER,
	split(value, ',' )[1] AS DATE,
	split(value, ',' )[2] AS OPEN,
    split(value, ',' )[4] AS HIGH,
	split(value, ',' )[5] AS LOW,
	split(value, ',' )[6] AS CLOSE,
	split(value, ',' )[7] AS VOL
	FROM ticker_refresh")

# ╔═╡ 8291a7f5-77f2-4d3e-86d0-4f8c23c9fe48
# Update the view to reflect the new columns and names

# ╔═╡ 1c6788d2-c2ea-419c-acfe-ea16af1d28e1
createOrReplaceTempView(wrangleData, "ticker_refresh")

# ╔═╡ 0b4b82a8-f21e-422c-b380-993138ed6c3e
#  Cast the data types. (Note: format for date is now "yyyy-mm-dd"; the CSV file was "yyyymmdd")

# ╔═╡ 27a8d2ed-8068-4088-93c2-954eef97d284
castTypes = sql(sprk, "SELECT TICKER, to_date(DATE, 'yyyy-mm-dd' ) AS DATE, cast(OPEN AS double) OPEN, cast(HIGH AS double) AS HIGH, cast(LOW as double) LOW, cast(CLOSE as double) as CLOSE, cast(VOL as integer) as VOL FROM ticker_refresh")

# ╔═╡ b8d36775-2f69-4866-80fa-a246af321c6a
# Replace the "stocks" view with the updated data from the Julia DataFrame.

# ╔═╡ c0bbb615-bef7-4153-9077-03681aae9c14
createOrReplaceTempView(castTypes, "stocks")

# ╔═╡ 50b99063-7958-4eac-995e-5a8aa4316955
# Show the row added to the Julia DataFrame is now in Spark. Query to find max date in Spark. (The expected output is "9999-1-30")

# ╔═╡ 526a5b0d-2a20-4274-8b69-6604509ca62f
getTopRecord = sql(sprk, "SELECT MAX(DATE) AS MaxDate FROM stocks")

# ╔═╡ cc5bc142-e1c2-47d5-85f4-d56fa59bde8c
# View the SparkSQL query results with the "toJuliaDF" function.

# ╔═╡ 3eaed1b0-6339-4ec9-9113-2d2ebfba2efc
maxDate = toJuliaDF(getTopRecord)

# ╔═╡ daae34e7-4ba0-4108-85c6-6fb8ce695327
# Congratulations! You now know the fundementals using the SparkSQL.jl Julia package.

# ╔═╡ Cell order:
# ╠═c031a91e-ef69-4b98-9b7a-66347f65f908
# ╠═b0c17e9b-21b9-44cc-924c-35f656a10793
# ╠═e3725bf4-a3dc-11eb-2ebb-bf94a858f39a
# ╠═863e724b-44a8-4dc6-b91f-4b596a167568
# ╠═b2203415-7169-41dc-8f80-bc843f5ed06c
# ╠═337c4b04-6ae5-495a-aa47-22bc7bdb66b4
# ╠═58199e82-0651-414c-863d-06cee4ce61ca
# ╠═689913c0-af7e-4616-aafa-560c3537d43d
# ╠═d5ce8110-6afa-4a31-8135-50e6d9fc938f
# ╠═8d1dd0b1-f422-4ed6-b540-f164d574a559
# ╠═8a78b16e-316e-4639-ba72-a4fd6c326822
# ╠═112e6688-220e-4a88-8fa5-711c0b45f1d7
# ╠═6b228ad4-60f8-4104-8c71-ac88d6444f27
# ╠═719cd252-573e-408c-8694-974e5a455860
# ╠═c489dd64-6abc-417f-a640-2fe235c8a94b
# ╠═aeb0d1df-70c9-47e0-a815-fffca74bda73
# ╠═d99f6039-6bc7-4d38-aa12-9cf0660a59ee
# ╠═9e2a359f-18e0-43e5-aa02-0f3884e93f24
# ╠═ddfc9503-e132-4832-869e-9db1ee79b273
# ╠═80e4f0ea-ded8-49f4-ac0a-181e23de530c
# ╠═f0886653-6f2a-4e0a-98fa-afa2b5e26417
# ╠═e8d16d90-e0cb-49d3-a62f-34091ed7a9c0
# ╠═e01e724f-198b-400a-832d-38186d587dc2
# ╠═76ec754a-f397-43c9-a34a-780aa0f6fa82
# ╠═cccca8da-6ee8-4ce1-9f46-fd13dd6c209b
# ╠═08d03e24-2e08-4d99-9d0a-c5bad13d3153
# ╠═97e46d43-ac63-44da-90ed-dc959a422147
# ╠═eb35660f-a0c8-4f76-9529-5a23cc627fb5
# ╠═f0f0bb02-7c25-4f23-98b7-5708f16015e7
# ╠═e3aee473-c606-4977-b616-bdaeda545c3d
# ╠═3ca39b89-016a-4126-a378-7b2678049eba
# ╠═92459661-2d84-44ff-a2fa-8a9cce682716
# ╠═2e250ef2-6f53-4519-96a9-b37ed8f099e7
# ╠═ae425c7d-69c5-4a2a-82d0-a98ec213d4a5
# ╠═8291a7f5-77f2-4d3e-86d0-4f8c23c9fe48
# ╠═1c6788d2-c2ea-419c-acfe-ea16af1d28e1
# ╠═0b4b82a8-f21e-422c-b380-993138ed6c3e
# ╠═27a8d2ed-8068-4088-93c2-954eef97d284
# ╠═b8d36775-2f69-4866-80fa-a246af321c6a
# ╠═c0bbb615-bef7-4153-9077-03681aae9c14
# ╠═50b99063-7958-4eac-995e-5a8aa4316955
# ╠═526a5b0d-2a20-4274-8b69-6604509ca62f
# ╠═cc5bc142-e1c2-47d5-85f4-d56fa59bde8c
# ╠═3eaed1b0-6339-4ec9-9113-2d2ebfba2efc
# ╠═daae34e7-4ba0-4108-85c6-6fb8ce695327

max_stars_repo_name: propelledanalytics/Tutorials

max_issues_repo_name: propelledanalytics/Tutorials

max_forks_repo_name: propelledanalytics/Tutorials


text:
using PastaQ
using Random
using ITensors
using Observers
import Optimisers

Random.seed!(1234)

# 1. Quantum process tomography of a unitary circuit

# Make the random circuit
N = 4
depth = 4
nshots = 10_000
circuit = randomcircuit(N, depth)

# Generate samples
data, U = getsamples(circuit, nshots; local_basis=["X", "Y", "Z"], process=true)
writesamples(data, U, "data/qpt_circuit.h5")

# Load target state and measurements. Each samples is built out
# of an input state (`first.(data)`) to the quantum channel, and the
# measurement output (`last.(data)`) after a local basis rotation.
data, Û = readsamples("data/qpt_circuit.h5")

# Split data into a train and test dataset
train_data, test_data = split_dataset(data; train_ratio=0.9)

# Set parameters
N = length(Û)     # Number of qubits
χ = maxlinkdim(Û) # Bond dimension of variational MPS

# Initialize the unitary MPO
U0 = randomprocess(Û; χ=χ)

opt = Optimisers.Descent(0.01)

F(U::MPO; kwargs...) = fidelity(U, Û; process=true)
obs = Observer(["F" => F])

# Initialize stochastic gradient descent optimizer
@show maxlinkdim(U0)

# Run process tomography
println("Run process tomography to learn noiseless circuit U")
U = tomography(
  train_data,
  U0;
  test_data=test_data,
  optimizer=opt,
  batchsize=500,
  epochs=5,
  (observer!)=obs,
  print_metrics=["F"],
)

@show maxlinkdim(U)
println()

# Noisy circuit
# Generate samples
data, Λ = getsamples(
  circuit,
  nshots;
  local_basis=["X", "Y", "Z"],
  process=true,
  noise=("amplitude_damping", (γ=0.01,)),
)
writesamples(data, Λ, "data/qpt_circuit_noisy.h5")

# Load data and target Choi matrix
data, Φ = readsamples("data/qpt_circuit_noisy.h5")

# Split data into a train and test dataset
train_data, test_data = split_dataset(data; train_ratio=0.9)

# Set up
N = length(Φ)
χ = 8
ξ = 2

# Initialize the Choi LPDO
Λ0 = randomprocess(Φ; mixed=true, χ=χ, ξ=ξ)

# Initialize stochastic gradient descent optimizer
opt = Optimisers.ADAM()

F(Λ::LPDO; kwargs...) = fidelity(Λ, Φ; process=true)
obs = Observer(["F" => F])

# Run process tomography
println("Run process tomography to learn noisy process Λ")
@disable_warn_order begin
  Λ = tomography(
    train_data,
    Λ0;
    test_data=test_data,
    optimizer=opt,
    batchsize=500,
    epochs=5,
    (observer!)=obs,
    print_metrics=["F"],
  )
end
@show maxlinkdim(Λ.X)
println()


text:
module TextMining

using ASCIIPlots, LightXML

import Base: getindex, setindex!, isempty, keys, values,
    copy, length, haskey, display

abstract FeatureSpace

include("feature_vector.jl")
include("cluster.jl")
include("data_set.jl")
include("clustering.jl")
include("distribution.jl")
include("text_processing.jl")

export  #types
  Cluster, DataSet, Distribution, FeatureSpace, FeatureVector,

  #functions
  #classification
  knn,
  #cluster
  centroid, distance,
  #clustering
  random_init, max_min_init, kmeans, elbow_method,
  #data set
  #distribution
  probability, features, entropy, info_gain, perplexity, set_smooth,
  remove_smoothing!, delta_smoothing!,
  #feature vector
  sanitize!, freq_list, find_common_type, add!, subtract!,
  multiply!, divide!, rationalize!, dist_cos, dist_zero,
  dist_zero_weighted, dist_taxicab, dist_euclidean, dist_infinite,
  #text processing
  clean, parse_xml, load_featurevector, load_cluster, load_dataset

end

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
module Epidemiology

using AlgebraicPetri
using Petri
using Catlab
using Catlab.Theories

export InfectiousDiseases, FunctorGenerators, F_epi, S, E, I, R, D, transmission, exposure, illness, recovery, death

ob = PetriCospanOb(1)
spontaneous_petri = PetriCospan([1], Petri.Model(1:2, [(Dict(1=>1), Dict(2=>1))]), [2])
transmission_petri = PetriCospan([1], Petri.Model(1:2, [(Dict(1=>1, 2=>1), Dict(2=>2))]), [2])
exposure_petri = PetriCospan([1, 2], Petri.Model(1:3, [(Dict(1=>1, 2=>1), Dict(3=>1, 2=>1))]), [3, 2])

""" InfectiousDiseases
"""
@present InfectiousDiseases(FreeBiproductCategory) begin
    S::Ob
    E::Ob
    I::Ob
    R::Ob
    D::Ob
    transmission::Hom(S⊗I, I)
    exposure::Hom(S⊗I, E⊗I)
    illness::Hom(E,I)
    recovery::Hom(I,R)
    death::Hom(I,D)
end

""" generators
"""
S,E,I,R,D,transmission,exposure,illness,recovery,death = generators(InfectiousDiseases);

""" FunctorGenerators
"""
const FunctorGenerators = Dict(S=>ob, E=>ob, I=>ob, R=>ob, D=>ob,
        transmission=>transmission_petri, exposure=>exposure_petri,
        illness=>spontaneous_petri, recovery=>spontaneous_petri, death=>spontaneous_petri)

""" F_epi
"""
F_epi(ex) = functor((PetriCospanOb, PetriCospan), ex, generators=FunctorGenerators)

end


text:
###############################################################################
#
# Calculate the product-limited photosynthetic rates
#
###############################################################################
"""
    product_limited_rate!(
                photo_set::C3ParaSet{FT},
                leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}
    product_limited_rate!(
                photo_set::C4ParaSet{FT},
                leaf::Leaf{FT}
    ) where {FT<:AbstractFloat}
    product_limited_rate!(
                photo_set::C4ParaSet{FT},
                leaf::Leaf{FT},
                envir::AirLayer{FT}
    ) where {FT<:AbstractFloat}

Calculate the product limited photosynthetic rate, given
- `photo_set` [`C3ParaSet`](@ref) or [`C4ParaSet`](@ref) type struct
- `leaf` [`Leaf`](@ref) type struct
- `envir` [`AirLayer`](@ref) type struct
"""
function product_limited_rate!(
            photo_set::C3Cytochrome{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    @unpack p_i, Vcmax, Γ_star = leaf;
    @unpack Eff_1, Eff_2 = photo_set;

    leaf.Ap = leaf.Vcmax / 2;
    leaf.J_P680_p = leaf.Ac * (Eff_1*p_i + Eff_2*Γ_star) / (p_i - Γ_star);
    leaf.J_P700_p = leaf.J_P680_c * leaf.η;

    return nothing
end

function product_limited_rate!(
            photo_set::C3ParaSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Ap = leaf.Vcmax / 2;

    return nothing
end

function product_limited_rate!(
            photo_set::C4ParaSet{FT},
            leaf::Leaf{FT}
) where {FT<:AbstractFloat}
    leaf.Ap = leaf.Vpmax * leaf.p_i / (leaf.p_i + leaf.Kpep);

    return nothing
end

function product_limited_rate!(
            photo_set::C4ParaSet{FT},
            leaf::Leaf{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    @unpack g_lc, Kpep, Rd, Vpmax = leaf;
    @unpack p_a, p_atm = envir;

    _a = Vpmax;
    _d = Kpep;
    _f = p_atm / g_lc * FT(1e-6);
    _p = p_a;

    _qa = _f;
    _qb = _f*Rd - _p - _d - _a*_f;
    _qc = _a*_p - Rd*(_p + _d);
    _an = lower_quadratic(_qa, _qb, _qc);

    leaf.Ap = _an + Rd;

    return nothing
end

max_stars_repo_name: climate-machine/Land


text:
"""
    grdview(cmd0::String="", arg1=[], arg2=[], arg3=[]; kwargs...)

Reads a 2-D grid file and produces a 3-D perspective plot by drawing a mesh, painting a
colored/grayshaded surface made up of polygons, or by scanline conversion of these polygons
to a raster image.

Full option list at [`grdview`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html)

- $(GMT.opt_J)
- $(GMT.opt_Jz)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **G** : **drapefile** : -- Str or GMTgrid or a Tuple with 3 GMTgrid types --

    Drape the image in drapefile on top of the relief provided by relief_file.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#g)
- **I** : **shade** : **intensity** : **intensfileintens** : -- Str or GMTgrid --

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#i)
- **N** : **plane** : -- Str or Int --

    Draws a plane at this z-level.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#n)
- $(GMT.opt_P)
- **Q** : **type** : -- Str or Int --

    Specify **m** for mesh plot, **s* for surface, **i** for image.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#q)
- **S** : **smooth** : -- Number --

    Smooth the contours before plotting.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#s)
- **T** : **no_interp** : -- Str --

    Plot image without any interpolation.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#t)
- **W** : **contour** : **mesh** : **facade** : -- Str --

    Draw contour, mesh or facade. Append pen attributes.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/grdview.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function grdview(cmd0::String="", arg1=[], arg2=[], arg3=[], arg4=[], arg5=[], arg6=[];
                 K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdview", cmd0, arg1)

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	cmd, opt_B, = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd, = parse_JZ(cmd, d)
	cmd  = parse_UVXY(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_n(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	cmd = add_opt(cmd, 'N', d, [:N :plane])
	cmd = add_opt(cmd, 'Q', d, [:Q :type])
	cmd = add_opt(cmd, 'S', d, [:S :smooth])
	cmd = add_opt(cmd, 'T', d, [:T :no_interp])
	cmd = add_opt(cmd, 'W', d, [:W])
	cmd = add_opt(cmd, "Wc", d, [:contour])
	cmd = add_opt(cmd, "Wm", d, [:mesh])
	cmd = add_opt(cmd, "Wf", d, [:facade])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)		# Find how data was transmitted

	if (isa(arg1, Array{<:Number}))
		arg1 = mat2grid(arg1)
	end

	N_used = got_fname == 0 ? 1 : 0				# To know whether a cpt will go to arg1 or arg2
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_used, arg1, arg2)

	for sym in [:I :shade :intensity :intensfile]
		if (haskey(d, sym))
			if (!isa(d[sym], GMTgrid))			# Uff, simple. Either a file name or a -A type modifier
				cmd = cmd * " -I" * arg2str(d[sym])
			else
				cmd, N_used = put_in_slot(cmd, d[sym], 'I', [arg1, arg2, arg3])
				if (N_used == 1)     arg1 = d[sym]
				elseif (N_used == 2) arg2 = d[sym]
				elseif (N_used == 3) arg3 = d[sym]
				end
			end
			break
		end
	end

	for sym in [:G :drapefile]
		if (haskey(d, sym))
			if (isa(d[sym], String))				# Uff, simple. Either a file name or a -A type modifier
				cmd = cmd * " -G" * d[sym]
			elseif (isa(d[sym], GMTgrid))			# A single drape grid (arg1-3 may be used already)
				cmd, N_used = put_in_slot(cmd, d[sym], 'G', [arg1, arg2, arg3, arg4])
				if (N_used == 1)     arg1 = d[sym]
				elseif (N_used == 2) arg2 = d[sym]
				elseif (N_used == 3) arg3 = d[sym]
				elseif (N_used == 4) arg4 = d[sym]
				end
			elseif (isa(d[sym], Tuple) && length(d[sym]) == 3)
				cmd, N_used = put_in_slot(cmd, d[sym][1], 'G', [arg1, arg2, arg3, arg4, arg5, arg6])
				cmd = cmd * " -G -G"	# Because the above only set one -G and we need 3
				if (N_used == 1)      arg1 = d[sym][1];	arg2 = d[sym][2];		arg3 = d[sym][3]
				elseif (N_used == 2)  arg2 = d[sym][1];	arg3 = d[sym][2];		arg4 = d[sym][3]
				elseif (N_used == 3)  arg3 = d[sym][1];	arg4 = d[sym][2];		arg5 = d[sym][3]
				elseif (N_used == 4)  arg4 = d[sym][1];	arg5 = d[sym][2];		arg6 = d[sym][3]
				end
			else
				error("Wrong way of setting the drape (G) option.")
			end
			break
		end
	end

	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "grdview", arg1, arg2, arg3, arg4, arg5, arg6)
end

# ---------------------------------------------------------------------------------------------------
grdview!(cmd0::String="", arg1=[], arg2=[], arg3=[], arg4=[], arg5=[], arg6=[];
        K=true, O=true, first=false, kw...) =
	grdview(cmd0, arg1, arg2, arg3, arg4, arg5, arg6; K=true, O=true, first=false, kw...)

grdview(arg1, cmd0::String="", arg2=[], arg3=[], arg4=[], arg5=[], arg6=[];
        K=false, O=false, first=true, kw...) =
	grdview(cmd0, arg1, arg2, arg3, arg4, arg5, arg6; K=K, O=O, first=first, kw...)

grdview!(arg1, cmd0::String="", arg2=[], arg3=[], arg4=[], arg5=[], arg6=[];
        K=true, O=true, first=false, kw...) =
	grdview(cmd0, arg1, arg2, arg3, arg4, arg5, arg6; K=true, O=true, first=false, kw...)