text:
_and() = throw(ArgumentError("at least one condition must be passed"))
_and(x::Bool) = x
_and(x::Bool, y::Bool...) = x && _and(y...)

function _and(x::Any...)
    loc = findfirst(x -> !(x isa Bool), x)
    # we know x has positive length and must contain non-boolean
    @assert !isnothing(loc)
    xv = x[loc]
    if ismissing(xv)
        throw(ArgumentError("missing was returned in condition number $loc " *
                            "but only true or false are allowed; pass " *
                            "skipmissing=true to skip missing values"))
    else
        throw(ArgumentError("value $xv was returned in condition number $loc " *
                            "but only true or false are allowed"))
    end
end

_and_missing() = throw(ArgumentError("at least one condition must be passed"))
_and_missing(x::Bool) = x
_and_missing(x::Bool, y::Union{Bool, Missing}...) = x && _and_missing(y...)
_and_missing(x::Missing, y::Union{Bool, Missing}...) = false

function _and_missing(x::Any...)
    loc = findfirst(x -> !(x isa Union{Bool, Missing}), x)
    # we know x has positive length and must contain non-boolean
    @assert !isnothing(loc)
    xv = x[loc]
    throw(ArgumentError("value $xv was returned in condition number $loc" *
                        "but only true, false, or missing are allowed"))
end

# we are guaranteed that ByRow returns a vector
# this workaround is needed for 0-argument ByRow
assert_bool_vec(fun::ByRow) = fun

function assert_bool_vec(@nospecialize(fun))
    return function(x...)
        val = fun(x...)
        if !(val isa AbstractVector)
            throw(ArgumentError("functions passed to `subset` must return an AbstractVector."))
        end
        return val
    end
end

# Note that _get_subset_conditions will have a large compilation time
# if more than 32 conditions are passed as `args`.
function _get_subset_conditions(df::Union{AbstractDataFrame, GroupedDataFrame},
                                (args,)::Ref{Any}, skipmissing::Bool)
    cs_vec = []
    for v in map(x -> broadcast_pair(df isa GroupedDataFrame ? parent(df) : df, x), args)
        if v isa AbstractVecOrMat{<:Pair}
            append!(cs_vec, v)
        elseif v isa MultiColumnIndex
            append!(cs_vec, names(df, v))
        else
            push!(cs_vec, v)
        end
    end

    # subset allows a transformation specification without a target column name or a column
    conditions = Any[if a isa ColumnIndex
                         a => Symbol(:x, i)
                     elseif a isa Pair{<:Any, <:Base.Callable}
                         first(a) => assert_bool_vec(last(a)) => Symbol(:x, i)
                     else
                         throw(ArgumentError("condition specifier $a is not supported by `subset`"))
                     end for (i, a) in enumerate(cs_vec)]
    isempty(conditions) && throw(ArgumentError("at least one condition must be passed"))

    if df isa AbstractDataFrame
        df_conditions = select(df, conditions..., copycols=!(df isa DataFrame))
    else
        df_conditions = select(df, conditions...,
                               copycols=!(parent(df) isa DataFrame), keepkeys=false)
    end

    @assert ncol(df_conditions) == length(conditions)

    if skipmissing
        cond = _and_missing.(eachcol(df_conditions)...)
    else
        cond = _and.(eachcol(df_conditions)...)
    end

    @assert eltype(cond) === Bool
    return cond
end

"""
    subset(df::AbstractDataFrame, args...; skipmissing::Bool=false, view::Bool=false)
    subset(gdf::GroupedDataFrame, args...; skipmissing::Bool=false, view::Bool=false,
           ungroup::Bool=true)

Return a copy of data frame `df` or parent of `gdf` containing only rows for
which all values produced by transformation(s) `args` for a given row are
`true`. All transformations must produce vectors containing `true` or `false`
(and optionally `missing` if `skipmissing=true`).

Each argument passed in `args` can be any specifier following the rules
described for [`select`](@ref) with the restriction that:
* specifying target column name is not allowed as `subset` does not create new
  columns;
* every passed transformation must return a scalar or a vector (returning
  `AbstractDataFrame`, `NamedTuple`, `DataFrameRow` or `AbstractMatrix`
  is not supported).

If `skipmissing=false` (the default) `args` are required to produce vectors
containing only `Bool` values. If `skipmissing=true`, additionally `missing` is
allowed and it is treated as `false` (i.e. rows for which one of the conditions
returns `missing` are skipped).

If `view=true` a `SubDataFrame` view  is returned instead of a `DataFrame`.

If `ungroup=false` the resulting data frame is re-grouped based on the same
grouping columns as `gdf` and a `GroupedDataFrame` is returned.

If a `GroupedDataFrame` is passed then it must include all groups present in the
`parent` data frame, like in [`select!`](@ref).

!!! note

    Note that as the `subset` function works in exactly the same way as other
    transformation functions defined in DataFrames.jl this is the preferred way to
    subset rows of a data frame or grouped data frame. In particular it uses a
    different set of rules for specifying transformations than [`filter`](@ref)
    which is implemented in DataFrames.jl to ensure support for the
    standard Julia API for collections.

See also: [`subset!`](@ref), [`filter`](@ref), [`select`](@ref)

# Examples

```jldoctest
julia> df = DataFrame(id=1:4, x=[true, false, true, false],
                      y=[true, true, false, false],
                      z=[true, true, missing, missing], v=[1, 2, 11, 12])
4×5 DataFrame
 Row │ id     x      y      z        v
     │ Int64  Bool   Bool   Bool?    Int64
─────┼─────────────────────────────────────
   1 │     1   true   true     true      1
   2 │     2  false   true     true      2
   3 │     3   true  false  missing     11
   4 │     4  false  false  missing     12

julia> subset(df, :x)
2×5 DataFrame
 Row │ id     x     y      z        v
     │ Int64  Bool  Bool   Bool?    Int64
─────┼────────────────────────────────────
   1 │     1  true   true     true      1
   2 │     3  true  false  missing     11

julia> subset(df, :v => x -> x .> 3)
2×5 DataFrame
 Row │ id     x      y      z        v
     │ Int64  Bool   Bool   Bool?    Int64
─────┼─────────────────────────────────────
   1 │     3   true  false  missing     11
   2 │     4  false  false  missing     12

julia> subset(df, :x, :y => ByRow(!))
1×5 DataFrame
 Row │ id     x     y      z        v
     │ Int64  Bool  Bool   Bool?    Int64
─────┼────────────────────────────────────
   1 │     3  true  false  missing     11

julia> subset(df, :x, :z, skipmissing=true)
1×5 DataFrame
 Row │ id     x     y     z      v
     │ Int64  Bool  Bool  Bool?  Int64
─────┼─────────────────────────────────
   1 │     1  true  true   true      1

julia> subset(df, :x, :z)
ERROR: ArgumentError: missing was returned in condition number 2 but only true or false are allowed; pass skipmissing=true to skip missing values

julia> subset(groupby(df, :y), :v => x -> x .> minimum(x))
2×5 DataFrame
 Row │ id     x      y      z        v
     │ Int64  Bool   Bool   Bool?    Int64
─────┼─────────────────────────────────────
   1 │     2  false   true     true      2
   2 │     4  false  false  missing     12
```
"""
function subset(df::AbstractDataFrame, @nospecialize(args...);
                skipmissing::Bool=false, view::Bool=false)
    row_selector = _get_subset_conditions(df, Ref{Any}(args), skipmissing)
    return view ? Base.view(df, row_selector, :) : df[row_selector, :]
end

function subset(gdf::GroupedDataFrame, @nospecialize(args...);
                skipmissing::Bool=false, view::Bool=false,
                        ungroup::Bool=true)
    row_selector = _get_subset_conditions(gdf, Ref{Any}(args), skipmissing)
    df = parent(gdf)
    res = view ? Base.view(df, row_selector, :) : df[row_selector, :]
    # TODO: in some cases it might be faster to groupby gdf.groups[row_selector]
    return ungroup ? res : groupby(res, groupcols(gdf))
end

"""
    subset!(df::AbstractDataFrame, args...; skipmissing::Bool=false)
    subset!(gdf::GroupedDataFrame{DataFrame}, args..., skipmissing::Bool=false,
            ungroup::Bool=true)

Update data frame `df` or the parent of `gdf` in place to contain only rows for
which all values produced by transformation(s) `args` for a given row is `true`.
All transformations must produce vectors containing `true` or `false` (and
optionally `missing` if `skipmissing=true`).

Each argument passed in `args` can be any specifier following the rules
described for [`select`](@ref) with the restriction that:
* specifying target column name is not allowed as `subset!` does not create new
  columns;
* every passed transformation must return a scalar or a vector (returning
  `AbstractDataFrame`, `NamedTuple`, `DataFrameRow` or `AbstractMatrix`
  is not supported).

If `skipmissing=false` (the default) `args` are required to produce vectors
containing only `Bool` values. If `skipmissing=true`, additionally `missing` is
allowed and it is treated as `false` (i.e. rows for which one of the conditions
returns `missing` are skipped).

If `ungroup=false` the resulting data frame is re-grouped based on the same
grouping columns as `gdf` and a `GroupedDataFrame` is returned.

If `GroupedDataFrame` is subsetted then it must include all groups present in
the `parent` data frame, like in [`select!`](@ref). In this case the passed
`GroupedDataFrame` is updated to have correct groups after its parent is
updated.

!!! note

    Note that as the `subset!` function works in exactly the same way as other
    transformation functions defined in DataFrames.jl this is the preferred way to
    subset rows of a data frame or grouped data frame. In particular it uses a
    different set of rules for specifying transformations than [`filter!`](@ref)
    which is implemented in DataFrames.jl to ensure support for the
    standard Julia API for collections.

See also: [`subset`](@ref), [`filter!`](@ref), [`select!`](@ref)

# Examples

```jldoctest
julia> df = DataFrame(id=1:4, x=[true, false, true, false], y=[true, true, false, false])
4×3 DataFrame
 Row │ id     x      y
     │ Int64  Bool   Bool
─────┼─────────────────────
   1 │     1   true   true
   2 │     2  false   true
   3 │     3   true  false
   4 │     4  false  false

julia> subset!(df, :x, :y => ByRow(!));

julia> df
1×3 DataFrame
 Row │ id     x     y
     │ Int64  Bool  Bool
─────┼────────────────────
   1 │     3  true  false

julia> df = DataFrame(id=1:4, y=[true, true, false, false], v=[1, 2, 11, 12]);

julia> subset!(groupby(df, :y), :v => x -> x .> minimum(x));

julia> df
2×3 DataFrame
 Row │ id     y      v
     │ Int64  Bool   Int64
─────┼─────────────────────
   1 │     2   true      2
   2 │     4  false     12

julia> df = DataFrame(id=1:4, x=[true, false, true, false],
                      z=[true, true, missing, missing], v=1:4)
4×4 DataFrame
 Row │ id     x      z        v
     │ Int64  Bool   Bool?    Int64
─────┼──────────────────────────────
   1 │     1   true     true      1
   2 │     2  false     true      2
   3 │     3   true  missing      3
   4 │     4  false  missing      4

julia> subset!(df, :x, :z)
ERROR: ArgumentError: missing was returned in condition number 2 but only true or false are allowed; pass skipmissing=true to skip missing values

julia> subset!(df, :x, :z, skipmissing=true);

julia> df
1×4 DataFrame
 Row │ id     x     z      v
     │ Int64  Bool  Bool?  Int64
─────┼───────────────────────────
   1 │     1  true   true      1

julia> df = DataFrame(id=1:4, x=[true, false, true, false], y=[true, true, false, false],
                      z=[true, true, missing, missing], v=[1, 2, 11, 12]);

julia> subset!(groupby(df, :y), :v => x -> x .> minimum(x));

julia> df
2×5 DataFrame
 Row │ id     x      y      z        v
     │ Int64  Bool   Bool   Bool?    Int64
─────┼─────────────────────────────────────
   1 │     2  false   true     true      2
   2 │     4  false  false  missing     12
```
"""
function subset!(df::AbstractDataFrame, @nospecialize(args...); skipmissing::Bool=false)
    row_selector = _get_subset_conditions(df, Ref{Any}(args), skipmissing)
    return deleteat!(df, findall(!, row_selector))
end

function subset!(gdf::GroupedDataFrame, @nospecialize(args...); skipmissing::Bool=false,
                 ungroup::Bool=true)
    ngroups = length(gdf)
    groups = gdf.groups
    lazy_lock = gdf.lazy_lock
    row_selector = _get_subset_conditions(gdf, Ref{Any}(args), skipmissing)
    df = parent(gdf)
    res = deleteat!(df, findall(!, row_selector))
    if nrow(res) == length(groups) # we have not removed any rows
        return ungroup ? res : gdf
    end
    newgroups = groups[row_selector]

    # TODO: add threading support
    seen = fill(false, ngroups)
    @inbounds for gix in newgroups
        @assert gix > 0 # having dropped groups in gdf is not allowed here
        seen[gix] = true
    end

    if sum(seen) < ngroups # subset! has dropped some groups
        remap = cumsum(seen)
        @inbounds for i in eachindex(newgroups)
            newgroups[i] = remap[newgroups[i]]
        end
        ngroups = remap[end]
    end

    # update GroupedDataFrame indices in a thread safe way
    Threads.lock(lazy_lock) do
        gdf.groups = newgroups
        gdf.idx = nothing
        gdf.starts = nothing
        gdf.ends = nothing
        gdf.ngroups = ngroups
        gdf.keymap = nothing
    end
    return ungroup ? res : gdf
end


text:
#=
Created on 09/01/2021 13:08:20
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Integrating for dummies. Compute the Riemann sum **without** making use of a for loop.

Compute the area of a circle.
=#

function riemannsum(f, a, b; n=100)
    dx = (b - a) / n
    return sum(f.(a:dx:b)) * dx
end

riemannsum(sin, 0, 2pi)

riemannsum(x->x*sin(x), 0, 2pi)

4riemannsum(x->(sqrt(1-x^2)), 0, 1, n=1000)

ext: jl


text:
n = 0.0

for i = 1:10
  global n
  n = n + 0.1
end

println(n)

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
@plottest begin
    # simulation parameters
    T = 1000

    # create three SRM0 neurons
    η₀ = 5.0
    τᵣ = 1.0
    vth = 1.0

    rng = StableRNG(123)

    # create population
    weights = Float32[ 0  5;
                    0  0]
    pop = Population(weights; cell = () -> SRM0(η₀, τᵣ),
                            synapse = Synapse.Alpha,
                            threshold = () -> Threshold.Ideal(vth),
                            learner = STDP(0.5, 0.5, size(weights, 1)))

    # create step input currents
    input = InputPopulation([ConstantRate(0.8; rng = rng)])

    # create network
    net = Network(Dict([:input => input, :pop => pop]))
    connect!(net, :input, :pop; weights = [1 0], synapse = Synapse.Alpha)

    # simulate
    w = Float64[]
    voltages = Dict([(i, Float64[]) for i in 1:2])
    cb = () -> begin
        push!(w, net[:pop].weights[1, 2])
        for id in 1:size(pop)
            push!(voltages[id], getvoltage(pop[id]))
        end
    end
    outputs = simulate!(net, T; cb = cb, dense = true)

    weight_plot = plot(1:T, w, label = "")
    title!("Synaptic Weights Over Simulation")
    xlabel!("Time (sec)")
    ylabel!("Weight")

    raster_plot = rasterplot(outputs[:pop])
    title!("Raster Plot")
    xlabel!("Time (sec)")

    plot(voltages[1], label = "Input")
    voltage_plot = plot!(voltages[2], label = "Neuron")
    title!("Membrane Voltages")
    xlabel!("Time (sec)")
    ylabel!("Voltage (V)")

    plot(weight_plot, raster_plot, voltage_plot, layout = grid(3, 1))
end joinpath(datadir, "Stdp-Test.png") !isci


text:
@testset "minimal_mulRevToPair_test" begin
    @test mul_rev_to_pair(emptyinterval(), interval(1.0, 2.0))[1] == emptyinterval() && mul_rev_to_pair(emptyinterval(), interval(1.0, 2.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(1.0, 2.0), emptyinterval())[1] == emptyinterval() && mul_rev_to_pair(interval(1.0, 2.0), emptyinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(emptyinterval(), emptyinterval())[1] == emptyinterval() && mul_rev_to_pair(emptyinterval(), emptyinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-2.1, -0.4))[1] == interval(0x1.999999999999AP-3, 0x1.5P+4) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-2.1, -0.4))[1] == interval(0x1.999999999999AP-3, Inf) && mul_rev_to_pair(interval(-2.0, 0.0), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-2.1, -0.4))[1] == interval(-Inf, -0x1.745D1745D1745P-2) && mul_rev_to_pair(interval(-2.0, 1.1), interval(-2.1, -0.4))[2] == interval(0x1.999999999999AP-3, Inf)
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-2.1, -0.4))[1] == interval(-Inf, -0x1.745D1745D1745P-2) && mul_rev_to_pair(interval(0.0, 1.1), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-2.1, -0.4))[1] == interval(-0x1.A400000000001P+7, -0x1.745D1745D1745P-2) && mul_rev_to_pair(interval(0.01, 1.1), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-2.1, -0.4))[1] == emptyinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-2.1, -0.4))[1] == interval(0.0, 0x1.5P+4) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-2.1, -0.4))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(-Inf, 0.0), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-2.1, -0.4))[1] == interval(-Inf, -0x1.745D1745D1745P-2) && mul_rev_to_pair(interval(-Inf, 1.1), interval(-2.1, -0.4))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-2.1, -0.4))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-2.0, Inf), interval(-2.1, -0.4))[2] == interval(0x1.999999999999AP-3, Inf)
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-2.1, -0.4))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(0.0, Inf), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-2.1, -0.4))[1] == interval(-0x1.A400000000001P+7, 0.0) && mul_rev_to_pair(interval(0.01, Inf), interval(-2.1, -0.4))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-2.1, -0.4))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(entireinterval(), interval(-2.1, -0.4))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-2.1, 0.0))[1] == interval(0.0, 0x1.5P+4) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-2.1, 0.0))[1] == interval(-0x1.A400000000001P+7, 0.0) && mul_rev_to_pair(interval(0.01, 1.1), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-2.1, 0.0))[1] == interval(0.0, 0x1.5P+4) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-2.1, 0.0))[1] == interval(-0x1.A400000000001P+7, 0.0) && mul_rev_to_pair(interval(0.01, Inf), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-2.1, 0.0))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(-2.1, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-2.1, 0.12))[1] == interval(-0x1.3333333333333P+0, 0x1.5P+4) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-2.1, 0.12))[1] == interval(-0x1.A400000000001P+7 , 0x1.8P+3) && mul_rev_to_pair(interval(0.01, 1.1), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-2.1, 0.12))[1] == interval(-0x1.3333333333333P+0, 0x1.5P+4) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-2.1, 0.12))[1] == interval(-0x1.A400000000001P+7 , 0x1.8P+3) && mul_rev_to_pair(interval(0.01, Inf), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-2.1, 0.12))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(-2.1, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(0.0, 0.12))[1] == interval(-0x1.3333333333333P+0, 0.0) && mul_rev_to_pair(interval(-2.0, -0.1), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(0.0, 0.12))[1] == interval(0.0, 0x1.8P+3) && mul_rev_to_pair(interval(0.01, 1.1), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(0.0, 0.12))[1] == interval(-0x1.3333333333333P+0, 0.0) && mul_rev_to_pair(interval(-Inf, -0.1), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(0.0, 0.12))[1] == interval(0.0, 0x1.8P+3) && mul_rev_to_pair(interval(0.01, Inf), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(0.0, 0.12))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(0.0, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(0.01, 0.12))[1] == interval(-0x1.3333333333333P+0, -0x1.47AE147AE147BP-8) && mul_rev_to_pair(interval(-2.0, -0.1), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(0.01, 0.12))[1] == interval(-Inf, -0x1.47AE147AE147BP-8) && mul_rev_to_pair(interval(-2.0, 0.0), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(0.01, 0.12))[1] == interval(-Inf, -0x1.47AE147AE147BP-8) && mul_rev_to_pair(interval(-2.0, 1.1), interval(0.01, 0.12))[2] == interval(0x1.29E4129E4129DP-7, Inf)
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(0.01, 0.12))[1] == interval(0x1.29E4129E4129DP-7, Inf) && mul_rev_to_pair(interval(0.0, 1.1), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(0.01, 0.12))[1] == interval(0x1.29E4129E4129DP-7, 0x1.8P+3) && mul_rev_to_pair(interval(0.01, 1.1), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(0.01, 0.12))[1] == emptyinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(0.01, 0.12))[1] == interval(-0x1.3333333333333P+0, 0.0) && mul_rev_to_pair(interval(-Inf, -0.1), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(0.01, 0.12))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-Inf, 0.0), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(0.01, 0.12))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-Inf, 1.1), interval(0.01, 0.12))[2] == interval(0x1.29E4129E4129DP-7, Inf)
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(0.01, 0.12))[1] == interval(-Inf, -0x1.47AE147AE147BP-8) && mul_rev_to_pair(interval(-2.0, Inf), interval(0.01, 0.12))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(0.0, Inf), interval(0.01, 0.12))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(0.0, Inf), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(0.01, 0.12))[1] == interval(0.0, 0x1.8P+3) && mul_rev_to_pair(interval(0.01, Inf), interval(0.01, 0.12))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(0.01, 0.12))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(entireinterval(), interval(0.01, 0.12))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(0.0, 0.0))[1] == interval(0.0, 0.0) && mul_rev_to_pair(interval(-2.0, -0.1), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(0.0, 0.0))[1] == interval(0.0, 0.0) && mul_rev_to_pair(interval(0.01, 1.1), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(0.0, 0.0))[1] == interval(0.0, 0.0) && mul_rev_to_pair(interval(-Inf, -0.1), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(0.0, 0.0))[1] == interval(0.0, 0.0) && mul_rev_to_pair(interval(0.01, Inf), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(0.0, 0.0))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(0.0, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-Inf, -0.1))[1] == interval(0x1.999999999999AP-5, Inf) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-Inf, -0.1))[1] == interval(0x1.999999999999AP-5 , Inf) && mul_rev_to_pair(interval(-2.0, 0.0), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-Inf, -0.1))[1] == interval(-Inf , -0x1.745D1745D1745P-4) && mul_rev_to_pair(interval(-2.0, 1.1), interval(-Inf, -0.1))[2] == interval(0x1.999999999999AP-5 , Inf)
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-Inf, -0.1))[1] == interval(-Inf, -0x1.745D1745D1745P-4) && mul_rev_to_pair(interval(0.0, 1.1), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-Inf, -0.1))[1] == interval(-Inf, -0x1.745D1745D1745P-4) && mul_rev_to_pair(interval(0.01, 1.1), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-Inf, -0.1))[1] == emptyinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-Inf, -0.1))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-Inf, -0.1))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(-Inf, 0.0), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-Inf, -0.1))[1] == interval(-Inf, -0x1.745D1745D1745P-4) && mul_rev_to_pair(interval(-Inf, 1.1), interval(-Inf, -0.1))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-Inf, -0.1))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-2.0, Inf), interval(-Inf, -0.1))[2] == interval(0x1.999999999999AP-5 , Inf)
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-Inf, -0.1))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(0.0, Inf), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-Inf, -0.1))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(0.01, Inf), interval(-Inf, -0.1))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-Inf, -0.1))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(entireinterval(), interval(-Inf, -0.1))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-Inf, 0.0))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-Inf, 0.0))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(0.01, 1.1), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-Inf, 0.0))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-Inf, 0.0))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(0.01, Inf), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-Inf, 0.0))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(-Inf, 0.0))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-Inf, 0.3))[1] == interval(-0x1.8P+1, Inf) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-Inf, 0.3))[1] == interval(-Inf, 0x1.EP+4) && mul_rev_to_pair(interval(0.01, 1.1), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-Inf, 0.3))[1] == interval(-0x1.8P+1, Inf) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-Inf, 0.3))[1] == interval(-Inf, 0x1.EP+4) && mul_rev_to_pair(interval(0.01, Inf), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-Inf, 0.3))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(-Inf, 0.3))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(-0.21, Inf))[1] == interval(-Inf , 0x1.0CCCCCCCCCCCDP+1) && mul_rev_to_pair(interval(-2.0, -0.1), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(-0.21, Inf))[1] == interval(-0x1.5P+4, Inf) && mul_rev_to_pair(interval(0.01, 1.1), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(-0.21, Inf))[1] == interval(-Inf, 0x1.0CCCCCCCCCCCDP+1) && mul_rev_to_pair(interval(-Inf, -0.1), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(-0.21, Inf))[1] == interval(-0x1.5P+4, Inf) && mul_rev_to_pair(interval(0.01, Inf), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(-0.21, Inf))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(-0.21, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(0.0, Inf))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-2.0, -0.1), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(0.0, Inf))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(0.01, 1.1), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(0.0, Inf))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-Inf, -0.1), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(0.0, Inf))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(0.01, Inf), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(0.0, Inf))[1] == entireinterval() && mul_rev_to_pair(entireinterval(), interval(0.0, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, -0.1), interval(0.04, Inf))[1] == interval(-Inf, -0x1.47AE147AE147BP-6) && mul_rev_to_pair(interval(-2.0, -0.1), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), interval(0.04, Inf))[1] == interval(-Inf, -0x1.47AE147AE147BP-6) && mul_rev_to_pair(interval(-2.0, 0.0), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), interval(0.04, Inf))[1] == interval(-Inf, -0x1.47AE147AE147BP-6) && mul_rev_to_pair(interval(-2.0, 1.1), interval(0.04, Inf))[2] == interval(0x1.29E4129E4129DP-5, Inf)
    @test mul_rev_to_pair(interval(0.0, 1.1), interval(0.04, Inf))[1] == interval(0x1.29E4129E4129DP-5, Inf) && mul_rev_to_pair(interval(0.0, 1.1), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), interval(0.04, Inf))[1] == interval(0x1.29E4129E4129DP-5, Inf) && mul_rev_to_pair(interval(0.01, 1.1), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), interval(0.04, Inf))[1] == emptyinterval() && mul_rev_to_pair(interval(0.0, 0.0), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), interval(0.04, Inf))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-Inf, -0.1), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), interval(0.04, Inf))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-Inf, 0.0), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), interval(0.04, Inf))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(interval(-Inf, 1.1), interval(0.04, Inf))[2] == interval(0x1.29E4129E4129DP-5, Inf)
    @test mul_rev_to_pair(interval(-2.0, Inf), interval(0.04, Inf))[1] == interval(-Inf, -0x1.47AE147AE147BP-6) && mul_rev_to_pair(interval(-2.0, Inf), interval(0.04, Inf))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(0.0, Inf), interval(0.04, Inf))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(0.0, Inf), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), interval(0.04, Inf))[1] == interval(0.0, Inf) && mul_rev_to_pair(interval(0.01, Inf), interval(0.04, Inf))[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), interval(0.04, Inf))[1] == interval(-Inf, 0.0) && mul_rev_to_pair(entireinterval(), interval(0.04, Inf))[2] == interval(0.0, Inf)
    @test mul_rev_to_pair(interval(-2.0, -0.1), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, -0.1), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 0.0), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 0.0), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, 1.1), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, 1.1), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 1.1), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 1.1), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, 1.1), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(0.01, 1.1), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, 0.0), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(0.0, 0.0), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, -0.1), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, -0.1), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 0.0), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 0.0), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-Inf, 1.1), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-Inf, 1.1), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(-2.0, Inf), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(-2.0, Inf), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.0, Inf), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(0.0, Inf), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(interval(0.01, Inf), entireinterval())[1] == entireinterval() && mul_rev_to_pair(interval(0.01, Inf), entireinterval())[2] == emptyinterval()
    @test mul_rev_to_pair(entireinterval(), entireinterval())[1] == entireinterval() && mul_rev_to_pair(entireinterval(), entireinterval())[2] == emptyinterval()
end


text:
# Predefined loss functions
# You may also implement your own loss by subtyping the abstract type Loss.
#
# Losses must have the following:
#   Fields:
#     `scale::Float64`
#           This field represents a scalar weight assigned to the loss function: w*l(u,a)
#     `domain::natural_Domain`
#           The "natural" domain that the loss function was meant to handle. E.g. BoolDomain for LogisticLoss,
#           RealDomain for QuadLoss, etc.

#   Other fields may be also be included to encode parameters of the loss function, encode the range or
#   set of possible values of the data, etc.
#
#   Methods:
#     `my_loss_type(args..., scale=1.0::Float64;
#                   domain=natural_Domain(args[range]...), kwargs...) ::my_loss_type`
#           Constructor for the loss type. The first few arguments are parameters for
#           which there isn't a rational default (a loss may not need any of these).
#           The last positional argument should be the scale, which should default to 1.
#           There must be a default domain which is a Domain, which may take arguments from
#           the list of positional arguments. Parameters besides the scale for which there are
#           reasonable defaults should be included as keyword arguments (there may be none).
#     `evaluate(l::my_loss_type, u::Float64, a::Number) ::Float64`
#           Evaluates the function l(u,a) where u is the approximation of a
#     `grad(l::my_loss_type, u::Float64, a::Number) ::Float64`
#           Evaluates the gradient of the loss at the given point (u,a)

#   In addition, loss functions should preferably implement methods:
#     `M_estimator(l::my_loss_type, a::AbstractArray) ::Float64`
#           Finds uₒ = argmin ∑l(u,aᵢ) which is the best single estimate of the array `a`
#           If `M_estimator` is not implemented, a live optimization procedure will be used when this function is
#           called in order to compute loss function scalings. The live optimization may be slow, so an analytic
#           implementation is preferable.
#     `impute(d::Domain, l::my_loss_type, u::Array{Float64})` (in impute_and_err.jl)
#           Finds a = argmin l(u,a), the most likely value for an observation given a parameter u

import Base: *, convert
import Optim: optimize, LBFGS
export Loss,
       DiffLoss, ClassificationLoss, SingleDimLoss, # categories of Losses
       QuadLoss, L1Loss, HuberLoss, QuantileLoss, # losses for predicting reals
       KullbackLeiblerLoss, #  loss for exponential family distributions
       PoissonLoss, # losses for predicting integers
       HingeLoss, WeightedHingeLoss, LogisticLoss, # losses for predicting booleans
       OrdinalHingeLoss, OrdisticLoss, MultinomialOrdinalLoss, BvSLoss, # losses for predicting ordinals
       MultinomialLoss, OvALoss, # losses for predicting nominals (categoricals)
       PeriodicLoss, # losses for predicting periodic variables
       evaluate, grad, M_estimator, # methods on losses
       avgerror, scale, mul!, *,
       embedding_dim, get_yidxs, datalevels, domain

abstract type Loss end
# a DiffLoss is one in which l(u,a) = f(u-a) AND argmin f(x) = 0
# for example, QuadLoss(u,a)=(u-a)² and we can write f(x)=x² and x=u-a
abstract type DiffLoss<:Loss end
# a ClassificationLoss is one in which observed values are true = 1 or false = 0 = -1 AND argmin_a L(u,a) = u>=0 ? true : false
abstract type ClassificationLoss<:Loss end
# Single Dimensional losses are DiffLosses or ClassificationLosses, which allow optimized evaluate and grad functions
const SingleDimLoss = Union{DiffLoss, ClassificationLoss}

mul!(l::Loss, newscale::Number) = (l.scale = newscale; l)
scale(l::Loss) = l.scale
*(newscale::Number, l::Loss) = (newl = copy(l); mul!(newl, newscale))
*(l::Loss, newscale::Number) = (newl = copy(l); mul!(newl, newscale))

domain(l::Loss) = l.domain

### embedding dimensions: mappings from losses/columns of A to columns of Y

# default number of columns
# number of columns is higher for multidimensional losses
embedding_dim(l::Loss) = 1
embedding_dim(l::Array{LossSubtype,1}) where LossSubtype<:Loss = sum(map(embedding_dim, l))

# find spans of loss functions (for multidimensional losses)
function get_yidxs(losses::Array{LossSubtype,1}) where LossSubtype<:Loss
    n = length(losses)
    ds = map(embedding_dim, losses)
    d = sum(ds)
    featurestartidxs = cumsum(append!([1], ds))
    # find which columns of Y map to which columns of A (for multidimensional losses)
    U = Union{UnitRange{Int}, Int}
    yidxs = Array{U}(undef, n)

    for f = 1:n
        if ds[f] == 1
            yidxs[f] = featurestartidxs[f]
        else
            yidxs[f] = featurestartidxs[f]:featurestartidxs[f]+ds[f]-1
        end
    end
    return yidxs
end

### promote integers to floats if given as the argument u
## causes ambiguity warnings
# evaluate(l::Loss, u::Number, a) = evaluate(l,convert(Float64,u),a)
# grad(l::Loss, u::Number, a) = grad(l,convert(Float64,u),a)
# evaluate{T<:Number}(l::Loss, u::Array{T,1}, a) = evaluate(l,convert(Array{Float64,1},u),a)
# grad{T<:Number}(l::Loss, u::Array{T,1}, a) = grad(l,convert(Array{Float64,1},u),a)

### -1,0,1::Int are translated to Booleans if loss is not defined on numbers
# convert(::Type{Bool}, x::Int) = x==1 ? true : (x==-1 || x==0) ? false : throw(InexactError("Bool method successfully overloaded by LowRankModels"))
myBool(x::Int) = x==1 ? true : (x==-1 || x==0) ? false : throw(InexactError())
evaluate(l::ClassificationLoss, u::Float64, a::Int) = evaluate(l,u,myBool(a))
grad(l::ClassificationLoss, u::Float64, a::Int) = grad(l,u,myBool(a))
M_estimator(l::ClassificationLoss, a::AbstractArray{Int,1}) = M_estimator(l,myBool(a))

### M-estimators

# The following is the M-estimator for loss functions that don't have one defined. It's also useful
# for checking that the analytic M_estimators are correct. To make sure this method is called instead
# of the loss-specific method (should only be done to test), simply pass the third paramter `test`.
# e.g. M_estimator(l,a) will call the implementation for l, but M_estimator(l,a,"test") will call the
# general-purpose optimizing M_estimator.
function M_estimator(l::Loss, a::AbstractArray; test="test")
    # the function to optimize over
    f = (u -> sum(map(ai->evaluate(l,u[1],ai), a))) # u is indexed because `optim` assumes input is a vector
    # the gradient of that function
    function g!(storage::Vector, u::Vector) # this is the format `optim` expects
        storage[1] = sum(map(ai->grad(l,u[1],ai), a))
    end
    m = optimize(f, g!, [median(a)], LBFGS()).minimum[1]
end

# Uses uₒ = argmin ∑l(u,aᵢ) to find (1/n)*∑l(uₒ,aᵢ) which is the
# average error incurred by using the estimate uₒ for every aᵢ
function avgerror(l::Loss, a::AbstractArray)
    b = collect(skipmissing(a))
    m = M_estimator(l,b)
    sum(map(ai->evaluate(l,m,ai),b))/length(b)
end

## Losses:

########################################## QUADRATIC ##########################################
# f: ℜxℜ -> ℜ
mutable struct QuadLoss<:DiffLoss
    scale::Float64
    domain::Domain
end
QuadLoss(scale=1.0::Float64; domain=RealDomain()) = QuadLoss(scale, domain)

evaluate(l::QuadLoss, u::Float64, a::Number) = l.scale*(u-a)^2

grad(l::QuadLoss, u::Float64, a::Number) = 2*(u-a)*l.scale

M_estimator(l::QuadLoss, a::AbstractArray) = mean(a)

########################################## L1 ##########################################
# f: ℜxℜ -> ℜ
mutable struct L1Loss<:DiffLoss
    scale::Float64
    domain::Domain
end
L1Loss(scale=1.0::Float64; domain=RealDomain()) = L1Loss(scale, domain)

evaluate(l::L1Loss, u::Float64, a::Number) = l.scale*abs(u-a)

grad(l::L1Loss, u::Float64, a::Number) = sign(u-a)*l.scale

M_estimator(l::L1Loss, a::AbstractArray) = median(a)

########################################## HUBER ##########################################
# f: ℜxℜ -> ℜ
mutable struct HuberLoss<:DiffLoss
    scale::Float64
    domain::Domain
    crossover::Float64 # where QuadLoss loss ends and linear loss begins; =1 for standard HuberLoss
end
HuberLoss(scale=1.0::Float64; domain=RealDomain(), crossover=1.0::Float64) = HuberLoss(scale, domain, crossover)

function evaluate(l::HuberLoss, u::Float64, a::Number)
    abs(u-a) > l.crossover ? (abs(u-a) - l.crossover + l.crossover^2)*l.scale : (u-a)^2*l.scale
end

grad(l::HuberLoss,u::Float64,a::Number) = abs(u-a)>l.crossover ? sign(u-a)*l.scale : (u-a)*l.scale

M_estimator(l::HuberLoss, a::AbstractArray) = median(a) # a heuristic, not the true estimator

########################################## QUANTILE ##########################################
# f: ℜxℜ -> ℜ
# define (u)_+ = max(u,0), (u)_- = max(-u,0) so (u)_+ + (u)_- = |u|
# f(u,a) = {    quantile (a - u)_+ + (1-quantile) (a - u)_-
# fits the `quantile`th quantile of the distribution
mutable struct QuantileLoss<:DiffLoss
    scale::Float64
    domain::Domain
    quantile::Float64 # fit the alphath quantile
end
QuantileLoss(scale=1.0::Float64; domain=RealDomain(), quantile=.5::Float64) = QuantileLoss(scale, domain, quantile)

function evaluate(l::QuantileLoss, u::Float64, a::Number)
    diff = a-u
    diff > 0  ?  l.scale * l.quantile * diff  :  - l.scale * (1-l.quantile) * diff
end

function grad(l::QuantileLoss,u::Float64,a::Number)
  diff = a-u
  diff > 0  ?  -l.scale * l.quantile  :  l.scale * (1-l.quantile)
end

M_estimator(l::QuantileLoss, a::AbstractArray) = quantile(a, l.quantile)

########################################## PERIODIC ##########################################
# f: ℜxℜ -> ℜ
# f(u,a) = w * (1 - cos((a-u)*(2*pi)/T))
# this measures how far away u and a are on a circle of circumference T.
mutable struct PeriodicLoss<:DiffLoss
    T::Float64 # the length of the period
    scale::Float64
    domain::Domain
end
PeriodicLoss(T, scale=1.0::Float64; domain=PeriodicDomain(T)) = PeriodicLoss(T, scale, domain)

evaluate(l::PeriodicLoss, u::Float64, a::Number) = l.scale*(1-cos((a-u)*(2*pi)/l.T))

grad(l::PeriodicLoss, u::Float64, a::Number) = -l.scale*((2*pi)/l.T)*sin((a-u)*(2*pi)/l.T)

function M_estimator(l::PeriodicLoss, a::AbstractArray{Float64})
    (l.T/(2*pi))*atan( sum(sin(2*pi*a/l.T)) / sum(cos(2*pi*a/l.T)) ) + l.T/2 # not kidding.
    # this is the estimator, and there is a form that works with weighted measurements (aka a prior on a)
    # see: http://www.tandfonline.com/doi/pdf/10.1080/17442507308833101 eq. 5.2
end

########################################## POISSON ##########################################
# f: ℜxℕ -> ℜ
# BEWARE:
# 1) this is a reparametrized poisson: we parametrize the mean as exp(u) so that u can take any real value and still produce a positive mean
# 2) THIS LOSS MAY CAUSE MODEL INSTABLITY AND DIFFICULTY FITTING.
mutable struct PoissonLoss<:Loss
    scale::Float64
    domain::Domain
end
PoissonLoss(max_count=2^31::Int; domain=CountDomain(max_count)::Domain) = PoissonLoss(1.0, domain)

function evaluate(l::PoissonLoss, u::Float64, a::Number)
    l.scale*(exp(u) - a*u + (a==0 ? 0 : a*(log(a)-1))) # log(a!) ~ a==0 ? 0 : a*(log(a)-1)
end

grad(l::PoissonLoss, u::Float64, a::Number) = l.scale*(exp(u) - a)

M_estimator(l::PoissonLoss, a::AbstractArray) = log(mean(a))

########################################## KULLBACK LEIBLER ##########################################
# f: ℜxℕ -> ℜ

mutable struct KullbackLeiblerLoss<:Loss
    scale::Float64
    domain::Domain
end

KullbackLeiblerLoss(scale=1.0::Float64; domain=RealDomain()) = KullbackLeiblerLoss(scale, domain)

# Following conventions found in:
# Elements of Information Theory, Cover and Thomas
function evaluate(l::KullbackLeiblerLoss, u::Float64, a::Number)
    if a == 0
        return 0.0
    end

    return l.scale * (a*log(a/u) - a + u)
end

grad(l::KullbackLeiblerLoss, u::Float64, a::Number) = l.scale*(1 - (a/(1e-8 + u)))

function M_estimator(l::KullbackLeiblerLoss, a::AbstractArray)
    N = length(a)
    a[a .== 0] .= 1.0

    exp.(sum(log.(a))/N)
end

########################################## ORDINAL HINGE ##########################################
# f: ℜx{min, min+1... max-1, max} -> ℜ
mutable struct OrdinalHingeLoss<:Loss
    min::Integer
    max::Integer
    scale::Float64
    domain::Domain
end
OrdinalHingeLoss(m1, m2, scale=1.0::Float64; domain=OrdinalDomain(m1,m2)) = OrdinalHingeLoss(m1,m2,scale,domain)
# this method should never be called directly but is needed to support copying
OrdinalHingeLoss() = OrdinalHingeLoss(1, 10, 1.0, OrdinalDomain(1,10))
OrdinalHingeLoss(m2) = OrdinalHingeLoss(1, m2, 1.0, OrdinalDomain(1, m2))

function evaluate(l::OrdinalHingeLoss, u::Float64, a::Number)
    #a = round(a)
    if u > l.max-1
        # number of levels higher than true level
        n = min(floor(u), l.max-1) - a
        loss = n*(n+1)/2 + (n+1)*(u-l.max+1)
    elseif u > a
        # number of levels higher than true level
        n = min(floor(u), l.max) - a
        loss = n*(n+1)/2 + (n+1)*(u-floor(u))
    elseif u > l.min+1
        # number of levels lower than true level
        n = a - max(ceil(u), l.min+1)
        loss = n*(n+1)/2 + (n+1)*(ceil(u)-u)
    else
        # number of levels higher than true level
        n = a - max(ceil(u), l.min+1)
        loss = n*(n+1)/2 + (n+1)*(l.min+1-u)
    end
    return l.scale*loss
end

function grad(l::OrdinalHingeLoss, u::Float64, a::Number)
    #a = round(a)
    if u > a
        # number of levels higher than true level
        n = min(ceil(u), l.max) - a
        g = n
    else
        # number of levels lower than true level
        n = a - max(floor(u), l.min)
        g = -n
    end
    return l.scale*g
end

M_estimator(l::OrdinalHingeLoss, a::AbstractArray) = median(a)

########################################## LOGISTIC ##########################################
# f: ℜx{-1,1}-> ℜ
mutable struct LogisticLoss<:ClassificationLoss
    scale::Float64
    domain::Domain
end
LogisticLoss(scale=1.0::Float64; domain=BoolDomain()) = LogisticLoss(scale, domain)

evaluate(l::LogisticLoss, u::Float64, a::Bool) = l.scale*log(1+exp(-(2a-1)*u))

grad(l::LogisticLoss, u::Float64, a::Bool) = (aa = 2a-1; -aa*l.scale/(1+exp(aa*u)))

function M_estimator(l::LogisticLoss, a::AbstractArray{Bool,1})
    d, N = sum(a), length(a)
    log(N + d) - log(N - d) # very satisfying
end

########################################## WEIGHTED HINGE ##########################################
# f: ℜx{-1,1} -> ℜ
# f(u,a) = {     w * max(1-a*u, 0) for a = -1
#        = { c * w * max(1-a*u, 0) for a =  1
mutable struct WeightedHingeLoss<:ClassificationLoss
    scale::Float64
    domain::Domain
    case_weight_ratio::Float64 # >1 for trues to have more confidence than falses, <1 for opposite
end
WeightedHingeLoss(scale=1.0; domain=BoolDomain(), case_weight_ratio=1.0) =
    WeightedHingeLoss(scale, domain, case_weight_ratio)
HingeLoss(scale=1.0::Float64; kwargs...) = WeightedHingeLoss(scale; kwargs...) # the standard HingeLoss is a special case of WeightedHingeLoss

function evaluate(l::WeightedHingeLoss, u::Float64, a::Bool)
    loss = l.scale*max(1-(2*a-1)*u, 0)
    if l.case_weight_ratio !==1. && a
        loss *= l.case_weight_ratio
    end
    return loss
end

function grad(l::WeightedHingeLoss, u::Float64, a::Bool)
    an = (2*a-1) # change to {-1,1}
    g = (an*u>=1 ? 0 : -an*l.scale)
    if l.case_weight_ratio !==1. && a
        g *= l.case_weight_ratio
    end
    return g
end

function M_estimator(l::WeightedHingeLoss, a::AbstractArray{Bool,1})
    r = length(a)/length(filter(x->x>0, a)) - 1
    if l.case_weight_ratio > r
        m = 1.0
    elseif l.case_weight_ratio == r
        m = 0.0
    else
        m = -1.0
    end
end

########################################## MULTINOMIAL ##########################################
# f: ℜx{1, 2, ..., max-1, max} -> ℜ
# f computes the (negative log likelihood of the) multinomial logit,
# often known as the softmax function
# f(u, a) = exp(u[a]) / (sum_{a'} exp(u[a']))
#         = 1         / (sum_{a'} exp(u[a'] - u[a]))
mutable struct MultinomialLoss<:Loss
    max::Integer
    scale::Float64
    domain::Domain
end
MultinomialLoss(m, scale=1.0::Float64; domain=CategoricalDomain(m)) = MultinomialLoss(m,scale,domain)
embedding_dim(l::MultinomialLoss) = l.max
datalevels(l::MultinomialLoss) = 1:l.max # levels are encoded as the numbers 1:l.max

# in Julia v0.4, argument u is a row vector (row slice of a matrix), which in julia is 2d
# function evaluate(l::MultinomialLoss, u::Array{Float64,2}, a::Int)
# this breaks compatibility with v0.4
function evaluate(l::MultinomialLoss, u::Array{Float64,1}, a::Int)
    sumexp = 0 # inverse likelihood of observation
    # computing soft max directly is numerically unstable
    # instead note logsumexp(a_j) = logsumexp(a_j - M) + M
    # and we'll pick a good big (but not too big) M
    M = maximum(u) - u[a] # prevents overflow
    for j in 1:length(u)
        sumexp += exp(u[j] - u[a] - M)
    end
    loss = log(sumexp) + M
    return l.scale*loss
end

# in Julia v0.4, argument u is a row vector (row slice of a matrix), which in julia is 2d
# function grad(l::MultinomialLoss, u::Array{Float64,2}, a::Int)
# this breaks compatibility with v0.4
function grad(l::MultinomialLoss, u::Array{Float64,1}, a::Int)
    g = zeros(size(u))
    # Using some nice algebra, you can show
    g[a] = -1
    # and g[b] = -1/sum_{a' \in S} exp(u[b] - u[a'])
    # the contribution of one observation to one entry of the gradient
    # is always between -1 and 0
    for j in 1:length(u)
        M = maximum(u) - u[j] # prevents overflow
        sumexp = 0
        for jp in 1:length(u)
            sumexp += exp(u[jp] - u[j] - M)
        end
        g[j] += exp(-M)/sumexp
    end
    return l.scale*g
end

## we'll compute it via a stochastic gradient method
## with fixed step size
function M_estimator(l::MultinomialLoss, a::AbstractArray)
    u = zeros(l.max)'
    for i = 1:length(a)
        ai = a[i]
        u -= .1*grad(l, u, ai)
    end
    return u
end

########################################## One vs All loss ##########################################
# f: ℜx{1, 2, ..., max-1, max} -> ℜ
mutable struct OvALoss<:Loss
    max::Integer
    bin_loss::Loss
    scale::Float64
    domain::Domain
end
OvALoss(m::Integer, scale::Float64=1.0; domain=CategoricalDomain(m), bin_loss::Loss=LogisticLoss(scale)) = OvALoss(m,bin_loss,scale,domain)
OvALoss() = OvALoss(1) # for copying correctly
embedding_dim(l::OvALoss) = l.max
datalevels(l::OvALoss) = 1:l.max # levels are encoded as the numbers 1:l.max

# in Julia v0.4, argument u is a row vector (row slice of a matrix), which in julia is 2d
# function evaluate(l::OvALoss, u::Array{Float64,2}, a::Int)
# this breaks compatibility with v0.4
function evaluate(l::OvALoss, u::Array{Float64,1}, a::Int)
    loss = 0
    for j in 1:length(u)
        loss += evaluate(l.bin_loss, u[j], a==j)
    end
    return l.scale*loss
end

# in Julia v0.4, argument u is a row vector (row slice of a matrix), which in julia is 2d
# function grad(l::OvALoss, u::Array{Float64,2}, a::Int)
# this breaks compatibility with v0.4
function grad(l::OvALoss, u::Array{Float64,1}, a::Int)
  g = zeros(length(u))
  for j in 1:length(u)
      g[j] = grad(l.bin_loss, u[j], a==j)
  end
  return l.scale*g
end

function M_estimator(l::OvALoss, a::AbstractArray)
    u = zeros(l.max)
    for j = 1:l.max
        u[j] = M_estimator(l.bin_loss, a==j)
    end
    return u
end

########################################## Bigger vs Smaller loss ##########################################
# f: ℜx{1, 2, ..., max-1} -> ℜ
mutable struct BvSLoss<:Loss
    max::Integer
    bin_loss::Loss
    scale::Float64
    domain::Domain
end
function BvSLoss(m::Integer, scale::Float64=1.0; domain=OrdinalDomain(1,m), bin_loss::Loss=LogisticLoss(scale))
  @assert(m >= 2, error("Number of levels of ordinal variable must be at least 2; got $m."))
  BvSLoss(m,bin_loss,scale,domain)
end
BvSLoss() = BvSLoss(10) # for copying correctly
embedding_dim(l::BvSLoss) = l.max-1
datalevels(l::BvSLoss) = 1:l.max # levels are encoded as the numbers 1:l.max

function evaluate(l::BvSLoss, u::Array{Float64,1}, a::Int)
    loss = 0
    for j in 1:length(u)
        loss += evaluate(l.bin_loss, u[j], a>j)
    end
    return l.scale*loss
end

function grad(l::BvSLoss, u::Array{Float64,1}, a::Int)
  g = zeros(length(u))
  for j in 1:length(u)
      g[j] = grad(l.bin_loss, u[j], a>j)
  end
  return l.scale*g
end

function M_estimator(l::BvSLoss, a::AbstractArray)
    u = zeros(l.max)
    for j = 1:l.max-1
        u[j] = M_estimator(l.bin_loss, a.>j)
    end
    return u
end

########################################## ORDERED LOGISTIC ##########################################
# f: ℜx{1, 2, ..., max-1, max} -> ℜ
# f computes the (negative log likelihood of the) multinomial logit,
# often known as the softmax function
# f(u, a) = exp(u[a]) / (sum_{a'} exp(u[a']))
mutable struct OrdisticLoss<:Loss
    max::Integer
    scale::Float64
    domain::Domain
end
OrdisticLoss(m::Int, scale=1.0::Float64; domain=OrdinalDomain(1,m)) = OrdisticLoss(m,scale,domain)
embedding_dim(l::OrdisticLoss) = l.max
datalevels(l::OrdisticLoss) = 1:l.max # levels are encoded as the numbers 1:l.max

function evaluate(l::OrdisticLoss, u::Array{Float64,1}, a::Int)
    diffusquared = u[a]^2 .- u.^2
    M = maximum(diffusquared)
    invlik = sum(exp, (diffusquared .- M))
    loss = M + log(invlik)
    return l.scale*loss
end

function grad(l::OrdisticLoss, u::Array{Float64,1}, a::Int)
    g = zeros(size(u))
    # Using some nice algebra, you can show
    g[a] = 2*u[a]
    sumexp = sum(map(j->exp(- u[j]^2), 1:length(u)))
    for j in 1:length(u)
        diffusquared = u[j]^2 .- u.^2
        M = maximum(diffusquared)
        invlik = sum(exp,(diffusquared .- M))
        g[j] -= 2 * u[j] * exp(- M) / invlik
    end
    return l.scale*g
end

## we'll compute it via a stochastic gradient method
## with fixed step size
function M_estimator(l::OrdisticLoss, a::AbstractArray)
    u = zeros(l.max)'
    for i = 1:length(a)
        ai = a[i]
        u -= .1*grad(l, u, ai)
    end
    return u
end

#################### Multinomial Ordinal Logit #####################
# l: ℜ^{max-1} x {1, 2, ..., max-1, max} -> ℜ
# l computes the (negative log likelihood of the) multinomial ordinal logit.
#
# the length of the first argument u is one less than
# the number of levels of the second argument a,
# since the entries of u correspond to the division between each level
# and the one above it.
#
# XXX warning XXX
# the documentation in the comment below this point is defunct
#
# To yield a sensible pdf, the entries of u should be increasing
# (b/c they're basically the -log of the cdf at the boundary between each level)
#
# The multinomial ordinal logit corresponds to a likelihood p with
# p(u, a > i) ~ exp(-u[i]), so
# p(u, a)     ~ exp(-u[1]) * ... * exp(-u[a-1]) * exp(u[a]) * ... * exp(u[end])
#             = exp(- u[1] - ... - u[a-1] + u[a] + ... + u[end])
# and normalizing,
# p(u, a)     = p(u, a) / sum_{a'} p(u, a')
#
# So l(u, a) = -log(p(u, a))
#            = u[1] + ... + u[a-1] - u[a] - ... - u[end] +
#              log(sum_{a'}(exp(u[1] + ... + u[a'-1] - u[a'] - ... - u[end])))
#
# Inspection of this loss function confirms that given u,
# the most probable value a is the index of the first
# positive entry of u

mutable struct MultinomialOrdinalLoss<:Loss
    max::Integer
    scale::Float64
    domain::Domain
end
MultinomialOrdinalLoss(m::Int, scale=1.0::Float64; domain=OrdinalDomain(1,m)) = MultinomialOrdinalLoss(m,scale,domain)
MultinomialOrdinalLoss() = MultinomialOrdinalLoss(10) # for copying
embedding_dim(l::MultinomialOrdinalLoss) = l.max - 1
datalevels(l::MultinomialOrdinalLoss) = 1:l.max # levels are encoded as the numbers 1:l.max

function enforce_MNLOrdRules!(u; TOL=1e-3)
  u[1] = min(-TOL, u[1])
  for j=2:length(u)
    u[j] = min(u[j], u[j-1]-TOL)
  end
  u
end
# argument u is a row vector (row slice of a matrix), which in julia is 2d
# todo: increase numerical stability
function evaluate(l::MultinomialOrdinalLoss, u::Array{Float64,1}, a::Int)
  enforce_MNLOrdRules!(u)
  if a == 1
    return -l.scale*log(exp(0) - exp(u[1])) # (log(1 - exp(u[a] - 1)))
  elseif a == l.max
    return -l.scale*u[a-1]
  else
    return -l.scale*log(exp(u[a-1]) - exp(u[a])) # (u[a-1] + log(1 - exp(u[a] - u[a-1])))
  end
end

# argument u is a row vector (row slice of a matrix), which in julia is 2d
function grad(l::MultinomialOrdinalLoss, u::Array{Float64,1}, a::Int)
  enforce_MNLOrdRules!(u)
  g = zeros(size(u))
  if a == 1
    g[1] = -exp(u[1])/(exp(0) - exp(u[1]))
    # g[1] = 1/(1 - exp(-u[1]))
  elseif a == l.max
    g[a-1] = 1
  else
    # d = exp(u[a] - u[a-1])
    # g[a] = d/(1-d)
    # g[a-1] = - g[a] - 1
    g[a] = -exp(u[a])/(exp(u[a-1]) - exp(u[a]))
    g[a-1] = exp(u[a-1])/(exp(u[a-1]) - exp(u[a]))
  end
  return -l.scale*g
end

## we'll compute it via a stochastic gradient method
## with fixed step size
## (we don't need a hyper accurate estimate for this)
function M_estimator(l::MultinomialOrdinalLoss, a::AbstractVector)
    u = zeros(l.max-1)'
    for i = 1:length(a)
        ai = a[i]
        u -= .1*grad(l, u, ai)
    end
    return u
end

### convenience methods for evaluating and computing gradients on vectorized arguments
function evaluate(l::Loss, u::Array{Float64,1}, a::AbstractVector)
  @assert size(u) == size(a)
  out = 0
  for i=1:length(a)
    out += evaluate(l, u[i], a[i])
  end
  return out
end

#Optimized vector evaluate on single-dimensional losses
function evaluate(l::SingleDimLoss, u::Vector{Float64}, a::AbstractVector)
  losseval = (x::Float64, y::Number) -> evaluate(l, x, y)
  mapped = fill!(similar(u),0.)
  map!(losseval, mapped, u, a)
  reduce(+, mapped)
end

# now for multidimensional losses
function evaluate(l::Loss, u::Array{Float64,2}, a::AbstractVector)
  # @show size(u,1)
  # @show size(a)
  @assert size(u,1) == length(a)
  out = 0
  for i=1:length(a)
    out += evaluate(l, u[i,:], a[i])
  end
  return out
end

function grad(l::Loss, u::Array{Float64,1}, a::AbstractVector)
  @assert size(u) == size(a)
  mygrad = zeros(size(u))
  for i=1:length(a)
    mygrad[i] = grad(l, u[i], a[i])
  end
  return mygrad
end

# Optimized vector grad on single-dimensional losses
function grad(l::SingleDimLoss, u::Vector{Float64}, a::AbstractVector)
  lossgrad = (x::Float64,y::Number) -> grad(l, x, y)
  mapped = fill!(similar(u),0.)
  map!(lossgrad, mapped, u, a)
end

# now for multidimensional losses
function grad(l::Loss, u::Array{Float64,2}, a::AbstractVector)
  @assert size(u,1) == length(a)
  mygrad = zeros(size(u))
  for i=1:length(a)
    mygrad[i,:] = grad(l, u[i,:], a[i])
  end
  return mygrad
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for geometry.jl

s = 3

vv1 = [  1f0,  2f0,   3f0]
vv2 = [1.5f0,  3f0, 4.5f0]
vv3 = [ 10f0, 20f0,  30f0]

sv1 = SVector{size(vv1)...}(vv1)
sv2 = SVector{size(vv2)...}(vv2)
sv3 = SVector{size(vv3)...}(vv3)

@testset "Vec" begin
    v1, v2, v3 = (sv1, vv2, vv3) .|> Vec

    @testset "constructor" begin
        @test v1 == sv1
        @test v2 == sv2
        @test Vec(vv1...) == sv1
    end

    @testset "operations" begin
        @test -v1 ≈ -sv1

        @test norm(v1) === norm(sv1)
        @test normalize(v1) === Vec(normalize(sv1))
        @test norm²(v1) ≈ norm(sv1)^2

        @test v1 ⋅ v2 === sv1 ⋅ sv2

        @test v1 + v2 === Vec(sv1 + sv2)
        @test v1 - v2 === Vec(sv1 - sv2)
        @test v1 × v2 === Vec(sv1 × sv2)

        @test s * v1 === v1 * s === Vec(s * sv1)

        @test Vec([15, 30, 45]) ≈ v2 * 10
    end
end

@testset "Normal" begin
    v1, v2, v3 = (sv1, vv2, vv3) .|> Normal

    @testset "constructor" begin
        @test v1 == sv1
        @test v2 == sv2
        @test Normal(vv1...) == sv1
    end

    @testset "operations" begin
        @test norm(v1) === norm(sv1)
        @test normalize(v1) === Normal{true}(normalize(sv1))
        @test norm²(v1) ≈ norm(sv1)^2

        @test v1 ⋅ v2 === sv1 ⋅ sv2

        @test v1 + v2 === Normal(sv1 + sv2)
        @test v1 - v2 === Normal(sv1 - sv2)
        @test v1 × v2 === Normal(sv1 × sv2)

        @test s * v1 === v1 * s === Normal(s * sv1)

        @test Normal([15, 30, 45]) ≈ v2 * 10
    end
end

@testset "Point" begin
    p1, p2 = (sv1, vv2) .|> Point
    v = Vec(vv3)

    @testset "constructor" begin
        @test p1.v == sv1
        @test p2.v == sv2
        @test_throws ArgumentError Point(ones(4))
    end

    @testset "operations" begin
        @test p1 ≈ Point(vv3 ./ 10)

        @test p1 - p2 == Vec(sv1 - sv2)

        @test p1 + v == Point(sv1 + sv3)
        @test p1 - v == Point(sv1 - sv3)

        @test Point([15, 30, 45]) ≈ p2 * 10
    end
end

@testset "ONB" begin
	pcg = PCG()

	@test begin
        for _ ∈ 1:100_000
            normal = Normal(rand(pcg, Float32, 3)) |> normalize

            e1, e2, e3 = create_onb_from_z(normal)

            # Verify that the z axis is aligned with the normal
            @assert e3 ≈ normal

            # Verify that the base is orthogonal
            atol = √(eps(eltype(normal))) # nonstandard approximation threshold
            @assert isapprox(e1 ⋅ e2, 0, atol=atol)
            @assert isapprox(e2 ⋅ e3, 0, atol=atol)
            @assert isapprox(e3 ⋅ e1, 0, atol=atol)

            # Verify that each component is normalized
            @assert norm²(e1) ≈ 1
            @assert norm²(e2) ≈ 1
            @assert norm²(e3) ≈ 1
        end
        true
    end
end


text:
const PyJuliaNumberValue_Type__ref = Ref(PyNULL)
PyJuliaNumberValue_Type() = begin
    ptr = PyJuliaNumberValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.NumberValue",
                base = base,
                as_number = (
                    bool = pyjlnumber_bool,
                    positive = pyjlnumber_positive,
                    negative = pyjlnumber_negative,
                    absolute = pyjlnumber_absolute,
                    power = pyjlnumber_power,
                    add = pyjlnumber_binop(+),
                    subtract = pyjlnumber_binop(-),
                    multiply = pyjlnumber_binop(*),
                    truedivide = pyjlnumber_binop(/),
                    divmod = pyjlnumber_binop((x, y) -> (fld(x, y), mod(x, y))),
                    floordivide = pyjlnumber_binop(fld),
                    remainder = pyjlnumber_binop(mod),
                    lshift = pyjlnumber_binop(<<),
                    rshift = pyjlnumber_binop(>>),
                    and = pyjlnumber_binop(&),
                    xor = pyjlnumber_binop(xor),
                    or = pyjlnumber_binop(|),
                ),
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyNumberABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaNumberValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaNumberValue_New(x::Number) = PyJuliaValue_New(PyJuliaNumberValue_Type(), x)
PyJuliaValue_From(x::Number) = PyJuliaNumberValue_New(x)

const PyJuliaComplexValue_Type__ref = Ref(PyNULL)
PyJuliaComplexValue_Type() = begin
    ptr = PyJuliaComplexValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaNumberValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.ComplexValue",
                base = base,
                getset = [
                    (name = "real", get = pyjlcomplex_real),
                    (name = "imag", get = pyjlcomplex_imag),
                ],
                methods = [
                    (
                        name = "conjugate",
                        flags = Py_METH_NOARGS,
                        meth = pyjlcomplex_conjugate,
                    ),
                    (
                        name = "__complex__",
                        flags = Py_METH_NOARGS,
                        meth = pyjlcomplex_complex,
                    ),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyComplexABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaComplexValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaComplexValue_New(x::Complex) = PyJuliaValue_New(PyJuliaComplexValue_Type(), x)
PyJuliaValue_From(x::Complex) = PyJuliaComplexValue_New(x)

const PyJuliaRealValue_Type__ref = Ref(PyNULL)
PyJuliaRealValue_Type() = begin
    ptr = PyJuliaRealValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaNumberValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.RealValue",
                base = base,
                as_number = (float = pyjlreal_float,),
                getset = [
                    (name = "real", get = pyjlreal_real),
                    (name = "imag", get = pyjlreal_imag),
                ],
                methods = [
                    (name = "conjugate", flags = Py_METH_NOARGS, meth = pyjlreal_conjugate),
                    (name = "__complex__", flags = Py_METH_NOARGS, meth = pyjlreal_complex),
                    (name = "__trunc__", flags = Py_METH_NOARGS, meth = pyjlreal_trunc),
                    (name = "__floor__", flags = Py_METH_NOARGS, meth = pyjlreal_floor),
                    (name = "__ceil__", flags = Py_METH_NOARGS, meth = pyjlreal_ceil),
                    (name = "__round__", flags = Py_METH_VARARGS, meth = pyjlreal_round),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyRealABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaRealValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaRealValue_New(x::Real) = PyJuliaValue_New(PyJuliaRealValue_Type(), x)
PyJuliaValue_From(x::Real) = PyJuliaRealValue_New(x)

const PyJuliaRationalValue_Type__ref = Ref(PyNULL)
PyJuliaRationalValue_Type() = begin
    ptr = PyJuliaRationalValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaRealValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.RationalValue",
                base = base,
                getset = [
                    (name = "numerator", get = pyjlrational_numerator),
                    (name = "denominator", get = pyjlrational_denominator),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyRationalABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaRationalValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaRationalValue_New(x::Rational) = PyJuliaValue_New(PyJuliaRationalValue_Type(), x)
PyJuliaValue_From(x::Rational) = PyJuliaRationalValue_New(x)

const PyJuliaIntegerValue_Type__ref = Ref(PyNULL)
PyJuliaIntegerValue_Type() = begin
    ptr = PyJuliaIntegerValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaRealValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.IntegerValue",
                base = base,
                as_number = (
                    int = pyjlinteger_int,
                    index = pyjlinteger_index,
                    invert = pyjlinteger_invert,
                ),
                getset = [
                    (name = "numerator", get = pyjlinteger_numerator),
                    (name = "denominator", get = pyjlinteger_denominator),
                ],
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyIntegralABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaIntegerValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaIntegerValue_New(x::Integer) = PyJuliaValue_New(PyJuliaIntegerValue_Type(), x)
PyJuliaValue_From(x::Integer) = PyJuliaIntegerValue_New(x)

pyjlnumber_bool(xo::PyPtr) =
    try
        iszero(PyJuliaValue_GetValue(xo)::Number) ? Cint(0) : Cint(1)
    catch err
        PyErr_SetJuliaError(err)
        Cint(-1)
    end

pyjlnumber_positive(xo::PyPtr) =
    try
        PyObject_From(+(PyJuliaValue_GetValue(xo)::Number))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlnumber_negative(xo::PyPtr) =
    try
        PyObject_From(-(PyJuliaValue_GetValue(xo)::Number))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlnumber_absolute(xo::PyPtr) =
    try
        PyObject_From(abs(PyJuliaValue_GetValue(xo)::Number))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

struct pyjlnumber_binop{F}
    f::F
end
(f::pyjlnumber_binop)(xo::PyPtr, yo::PyPtr) = begin
    r = PyObject_TryConvert(xo, Number)
    r == -1 && return PyNULL
    r == 0 && return PyNotImplemented_New()
    x = takeresult(Number)
    r = PyObject_TryConvert(yo, Number)
    r == -1 && return PyNULL
    r == 0 && return PyNotImplemented_New()
    y = takeresult(Number)
    try
        PyObject_From(f.f(x, y))
    catch err
        if err isa MethodError && err.f === f.f
            PyNotImplemented_New()
        else
            PyErr_SetJuliaError(err)
            PyNULL
        end
    end
end

pyjlnumber_power(xo::PyPtr, yo::PyPtr, zo::PyPtr) = begin
    r = PyObject_TryConvert(xo, Number)
    r == -1 && return PyNULL
    r == 0 && return PyNotImplemented_New()
    x = takeresult(Number)
    r = PyObject_TryConvert(yo, Number)
    r == -1 && return PyNULL
    r == 0 && return PyNotImplemented_New()
    y = takeresult(Number)
    if PyNone_Check(zo)
        try
            PyObject_From(x^y)
        catch err
            if err isa MethodError && err.f === ^
                PyNotImplemented_New()
            else
                PyErr_SetJuliaError(err)
                PyNULL
            end
        end
    else
        r = PyObject_TryConvert(zo, Number)
        r == -1 && return PyNULL
        r == 0 && return PyNotImplemented_New()
        z = takeresult(Number)
        try
            PyObject_From(powermod(x, y, z))
        catch err
            if err isa MethodError && err.f === powermod
                PyNotImplemented_New()
            else
                PyErr_SetJuliaError(err)
                PyNULL
            end
        end
    end
end

pyjlcomplex_real(xo::PyPtr, ::Ptr{Cvoid}) =
    try
        PyObject_From(real(PyJuliaValue_GetValue(xo)::Complex))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlcomplex_imag(xo::PyPtr, ::Ptr{Cvoid}) =
    try
        PyObject_From(imag(PyJuliaValue_GetValue(xo)::Complex))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlcomplex_conjugate(xo::PyPtr, ::PyPtr) =
    try
        PyObject_From(conj(PyJuliaValue_GetValue(xo)::Complex))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlcomplex_complex(xo::PyPtr, ::PyPtr) =
    try
        PyComplex_From(convert(Complex{Float64}, PyJuliaValue_GetValue(xo)::Complex))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlreal_real(xo::PyPtr, ::Ptr{Cvoid}) = (Py_IncRef(xo); xo)

pyjlreal_imag(xo::PyPtr, ::Ptr{Cvoid}) = PyLong_From(0)

pyjlreal_conjugate(xo::PyPtr, ::PyPtr) = (Py_IncRef(xo); xo)

pyjlreal_complex(xo::PyPtr, ::PyPtr) =
    try
        PyComplex_From(convert(Float64, PyJuliaValue_GetValue(xo)::Real))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlreal_float(xo::PyPtr) =
    try
        PyFloat_From(convert(Float64, PyJuliaValue_GetValue(xo)::Real))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlreal_trunc(xo::PyPtr, ::PyPtr) =
    try
        PyObject_From(trunc(Integer, PyJuliaValue_GetValue(xo)::Real))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlreal_floor(xo::PyPtr, ::PyPtr) =
    try
        PyObject_From(floor(Integer, PyJuliaValue_GetValue(xo)::Real))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlreal_ceil(xo::PyPtr, ::PyPtr) =
    try
        PyObject_From(ceil(Integer, PyJuliaValue_GetValue(xo)::Real))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlreal_round(xo::PyPtr, args::PyPtr) = begin
    ism1(PyArg_CheckNumArgsLe("round", args, 1)) && return PyNULL
    ism1(PyArg_GetArg(Union{Int,Nothing}, "round", args, 0, nothing)) && return PyNULL
    ndigits = takeresult(Union{Int,Nothing})
    x = PyJuliaValue_GetValue(xo)::Real
    try
        if ndigits === nothing
            PyObject_From(round(Integer, x))
        else
            PyObject_From(round(x; digits = ndigits))
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlrational_numerator(xo::PyPtr, ::Ptr{Cvoid}) =
    try
        PyObject_From(numerator(PyJuliaValue_GetValue(xo)::Rational))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlrational_denominator(xo::PyPtr, ::Ptr{Cvoid}) =
    try
        PyObject_From(denominator(PyJuliaValue_GetValue(xo)::Rational))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end

pyjlinteger_numerator(xo::PyPtr, ::Ptr{Cvoid}) = (Py_IncRef(xo); xo)

pyjlinteger_denominator(xo::PyPtr, ::Ptr{Cvoid}) = PyLong_From(1)

pyjlinteger_int(xo::PyPtr) = PyLong_From(PyJuliaValue_GetValue(xo)::Integer)

pyjlinteger_index(xo::PyPtr) = PyLong_From(PyJuliaValue_GetValue(xo)::Integer)

pyjlinteger_invert(xo::PyPtr) =
    try
        PyObject_From(~(PyJuliaValue_GetValue(xo)::Integer))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end


text:
# Case insensitive Dict - a simple wrapper over Dict

struct CIDict{K, T}

    dct::Dict{K, T}

    # type checking
    check(K) = K <: AbstractString || K <: Symbol ||
        throw(ArgumentError("Key must be Symbol or String type"))

    # constructors
    CIDict{K, T}() where {K,T} = (check(K); new(Dict{K,T}()))
    CIDict{K, T}(d::Dict{K,T}) where {K,T} = begin
        check(K)
        d2 = Dict{K,T}()
        for k in keys(d)
            d2[lcase(k)] = d[k]
        end
        new(d2)
    end
end

lcase(s::Symbol) = Symbol(lowercase(String(s)))
lcase(s::AbstractString) = lowercase(s)

Base.getindex(d::CIDict, s::Symbol) = d.dct[lcase(s)]
Base.getindex(d::CIDict, s::String) = d.dct[lcase(s)]

Base.setindex!(d::CIDict, v, s::Symbol) = d.dct[lcase(s)] = v
Base.setindex!(d::CIDict, v, s::String) = d.dct[lcase(s)] = v

Base.haskey(d::CIDict, s::Symbol) = haskey(d.dct, lcase(s))
Base.haskey(d::CIDict, s::String) = haskey(d.dct, lcase(s))

Base.keys(d::CIDict) = keys(d.dct)
Base.values(d::CIDict) = values(d.dct)

Base.iterate(d::CIDict) = Base.iterate(d.dct)
Base.iterate(d::CIDict, state) = Base.iterate(d.dct, state)

Base.length(d::CIDict) = length(d.dct)

issym(x) = typeof(x) == Symbol

function Base.show(io::IO, d::SASLib.CIDict)
    print(io, "CIDict(")
    for (i, (k,v)) in enumerate(d.dct)
        i > 1 && print(io, ", ")
        print(io, issym(k) ? ":" : "", k, " => ", v)
    end
    print(io, ")")
end


text:
# a contiguous matrix with some ell information slapped on for safety
struct BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA} <: AbstractArray{T,2}
    parent::AA
    m_ells::NTuple{M_BLOCKS, UnitRange{Int64}}
    n_ells::NTuple{N_BLOCKS, UnitRange{Int64}}
end
@forward BlockSpectralMatrix.parent (Base.getindex, Base.setindex, Base.setindex!,
    Base.size, Base.iterate, Base.axes, Base.show, Base.in, Base.strides, Base.elsize,
    Base.dataids)

Base.parent(A::BlockSpectralMatrix) = A.parent

function Base.hvcat(rows::Tuple{Vararg{Int}}, values::SA...)  where {T, SA<:SpectralArray{T}}
    𝐀 = hvcat(rows, map(A -> parent(A), values)...)
    m_blocks = length(rows)
    n_blocks = first(rows)

    m_ells = map(A->axes(A,1), values[1:m_blocks])
    n_ells = map(A->axes(A,2), values[1:n_blocks:end])
    return BlockSpectralMatrix{T, m_blocks, n_blocks, typeof(𝐀)}(𝐀, m_ells, n_ells)
end

function Base.vcat(values::SA...)  where {T, SA<:SpectralArray{T}}
    𝐀 = vcat(map(A -> parent(A), values)...)
    m_blocks = length(values)
    n_blocks = 1
    m_ells = map(A->axes(A,1), values)
    n_ells = (axes(first(values),2),)
    return BlockSpectralMatrix{T, m_blocks, n_blocks, typeof(𝐀)}(
        𝐀, m_ells, n_ells)
end

function Base.hcat(values::SA...)  where {T, SA<:SpectralArray{T}}
    𝐀 = hcat(map(A -> parent(A), values)...)
    m_blocks = 1
    n_blocks = length(values)
    m_ells = (axes(first(values),1),)
    n_ells = map(A->axes(A,2), values)
    return BlockSpectralMatrix{T, m_blocks, n_blocks, typeof(𝐀)}(
        𝐀, m_ells, n_ells)
end

function Base.show(io::IO, m::MIME"text/plain",
        ba::BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS}) where {T,M_BLOCKS,N_BLOCKS}
    print(io, "BlockSpectralMatrix (column blocks=$(M_BLOCKS), row blocks=$(N_BLOCKS))\n")
    print(io, "row ℓ:    ")
    Base.show(io, m, ba.m_ells)
    print(io, "\ncolumn ℓ: ")
    Base.show(io, m, ba.n_ells)

    print(io, "\nparent: ")

    nrows, ncols = displaysize(io)
    Base.show(
        IOContext(io, :compact => true, :displaysize => (min(10,nrows),ncols)),
        m, parent(ba))
end

struct BlockSpectralFactorization{T,M_BLOCKS,N_BLOCKS,AA,
        F<:Factorization{T}} <: Factorization{T}
    parent::F
    m_ells::NTuple{M_BLOCKS, UnitRange{Int64}}
    n_ells::NTuple{N_BLOCKS, UnitRange{Int64}}
end

function LinearAlgebra.inv(
        F::BlockSpectralFactorization{T,M_BLOCKS,N_BLOCKS,AA}) where {T,M_BLOCKS,N_BLOCKS,AA}
    return BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA}(inv(F.parent), F.m_ells, F.n_ells)
end

Base.parent(F::BlockSpectralFactorization) = F.parent

function BlockSpectralMatrix(A::SpectralArray{T,2,AA}) where {T,AA}
    return  BlockSpectralMatrix{T,1,1,AA}(parent(A),
        (UnitRange(axes(A,1)),), (UnitRange(axes(A,2)),))
end

function BlockSpectralMatrix(A::SpectralArray{T,1,AA}) where {T,AA}
    return  BlockSpectralMatrix{T,1,1,AA}(parent(A),
        (UnitRange(axes(A,1)),), (UnitRange(axes(A,2)),))
end

function LinearAlgebra.lu(A::SpectralArray{T,2,AA}) where {T,AA}
    F = lu(parent(A))
    return BlockSpectralFactorization{T,1,1,AA,typeof(F)}(F,
        (UnitRange(axes(A,1)),), (UnitRange(axes(A,2)),))
end

function LinearAlgebra.lu(A::BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA}) where {T,M_BLOCKS,N_BLOCKS,AA}
    F = lu(parent(A))
    return BlockSpectralFactorization{T,M_BLOCKS,N_BLOCKS,AA,typeof(F)}(F,
        A.m_ells, A.n_ells)
end

function (\)(F::BlockSpectralFactorization, B::SpectralVector{T}) where T
    x = parent(F) \ parent(B)
    return SpectralVector(x, B.parent.offsets)
end

function (\)(F::BlockSpectralFactorization,
        B::BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA}) where {T,M_BLOCKS,N_BLOCKS,AA}
    x = parent(F) \ parent(B)
    return BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA}(x, B.m_ells, B.n_ells)
end

function (\)(A::BlockSpectralMatrix,
        B::BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA}) where {T,M_BLOCKS,N_BLOCKS,AA}
    # check dims match
    for blockindex = 1:N_BLOCKS
        @assert A.n_ells[blockindex] == B.m_ells[blockindex]
    end
    F = lu(A)
    x = parent(F) \ parent(B)
    return BlockSpectralMatrix{T,M_BLOCKS,N_BLOCKS,AA}(x, B.m_ells, B.n_ells)
end

function (\)(A::SpectralArray{T,2}, B::SpectralVector{T}) where T
    @assert firstindex(A,2) == firstindex(B,1)
    F = lu(A)
    x = parent(F) \ parent(B)
    return SpectralVector(x, B.parent.offsets)
end

"""
    getblock(A::BlockSpectralMatrix{T,M_BLOCKS,1}, i) where {T,M_BLOCKS}

Extract a sub-block from a `BlockSpectralMatrix`.

# Arguments:
- `A::BlockSpectralMatrix{T,M_BLOCKS,1}`: array to extract from
- `i::Int`: index of the sub-block (1-indexed)

# Returns:
- `Array{T,2}`: sub-blocks
"""
function getblock(A::BlockSpectralMatrix{T,M_BLOCKS,1},
                  i::Int) where {T,M_BLOCKS}
    @assert i ≤ M_BLOCKS
    offset = 0
    for block_i in 1:(i-1)
        offset += length(A.m_ells[block_i])
    end
    return SpectralVector(parent(A)[(1+offset):(length(A.m_ells[i])+offset)], A.m_ells[i])
end

# based on UnPack's macro
@doc raw"""
    @spectra [expr=BlockSpectralMatrix]

Unpack a block vector. This is equivalent to calling `getblock` for all the
sub-blocks and putting them in a Tuple.

# Example
```julia
# compute stacked EE,BB mode-coupling matrix from mask alm
M_EE_BB = mcm(:EE_BB, map2alm(mask1), map2alm(mask2))
# apply the 2×2 block mode-coupling matrix to the stacked EE and BB spectra
@spectra Cl_EE, Cl_BB = M_EE_BB \ [pCl_EE; pCl_BB]
```
"""
macro spectra(args)
    args.head!=:(=) && error("Expression needs to be of form `a, b = c`")
    items, suitcase = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitcase_instance = gensym()
    kd = [:( $key = PowerSpectra.getblock($suitcase_instance, $i) ) for
        (i, key) in enumerate(items)]
    kdblock = Expr(:block, kd...)
    expr = quote
        local $suitcase_instance = $suitcase # handles if suitcase is not a variable but an expression
        $kdblock
        $suitcase_instance # return RHS of `=` as standard in Julia
    end
    esc(expr)
end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov