text:
require("bspline.jl")

function exercise_zero()
    controlPoints = [
                     [-0.2 2];
                     [-0.3 6.2];
                     [-1.2 4.8];
                     [-2.8 8.8];
                     [-0.7 14];
                     [1.4 14.7];
                     [3.6 10.2];
                     [3.2 5.1];
                     [1.5 6.2];
                     [1.4 2];
                     ]

    ## continuityVector = ones(length(controlPoints[:,1]))
    k = 4
    extendedPartition = extendedPartitionAttempt(k, ones(6), Clumped())
    M = sum(ones(6))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition, k, M, controlPoints, matrix, fittingPoints)

    writeArrayForGnuplot(controlPoints,
                         "exercise-zero-clumped-control-poly.coordinates")

    writeArrayForGnuplot(bspline,
                         "exercise-zero-clumped-bspline.coordinates")

    extendedPartition = extendedPartitionAttempt(k, ones(6), Uniformed())
    M = sum(ones(6))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition, k, M, controlPoints, matrix, fittingPoints)

    writeArrayForGnuplot(controlPoints,
                         "exercise-zero-uniformed-control-poly.coordinates")

    writeArrayForGnuplot(bspline,
                         "exercise-zero-uniformed-bspline.coordinates")

    closedControlPoints = [controlPoints; controlPoints[1,:];]
    closedBSpline = drawBSpline(k, closedControlPoints, 200)
    writeArrayForGnuplot(closedControlPoints,
                         "exercise-zero-closed-control-poly.coordinates")

    writeArrayForGnuplot(closedBSpline,
                         "exercise-zero-closed-bspline.coordinates")

end

function exercise_one()
    controlPoints = [
                     [-2.5 -1];
                     [-3 2];
                     [0.3 2.0];
                     [-0.3 -2.0];
                     [4 -2.0];
                     [3 1.0];
                     ]

    writeArrayForGnuplot(controlPoints,
                         "exercise-one-control-poly.coordinates")

    fittingPoints = 200

    k = 2
    continuityVector = ones(4)
    k_two_extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, k_two_extendedPartition)
    k_two_bspline = drawCurve(k_two_extendedPartition,
                              k,
                              M,
                              controlPoints,
                              matrix,
                              fittingPoints)
    writeArrayForGnuplot(k_two_bspline,
                         "exercise-one-clumped-bspline-k-two.coordinates")

    k = 3
    continuityVector = ones(3)
    k_three_extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, k_three_extendedPartition)
    k_three_bspline = drawCurve(k_three_extendedPartition,
                              k,
                              M,
                              controlPoints,
                              matrix,
                              fittingPoints)
    writeArrayForGnuplot(k_three_bspline,
                         "exercise-one-clumped-bspline-k-three.coordinates")

    k = 4
    continuityVector = ones(2)
    k_four_extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, k_four_extendedPartition)
    k_four_bspline = drawCurve(k_four_extendedPartition,
                              k,
                              M,
                              controlPoints,
                              matrix,
                              fittingPoints)
    writeArrayForGnuplot(k_four_bspline,
                         "exercise-one-clumped-bspline-k-four.coordinates")

    k = 5
    continuityVector = ones(1)
    k_five_extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, k_five_extendedPartition)
    k_five_bspline = drawCurve(k_five_extendedPartition,
                              k,
                              M,
                              controlPoints,
                              matrix,
                              fittingPoints)
    writeArrayForGnuplot(k_five_bspline,
                         "exercise-one-clumped-bspline-k-five.coordinates")

    k = 6
    continuityVector = ones(0)
    k_six_extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, k_six_extendedPartition)
    k_six_bspline = drawCurve(k_six_extendedPartition,
                              k,
                              M,
                              controlPoints,
                              matrix,
                              fittingPoints)
    writeArrayForGnuplot(k_six_bspline,
                         "exercise-one-clumped-bspline-k-six.coordinates")

end

function exercise_two()
    controlPoints = [
                     [-0.7 -0.4];
                     [-1 -1];
                     [-1.3 -0.3];
                     [-1 1.2];
                     [-0.3 1.4];
                     [-0.15 0];
                     [0.15 0];
                     [0.3 1.6];
                     [2 0.6];
                     [2.6 0.2];
                     ]

    writeArrayForGnuplot(controlPoints,
                         "exercise-two-control-poly.coordinates")

    k = 4
    extendedPartition = expandPartition(0:10, [4 [1 for _ in 1:10]'])
    extendedPartition = normalizeRespect(extendedPartition, Maximum(), k)
    M = sum(ones(6))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)

    writeArrayForGnuplot(bspline,
                         "exercise-two-bspline-first.coordinates")

    k = 4
    extendedPartition = expandPartition(0:4, [4 2 2 2 4])
    extendedPartition = normalizeRespect(extendedPartition, Maximum())
    M = sum(ones(6))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition, k, M, controlPoints, matrix, fittingPoints)

    writeArrayForGnuplot(bspline,
                         "exercise-two-bspline-second.coordinates")

    k = 6
    extendedPartition = expandPartition(0:5, [6 1 1 1 1 6])
    extendedPartition = normalizeRespect(extendedPartition, Maximum())
    M = sum(ones(4))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition, k, M, controlPoints, matrix, fittingPoints)

    writeArrayForGnuplot(bspline,
                         "exercise-two-bspline-three.coordinates")

    k = 6
    extendedPartition = expandPartition(0:4, [6 1 1 2 6])
    extendedPartition = normalizeRespect(extendedPartition, Maximum())
    M = sum(ones(4))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition, k, M, controlPoints, matrix, fittingPoints)

    writeArrayForGnuplot(bspline,
                         "exercise-two-bspline-four.coordinates")

    k = 8
    extendedPartition = expandPartition(0:3, [8 1 1 8])
    extendedPartition = normalizeRespect(extendedPartition, Maximum())
    M = sum(ones(2))
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    fittingPoints = 200
    bspline = drawCurve(extendedPartition, k, M, controlPoints, matrix, fittingPoints)

    writeArrayForGnuplot(bspline,
                         "exercise-two-bspline-five.coordinates")

end

function exercise_three()
    controlPoints = [
                     [-3.0 2];
                     [-1.0 -2];
                     [1 1.0];
                     [1 1.0];
                     [3.0 -2];
                     [5 3.0];
                     ]

    writeArrayForGnuplot(controlPoints,
                         "exercise-three-control-poly.coordinates")

    fittingPoints = 200

    k = 3
    continuityVector = ones(3)
    extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)
    writeArrayForGnuplot(bspline,
                         "exercise-three-bspline.coordinates")
end

function exercise_four()
    controlPoints = [
                     [-2.0 2];
                     [-3.0 4];
                     [0 6];
                     [0 6];
                     [3 4];
                     [2.0 2];
                     ]

    writeArrayForGnuplot(controlPoints,
                         "exercise-four-control-poly.coordinates")

    fittingPoints = 200

    k = 4
    ## continuityVector = ones(3)
    ## extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    extendedPartition = [0 0 0 0 1/4 3/4 1 1 1 1]
    M = 2
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)
    writeArrayForGnuplot(bspline,
                         "exercise-four-bspline-first-partition.coordinates")

    extendedPartition = [0 0 0 0 1/2 1/2 1 1 1 1]
    M = 2
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)
    writeArrayForGnuplot(bspline,
                         "exercise-four-bspline-second-partition.coordinates")

end

function exercise_five()
    controlPoints = [
                     [-2.0 1];
                     [-1.0 -1];
                     [0 1];
                     [1 -1];
                     [2.0 1];
                     ]

    writeArrayForGnuplot(controlPoints,
                         "exercise-five-control-poly.coordinates")

    fittingPoints = 200

    k = 4
    continuityVector = ones(1)
    extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)
    writeArrayForGnuplot(bspline,
                         "exercise-five-bspline-single.coordinates")

    controlPoints = [
                     [-2.0 1];
                     [-1.0 -1];
                     [0 1];
                     [0 1];
                     [1 -1];
                     [2.0 1];
                     ]

    fittingPoints = 200

    k = 4
    continuityVector = ones(2)
    extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)
    writeArrayForGnuplot(bspline,
                         "exercise-five-bspline-doubled.coordinates")

    controlPoints = [
                     [-2.0 1];
                     [-1.0 -1];
                     [0 1];
                     [0 1];
                     [0 1];
                     [1 -1];
                     [2.0 1];
                     ]

    fittingPoints = 200

    k = 4
    continuityVector = ones(3)
    extendedPartition = extendedPartitionAttempt(k, continuityVector, Clumped())
    M = sum(continuityVector)
    matrix = buildFunctionsMatrix(k, M, extendedPartition)
    bspline = drawCurve(extendedPartition,
                        k,
                        M,
                        controlPoints,
                        matrix,
                        fittingPoints)
    writeArrayForGnuplot(bspline,
                         "exercise-five-bspline-tripled.coordinates")

end

function exercise_six()

    k = 4

    firstClosedControlPoints = [
                                [1.0 4];
                                [-0.8 4];
                                [-1.9 8.2];
                                [0.3 13.2];
                                [2.8 4];
                                [1.0 4];
                                ]

    firstClosedBSpline = drawBSpline(k, firstClosedControlPoints, 200)
    writeArrayForGnuplot(firstClosedControlPoints,
                         "exercise-six-first-closed-control-poly.coordinates")

    writeArrayForGnuplot(firstClosedBSpline,
                         "exercise-six-first-closed-bspline.coordinates")

    secondClosedControlPoints = [
                                 [0.4 6];
                                 [-1.8 3];
                                 [-2 14.8];
                                 [0.4 11];
                                 [3.0 14];
                                 [3.4 2.2];
                                 [0.4 6];
                                 ]

    secondClosedBSpline = drawBSpline(k, secondClosedControlPoints, 200)
    writeArrayForGnuplot(secondClosedControlPoints,
                         "exercise-six-second-closed-control-poly.coordinates")

    writeArrayForGnuplot(secondClosedBSpline,
                         "exercise-six-second-closed-bspline.coordinates")

end

function exercise_knotsInsertion()
    controlPoints = [
                     [0.0 0.0];
                     [1.0 2.0];
                     [2.0 2.5];
                     [3.0 2.2];
                     [4.0 0.5];
                     ]

    partition = [-1; 0.5; 1.8; 2.8; 3.2; 4;]

    k = 3
    refined_controlPoints =
        knotsInsertion_head(controlPoints, partition, k, 1)

    writeArrayForGnuplot(controlPoints,
                         "exercise-knots-insertion-original-control-points.coordinates")

    writeArrayForGnuplot(refined_controlPoints,
                         "exercise-knots-insertion.coordinates")

    ## deBoorBSpline = drawBSpline(k, controlPoints, 200)
    ## writeArrayForGnuplot(deBoorBSpline,
    ##                      "exercise-knots-insertion-deBoor.coordinates")

end

max_stars_repo_name: massimo-nocentini/cagd

max_issues_repo_name: massimo-nocentini/cagd

max_forks_repo_name: massimo-nocentini/cagd


text:
abstract type AbstractActivityMethod end

struct IntervalArithmetic <: AbstractActivityMethod end

struct LinearProgrammingRelaxation <: AbstractActivityMethod
    opt_factory::Any
end

struct FullFormulation <: AbstractActivityMethod
    opt_factory::Any
end

struct PartiallyEnumerateDisjunction{T<:Real} <: AbstractActivityMethod
    disjunction_ci::DisjunctionCI{T}
    fixings::Vector{Int}
    sub_method::AbstractActivityMethod
end

function minimum_activity(
    method::AbstractActivityMethod,
    model::MOI.ModelLike,
    x,
)
    return -maximum_activity(method, model, -x)
end

function maximum_activity(
    method::AbstractActivityMethod,
    model::MOI.ModelLike,
    x::T,
) where {T<:Real}
    return maximum_activity(
        method,
        model,
        convert(MOI.ScalarAffineFunction{T}, x),
    )
end

function maximum_activity(
    method::IntervalArithmetic,
    model::MOI.ModelLike,
    aff::MOI.ScalarAffineFunction{T},
)::T where {T<:Real}
    val = aff.constant
    for term in aff.terms
        lb, ub = MOIU.get_bounds(model, T, term.variable_index)
        coeff = term.coefficient
        if coeff > 0
            val += coeff * ub
        elseif coeff < 0
            val += coeff * lb
        end
    end
    return val
end

function _is_linear_constraint(::MOI.ConstraintIndex{F,S}) where {F,S}
    if !(
        F <: Union{
            MOI.SingleVariable,
            MOI.ScalarAffineFunction,
            MOI.VectorOfVariables,
            MOI.VectorAffineFunction,
        }
    )
        return false
    end
    if !(
        S <: Union{
            MOI.EqualTo,
            MOI.LessThan,
            MOI.GreaterThan,
            MOI.Interval,
            MOI.Zeros,
            MOI.Nonnegatives,
            MOI.Nonpositives,
        }
    )
        return false
    end
    return true
end

_get_constraint_filter(::LinearProgrammingRelaxation) = _is_linear_constraint
_get_constraint_filter(::FullFormulation) = nothing

function maximum_activity(
    method::Union{LinearProgrammingRelaxation,FullFormulation},
    model::MOI.ModelLike,
    aff::MOI.ScalarAffineFunction{T},
)::T where {T<:Real}
    # TODO: This extra copy is to handle optimizers that implement their own
    #       MOI.copy_to, and so crash if you try to call MOI.automatic_copy_to.
    #       We really should implement our own copy_to variant.
    placeholder_opt = MOIU.Model{T}()
    MOIU.automatic_copy_to(
        placeholder_opt,
        model,
        copy_names = false,
        filter_constraints = _get_constraint_filter(method),
    )
    MOI.set(
        placeholder_opt,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        aff,
    )
    MOI.set(placeholder_opt, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    opt = method.opt_factory()
    # TODO: What happens to ZeroOne constraints? Should get relaxed to [0,1]?
    MOI.copy_to(opt, placeholder_opt, copy_names = false)
    MOI.optimize!(opt)
    term = MOI.get(opt, MOI.TerminationStatus())
    if term == MOI.INFEASIBLE
        return T(-Inf)
    elseif term == MOI.DUAL_INFEASIBLE
        return T(Inf)
    elseif term == MOI.OPTIMAL
        # TODO: Remove this assertion
        @assert MOI.get(opt, MOI.ObjectiveValue()) <=
                maximum_activity(IntervalArithmetic(), model, aff)
        return T(MOI.get(opt, MOI.ObjectiveValue()))
    else
        @info "Unusual solution status $term; returning valid dual bound."
        return T(MOI.get(opt, MOI.ObjectiveBound()))
    end
end

function maximum_activity(
    method::PartiallyEnumerateDisjunction,
    model::MOI.ModelLike,
    aff::MOI.ScalarAffineFunction{T},
)::T where {T<:Real}
    vector_f = MOI.get(model, MOI.ConstraintFunction(), method.disjunction_ci)
    vector_s = MOI.get(model, MOI.ConstraintSet(), method.disjunction_ci)
    max_val = T(-Inf)
    for i in method.fixings
        # TODO: Copy once and change constraint right-hand sides. However, this
        # only works for non-ranged constraints whose sense is the same for
        # each alternative.
        constrained_model = MOIU.Model{T}()
        MOI.copy_to(constrained_model, model, copy_names = false)
        for (j, scalar_f) in enumerate(MOIU.scalarize(vector_f))
            lb = vector_s.lbs[j, i]
            ub = vector_s.ubs[j, i]
            scalar_s = MOI.Interval(lb, ub)
            MOI.add_constraint(constrained_model, scalar_f, scalar_s)
        end
        max_val = max(
            max_val,
            maximum_activity(constrained_model, aff, method.sub_method),
        )
    end
    return max_val
end

ext: jl


text:
module SparseMatrixVBCs

using SparseArrays
using ChainPartitioners
using SIMD
using Hwloc
using BenchmarkTools
using Statistics
using LinearAlgebra
using MemoizedMethods
using DiskCaches
using Scratch
using SHA
using Base.Threads

using CpuId
const DEFAULT_SIMD_SIZE = CpuId.simdbytes()

export SparseMatrix1DVBC
export model_SparseMatrix1DVBC_blocks
export model_SparseMatrix1DVBC_memory
export model_SparseMatrix1DVBC_TrSpMV_time
export model_SparseMatrix1DVBC_TrSpMV_fancy
export SparseMatrixVBC
export model_SparseMatrixVBC_blocks
export model_SparseMatrixVBC_memory
export model_SparseMatrixVBC_TrSpMV_time
export model_SparseMatrixVBC_TrSpMV_fancy
export TrSpMV!

"""
    SparseMatrix1DVBC{W, Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
Matrix type for storing sparse matrices in the
One-Dimensional Variable Block Column format.
"""
struct SparseMatrix1DVBC{W, Tv, Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Int
    n::Int
    Φ::SplitPartition{Ti}
    pos::Vector{Ti}
    idx::Vector{Ti}
    ofs::Vector{Ti}
    val::Vector{Tv}
    function SparseMatrix1DVBC{W, Tv, Ti}(m::Integer, n::Integer, Φ::SplitPartition{Ti}, pos::Vector{Ti}, idx::Vector{Ti}, ofs::Vector{Ti}, val::Vector{Tv}) where {W, Tv, Ti<:Integer}
        @noinline throwsz(str, lbl, K) =
            throw(ArgumentError("number of $str ($lbl) must be ≥ 0, got $K"))
        m < 0 && throwsz("rows", 'm', m)
        n < 0 && throwsz("columns", 'n', n)
        W isa Int || throw(ArgumentError("W must be an Int"))
        W > 0 || throw(ArgumentError("W must be > 0"))
        new(m, n, Φ, pos, idx, ofs, val)
    end
end

Base.size(A::SparseMatrix1DVBC) = (A.m, A.n)

"""
    SparseMatrixVBC{W, Tv, Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
Matrix type for storing sparse matrices in the
Variable Block Column format.
"""
struct SparseMatrixVBC{U, W, Tv, Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Int
    n::Int
    Π::SplitPartition{Ti}
    Φ::SplitPartition{Ti}
    pos::Vector{Ti}
    idx::Vector{Ti}
    ofs::Vector{Ti}
    val::Vector{Tv}
    function SparseMatrixVBC{U, W, Tv, Ti}(m::Integer, n::Integer, Π::SplitPartition{Ti}, Φ::SplitPartition{Ti}, pos::Vector{Ti}, idx::Vector{Ti}, ofs::Vector{Ti}, val::Vector{Tv}) where {U, W, Tv, Ti<:Integer}
        @noinline throwsz(str, lbl, x) =
            throw(ArgumentError("number of $str ($lbl) must be ≥ 0, got $x"))
        m < 0 && throwsz("rows", 'm', m)
        n < 0 && throwsz("columns", 'n', n)
        U isa Int || throw(ArgumentError("U must be an Int"))
        W isa Int || throw(ArgumentError("W must be an Int"))
        U > 0 || throw(ArgumentError("U must be > 0"))
        W > 0 || throw(ArgumentError("W must be > 0"))
        new(m, n, Π, Φ, pos, idx, ofs, val)
    end
end

Base.size(A::SparseMatrixVBC) = (A.m, A.n)

include("util.jl")
include("TrSpMV.jl")
include("costs.jl")
include("constructors_1DVBC.jl")
include("constructors_VBC.jl")
include("hacks.jl")

end # module


text:
# comparing iSARA policy to canonical acquisition functions
using Plots
using HDF5
using Statistics
function plot_kernel_benchmark(f)
    # error = read(f, "mean_squared_error")
    # error = read(f, "rmse")
    error = read(f, "r2")
    μ = mean(error, dims = (2, 5)) # average over independent experiments and stripes
    σ = std(error, dims = (2, 5))
    nsample, nexp, nlength, nkernel = read(f, "nsample"), read(f, "nexp"), read(f, "nlength"), read(f, "nkernel")
    μ = reshape(μ, nsample, nkernel, nlength)
    σ = reshape(σ, nsample, nkernel, nlength)
    # for each kernel, use length scale with best performance
    lind = zeros(Int, nkernel)
    for i in 1:nkernel
        A = @view μ[end, i, :]
        # lind[i] = argmin(A)
        lind[i] = argmax(A)
    end
    A = zeros(nsample, nkernel)
    B = zeros(nsample, nkernel)
    for i in 1:nkernel
        @. A[:, i] = μ[:, i, lind[i]]
        @. B[:, i] = σ[:, i, lind[i]]
    end
    n = prod(size.((error,), (2, 5)))
    confidence = maximum(B) / sqrt(n)
    println("confidence = $confidence")
    lengthscales = read(f, "lengthscales")
    println("lengthscales = $(lengthscales[lind])")
    names = read(f, "kernel names")
    plot(xlabel = "# of points", ylabel = "R2") #i, yscale = :log10)
    plot!(A, label = permutedims(names), linewidth = 2) #, linestyle = :auto)
    gui()
end

path = "SARA/NatCom2020/results/"
file = "inner_kernel_benchmark_short.h5"
# file = "inner_kernel_benchmark_short_endpoints.h5"
f = h5open(path * file)

plotly()
# plotlyjs()
# pyplot()
# Plots.scalefontsizes(1.5)
plot_kernel_benchmark(f)
# savefig("inner_kernel_plot.pdf")


text:
#=

   From https://reference.wolfram.com/language/ref/HypergeometricDistribution.html
   """
   Suppose an urn has 100 elements, of which 40 are special.
   ...
   Compute the probability that there are more than 25 special elements in a
   draw of 50 elements.
     Answer: 0.0120902

   Compute the expected number of special elements in a draw of 50 elements.
     Answer: 20
   """

   mean(sumSpecialOf50): 19.9992

   Distributions of variable sumSpecialOf50 (num:0)
   19.00000 =>    1210  (0.121000)
   20.00000 =>    1139  (0.113900)
   21.00000 =>    1118  (0.111800)
   18.00000 =>     972  (0.097200)
   22.00000 =>     902  (0.090200)
   17.00000 =>     794  (0.079400)
   23.00000 =>     732  (0.073200)
   24.00000 =>     619  (0.061900)
   16.00000 =>     596  (0.059600)
   15.00000 =>     429  (0.042900)
   25.00000 =>     384  (0.038400)
   26.00000 =>     256  (0.025600)
   14.00000 =>     255  (0.025500)
   13.00000 =>     151  (0.015100)
   27.00000 =>     151  (0.015100)
   28.00000 =>     102  (0.010200)
   12.00000 =>      72  (0.007200)
   29.00000 =>      43  (0.004300)
   11.00000 =>      28  (0.002800)
   30.00000 =>      21  (0.002100)
   10.00000 =>      13  (0.001300)
   31.00000 =>       7  (0.000700)
   35.00000 =>       1  (0.000100)
   8.00000 =>       1  (0.000100)
   32.00000 =>       1  (0.000100)
   9.00000 =>       1  (0.000100)
   7.00000 =>       1  (0.000100)
   34.00000 =>       1  (0.000100)

   moreThan25SpecialOf50): 0.0583

   Cf ~/blog/urn_model1.blog
      ~/webppl/urn_model1.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function urn_model1(n=100,NumDraws=50,NumSpecial=40)

    nonspecial = 1
    special = 2

    element ~ filldist(Categorical(simplex([n-NumSpecial,NumSpecial])),n) # [nonspecial,special]

    # We have (/observe) exactly 40 special elements (no random there!)
    numSpecial ~ Dirac(sum([element[i] == special for i in 1:n]))

    # Compute the expected number of special elements in a draw of 50 elements.
    sumSpecialOf50 ~ Dirac(sum([element[i] == special ? 1 : 0 for i in 1:NumDraws]))

    # What's the probability that there are more than 25 special elements in a draw of 50 elements
    moreThan25SpecialOf50 ~ Dirac(sumSpecialOf50 > 25)

end

model = urn_model1()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# display(chns)
# display(plot(chns))

println("mean(sumSpecialOf50): ", mean_val(chns,:sumSpecialOf50))
show_var_dist_pct(chns, :sumSpecialOf50)
println("moreThan25SpecialOf50): ", mean_val(chns,:moreThan25SpecialOf50))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using LinearAlgebra
"Sigmoid function"
sigm(u) = 1/(1 + exp(-u))

"Convert an n-vector to an nx1 array"
matrix(x::Array{T,1}) where {T<:Any} = reshape(x, length(x),1)
matrix(x::Array{T,2}) where {T<:Any} = x
matrix(x::AbstractRange) = reshape(collect(x), length(x), 1)
matrix(X::Array{Any,2}) = X
matrix(x::Adjoint) = reshape(transpose(x), 1, length(x))

"Apply f to each row of a matrix. f should map vectors to vectors"
function rowwise(f, X::Array{T,2}) where {T<:Any}
    y1 = f(X[1,:])
    n = size(X,1)
    d = length(y1)
    Y = Array{Any}(n,d)
    Y[1,:] = y1
    for i=2:n
        Y[i,:] = f(X[i,:])
    end
    return Y
end

"compute the rms of a matrix or a vector"
rms(x::Array{Float64,1}) = sqrt.( (1/length(x)) * (x'*x))
rms(x::Array{Float64,2}) = rms(x[:])

"`findvalue(s, lst)`: Find the first index of `s` in a list `lst`"
findvalue(s, lst) = findfirst(isequal(s), lst)

number(x::Number) = convert(Float64,x)
number(x::String) = parse(Float64,x)


text:
using Merly

global u
u="hello"

server = Merly.app()
server.use("CORS")

@page "/" "Hello World!"
@page "/hola/:usr" "<b>Hello {{usr}}!</b>"

@route GET "/get/:data" begin
  "get this back: {{data}}"
end

@route POST "/post" begin
  "I did something!"
end

@route POST|PUT|DELETE "/" begin
  println("params: ",q.params)
  println("query: ",q.query)
  println("body: ",q.body)

  r.headers["Content-Type"]="text/plain"

  "I did something!"
end

Get("/data", (q,r)->(begin
  r.headers["Content-Type"]="text/plain"
  "$u data"
end))

Post("/data", (q,r)->(begin
  println("params: ",q.params)
  println("query: ",q.query)
  println("body: ",q.body)
  r.headers["Content-Type"]="text/plain"
  global u="bye"
  "I did something!"
end))

server.start("localhost", 8080)


text:
using Jive # @If
@If VERSION >= v"1.1.0-DEV.764" module test_julia_locals

using Test

output = let x = 1
    sprint(show, Base.@locals)
end

if VERSION >= v"1.6.0-DEV.849"
    @test output == "Dict{Symbol, Any}(:x => 1)"
end

end # @If VERSION >= v"1.1.0-DEV.764" module test_julia_locals


text:
@testset "fileio" begin
    tmp_file = tempname() * ".sixel"
    for img in (
        repeat(Gray.(0:0.1:0.9), inner=(10, 50)),
        repeat(RGB.(0:0.1:0.9), inner=(10, 50))
    )
        save(tmp_file, img)
        img_readback = load(tmp_file)
        @test eltype(img_readback) == RGB{N0f8}
        @test img_readback isa IndirectArray

        @test assess_psnr(img, eltype(img).(img_readback)) > 30
    end
end


text:
function get_fit_window(counts, risefactor, tailfactor)

    maxval = maximum(counts)
    ind1::Int = findfirst(x -> x > maxval*risefactor, counts)
    ind2::Int =  findlast(x -> x >maxval*tailfactor, counts)

    return ind1, ind2
end

function pads(u::Array{Float64,1})
    outsize = 2 .* length(u) .- 1
    out = zeros(eltype(u), nextfastfft(outsize))
    copyto!(out, u)
end

function compute_plans(IRFpad)
    pl = plan_rfft(IRFpad)
    fftIRF = pl * IRFpad
    pli = plan_irfft(fftIRF, length(IRFpad))

    return pl, fftIRF, pli
end

function _conv_DT(t, β::Array{Float64,1}, data::fitDTOF, out, fftIRF, pl, pli, mod::Val{TPSF_DA_semiinf_refl};
    RtDT = Array{Float64}(undef, length(data.IRF))
    )

        RtDT = TPSF_DA_semiinf_refl(data.t, β, data.ρ, data.ndet, data.nmed)
        RtDT = RtDT./maximum(RtDT)
        copyto!(out, RtDT)

        fftRtDT = pl * out
        fftRtDT .*= fftIRF
        out = pli * fftRtDT

        out = out./maximum(out)

	    tidx = findfirst(x -> x == t[1], data.t)

	    out = out[tidx:tidx+length(t)-1]
	return log.(out)
end

function _conv_DT(t, β::Array{Float64,1}, data::fitDTOF, out, fftIRF, pl, pli, mod::Val{TPSF_DA_slab_refl};
    RtDT = Array{Float64}(undef, length(data.IRF))
    )

        RtDT = TPSF_DA_slab_refl(data.t, β, data.ρ, data.ndet, data.nmed, data.s)
        RtDT = RtDT./maximum(RtDT)
        copyto!(out, RtDT)

        fftRtDT = pl * out
        fftRtDT .*= fftIRF
        out = pli * fftRtDT

        out = out./maximum(out)

	    tidx = findfirst(x -> x == t[1], data.t)

	    out = out[tidx:tidx+length(t)-1]
	return log.(out)
end

function _conv_DT(t, β::Array{Float64,1}, data::fitDTOF, out, fftIRF, pl, pli, mod::Val{TPSF_DA_paralpip_refl};
    RtDT = Array{Float64}(undef, length(data.IRF))
    )

        RtDT = TPSF_DA_paralpip_refl(data.t, β, data.ndet, data.nmed, data.rd, data.rs, data.L)
        RtDT = RtDT./maximum(RtDT)
        copyto!(out, RtDT)

        fftRtDT = pl * out
        fftRtDT .*= fftIRF
        out = pli * fftRtDT

        out = out./maximum(out)

	    tidx = findfirst(x -> x == t[1], data.t)

	    out = out[tidx:tidx+length(t)-1]
	return log.(out)
end

function getfit(data::fitDTOF, model_params::DTOF_fitparams)
     ## get FFT padding and plans
     out = pads(data.IRF)

     pl, fftIRF, pli = compute_plans(out)

    #get fit windows
    ind1, ind2 = get_fit_window(data.DTOF, model_params.risefactor, model_params.tailfactor)

    fit = curve_fit((t, β) -> _conv_DT(t, β, data, out, fftIRF, pl, pli, Val(model_params.model)), data.t[ind1:ind2], log.(data.DTOF[ind1:ind2]),
        model_params.initparams, lower=model_params.lb, upper=model_params.ub)#; autodiff=:finiteforward)

    return fit, data.t[ind1:ind2], log.(data.DTOF[ind1:ind2]), _conv_DT(data.t[ind1:ind2], fit.param, data, out, fftIRF, pl, pli, Val(model_params.model))

end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
import MathOptInterface
const MOI = MathOptInterface

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Ptr{Cvoid}
    objective_sense::MOI.OptimizationSense
    variable_map::Dict{MOI.VariableIndex, String}
    objective_constant::Float64

    """
        Optimizer()

    Create a new Optimizer object.
    """
    function Optimizer()
        ptr = Highs_create()
        if ptr == C_NULL
            error("Unable to create an internal model via the C API.")
        end
        model = new(
            ptr,
            MOI.FEASIBILITY_SENSE,
            Dict{MOI.VariableIndex, String}(),
            0.0,
        )
        finalizer(Highs_destroy, model)
        return model
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, model::Optimizer) = model
Base.unsafe_convert(::Type{Ptr{Cvoid}}, model::Optimizer) = model.inner

function _check_ret(model::Optimizer, ret::Cint)
    # TODO: errors for invalid return codes!
    return nothing
end

function MOI.empty!(o::Optimizer)
    Highs_destroy(o)
    o.inner = Highs_create()
    o.objective_sense = MOI.FEASIBILITY_SENSE
    empty!(o.variable_map)
    o.objective_constant = 0.0
    return
end

function MOI.is_empty(o::Optimizer)
    return Highs_getNumRows(o) == 0 &&
      Highs_getNumCols(o) == 0 &&
      o.objective_sense == MOI.FEASIBILITY_SENSE &&
      o.objective_constant == 0.0
end

function MOI.optimize!(o::Optimizer)
    if o.objective_sense == MOI.FEASIBILITY_SENSE
        obj_func = MOI.get(o, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        if any(!=(0.0), obj_func.terms)
            error("Feasibility sense with non-constant objective function, set the sense to min/max, the objective to 0 or reset the sense to feasibility to erase the objective")
        end
    end
    Highs_run(o)
    return
end

function MOI.add_variable(o::Optimizer)
    _ = Highs_addCol(o, 0.0, -Inf, Inf, Cint(0), Cint[], Cint[])
    col_idx = MOI.get(o, MOI.NumberOfVariables()) - 1
    return MOI.VariableIndex(col_idx)
end

function MOI.add_constrained_variable(o::Optimizer, set::S) where {S <: MOI.Interval}
    _ = Highs_addCol(o, 0.0, Cdouble(set.lower), Cdouble(set.upper), Cint(0), Cint[], Cint[])
    col_idx = MOI.get(o, MOI.NumberOfVariables()) - 1
    return (MOI.VariableIndex(col_idx), MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}(col_idx))
end

MOI.supports_constraint(::Optimizer, ::MOI.SingleVariable, ::MOI.Interval) = true

function MOI.add_constraint(o::Optimizer, sg::MOI.SingleVariable, set::MOI.Interval)
    var_idx = Cint(sg.variable.value)
    _ = Highs_changeColBounds(o, var_idx, Cdouble(set.lower), Cdouble(set.upper))
    return
end

function MOI.get(o::Optimizer, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}})
    nrows = Highs_getNumRows(o)
    return Int(nrows)
end

MOI.supports_constraint(::Optimizer, ::MOI.ScalarAffineFunction, ::MOI.Interval) = true

function MOI.add_constraint(o::Optimizer, func::MOI.ScalarAffineFunction, set::MOI.Interval)
    number_nonzeros = length(func.terms)
    coefficients = Vector{Cdouble}(undef, number_nonzeros)
    col_indices = Vector{Cint}(undef, number_nonzeros)
    for j in Base.OneTo(number_nonzeros)
        term = func.terms[j]
        coefficients[j] = term.coefficient
        col_indices[j] = Cint(term.variable_index.value)
    end
    coefficients_ptr = pointer(coefficients)
    col_indices_ptr = pointer(col_indices)
    lower = convert(Cdouble, set.lower - func.constant)
    upper = convert(Cdouble, set.upper - func.constant)
    row_idx = Highs_addRow(o, lower, upper, Cint(number_nonzeros), col_indices_ptr, coefficients_ptr)
    return MOI.ConstraintIndex{typeof(func), typeof(set)}(row_idx)
end

MOI.supports(::Optimizer, ::MOI.SolverName) = true
MOI.get(::Optimizer, ::MOI.SolverName) = "HiGHS"

MOI.supports(::Optimizer, param::MOI.RawParameter) = true

# setting HiGHS options
function MOI.set(model::Optimizer, param::MOI.RawParameter, value::Integer)
    ret = Highs_setHighsIntOptionValue(model, param.name, Cint(value))
    return _check_ret(model, ret)
end

function MOI.set(model::Optimizer, param::MOI.RawParameter, value::Bool)
    ret = Highs_setHighsBoolOptionValue(model, param.name, Cint(value))
    return _check_ret(model, ret)
end

function MOI.set(model::Optimizer, param::MOI.RawParameter, value::AbstractFloat)
    ret = Highs_setHighsDoubleOptionValue(model, param.name, Cdouble(value))
    return _check_ret(model, ret)
end

function MOI.set(model::Optimizer, param::MOI.RawParameter, value::String)
    ret = Highs_setHighsStringOptionValue(model, param.name, value)
    return _check_ret(model, ret)
end

function _get_option(model::Optimizer, option::String, ::Type{Cint})
    value = Ref{Cint}(0)
    Highs_getHighsIntOptionValue(model, option, value)
    return value[]
end

function _get_option(model::Optimizer, option::String, ::Type{Bool})
    value = Ref{Cint}(0)
    Highs_getHighsBoolOptionValue(model, option, value)
    return Bool(value[])
end

function _get_option(model::Optimizer, option::String, ::Type{Cdouble})
    value = Ref{Cdouble}()
    Highs_getHighsDoubleOptionValue(model, option, value)
    return value[]
end

function _get_option(model::Optimizer, option::String, ::Type{String})
    buffer = Vector{Cchar}(undef, 100)
    bufferP = pointer(buffer)
    GC.@preserve buffer begin
        Highs_getHighsStringOptionValue(model, option, bufferP)
        return unsafe_string(bufferP)
    end
end

const _OPTIONS = Dict{String,DataType}(
    "presolve" => String,
    "solver" => String,
    "parallel" => String,
    "simplex_strategy" => Cint,
    "simplex_iteration_limit" => Cint,
    "highs_min_threads" => Cint,
    "message_level" => Cint,
    "time_limit" => Cdouble,
)

function MOI.get(o::Optimizer, param::MOI.RawParameter)
    param_type = get(_OPTIONS, param.name, nothing)
    if param_type === nothing
        throw(ArgumentError("Parameter $(param.name) is not supported"))
    end
    return _get_option(o, param.name, param_type)
end

const _SUPPORTED_MODEL_ATTRIBUTES = Union{
    MOI.ObjectiveSense,
    MOI.NumberOfVariables,
    MOI.ListOfVariableIndices,
    MOI.ListOfConstraintIndices,
    MOI.NumberOfConstraints, # TODO single variables
    MOI.ListOfConstraints, # TODO single variables
    MOI.ObjectiveFunctionType,
    MOI.ObjectiveValue,
    # MOI.DualObjectiveValue, # TODO
    # MOI.SolveTime,  # TODO
    MOI.SimplexIterations,
    MOI.BarrierIterations,
    MOI.RawSolver,
    # MOI.RawStatusString,  # TODO
    MOI.ResultCount,
    # MOI.Silent, # TODO
    MOI.TerminationStatus,
    # MOI.PrimalStatus,
    # MOI.DualStatus
}
@enum HighsModelStatus begin
    NOTSET = 0
    LOAD_ERROR
    MODEL_ERROR
    PRESOLVE_ERROR
    SOLVE_ERROR
    POSTSOLVE_ERROR
    MODEL_EMPTY
    PRIMAL_INFEASIBLE
    PRIMAL_UNBOUNDED
    OPTIMAL
    REACHED_DUAL_OBJECTIVE_VALUE_UPPER_BOUND
    REACHED_TIME_LIMIT
    REACHED_ITERATION_LIMIT
    PRIMAL_DUAL_INFEASIBLE
    DUAL_INFEASIBLE
end
MOI.supports(::Optimizer, ::_SUPPORTED_MODEL_ATTRIBUTES) = true

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

MOI.get(o::Optimizer, ::MOI.RawSolver) = o

function MOI.get(o::Optimizer, ::MOI.ResultCount)
    status = HighsModelStatus(Highs_getModelStatus(o, Cint(0)))
    return status == OPTIMAL ? 1 : 0
end
function MOI.get(o::Optimizer, ::MOI.ObjectiveSense)
    return o.objective_sense
end
function MOI.get(o::Optimizer, ::MOI.TerminationStatus)
    status = HighsModelStatus(Highs_getModelStatus(o.inner, Cint(0)))
    if status == OPTIMAL
        return MOI.OPTIMAL
    elseif status == PRIMAL_INFEASIBLE
        return MOI.INFEASIBLE
    elseif status == DUAL_INFEASIBLE
        return MOI.DUAL_INFEASIBLE
    elseif status == REACHED_ITERATION_LIMIT
        return MOI.ITERATION_LIMIT
    elseif status == REACHED_TIME_LIMIT
        return MOI.TIME_LIMIT
    elseif status == MODEL_ERROR
        return MOI.INVALID_MODEL
    else
        return MOI.OTHER_ERROR
    end
end

function MOI.set(o::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    sense_code = sense == MOI.MAX_SENSE ? Cint(-1) : Cint(1)
    _ = Highs_changeObjectiveSense(o, sense_code)
    o.objective_sense = sense
    # if feasibility sense set, erase the function
    if sense == MOI.FEASIBILITY_SENSE
        MOI.set(o, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction{Float64}([], 0))
    end
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.NumberOfVariables)
    return Int(Highs_getNumCols(o))
end

function MOI.get(o::Optimizer, ::MOI.ListOfVariableIndices)
    ncols = MOI.get(o, MOI.NumberOfVariables())
    return [MOI.VariableIndex(j) for j in 0:(ncols-1)]
end

MOI.get(::Optimizer, ::MOI.ObjectiveFunctionType) = MOI.ScalarAffineFunction{Float64}

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true

function MOI.set(o::Optimizer, ::MOI.ObjectiveFunction{F}, func::F) where {F <: MOI.ScalarAffineFunction{Float64}}
    total_ncols = MOI.get(o, MOI.NumberOfVariables())
    coefficients = zeros(Cdouble, total_ncols)
    for term in func.terms
        j = term.variable_index.value
        coefficients[j+1] = Cdouble(term.coefficient)
    end
    coefficients_ptr = pointer(coefficients)
    mask = pointer(ones(Cint, total_ncols))
    Highs_changeColsCostByMask(o, mask, coefficients_ptr)
    o.objective_constant = MOI.constant(func)
    return
end

function MOI.get(o::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}})
    ncols = MOI.get(o, MOI.NumberOfVariables())
    num_cols = Ref{Cint}(0)
    costs = Vector{Cdouble}(undef, ncols)
    if ncols > 0
        _ = Highs_getColsByRange(
            o,
            Cint(0), Cint(ncols-1), # column range
            num_cols, costs,
            C_NULL, C_NULL, # lower, upper
            Ref{Cint}(0), C_NULL, C_NULL, C_NULL # coefficients
        )
        num_cols[] == ncols || error("Unexpected number of columns, inconsistent HiGHS state")
    end
    terms = Vector{MOI.ScalarAffineTerm{Float64}}()
    for (j, cost) in enumerate(costs)
        if cost != 0.0
            var_idx = MOI.VariableIndex(j-1)
            push!(terms, MOI.ScalarAffineTerm(cost, var_idx))
        end
    end
    return MOI.ScalarAffineFunction{Float64}(terms, o.objective_constant)
end

function MOI.get(o::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(o, attr)
    value = Ref{Cdouble}()
    Highs_getHighsDoubleInfoValue(o, "objective_function_value", value);
    return o.objective_constant + value[]
end

function MOI.get(o::Optimizer, ::MOI.SimplexIterations)
    simplex_iteration_count = Ref{Cint}(0)
    Highs_getHighsIntInfoValue(o, "simplex_iteration_count", simplex_iteration_count)
    return Int(simplex_iteration_count[])
end

function MOI.get(o::Optimizer, ::MOI.BarrierIterations)
    return 0
end

function MOI.get(o::Optimizer, ::MOI.VariableName, v::MOI.VariableIndex)
    return get(o.variable_map, v, "")
end

function MOI.set(o::Optimizer, ::MOI.VariableName, v::MOI.VariableIndex, name::String)
    o.variable_map[v] = name
    return
end

function MOI.get(o::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    for (vi, vname) in o.variable_map
        if vname == name
            return vi
        end
    end
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(o, MOI.RawParameter("time_limit"))
end

function MOI.set(o::Optimizer, ::MOI.TimeLimitSec, value::Real)
    return MOI.set(o, MOI.RawParameter("time_limit"), Cdouble(value))
end

function MOI.set(o::Optimizer, ::MOI.TimeLimitSec, ::Nothing)
    # TODO(odow): handle default time limit?
    return
end


text:
module RoboDojo

using LinearAlgebra
using Symbolics
using SparseArrays
using BenchmarkTools
using MeshCat
using GeometryBasics
using Colors
using CoordinateTransformations
using Rotations
using Parameters
using JLD2
using Scratch

# Utilities
include("utils.jl")

# Solver
include("solver/lu.jl")
include("solver/cones.jl")
include("solver/indices.jl")
include("solver/interior_point.jl")

# Simulator
include(joinpath("../src/simulator/model.jl"))
include(joinpath("../src/simulator/dimensions.jl"))
include(joinpath("../src/simulator/indices.jl"))
include(joinpath("../src/simulator/disturbances.jl"))
include(joinpath("../src/simulator/trajectory.jl"))
include(joinpath("../src/simulator/policy.jl"))
include(joinpath("../src/simulator/simulator.jl"))
include(joinpath("../src/simulator/residual.jl"))
include(joinpath("../src/simulator/codegen.jl"))

# Robots
include(joinpath("../src/robots/visual_utils.jl"))

include(joinpath("../src/robots/hopper/model.jl"))
include(joinpath("../src/robots/hopper/simulator.jl"))
include(joinpath("../src/robots/hopper/visuals.jl"))

include(joinpath("../src/robots/biped/model.jl"))
include(joinpath("../src/robots/biped/visuals.jl"))

include(joinpath("../src/robots/quadruped/model.jl"))
include(joinpath("../src/robots/quadruped/visuals.jl"))

include(joinpath("../src/robots/box/model.jl"))
include(joinpath("../src/robots/box/visuals.jl"))

include(joinpath("../src/robots/integrator.jl"))
include(joinpath("../src/robots/codegen.jl"))

include(joinpath("../src/robots/load.jl"))

# load simulation environments

# RESIDUAL_EXPR = Dict{String, Any}()
# residual_expr(model::Model) = RESIDUAL_EXPR[String(name(model)) * "_r"]
# jacobian_var_expr(model::Model) = RESIDUAL_EXPR[String(name(model)) * "_rz"]
# jacobian_data_expr(model::Model) = RESIDUAL_EXPR[String(name(model)) * "_rθ"]

# model = hopper
# path = joinpath(@__DIR__, "robots", String(name(model)), "expr/expr.jld2")
# @load path r_model rz_model rθ_model
# RESIDUAL_EXPR[String(name(model)) * "_r"] = eval(r_model)
# RESIDUAL_EXPR[String(name(model)) * "_rz"] = eval(rz_model)
# RESIDUAL_EXPR[String(name(model)) * "_rθ"] = eval(rθ_model)

# model = biped
# path = joinpath(@__DIR__, "robots", String(name(model)), "expr/expr.jld2")
# @load path r_model rz_model rθ_model
# RESIDUAL_EXPR[String(name(model)) * "_r"] = eval(r_model)
# RESIDUAL_EXPR[String(name(model)) * "_rz"] = eval(rz_model)
# RESIDUAL_EXPR[String(name(model)) * "_rθ"] = eval(rθ_model)

# model = quadruped
# path = joinpath(@__DIR__, "robots", String(name(model)), "expr/expr.jld2")
# @load path r_model rz_model rθ_model
# RESIDUAL_EXPR[String(name(model)) * "_r"] = eval(r_model)
# RESIDUAL_EXPR[String(name(model)) * "_rz"] = eval(rz_model)
# RESIDUAL_EXPR[String(name(model)) * "_rθ"] = eval(rθ_model)

# Policy

export
    hopper, biped, quadruped,
    nominal_configuration,
    step!, simulate!, Simulator,
    visualize!, Visualizer,
    friction_coefficients,
    lagrangian,
    num_var, num_data,
    residual

end # module


text:
module ISCLib

using Stats, NIfTI, Logging, ExcelReaders, DataFrames

include("consts.jl")

export NSUBJ, SHAPE, SIZE, get_nii,
      get_niis, cov_n!, LEN,
      get_mean, triangular_index, ISC_res,
      get_indexes, NEXP

"""
# indexes for the end of response for each question
# i - Subject_id
"""

function get_indexes(subj::Int, t::Array{String,1} = ["S1","S2"], rb::Int64=1)
    data = readxl(DataFrame,joinpath(PATH,"logs","$subj.xlsx"), "Лист1!A1:G$(NEXP+1)");
    data = data[reduce((x,y)-> x .| y,(Array(data[:Stimul_1].== w) for w in t)),:]
    ind_end = data[:Response]
    # map(Int,(ceil(hcat(ind_end-minimum(data[Symbol("Stim dur")]),ind_end))*2)) #ceil returns float64
    map(Int,(ceil(hcat(ind_end-6,ind_end+2))*2))[1:end-1-rb,:]
end

function get_nii(subj::Int, ind::Int)
    fnames = readdir(PATH_Fmt(subj))[2:end] # first elem is  mean
    return niread(joinpath(PATH_Fmt(subj),fnames[ind]),mmap= true)
end

function get_niis(subj::Int, rng::AbstractArray=1:LEN)
    fnames = readdir(PATH_Fmt(subj))[2:end] # first elem is  mean
    return [niread(joinpath(PATH_Fmt(subj),fnames[ind]),mmap = true) for ind in rng]
end

@inline function rema!{T<: Real}(shadow::Array{T}, variable::Array{T}, decay = 0.9)
    shadow .-= (1 - decay) * (shadow - variable)
end

function get_mean(subj::Int)
  debug("mean",subj)
  if isfile(joinpath(PATH_Fmt(subj),"mean.nii"))
    return niread(joinpath(PATH_Fmt(subj),"mean.nii"))
  else
    seq = get_niis(subj,1:LEN)
    m = mean(seq)
    ni = NIfTI.NIVolume(seq[1].header, seq[1].extensions, m)
    niwrite(joinpath(PATH_Fmt(subj),"mean.nii"), ni)
  end
  m
end

function triangular_index(i::Int,j::Int)
  Int(round(j * (j - 3) / 2 + i + 1))
end

function reverse_tri_index(k::Int)
  j = round(floor(-.5 + .5 * sqrt(1 + 8 * (k - 1))) + 2);
  i = round(j * (3 - j) / 2 + k - 1);
  Int(i),Int(j)
end

function cov_n!{T<:Real}(covs::Array{T,2},ds::Array{T,2},seq::AbstractArray,means::AbstractArray,lim::Int)
  info("cov_n called")
  srol = [seq[i][1].raw - means[i] for i in 1:NSUBJ ]
  for i in 2:lim
      @simd for subj in 1:NSUBJ
          rema!(srol[subj],seq[subj][i].raw - means[subj])
          # in-place Фильтр Брауна по 1 аргументу ( x = δ x + (1-δ) χ, χ - new x )
          @inbounds ds[:,subj] .+= reshape(srol[subj].^2, SIZE)
      end
      @simd for k in 2:NSUBJ for j in 1:k-1
              # info(i,j,k,"  ",triangular_index(j,k))
              @inbounds covs[:,triangular_index(j,k)] .+= reshape(srol[j] .* srol[k], SIZE)
      end end
  end
end

function ISC_res{T<:Real}(covs::Array{T,2},ds::Array{T,2})
    res = zeros(SIZE)
    for i in 1:SIZE
      for k in 2:NSUBJ for j in 1:k-1
          res[i] += covs[i,triangular_index(j,k)]/sqrt(ds[i,j]*ds[i,k])
      end end
      res[i] /= (NSUBJ*(NSUBJ-1))/2
    end
    res
end

end


text:
# heavily modeled after CDK

# a list of functions that any widget would have
# At a minimum, a new widget MUST have their own drawTw<widgetname> function
# All other methods with their TwObj dummy default behavior, which may
# or may not do what you want.
# The functions are automatically discovered as long as they are either defined in
# TermWin, or exported to Main.

if VERSION < v"0.4.0-dev+2275"
    immutable Val{T}
    end
    export Val
end

type TwWindow
    parent::WeakRef # to another TwObj, or nothing
    yloc::Int # 0-based
    xloc::Int
    height::Int # this is to help do box/erase
    width::Int
end

type TwObj{T,S}
    screen::WeakRef # the parent screen
    screenIndex::Int
    window::Union{ Void, Ptr{Void}, TwWindow }
    panel::Union{ Void, Ptr{Void} } # when window is a TwWindow, this is nothing
    height::Int
    width::Int
    xpos::Int
    ypos::Int
    box::Bool
    borderSizeV::Int # how thick is the border at top and botom
    borderSizeH::Int # how thick is the border at left and right edges
    acceptsFocus::Bool
    hasFocus::Bool
    grabUnusedKey::Bool
    isVisible::Bool
    data::T
    value::Any # the logical "content" that this object contains (return value if editable)
    title::UTF8String
    listeners::Dict{ Symbol, Array } # event=>array of registered listeners. each listener is of the type (o, ev)->Void
    function TwObj( data::T )
        log( "TwObj datatype=" * string( T ) * " TwObjSubtype="*string(S) )
        x = new( WeakRef(), 0,
            nothing,
            nothing,
            0, 0, 0, 0,
            false, 0, 0,
            true, true, false, true, data, nothing, utf8(""), Dict{Symbol, Array{Function,1} }() )
        finalizer( x, y->begin
            global rootwin
            if y.panel != nothing
                del_panel( y.panel )
                y.panel = nothing
            end
            if typeof( y.window ) <: Ptr && y.window != rootwin
                delwin( y.window )
                y.window = nothing
            end
            if y.screen.value != nothing
                unregisterTwObj( y.screen.value, y )
                y.screen = WeakRef()
            end
            y.listeners = Dict{Symbol,Array{Function,1}}()
        end )
        x
    end
end

# bookkeeping data for a screen
type TwScreenData
    objects::Array{TwObj, 1 }
    focus::Int
    TwScreenData() = new( TwObj[], 0 )
end

type TwListData
    horizontal::Bool
    widgets::Array{TwObj,1} # this is static.
    focus::Int # which of the widgets has the focus
    canvasheight::Int
    canvaswidth::Int
    pad::Union{ Void, Ptr{Void} } # nothing, or Ptr{Void} to the WINDOW from calling newpad()
    canvaslocx::Int # 0-based, view's location on canvas
    canvaslocy::Int # 0-based
    showLineInfo::Bool
    navigationmode::Bool
    bottomText::UTF8String
    function TwListData()
        ret = new( false, TwObj[], 0, 0, 0, nothing, 0, 0, false, false, utf8("") )
        finalizer( ret, y->begin
            if y.pad != nothing
                delwin( y.pad )
            end
        end)
        ret
    end
end

typealias TwScreen TwObj{TwScreenData}

function TwObj{T,S}( d::T, ::Type{Val{S}} )
    return( TwObj{T,S}(d) )
end
import Base.show

function Base.show( io::IO, o::TwObj{TwListData} )
    if o.data.horizontal
        print( io, "HList(")
    else
        print( io, "VList(")
    end
    for w in o.data.widgets
        strs = split( string(w), "\n" )
        for s in strs
            print( io, "\n  ")
            print( io, s )
        end
    end
    print( io, ")")
end

function Base.show( io::IO, o::TwObj )
    print( io, "TwObj("*string(objtype(o))*"="*string(o.value)*")")
end

draw( p::TwObj ) = error( string( p ) * " draw is undefined.")
objtype{T,S}( _::TwObj{T,S} ) = S


text:
###############################################################################
#
# Plot ellipse on axis
#
###############################################################################
"""
    plot_ellipse!(
                ax::PyObject,
                xy::Tuple{Number,Number};
                width::Number = 10,
                height::Number = 10,
                angle::Number = 0,
                color::String = "black",
                edgecolor::String = color,
                facecolor::String = color,
                alpha::Number = 0.5
    )

Plot an ellipse on axis, given
- `ax` Axis to plot on
- `xy` Center of the ellipse
- `width` Width of the ellipse
- `height` Height of the ellipse
- `angle` Rotation angle of the ellipse
- `color` Color of the ellipse
- `edgecolor` Edgecolor of the ellipse
- `facecolor` Face color of the ellipse
- `alpha` Transparency of the ellipse
"""
function plot_ellipse!(
            ax::PyObject,
            xy::Tuple{Number,Number};
            width::Number = 10,
            height::Number = 10,
            angle::Number = 0,
            color::String = "black",
            edgecolor::String = color,
            facecolor::String = color,
            alpha::Number = 0.5
)
    patch = pyimport("matplotlib.patches");
    outl  = patch.Ellipse(xy=xy,
                          width=width,
                          height=height,
                          edgecolor=edgecolor,
                          facecolor=facecolor,
                          angle=angle,
                          clip_on=false,
                          alpha=alpha);
    ax.add_artist(outl);

    return nothing
end


text:
# # TODO -- find a better home for Domains (or delete them...)
# abstract Domain

# immutable RealDomain <: Domain end
# Base.in(x,::Type{RealDomain}) = typeof(x) <: Real

# immutable NonNegRealDomain <: Domain end
# Base.in(x,::Type{NonNegRealDomain}) = (typeof(x) <: Real && x >= 0)

# domain(::GaussianTransformation) = Real
# domain(::PoissonTransformation) = NonNegRealDomain

## GAUSSIAN ##
"""
    tfm = GaussianTransformation(σ)

Adds zero-mean gaussian noise to an input variable x with standard deviation σ.

X -> ξ,  ξ ∼ N(X,σ)
"""
immutable GaussianTransformation{T<:Real} <: StochasticTransformation
    σ::T
end

# constructors
"""
Adds zero-mean gaussian noise to an input variable x with standard deviation σ.

X -> ξ,  ξ ∼ N(X,σ)
"""
GaussianTransformation(d::Normal) = GaussianTransformation(d.σ)
GaussianTransformation() = GaussianTransformation(1.0)
Transformation(d::Normal) = GaussianTransformation(d.σ)
StochasticTransformation(d::Normal) = GaussianTransformation(d.σ)

# functions
Base.rand(tfm::GaussianTransformation, x) = x + randn()*tfm.σ
transform(tfm::GaussianTransformation, x) = Normal(x,tfm.σ)

## POISSON ##
"""
Transform a non-negative variable, x, into poisson random variable with mean x.

x -> ξ,  ξ ∼ Poisson(x)
"""
immutable PoissonTransformation <: StochasticTransformation end

# constructors
Transformation(d::Poisson) = PoissonTransformation()
StochasticTransformation(d::Poisson) = PoissonTransformation()

# functions
Base.rand(::PoissonTransformation, x) = rand(Poisson(x))
transform(::PoissonTransformation, x) = Poisson(x)

## BERNOULLI ##
"""
Transform an input 0 ≤ x ≤ 1 into a Bernoulli random variable with mean x.

x -> ξ,  ξ ∼ Bernoulli(x)
"""
immutable BernoulliTransformation <: StochasticTransformation end

# constructors
Transformation(d::Bernoulli) = BernoulliTransformation()
StochasticTransformation(d::Bernoulli) = BernoulliTransformation()

# functions
Base.rand(::BernoulliTransformation, x) = rand(Bernoulli(x))
transform(::BernoulliTransformation, x) = Bernoulli(x)

## NEGATIVE BINOMIAL ##
"""
Transform a non-negative variable, x, into a negative binomial random variable
with mean x and (fixed) variance σ²

x -> ξ,  ξ ∼ NegativeBinomial,  E{ξ} = x
"""
immutable NegativeBinomialTransformation{T<:Real} <: StochasticTransformation
    σ²::T
end

# constructors
Transformation(d::NegativeBinomial) = NegativeBinomialTransformation(var(d))
StochasticTransformation(d::NegativeBinomial) = NegativeBinomialTransformation(var(d))
NegativeBinomialTransformation(d::NegativeBinomial) = NegativeBinomialTransformation(var(d))
NegativeBinomialTransformation() = NegativeBinomialTransformation(1.0)

# functions
function specify_moments(::Type{NegativeBinomial},μ,σ²)
    p = μ / σ²
    r = μ*p / (1.0-p)
    return NegativeBinomial(r,p)
end
Base.rand(tfm::NegativeBinomialTransformation, x) = rand(specify_moments(NegBinomial,x,tfm.σ²))
transform(tfm::NegativeBinomialTransformation, x) = specify_moments(NegBinomial,x,tfm.σ²)


text:
"""
    loadjdf(indir, verbose = true)

    loadjdf(indir, cols = Vector{Symbol}, verbose = true)

Load a `DataFrame` from JDF saved at `outdir`. On Julia > v1.3, a multithreaded
version is used.
"""
loadjdf(indir; cols = Symbol[], verbose = false) = begin
	if VERSION < v"1.3.0-rc1.0"
		return sloadjdf(indir, cols = cols, verbose = verbose)
	end

	if verbose
		println("loading $indir in parallel")
	end

	metadatas = open(joinpath(indir,"metadata.jls")) do io
		 deserialize(io)
	end

	if length(cols) == 0
		cols = metadatas.names
	else
		scmn = setdiff(cols, metadatas.names)
		if length(scmn) > 0
			throw("columns $(reduce((x,y) -> string(x) * ", " * string(y), scmn)) are not available, please ensure you have spelt them correctly")
		end
	end

    df = DataFrame()

	# get the maximum number of bytes needs to read
	# bytes_needed = maximum(get_bytes.(metadatas.metadatas))

	results = Vector{Any}(undef, length(cols))

	# rate limit channel
	c1 = Channel{Bool}(Threads.nthreads())
	atexit(()->close(c1))

	i = 1
    for (name, metadata) in zip(metadatas.names, metadatas.metadatas)
		if name in cols
			put!(c1, true)
			results[i] = @spawn begin
				io = BufferedInputStream(open(joinpath(indir,string(name)), "r"))
				new_result = column_loader(metadata.type, io, metadata)
				close(io)
				(name = name, task = new_result)
			end
			take!(c1)
			i+=1
		end
    end

	# run this serially
	for result in results
		if verbose
			println("Extracting $(result.name)")
		end

		new_result = fetch(result).task
		colname = fetch(result).name
		if new_result == nothing
			df[!, colname] = Vector{Missing}(missing, metadatas.rows)
		else
			df[!, colname] = new_result
		end
	end
 	df
end

loadjdf(jdf::JDFFile; args...) = loadjdf(jdf.path; args...)
sloadjdf(jdf::JDFFile; args...) = sloadjdf(jdf.path; args...)

# load the data from file with a schema
sloadjdf(indir; cols = Symbol[], verbose = false) = begin
	metadatas = open(joinpath(indir,"metadata.jls")) do io
		deserialize(io)
	end

	if length(cols) == 0
		cols = metadatas.names
	else
		scmn = setdiff(cols, metadatas.names)
		if length(scmn) > 0
			throw("columns $(reduce((x,y) -> string(x) * ", " * string(y), scmn)) are not available, please ensure you have spelt them correctly")
		end
	end

    df = DataFrame()

	# get the maximum number of bytes needs to read
	# bytes_needed = maximum(get_bytes.(metadatas.metadatas))

	# rate limit channel
	#results = Vector{Any}(undef, length(metadatas.names))
	results = Vector{Any}(undef, length(cols))

	i = 1
    for (name, metadata) in zip(metadatas.names, metadatas.metadatas)
		if name in cols
			if verbose
				println("Loading $name")
			end
			results[i] = begin
				io = BufferedInputStream(open(joinpath(indir,string(name)), "r"))
				new_result = column_loader(metadata.type, io, metadata)
				close(io)
				(name = name, task = new_result)
			end
			i+=1
		end
    end

	# run this serially
	for result in results
		if verbose
			println("Extracting $(result.name)")
			println(result.task)
		end

		new_result = result.task
		colname = result.name
		if new_result == nothing
			df[!, colname] = Vector{Missing}(missing, metadatas.rows)
		else
			df[!, colname] = new_result
		end
	end
 	df
end


text:
using Base.Iterators

function plot_along_param_curve!(curve, contact_str, step; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0, frame=false)
	ts = collect(0:step:1)
	points = curve.(ts)
	plt = nothing
	for (x,y,z) in points
		basis = contact_str.distribution(x,y,z)
		basis_a = normalize(basis[1,:])
		basis_b = normalize(basis[2,:])
		plt = plot_plane([x,y,z],basis_a,basis_b, :grey, ε, frame)
	end
	scatter3d!(plt,(x -> x[1]).(points),(x -> x[2]).(points),(x -> x[3]).(points), color=midpoint_marker_color, markersize = markersize)
	return plt
end

function plot_along_param_curve(curve, contact_str, step; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	surface()
	plot_along_param_curve!(curve, contact_str, step; midpoint_marker_color, markersize, ε)
end

function plot_along_param_curve!(curve, contact_form :: OneForm, step; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	plot_along_param_curve!(curve, structuralize_form(contact_form), step; midpoint_marker_color, markersize, ε)
end

function plot_along_param_curve(curve, contact_form :: OneForm, step; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	plot_along_param_curve(curve, structuralize_form(contact_form), step; midpoint_marker_color, markersize, ε)
end

function plot_along_xy_grid!(xs, ys, z, contact_str; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	grid = vcat(collect(product(xs,ys))...)
	plt = nothing
	for (x,y) in grid
		basis = contact_str.distribution(x,y,z)
		basis_a = normalize(basis[1,:])
		basis_b = normalize(basis[2,:])
		plt = plot_plane([x,y,z],basis_a,basis_b, :grey, ε)
	end
	scatter3d!((x -> x[1]).(grid),(x -> x[2]).(grid),repeat([z], size(grid)...), color=midpoint_marker_color, markersize=markersize)
	return plt
end

function plot_along_xy_grid(xs, ys, z, contact_str; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	surface()
	plot_along_xy_grid!(xs, ys, z, contact_str; midpoint_marker_color, markersize, ε)
end

function plot_along_xy_grid!(xs, ys, z, contact_form::OneForm; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	plot_along_xy_grid!(xs, ys, z,structuralize_form(contact_form); midpoint_marker_color, markersize, ε)
end

function plot_along_xy_grid(xs, ys, z, contact_form::OneForm; midpoint_marker_color = :red, markersize = default(:markersize), ε = 1.0)
	plot_along_xy_grid(xs, ys, z,structuralize_form(contact_form); midpoint_marker_color, markersize, ε)
end


text:
module ExceptionUnwrappingTest

using Test
using ExceptionUnwrapping
using ExceptionUnwrapping: UnwrappedExceptionNotFound

# TaskFailedException is available in Julia 1.3+
if VERSION >= v"1.3.0-"
    @testset "Wrapped TaskFailedException" begin
        try
            fetch(@async fetch(@async error("hi")))
        catch e
            @assert (e isa ErrorException) == false
            @assert e isa TaskFailedException

            @test is_wrapped_exception(e)
            @test unwrap_exception(e) isa TaskFailedException
            @test unwrap_exception_to_root(e) isa ErrorException
            @test has_wrapped_exception(e, ErrorException)
            @test has_wrapped_exception(e, TaskFailedException)
            @test has_wrapped_exception(e, ArgumentError) == false
            @test unwrap_exception_until(e, ErrorException) isa ErrorException
            @test_throws UnwrappedExceptionNotFound{ArgumentError} unwrap_exception_until(e, ArgumentError) isa ErrorException
        end
    end
end

struct MyWrappedException{T}
    wrapped_exception::T
end
# Implement the ExceptionUnwrapping API:
ExceptionUnwrapping.unwrap_exception(e::MyWrappedException) = e.wrapped_exception

@testset "Custom Exception Types" begin
    e1 = ErrorException("1")
    e2 = MyWrappedException(ErrorException("1"))
    try
        throw(e2)
    catch e
        @assert e === e2
        @assert (e isa ErrorException) == false
        @assert e isa MyWrappedException

        @test is_wrapped_exception(e2)
        @test !is_wrapped_exception(e1)
        @test unwrap_exception(e) === e1
        @test unwrap_exception_to_root(e) === e1
        @test has_wrapped_exception(e, ErrorException)
        @test has_wrapped_exception(e, MyWrappedException)
        @test has_wrapped_exception(e, ArgumentError) == false
        @test unwrap_exception_until(e, ErrorException) === e1
        @test_throws UnwrappedExceptionNotFound{ArgumentError} unwrap_exception_until(e, ArgumentError) isa ErrorException
    end
end

struct MyWrappedException2
    exc::Any
end
# Implement the ExceptionUnwrapping API:
ExceptionUnwrapping.unwrap_exception(e::MyWrappedException2) = e.exc

@testset "Multiple Layers" begin
    e1 = ErrorException("1")
    e2 = MyWrappedException(ErrorException("1"))
    e3 = MyWrappedException2(MyWrappedException(ErrorException("1")))
    try
        throw(e3)
    catch e
        @assert e === e3

        @test is_wrapped_exception(e3)
        @test is_wrapped_exception(e2)
        @test !is_wrapped_exception(e1)
        @test unwrap_exception(e) === e2
        @test unwrap_exception(unwrap_exception(e)) === e1
        @test unwrap_exception_to_root(e) === e1
        @test has_wrapped_exception(e, ErrorException)
        @test has_wrapped_exception(e, MyWrappedException)
        @test has_wrapped_exception(e, MyWrappedException2)
        @test unwrap_exception_until(e, MyWrappedException) === e2
        @test unwrap_exception_until(e, ErrorException) === e1
    end
end

end # module


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
using Test
import Random

@testset "FMI.jl" begin
    if Sys.iswindows()
        @info "Automated testing for Windows is supported."
        @testset "FMU functions" begin
            include("getterSetterTest_fmu.jl")
            include("independentFunctionsTest_fmu.jl")
            include("stateTest_fmu.jl")
        end
        @testset "fmi2Component functions" begin
            include("getterSetterTest_comp.jl")
            include("independentFunctionsTest_comp.jl")
            include("stateTest_comp.jl")
        end
        @testset "Simulation Tests" begin
            include("test_setter_getter.jl")
            include("test_sim_cs.jl")
            include("test_sim_me.jl")
        end
    elseif Sys.islinux()
        @warn "Test-sets are currrently using Windows-FMUs, automated testing for Linux is currently not supported."
    elseif Sys.isapple()
        @warn "Test-sets are currrently using Windows-FMUs, automated testing for macOS is currently not supported."
    end
end


text:
"""
    timescale(ep)

Return the time scale of epoch `ep`.

# Example

```jldoctest; setup = :(using AstroTime)
julia> ep = TTEpoch(2000, 1, 1)
2000-01-01T00:00:00.000 TT

julia> timescale(ep)
TT
```
"""
timescale(ep::Epoch) = ep.scale

"""
    DateTime(ep::Epoch)

Convert the epoch `ep` to an `AstroDates.DateTime`.
"""
function DateTime(ep::Epoch)
    if !isfinite(ep.fraction)
        if ep.fraction < 0
            return DateTime(AstroDates.MIN_EPOCH, AstroDates.H00)
        else
            return DateTime(AstroDates.MAX_EPOCH, Time(23, 59, 59.999))
        end
    end

    sec = ep.second + 43200
    time = sec % 86400
    if time < 0
        time += 86400
    end
    date = (sec - time) ÷ 86400

    date_comp = Date(AstroDates.J2000_EPOCH, date)
    time_comp = Time(time, ep.fraction)

    return DateTime(date_comp, time_comp)
end

"""
    year(ep::Epoch)

Get the year of the epoch `ep`.
"""
year(ep::Epoch) = year(DateTime(ep))

"""
    month(ep::Epoch)

Get the month of the epoch `ep`.
"""
month(ep::Epoch) = month(DateTime(ep))

"""
    day(ep::Epoch)

Get the day of the epoch `ep`.
"""
day(ep::Epoch) = day(DateTime(ep))

"""
    yearmonthday(ep::Epoch)

Get the year, month, and day of the epoch `ep` as a tuple.
"""
yearmonthday(ep::Epoch) = yearmonthday(DateTime(ep))

"""
    dayofyear(ep::Epoch)

Get the day of the year of the epoch `ep`.
"""
dayofyear(ep::Epoch) = dayofyear(DateTime(ep))

"""
    hour(ep::Epoch)

Get the hour of the epoch `ep`.
"""
hour(ep::Epoch) = hour(DateTime(ep))

"""
    minute(ep::Epoch)

Get the minute of the epoch `ep`.
"""
minute(ep::Epoch) = minute(DateTime(ep))

"""
    second(type, ep::Epoch)

Get the second of the epoch `ep` as a `type`.
"""
second(typ, ep::Epoch) = second(typ, DateTime(ep))

"""
    second(ep::Epoch) -> Int

Get the second of the epoch `ep` as an `Int`.
"""
second(ep::Epoch) = second(Int, DateTime(ep))

"""
    millisecond(ep::Epoch)

Get the number of milliseconds of the epoch `ep`.
"""
millisecond(ep::Epoch) = millisecond(DateTime(ep))

"""
    microsecond(ep::Epoch)

Get the number of microseconds of the epoch `ep`.
"""
microsecond(ep::Epoch) = microsecond(DateTime(ep))

"""
    nanosecond(ep::Epoch)

Get the number of nanoseconds of the epoch `ep`.
"""
nanosecond(ep::Epoch) = nanosecond(DateTime(ep))

"""
    subsecond(ep::Epoch, n)

Get the number of fractional seconds with the unit ``s * \\frac{1}{10^n}``,
e.g. `subsecond(ep, 3)` for milliseconds, of the epoch `ep`. `n` must be
divisible by 3.
"""
AstroDates.subsecond(ep::Epoch, n) = subsecond(DateTime(ep), n)

"""
    Time(ep::Epoch)

Get the `Time` of the epoch `ep`.
"""
Time(ep::Epoch) = Time(DateTime(ep))

"""
    Date(ep::Epoch)

Get the `Date` of the epoch `ep`.
"""
Date(ep::Epoch) = Date(DateTime(ep))

"""
    fractionofday(ep::Epoch)

Get the time of the day of the epoch `ep` as a fraction.
"""
fractionofday(ep::Epoch) = fractionofday(Time(ep))

"""
    fractionofsecond(ep::Epoch)

Get the fraction of the current second of the epoch `ep`.
"""
AstroDates.fractionofsecond(ep::Epoch) = ep.fraction

"""
    Dates.DateTime(ep::Epoch)

Convert the epoch `ep` to a `Dates.DateTime`.
"""
Dates.DateTime(ep::Epoch) = Dates.DateTime(DateTime(ep))

"""
    Dates.Date(ep::Epoch)

Convert the date of epoch `ep` to a `Dates.Date`.
"""
Dates.Date(ep::Epoch) = Dates.Date(Date(ep))

"""
    Dates.Time(ep::Epoch)

Convert the time of epoch `ep` to a `Dates.Time`.
"""
Dates.Time(ep::Epoch) = Dates.Time(Time(ep))

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
struct VariableStd <: DataItem
    name :: String
    value :: Vector{<:Real}
    offset :: Real
    scale :: Real
end

"""
    VariableStd(name::String, z::Vector{<:Real})::VariableStd

Construct a VariableStd type from the standardized vector `z` with the name `name`.

"""
function VariableStd(name::String, z::Vector{<:Real})::VariableStd

    m = mean(z)
    s = std(z)

    @assert isapprox(m, 0.0, atol = sqrt(eps())) "the mean should be equal to zero. Use the type Variable instead."
    @assert isapprox(s, 1.0, atol = sqrt(eps())) "the standard deviation should be equal to one. Use the type Variable instead."

    return VariableStd(name,z,0,1)
end

"""
    reconstruct(vstd::VariableStd)

Reconstruct the original Variable from the standardized one.

See also

See also [`Variable`](@ref), [`VariableStd`](@ref)  and [`standardize`](@ref).
"""
function reconstruct(vstd::VariableStd)::Variable

    z = vstd.value

    x = vstd.scale*z .+ vstd.offset

    v = Variable(vstd.name, x)

    return v

end


text:
#! format: off
get_variable_multiplier(::SystemBalanceSlackUp, ::Type{<: Union{PSY.Bus, PSY.System}}, _) = 1.0
get_variable_multiplier(::SystemBalanceSlackDown, ::Type{<: Union{PSY.Bus, PSY.System}}, _) = -1.0
#! format: on

function add_variables!(
    container::OptimizationContainer,
    ::Type{T},
    ::PSY.System,
    ::Type{U},
) where {
    T <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    U <: Union{CopperPlatePowerModel, StandardPTDFModel},
}
    time_steps = get_time_steps(container)
    variable = add_variable_container!(container, T(), PSY.System, time_steps)

    for t in time_steps
        variable[t] = JuMP.@variable(
            container.JuMPmodel,
            base_name = "slack_{$(T), $t}",
            lower_bound = 0.0
        )
    end
    return
end

function add_variables!(
    container::OptimizationContainer,
    ::Type{T},
    sys::PSY.System,
    ::Type{U},
) where {
    T <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    U <: PM.AbstractActivePowerModel,
}
    time_steps = get_time_steps(container)
    bus_numbers = PSY.get_number.(PSY.get_components(PSY.Bus, sys))
    variable = add_variable_container!(container, T(), PSY.Bus, bus_numbers, time_steps)

    for t in time_steps, n in bus_numbers
        variable[n, t] = JuMP.@variable(
            container.JuMPmodel,
            base_name = "slack_{$(T), $n, $t}",
            lower_bound = 0.0
        )
    end
    return
end

function add_variables!(
    container::OptimizationContainer,
    ::Type{T},
    sys::PSY.System,
    ::Type{U},
) where {
    T <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
    U <: PM.AbstractPowerModel,
}
    time_steps = get_time_steps(container)
    bus_numbers = PSY.get_number.(PSY.get_components(PSY.Bus, sys))
    variable_active =
        add_variable_container!(container, T(), PSY.Bus, "P", bus_numbers, time_steps)
    variable_reactive =
        add_variable_container!(container, T(), PSY.Bus, "Q", bus_numbers, time_steps)

    for t in time_steps, n in bus_numbers
        variable_active[n, t] = JuMP.@variable(
            container.JuMPmodel,
            base_name = "slack_{p, $(T), $n, $t}",
            lower_bound = 0.0
        )
        variable_reactive[n, t] = JuMP.@variable(
            container.JuMPmodel,
            base_name = "slack_{q, $(T), $n, $t}",
            lower_bound = 0.0
        )
    end
    return
end

function objective_function!(
    container::OptimizationContainer,
    ::Type{PSY.System},
    model::NetworkModel{T},
    S::Type{T},
) where {T <: Union{CopperPlatePowerModel, StandardPTDFModel}}
    variable_up = get_variable(container, SystemBalanceSlackUp(), PSY.System)
    variable_dn = get_variable(container, SystemBalanceSlackDown(), PSY.System)

    for t in get_time_steps(container)
        add_to_objective_invariant_expression!(
            container,
            (variable_dn[t] + variable_up[t]) * BALANCE_SLACK_COST,
        )
    end
    return
end

function objective_function!(
    container::OptimizationContainer,
    ::Type{PSY.Bus},
    model::NetworkModel{T},
    S::Type{T},
) where {T <: PM.AbstractActivePowerModel}
    variable_up = get_variable(container, SystemBalanceSlackUp(), PSY.Bus)
    variable_dn = get_variable(container, SystemBalanceSlackDown(), PSY.Bus)
    bus_numbers = axes(variable_up)[1]
    @assert_op bus_numbers == axes(variable_dn)[1]
    for t in get_time_steps(container), n in bus_numbers
        add_to_objective_invariant_expression!(
            container,
            (variable_dn[n, t] + variable_up[n, t]) * BALANCE_SLACK_COST,
        )
    end
    return
end

function objective_function!(
    container::OptimizationContainer,
    ::Type{PSY.Bus},
    model::NetworkModel{T},
    S::Type{T},
) where {T <: PM.AbstractPowerModel}
    variable_p_up = get_variable(container, SystemBalanceSlackUp(), PSY.Bus, "P")
    variable_p_dn = get_variable(container, SystemBalanceSlackDown(), PSY.Bus, "P")
    variable_q_up = get_variable(container, SystemBalanceSlackUp(), PSY.Bus, "Q")
    variable_q_dn = get_variable(container, SystemBalanceSlackDown(), PSY.Bus, "Q")
    bus_numbers = axes(variable_p_up)[1]
    @assert_op bus_numbers == axes(variable_q_dn)[1]
    for t in get_time_steps(container), n in bus_numbers
        add_to_objective_invariant_expression!(
            container,
            (
                variable_p_dn[n, t] +
                variable_p_up[n, t] +
                variable_q_dn[n, t] +
                variable_q_up[n, t]
            ) * BALANCE_SLACK_COST,
        )
    end
    return
end


text:
module OpenTelemetryInstrumentationHTTP

using OpenTelemetryAPI

using HTTP
using Sockets
using URIs

using TOML

const PKG_VERSION =
    VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
const INSTRUMENTATION_INFO = InstrumentationInfo(
    name = string(@__MODULE__),
    version = PKG_VERSION,
    schema_url = "https://oolong.dev/OpenTelemetry.jl/dev/OpenTelemetryInstrumentationHTTP/",
)

abstract type OpenTelemetryLayer{Next<:HTTP.Layer} <: HTTP.Layer{Next} end

const HTTP_TRACER = Ref{Tracer}()
const HTTP_METRICS = Ref{Counter{UInt}}()

function HTTP.request(
    ::Type{OpenTelemetryLayer{Next}},
    method,
    url::URI,
    headers,
    body;
    kw...,
)::HTTP.Response where {Next}
    with_span(
        "http $method",
        HTTP_TRACER[];
        kind = SPAN_KIND_CLIENT,
        attributes = Dict{String,TAttrVal}(
            "http.method" => method,
            "http.url" => string(
                URI(;
                    scheme = url.scheme,
                    userinfo = "",  # !!! ignore credentials
                    host = url.host,
                    port = url.port,
                    path = url.path,
                    query = url.query,
                    fragment = url.fragment,
                ),
            ),
            "http.target" => url.path * "?" * url.query,
            "http.host" => url.host |> string,
            "http.scheme" => url.scheme |> string,
            "net.peer.name" => url.host |> string,
            "net.peer.port" => parse(Int, url.port),
        ),
    ) do
        resp = HTTP.request(Next, method, url, inject!(headers), body; kw...)
        s = current_span()
        s["http.status_code"] = resp.status |> Int
        s["http.flavor"] = string(resp.version)
        s["http.response_content_length"] = length(resp.body)
        HTTP_METRICS[](; scheme = string(url.scheme), status = Int(resp.status))
        resp
    end
end

"""
## Metrics

| Meter Name | Instrument Name | Instrument Type | Unit | Dimensions        | Description         |
|:---------- |:--------------- |:--------------- |:---- |:----------------- |:------------------- |
| `http`     | `request`       | `Counter{UInt}` |      | `route`, `status` | Number of requests. |

## Spans

We try to follow [semantic conventions for HTTP
spans](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/trace/semantic_conventions/http.md)
here.

The following attributes are added on span creation:

  - `http.method`
  - `http.target`
  - `http.flavor`
  - `http.request_content_length`
  - `http_route`

The following attributes are added when response is received:

  - `http.status_code`
  - `http.response_content_length`
"""
function init(;
    tracer_provider = global_tracer_provider(),
    meter_provider = global_meter_provider(),
)
    HTTP_METRICS[] = Counter{UInt}(
        "http",
        Meter(
            "request";
            provider = meter_provider,
            instrumentation_info = INSTRUMENTATION_INFO,
        );
        unit = "",
        description = "Number of requests received.",
    )
    HTTP_TRACER[] =
        Tracer(; provider = tracer_provider, instrumentation_info = INSTRUMENTATION_INFO)
    top = HTTP.top_layer(stack())
    insert_default!(top, OpenTelemetryLayer)
end

function HTTP.serve(
    f,
    host::Union{IPAddr,String} = Sockets.localhost,
    port::Integer = 8081;
    stream::Bool = false,
    kw...,
)
    handler = if f isa HTTP.Handlers.Handler
        f
    elseif stream
        HTTP.StreamHandlerFunction(f)
    else
        HTTP.RequestHandlerFunction(f)
    end
    return HTTP.Servers.listen(x -> otel_handle(handler, x), host, port; kw...)
end

otel_handle(handler, x) = HTTP.handle(handler, x)

function otel_handle(r::HTTP.Router, http::HTTP.Stream)
    req::HTTP.Request = http.message
    req.body = read(http)
    closeread(http)

    # directly taken from https://github.com/JuliaWeb/HTTP.jl/blob/46687c8e594cbde40cdcfa63fe66123ebad1ea5c/src/Handlers.jl#L459-L468
    #
    # BEGIN
    # get the url/path of the request
    m = Val(Symbol(req.method))
    # get scheme, host, split path into strings and get Vals
    uri = URI(req.target)
    s = get(HTTP.Handlers.SCHEMES, uri.scheme, HTTP.Handlers.EMPTYVAL)
    h = Val(Symbol(uri.host))
    p = uri.path
    segments = split(p, '/'; keepempty = false)
    # dispatch to the most specific handler, given the path
    vals = (get(r.segments, s, Val(Symbol(s))) for s in segments)
    # END

    m_handler = which(
        HTTP.Handlers.gethandler,
        (typeof(r), typeof(m), typeof(s), typeof(h), map(typeof, vals)...),
    )

    handler = HTTP.Handlers.gethandler(r, m, s, h, vals...)

    name = if handler === r.default
        "DEFAULT ROUTE"
    else
        "/" * join(map(type2segment, m_handler.sig.parameters[6:end]), "/")
    end

    with_context(extract(req.headers)) do
        with_span(
            name,
            HTTP_TRACER[];
            kind = SPAN_KIND_SERVER,
            attributes = Dict{String,TAttrVal}(
                "http.method" => req.method,
                "http.target" => req.target,
                "http.flavor" => string(req.version),
                "http.request_content_length" => length(req.body),
                "http.route" => name,
            ),
        ) do
            s = current_span()
            resp = HTTP.Handlers.handle(handler, req)
            req.response::HTTP.Response = resp
            req.response.request = req
            startwrite(http)
            write(http, req.response.body)
            s["http.status_code"] = string(resp.status)
            s["http.response_content_length"] = length(resp.body)
            HTTP_METRICS[](; route = name, status = Int(resp.status))
        end
    end
    nothing
end

type2segment(x) = ""
type2segment(x::Type{Val{S}}) where {S} = S
type2segment(x::Type{Any}) = "*"

function __init__()
    init()
end

end # module


text:
# p4.jl - periodic spectral differentiation

# Set up grid and differentiation matrix:
N = 24; h = 2*pi/N; x = h*(1:N);
column = [0; @. .5*(-1)^(1:N-1)*cot((1:N-1)*h/2)];
D = toeplitz(column,column[[1;N:-1:2]]);

# Differentiation of a hat function:
v = @. max(0,1-abs(x-pi)/2); clf();
subplot(2,2,1);  plot(x,v,".-",markersize=6);
axis([0,2pi,-.5,1.5]); grid(true); title("function");
subplot(2,2,2), plot(x,D*v,".-",markersize=6);
axis([0,2pi,-1,1]); grid(true); title("spectral derivative");

# Differentiation of exp(sin(x)):
v = @. exp(sin(x)); vprime = @. cos(x)*v;
subplot(2,2,3), plot(x,v,".-",markersize=6);
axis([0,2pi,0,3]), grid(true);
subplot(2,2,4), plot(x,D*v,".-",markersize=6);
axis([0,2pi,-2,2]), grid(true);
error = signif(norm(D*v-vprime,Inf),5);
text(2.2,1.4,"max error = $error",fontsize=8);


text:
module Master

include("./constraints/geometric/variables.jl")
include("./constraints/geometric/sentinel.jl")
include("./constraints/geometric/rotations_currence.jl")
include("./constraints/geometric/non_overlapping.jl")
include("./constraints/geometric/inside_pallet.jl")
include("./constraints/static_stability/partial_base/0_partial_base.jl")
include("./constraints/static_stability/partial_base_2/0_partial_base.jl")
include("./constraints/static_stability/on_top/0_on_top.jl")
include("./constraints/static_stability/ss_mult_boxes/0_ss_mult_boxes.jl")
include("./constraints/dynamic_stability/guillotine/0_guillotine.jl")
include("./constraints/dynamic_stability/guillotine/callback.jl")
include("./constraints/dynamic_stability/middle/0_middle.jl")
include("./constraints/objective/constraints/objective.jl")
include("./constraints/geometric/constants.jl")
include("./constraints/dynamic_stability/mult_boxes/0_mult_boxes.jl")
include("./constraints/load_balance/0_load_balancing.jl")
include("./constraints/robot/robot_pattern/0_robot_pattern.jl")
include("./constraints/robot/simple/0_robot_simple.jl")
include("./constraints/robot/simple/0_robot_simple_fixed.jl")
include("./constraints/load_bearing/partial_base/0_partial_base_lb.jl")
include("./constraints/load_bearing/simple/0_pressure_lb.jl")
include("./constraints/load_bearing/simple/0_number_lb.jl")
include("./results/write_to_file.jl")
include("./results/write_to_json.jl")
using JuMP, Gurobi
using .Variables
using .Sentinel
using .RotationsCurrence
using .NonOverlapping
using .InsidePallet
using .PartialBase
using .PartialBase2
using .OnTop
using .CMultipleBoxes
using .Guillotine
using .GuillotineCallback
using .Middle
using .Objective
using .Constants
using .MultBoxes
using .LoadBalancing
using .RobotPattern
using .RobotSimple
using .RobotSimpleFixed
using .LBPartialBase
using .LBSimple
using .LBNumber
using .WriteToFile
using .WriteToJSON

function init_model()
    return direct_model(Gurobi.Optimizer())
end

export main

function main(mode::Int64)
    m = init_model()
    m = init_variables(m)
    m = init_sentinel(m)
    m = cons_rotations_currence(m)
    m = cons_non_overlapping(m)
    m = cons_inside_pallet(m)

    if (mode == 1)
        # LB

        m = cons_lb_load_balancing(m)
    elseif (mode == 2)
        # LBP

        m = cons_lb_simple_load_bearing(m)
    elseif (mode == 3)
        # Robot

        m = cons_rp_robot_simple_fixed(m)
    elseif (mode == 4)
        #PBS1

        m = cons_ss_partial_base(m)
    elseif (mode == 5)
        #PBS2

        m = cons_ss_partial_base_two(m)
    elseif (mode == 6)
        # LBP OnTop

        m = cons_ss_on_top(m)
        m = cons_lb_simple_load_bearing(m)
    elseif (mode == 7)
        # OnTop

        m = cons_ss_on_top(m)
    elseif (mode == 8)
        # Guillotine

        m = cons_ds_guillotine(m)
    elseif (mode == 9)
        # Guillotine OnTop

        m = cons_ss_on_top(m)
        m = cons_ds_guillotine(m)
    elseif (mode == 10)
        # Guillotine OnTop LB LBP

        m = cons_ss_on_top(m)
        m = cons_ds_guillotine(m)
        m = cons_lb_load_balancing(m)
        m = cons_lb_simple_load_bearing(m)
    elseif (mode == 11)
        # PBS2 LB LBP

        m = cons_ss_partial_base_two(m)
        m = cons_lb_load_balancing(m)
        m = cons_lb_simple_load_bearing(m)
    elseif (mode == 12)
        # PBS2 Guillotine LB LBP

        m = cons_ss_partial_base_two(m)
        m = cons_ds_guillotine(m)
        m = cons_lb_load_balancing(m)
        m = cons_lb_simple_load_bearing(m)
    elseif (mode == 13)
        #PBS1 LB LBP MIDDLE MB

        m = cons_ss_partial_base(m)
        m = cons_lb_load_balancing(m)
        m = cons_lb_simple_load_bearing(m)
        m = cons_ds_mult_boxes(m)
        m = cons_ds_middle(m)
    elseif (mode == 14)
        # MB OnTop

        m = cons_ss_on_top(m)
        m = cons_ds_mult_boxes(m)
    elseif (mode == 15)
        # MIDDLE OnTop

        m = cons_ss_on_top(m)
        m = cons_ds_mult_boxes(m)
    elseif (mode == 16)
        #MBSS ONTOP

        m = cons_ss_partial_base_two(m)
        m = cons_ss_mult_boxes(m)
    end

    if (mode == 8 || mode == 9 || mode == 10 || mode == 12)
        # Guillotine

        m = objective(m, 1)
    elseif (mode == 13)
        # MIDDLE MB

        m = objective(m, 2)
    elseif (mode == 14)
        # MB

        m = objective(m, 3)
    elseif (mode == 13)
        # MIDDLE

        m = objective(m, 4)
    else
        # NONE

        m = objective(m, 0)
    end

    function callback_function(cb_data)
        #guillotine_callback(cb_data, m)
        return
    end

    MOI.set(m, MOI.RawParameter("TimeLimit"), 2400)

    try
        optimize!(m)
    catch err
        print("\n ERROR \n")
    end

    try
        print("\n MaxHEIGHT \n")
        print(value(m[:maxHeight]))
        print('\n')
    catch err
        print("\n ERROR \n")
    end

    WriteToFile.append_to_file(m, mode)
    """
    if (has_values(m))
        WriteToJSON.write_to_json_simple(m)
    end
    """

end

#SPECIFY MODE
for i in [5, 11, 12]
    main(i)
end

end


text:
using DoubleDouble
using Base.Test

x = sqrt(2.0)
bx = big(x)
sx = Single(x)
dx = Double(x)

y = 0.1
by = big(y)
sy = Single(y)
dy = Double(y)

@test x == sx == dx
@test y == sy == dy

dxy = dx*dy
bxy = bx*by
@test sx*sy == dxy
@test x*y == Float64(dxy)
@test dxy == Double(bxy)

@test x+y == Float64(dx+dy)
@test dx+dy == Double(bx+by)

@test x-y == Float64(dx-dy)
@test dx-dy == Double(bx-by)

@test x/y == Float64(dx/dy)
@test dx/dy == Double(bx/by)

@test sqrt(y) == Float64(sqrt(dy))
@test sqrt(dy) == Double(sqrt(by))

#@test rem(dxy,1.0) == Double(rem(bxy,1.0))

## New
@test Double(pi) == Double{Float64}(3.141592653589793, 1.2246467991473532e-16)
@test Double(3.5) == Double{Float64}(3.5, 0.0)
@test Double(3.5) == Double{Float64}(3.5, 0.0)

a = Double(big"3.1")
@test a == Double(3.1, -8.881784197001253e-17)
@test BigFloat(a) == big"3.099999999999999999999999999999995069619342368676216176696466982586064542459781"

@test Double(3) == Double(3.0, 0.0)
@test Double(big(3)) == Double(3.0, 0.0)


text:
# This file includes and example file for 2D interpolation

using FractalTools
using GeometryBasics
using CairoMakie

# Generate data
f(x, y) = x^2 + y^2 + 1
ngon = Triangle(
    Point(BigFloat(0.), BigFloat(0.)),
    Point(BigFloat(1.), BigFloat(0.)),
    Point(BigFloat(0.5), BigFloat(1.)))
npts = 100
pts = getdata(f, ngon, npts)
interp = interpolate(pts, Interp2D(0.01))

function err(x, y)
    fval = f(x, y)
    ival = interp(x, y)
    abs(fval - ival) / abs(fval) * 100
end

tpts = getdata(ngon, npts)

fig, ax, plt = trisurf(tpts, f, meshcolor3=:red)
trisurf!(ax, tpts, interp, meshcolor3=:blue)
fig, ax, plt = trisurf(tpts, err, meshcolor3=:red)


text:
function _end(obj::RealVectorChangingVisitor)
    return jcall(obj, "end", jdouble, ())
end

function start(obj::RealVectorChangingVisitor, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "start", void, (jint, jint, jint), arg0, arg1, arg2)
end

function visit(obj::RealVectorChangingVisitor, arg0::jint, arg1::jdouble)
    return jcall(obj, "visit", jdouble, (jint, jdouble), arg0, arg1)
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _altaz2hadec{T<:AbstractFloat}(alt::T, az::T, lat::T)
    # Convert to radians.
    alt_r = deg2rad(alt)
    az_r = deg2rad(az)
    lat_r = deg2rad(lat)
    # Find local hour angle (in degrees, from 0. to 360.).
    ha = rad2deg(atan2(-sin(az_r)*cos(alt_r),
                       -cos(az_r)*sin(lat_r)*cos(alt_r) +
                       sin(alt_r)*cos(lat_r)))
    ha < 0.0 && (ha += 360.0)
    ha = ha.%360.0
    # Find declination (positive if north of Celestial Equator, negative if
    # south)
    sindec = sin(lat_r)*sin(alt_r) + cos(lat_r)*cos(alt_r)*cos(az_r)
    dec = rad2deg(asin(sindec))  # convert dec to degrees
    return ha, dec
end

"""
    altaz2hadec(alt, az, lat) -> ha, dec

### Purpose ###

Convert Horizon (Alt-Az) coordinates to Hour Angle and Declination.

### Explanation ###

Can deal with the NCP singularity.  Intended mainly to be used by program
`hor2eq`.

### Arguments ###

Input coordinates may be either a scalar or an array, of the same dimension.

* `alt`: local apparent altitude, in degrees, scalar or array.
* `az`: the local apparent azimuth, in degrees, scalar or vector, measured
 *east* of *north*!!!  If you have measured azimuth west-of-south (like the book
 Meeus does), convert it to east of north via: `az = (az + 180) % 360`.
* `lat`: the local geodetic latitude, in degrees, scalar or array.

`alt` and `az` can be given as a 2-tuple `(alt, az)`.

### Output ###

2-tuple `(ha, dec)`

* `ha`: the local apparent hour angle, in degrees.  The hour angle is the time
  that right ascension of 0 hours crosses the local meridian.  It is
  unambiguously defined.
* `dec`: the local apparent declination, in degrees.

The output coordinates are always floating points and have the same type (scalar
or array) as the input coordinates.

### Example ###

Arcturus is observed at an apparent altitude of 59d,05m,10s and an azimuth
(measured east of north) of 133d,18m,29s while at the latitude of +43.07833
degrees.  What are the local hour angle and declination of this object?

``` julia
ha, dec = altaz2hadec(ten(59,05,10), ten(133,18,29), 43.07833)
# => (336.6828582472844,19.182450965120402)
```
The widely available XEPHEM code gets:

    Hour Angle = 336.683
    Declination = 19.1824

### Notes ###

`hadec2altaz` converts Hour Angle and Declination to Horizon (Alt-Az)
coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
altaz2hadec(alt::Real, az::Real, lat::Real) =
    _altaz2hadec(promote(float(alt), float(az), float(lat))...)

altaz2hadec(altaz::Tuple{Real, Real}, lat::Real) = altaz2hadec(altaz..., lat)

function altaz2hadec{R1<:Real, R2<:Real, R3<:Real}(alt::AbstractArray{R1},
                                                   az::AbstractArray{R2},
                                                   lat::AbstractArray{R3})
    @assert length(alt) == length(az) == length(lat)
    typealt = typeof(float(one(R1)))
    ha  = similar(alt, typealt)
    dec = similar(alt, typealt)
    for i in eachindex(alt)
        ha[i], dec[i] = altaz2hadec(alt[i], az[i], lat[i])
    end
    return ha, dec
end


text:
abstract type AbstractWorkbook end
abstract type AbstractWorksheet{W <: AbstractWorkbook} end

struct WorksheetRow{W <: AbstractWorksheet}
    parent::W
    row_data::st_row_data
    cell_data::Dict{UInt16,st_cell_data}

    function WorksheetRow(parent::W, row_data::st_row_data) where {W <: AbstractWorksheet}
        return new{W}(parent, row_data, Dict{UInt16,st_cell_data}())
    end
end

mutable struct Worksheet{W <: AbstractWorkbook} <: AbstractWorksheet{W}
    parent::W
    sheet_index::Int
    handle::Ptr{xlsWorkSheet}
    rows::st_row
    worksheet_rows::Dict{UInt16,WorksheetRow}

    function Worksheet(parent::W, sheet_index::Int, handle::Ptr{xlsWorkSheet}, rows::st_row) where {W <: AbstractWorkbook}
        new_ws = new{W}(parent, sheet_index, handle, rows, Dict{UInt16,WorksheetRow}())
        finalizer(close, new_ws)
        return new_ws
    end
end

struct WorksheetInfo
    name::String
    isvisible::Bool
end

mutable struct Workbook <: AbstractWorkbook
    handle::Ptr{xlsWorkBook}
    is1904::Bool
    charset::String
    sheets_info::Vector{WorksheetInfo}
    sheetname_index::Dict{String,Int}
    sheets::Dict{Int,Worksheet}

    function Workbook(handle::Ptr{xlsWorkBook}, is1904::Bool, charset::String, sheets_info::Vector{WorksheetInfo}, sheetname_index::Dict{String,Int}, sheets::Dict{Int,Worksheet})
        new_wb = new(handle, is1904, charset, sheets_info, sheetname_index, sheets)
        finalizer(close, new_wb)
        return new_wb
    end
end


text:
println("Dj) Test of the MultiPhotonDeExcitation module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-MultiPhotonDeExcitation.sum")
wa = Atomic.Computation("xx",  Nuclear.Model(26.),
                        initialConfigs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                        intermediateConfigs=[Configuration("1s 2s^2 2p"), Configuration("1s 2s 2p^2") ],
                        finalConfigs  =[Configuration("1s^2 2s^2"), Configuration("1s^2 2s 2p") ],
                        process = JAC.MultiPhotonDE,
                        processSettings=MultiPhotonDeExcitation.Settings()  )
#

wb = perform(wa)
setDefaults("print summary: close", "")


text:
Ω = Circle(-1,1) × Circle(-1,1) × Circle(-1,1)
dims = ndims(Ω)
ex, ey, ez = (3,4,5) .* 10
(npx, npy, npz) = (1,2,3) # needs to be at least these polynomial orders due to last test
ClimateMachine.gpu_allowscalar(true)
# Define Grid: might want to loop over element sizes and polynomial orders
grid = DiscontinuousSpectralElementGrid(Ω, elements = (ex, ey, ez), polynomialorder = (npx, npy, npz), array = ArrayType)
gridhelper = GridHelper(grid)
addup(xC, tol) = sum(abs.(xC[1] .- xC) .≤ tol)
npx, npy, npz = polynomialorders(grid)
x, y, z = coordinates(grid)
xC, yC, zC = cellcenters(grid)

ne = size(x)[2]
ex = round(Int64, ne / addup(xC, 10^4 * eps(maximum(abs.(x)))))
ey = round(Int64, ne / addup(yC, 10^4 * eps(maximum(abs.(y)))))
ez = round(Int64, ne / addup(zC, 10^4 * eps(maximum(abs.(z)))))
check = ne == ex * ey * ez
check ? true : error("improper counting")
p = getperm(xC, yC, zC, ex, ey, ez)
lin = reshape(collect(1:length(xC)), (ex, ey, ez))

ne = size(x)[2]
rx, wx = GaussQuadrature.legendre(npx+1, both)
ωx = baryweights(rx)
ry, wy = GaussQuadrature.legendre(npy+1, both)
ωy = baryweights(ry)
rz, wz = GaussQuadrature.legendre(npz+1, both)
ωz = baryweights(rz)
linlocal = reshape(collect(1:(npx+1)*(npy+1)*(npz+1)), (npx+1, npy+1, npz+1))

# computation

# now test
@testset "Check getvalue Interpolation" begin
    f = copy(z)
    fl = reshape(f, (npx+1, npy+1, npz+1, ex*ey*ez))
    location = (0.0, 0.0, -0.3)
    @test location[3] ≈ getvalue(fl, xC, yC, zC, location, p, lin, linlocal, x, y, z, rx, ry, rz, ωx, ωy, ωz)
    @. f = x
    fl = reshape(f, (npx+1, npy+1, npz+1, ex*ey*ez))
    location = (0.1, 0.0, -0.3)
    @test location[1] ≈ getvalue(fl, xC, yC, zC, location, p, lin, linlocal, x, y, z, rx, ry, rz, ωx, ωy, ωz)
    @. f = y
    fl = reshape(f, (npx+1, npy+1, npz+1, ex*ey*ez))
    location = (0.1, -0.9, -0.3)
    @test location[2] ≈ getvalue(fl, xC, yC, zC, location, p, lin, linlocal, x, y, z, rx, ry, rz, ωx, ωy, ωz)
end
##
@testset "Check GridHelper Interpolation" begin
    f = copy(z)
    location = (0.0, 0.0, -0.3)
    @test location[3] ≈ getvalue(f, location, gridhelper)
    @. f = x
    location = (0.1, 0.0, -0.3)
    @test location[1] ≈ getvalue(f, location, gridhelper)
    @. f = y
    location = (0.1, -0.9, -0.3)
    @test location[2] ≈ getvalue(f, location, gridhelper)
    @. f = x + y^2 + z^3
    location = (0.1, -0.9, -0.3)
    value = location[1] + location[2]^2 + location[3]^3
    @test value ≈ getvalue(f, location, gridhelper)
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks