text:
#############
#  PDEOperatorSchur
# represent a splitting rank 2 operator that's been discretized in the Y direction using
# the Schur decomposition
#
#  S is an AbstractOperatorSchur that encodes the schur decomposition
#############

abstract AbstractPDEOperatorSchur

immutable PDEOperatorSchur{OS<:AbstractOperatorSchur,LT<:Number,MT<:Number,FT<:Operator} <: AbstractPDEOperatorSchur
    Bx::Vector{FT}
    Lx::Operator{LT}
    Mx::Operator{MT}
    S::OS

    indsBx::Vector{Int}
    indsBy::Vector{Int}

    Rdiags::Vector{QROperator}
end

#TODO: Functional type?
Base.eltype{OS,LT,MT,FT}(PO::PDEOperatorSchur{OS,LT,MT,FT})=promote_type(eltype(PO.S),LT,MT)

function PDEOperatorSchur{LT<:Number,MT<:Number,BT<:Number,ST<:Number}(Bx,Lx::Operator{LT},Mx::Operator{MT},S::AbstractOperatorSchur{BT,ST},indsBx,indsBy)
    Xops=promotespaces([Lx,Mx])
    Lx=cache(Xops[1]);Mx=cache(Xops[2])

    ny=size(S,1)
    nbcs=numbcs(S)
    RT=promote_type(LT,MT,BT,ST)
    Rdiags=Array(QROperator,ny)
    resizedata!(Lx,ny,:);resizedata!(Mx,ny,:)

    for k=1:ny-nbcs
        ##TODO: Do block case
        Rdiags[k]=qrfact([Bx;PlusOperator(Operator{RT}[getdiagonal(S,k,1)*Lx,getdiagonal(S,k,2)*Mx])])
        resizedata!(Rdiags[k],:,ny)
    end

    if isempty(Bx)
        funcs=Array(CachedOperator{Float64},0)
    else
        funcs=Array(CachedOperator{mapreduce(eltype,promote_type,Bx)},length(Bx))
        for k=1:length(Bx)
            funcs[k]=cache(Bx[k])
            resizedata!(funcs[k],:,ny)
        end
    end

    PDEOperatorSchur(funcs,Lx,Mx,S,indsBx,indsBy,Rdiags)
end

PDEOperatorSchur(Bx::Vector,Lx::Operator,Mx::Operator,S::AbstractOperatorSchur) =
    PDEOperatorSchur(Bx,Lx,Mx,S,[1:length(Bx)],length(Bx)+[1:numbcs(S)])
PDEOperatorSchur(Bx::Vector,Lx::Operator,Mx::UniformScaling,S::AbstractOperatorSchur) =
    PDEOperatorSchur(Bx,Lx,ConstantOperator(Mx.λ),S)
PDEOperatorSchur(Bx::Vector,Lx::UniformScaling,Mx::Operator,S::AbstractOperatorSchur) =
    PDEOperatorSchur(Bx,ConstantOperator(Lx.λ),Mx,S)

function PDEOperatorSchur(Bx,By,A::Operator,ny::Integer,indsBx,indsBy)
    @assert iskronsumop(A)
    ops=sumops(A)
    @assert all(isproductop,ops)
    @assert length(ops)==2

    PDEOperatorSchur(Bx,dekron(ops[1],1),dekron(ops[2],1),schurfact(By,[dekron(ops[1],2),dekron(ops[2],2)],ny),indsBx,indsBy)
end

bcinds(S::PDEOperatorSchur,k)=k==1?S.indsBx:S.indsBy
numbcs(S::AbstractPDEOperatorSchur,k)=length(bcinds(S,k))
Base.length(S::PDEOperatorSchur)=size(S.S,1)

function PDEOperatorSchur(A::Vector,ny::Integer)
    indsBx,Bx=findfunctionals(A,1)
    indsBy,By=findfunctionals(A,2)

    LL=A[end]

    PDEOperatorSchur(Bx,By,LL,ny,indsBx,indsBy)
end

for op in (:domainspace,:rangespace)
    @eval begin
        $op(P::PDEOperatorSchur,k::Integer)=k==1?$op(P.Lx):$op(P.S)
        $op(L::PDEOperatorSchur)=$op(L,1)⊗$op(L,2)
    end
end

domain(P::PDEOperatorSchur,k::Integer)=k==1?domain(P.Lx):domain(P.S)
domain(P::PDEOperatorSchur)=domain(domainspace(P))

## Constructuor

Base.schurfact{LT<:Number,MT<:Number,BT<:Number,ST<:Number}(Bx,
            Lx::Operator{LT},Mx::Operator{MT},
            S::AbstractOperatorSchur{BT,ST},
            indsBx,indsBy) =
    PDEOperatorSchur(Bx,Lx,Mx,S,indsBx,indsBy)

function Base.schurfact(Bx,By,A::Operator,ny::Integer,indsBx,indsBy)
    @assert iskronsumop(A)
    opsx,opsy=simplifydekron(A)
    @assert length(opsx)==length(opsy)==2

    schurfact(Bx,
              opsx[1],opsx[2],
              schurfact(By,opsy,ny),
              indsBx,indsBy)
end

function Base.schurfact{BT<:Operator}(A::Vector{BT},S::ProductDomain,ny::Integer)
    indsBx,Bx=findfunctionals(A,1)
    indsBy,By=findfunctionals(A,2)

    LL=A[end]

    schurfact(Bx,By,LL,ny,indsBx,indsBy)
end

Base.schurfact{BT<:Operator}(A::Vector{BT},n::Integer) = schurfact(A,domain(A[end]),n)
Base.schurfact(A::Operator,n::Integer) = schurfact([A],n)

## discretize

function discretize{OT<:Operator}(A::Vector{OT},S...)
    if iskronsumop(A[end])&&length(simplifydekron(A[end])[1])==2
        schurfact(A,S...)
    elseif isdiagop(A[end],2)
        schurfact(A,S...)
    else
        kronfact(A,S...)
    end
end
discretize(A::Operator,n::Integer)=discretize([A],n)


text:
using IndexableIterators
using Test

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

ext: jl


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_bboxes(fname)
    pagewidth, pageheight = 1200, 1400

    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setline(40)
    sethue("blue")
    box(BoundingBox(), :stroke)
    sethue("red")
    box(BoundingBox() * 0.9, :stroke)

    sethue("green")
    box(BoundingBox(centered=false) * 0.9, :stroke)

    sethue("purple")
    setopacity(0.7)
    fontsize(100)
    box(BoundingBox("purple"), :fill)
    sethue("white")
    text("purple")

    setline(10)
    sethue("black")
    bb = BoundingBox()
    circle(midpoint(bb), 30, :stroke)
    for i in 0:0.05:1
        randomhue()
        circle(between(bb, i), 30, :fill)
    end

    # combine and intersections
    box1 = BoundingBox()/3 - (200, 200)
    box2 = BoundingBox()/4 - (300, 300)
    box3 = box1 + box2
    box4 = intersectboundingboxes(box1, box2)
    sethue("black")
    setline(2)
    box.([box1, box2, box3, box4], :stroke)

    sethue("pink")
    box(box1, :fill)

    sethue("yellow")
    box(box2, :fill)

    sethue("purple")
    box(box3, :fill)

    sethue("green")
    box(box4, :fill)

    @test isapprox(boxaspectratio(box1), 1.166, atol = 0.01)
    @test isapprox(boxdiagonal(box1), 614.636, atol = 0.01)

    # get vertices
    bv = box(box1, vertices=true)
    @test isapprox(bv[1].y, 33.333, atol = 0.1)

    setline(1)
    # contains
    for i in 1:1000
        pt = Point(rand(-1000:1000), rand(-1000:1000))
        if isinside(pt, box4)
            sethue("red")
            circle(pt, 15, :fill)
            sethue("black")
            circle(pt, 1, :fill)
        elseif isinside(pt, box2)
            sethue("white")
            circle(pt, 20, :fill)
            sethue("black")
            circle(pt, 1, :fill)
        elseif isinside(pt, box1)
            sethue("black")
            circle(pt, 5, :stroke)
        end
    end

    # poly
    s = star(O, 300, 7, 0.5, 0, vertices=true)
    box5 = BoundingBox(s)
    poly(s, :stroke, close=true)
    prettypoly(box5, :stroke, close=true, vertexlabels=(n, t) -> label(string(n)))
    @test finish() == true
    println("...finished bounding box test, saved in $(fname)")
end

Random.seed!(42)

fname = "bounding-box-test.png"

test_bboxes(fname)


text:
using StatsBase
using MultipleTesting
using IndependentHypothesisWeighting
using Test

pvs = rand(100)
for method in [Bonferroni(), BenjaminiHochberg()]
	unit_priority_ws = PriorityWeights(ones(100))
	unit_ws = UnitWeights{Float64}(100)
	adjust_unwt = adjust(pvs, method)
	adjust_wt = adjust(pvs, unit_priority_ws, method)
	adjust_unit_wt = adjust(pvs, unit_ws, method)

	@test adjust_unwt == adjust_wt
	@test adjust_unwt == adjust_unit_wt

	half_ws = PriorityWeights([2 .* ones(50); zeros(50)])
	adjust_half_wt = adjust(pvs, half_ws, method)
	adjust_half_unwt = adjust(pvs[1:50], method)
	@test  adjust_half_wt[1:50] == adjust_half_unwt
end


text:
function plotBaseEdificio3d_v2(fpe, x, alturaPiso, ps_predio,
                             ps_volteor, matConexionVertices_ss, vecVertices_ss,
                            ps_restSombra, matConexionVertices_cs, vecVertices_cs,
                            ps_publico, ps_calles, ps_base, ps_baseSeparada)

    f_predio = fpe.predio
    f_volTeor = fpe.volTeor
    f_restSombra = fpe.restSombra
    f_edif = fpe.edif

    f_sombraVolTeor_p = fpe.sombraVolTeor_p
    f_sombraVolTeor_o = fpe.sombraVolTeor_o
    f_sombraVolTeor_s = fpe.sombraVolTeor_s

    f_sombraEdif_p = fpe.sombraEdif_p
    f_sombraEdif_o = fpe.sombraEdif_o
    f_sombraEdif_s = fpe.sombraEdif_s

    alt = x[1]

    if f_predio
        # Grafica Predio
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_predio, 0, nothing, nothing, nothing, "green", .3)
    end

    if f_volTeor
        # Grafica Volumen Teórico
        V_volteor = ps_volteor.Vertices[1]
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v5(ps_volteor, matConexionVertices_ss, vecVertices_ss, fig, ax, ax_mat)
    end

    if f_restSombra
        # Grafica Volumen Teórico
        V_restSombra = ps_restSombra.Vertices[1]
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v5(ps_restSombra, matConexionVertices_cs, vecVertices_cs, fig, ax, ax_mat, "gray", 0.1)
    end

    #ps_calles = polyShape.polyUnion_v2(ps_calles, ps_calles)
    fig, ax = polyShape.plotPolyshape3d_v1(ps_calles, 0, fig, ax, "grey", 0.25)

    if f_edif
        # Grafica cabida óptima
        numPisos = Int(round(alt/alturaPiso))
        for k = 1:numPisos
            for j=1:ps_base.NumRegions
                V_base_j = ps_base.Vertices[j]
                numVerticesBase_j = size(V_base_j,1)
                V_kj = [V_base_j alturaPiso * (k - 1) * ones(numVerticesBase_j,1);
                        V_base_j alturaPiso * k * ones(numVerticesBase_j,1)]
                fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(PolyShape([V_kj],1), alturaPiso * (k - 1), fig, ax, ax_mat, "teal", 1)
            end
            fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_base, alturaPiso * k, fig, ax, ax_mat, "teal", 1)
        end
    end

    ps_SombraVolTeor_p, ps_SombraVolTeor_o, ps_SombraVolTeor_s = generaSombraTeor_v3(ps_volteor, matConexionVertices_ss, vecVertices_ss, ps_publico, ps_calles)
    if f_sombraVolTeor_p
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_SombraVolTeor_p, 0, fig, ax, ax_mat, "gold", 0.3)
    end
    if f_sombraVolTeor_o
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_SombraVolTeor_o, 0, fig, ax, ax_mat, "gold", 0.3)
    end
    if f_sombraVolTeor_s
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_SombraVolTeor_s, 0, fig, ax, ax_mat, "gold", 0.3)
    end

    ps_sombraEdif_p, ps_sombraEdif_o, ps_sombraEdif_s = generaSombraEdificio(ps_baseSeparada, alt, ps_publico, ps_calles)
    if f_sombraEdif_p
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_sombraEdif_p, 0, fig, ax, ax_mat, "red", 0.25)
    end
    if f_sombraEdif_o
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_sombraEdif_o, 0, fig, ax, ax_mat, "red", 0.25)
    end
    if f_sombraEdif_s
        fig, ax, ax_mat = polyShape.plotPolyshape3d_v3(ps_sombraEdif_s, 0, fig, ax, ax_mat, "red", 0.25)
    end

    return fig
end


text:
"""
JuMP constraints
"""
sense_to_set(::Function, ::Val{:!=}) = NotEqualTo(0.0)
sense_to_set(::Function, ::Val{:<}) = Strictly(MOI.LessThan(0.0))
sense_to_set(::Function, ::Val{:>}) = Strictly(MOI.GreaterThan(0.0))

MOIU.shift_constant(set::NotEqualTo, value) = NotEqualTo(set.value + value)

include("indicator.jl")
include("reified.jl")
include("and.jl")
include("or.jl")

"""
MOI constraints
"""

"""
Linear constraints
"""
MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{MOI.EqualTo{T}},
) where {T<:Real} = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{MOI.LessThan{T}},
) where {T<:Real} = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{EqualSetInternal},
) = true
MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{AllDifferentSetInternal},
) = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{TableSetInternal},
) = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{NotEqualTo{T}},
) where {T<:Real} = true

# Don't directly support StrictlyGreaterThan => use a bridge
MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{Strictly{T, MOI.LessThan{T}}},
) where {T<:Real} = true

function MOI.supports_constraint(
    ::Optimizer,
    func::Type{VAF{T}},
    set::Type{OS},
) where {A,T<:Real,IS<:MOI.AbstractScalarSet,OS<:MOI.IndicatorSet{A,IS}}
    if IS <: MOI.GreaterThan || IS <: Strictly{T, MOI.GreaterThan{T}}
        return false
    end
    return A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{VAF{T}},
    set::Type{OS},
) where {A,T<:Real,IS,OS<:CS.IndicatorSet{A,IS}}
    if IS <: BoolSet
        return is_boolset_supported(optimizer, IS)
    end
    return A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{MOI.VectorOfVariables},
    set::Type{OS},
) where {A,IS,OS<:CS.IndicatorSet{A,IS}}
    !(A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO) && return false
    return MOI.supports_constraint(optimizer, func, IS)
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{MOI.VectorOfVariables},
    set::Type{OS}
) where {A,IS,OS<:CS.ReifiedSet{A,IS}}
    !(A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO) && return false
    return MOI.supports_constraint(optimizer, func, IS)
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{VAF{T}},
    set::Type{OS},
) where {A,T<:Real,IS,OS<:CS.ReifiedSet{A,IS}}
    if IS <: MOI.GreaterThan || IS <: Strictly{T, MOI.GreaterThan{T}}
        return false
    end
    if IS <: BoolSet
        return is_boolset_supported(optimizer, IS)
    end
    return A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{VAF{T}},
    set::Type{BS},
) where {T,F1,F2,F1dim,F2dim,S1,S2,BS<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
    return is_boolset_supported(optimizer, set)
end

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{GeqSetInternal},
) = true

"""
    Return whether the two constraint inside the `BoolSet` are supported directly by the solver
"""
function is_boolset_supported(optimizer::Optimizer, ::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}) where {F1, F2, F1dim, F2dim, S1, S2}
    is_supported = MOI.supports_constraint(optimizer, F1, S1) && MOI.supports_constraint(optimizer, F2, S2)
    return is_supported
end

function check_inbounds(model::Optimizer, aff::SAF{T}) where {T<:Real}
    for term in aff.terms
        check_inbounds(model, term.variable_index)
    end
    return
end

function check_inbounds(model::Optimizer, vov::MOI.VectorOfVariables)
    for var in vov.variables
        check_inbounds(model, var)
    end
    return
end

"""
    add_constraint!(model::Optimizer, constraint::Constraint)

Calls `add_constraint!` on the `CoM` without setting pvals.
The variable bounds might not be ready at this stage.
"""
function add_constraint!(model::Optimizer, constraint::Constraint)
    add_constraint!(model.inner, constraint; set_pvals = false)
end

"""
    create_interals(com::CoM, vars, set)

Create ConstraintInternals for a vector of variables constraint
"""
function create_interals(com::CoM, vars, set)
    internals = create_interals(vars, set)
    internals.idx = length(com.constraints) + 1
    return internals
end

"""
    create_interals(vars::MOI.VectorOfVariables, set)

Create ConstraintInternals for a vector of variables constraint
"""
function create_interals(vars::MOI.VectorOfVariables, set)
    ConstraintInternals(
        0,
        vars,
        set,
        Int[v.value for v in vars.variables],
    )
end

"""
    create_interals(func::VAF{T}, set) where {T}

Create ConstraintInternals for a vector of variables constraint
"""
function create_interals(func::VAF{T}, set) where {T}
    ConstraintInternals(
        0,
        func,
        set,
        get_indices(func)
    )
end

"""
    get_anti_constraint(mode, constraint)

Return the anti constraint if it exists and `nothing` otherwise.
The anti constraint is the constraint that expresses the opposite i.e
input: 2x + 7 <= 5 => 2x + 7 > 5
 - it will actually output only less than constraints not great than as it's not supported
input 5x == 2 => 5x != 2

Currently it's only implemented for linear constraints
"""
function get_anti_constraint(model, constraint::Constraint)
    return nothing
end

function get_anti_constraint(model, constraint::LinearConstraint{T}) where T
    set = constraint.set
    anti_fct = nothing
    anti_set = nothing
    if constraint.set isa MOI.LessThan
        anti_fct = MOIU.operate(-, T, constraint.fct)
        anti_set = Strictly(MOI.LessThan(-set.upper))
    elseif constraint.set isa Strictly{T, MOI.LessThan{T}}
        anti_fct = MOIU.operate(-, T, constraint.fct)
        anti_set = MOI.LessThan(-set.set.upper)
    elseif constraint.set isa MOI.EqualTo
        anti_fct = copy(constraint.fct)
        anti_set = NotEqualTo(set.value)
    elseif constraint.set isa NotEqualTo
        anti_fct = copy(constraint.fct)
        anti_set = MOI.EqualTo(set.value)
    end
    anti_lc = new_linear_constraint(model, anti_fct, anti_set)
    anti_lc.idx = 0
    return anti_lc
end

function get_anti_constraint(model, constraint::BoolConstraint)
    lhs_anti_constraint = get_anti_constraint(model, constraint.lhs)
    rhs_anti_constraint = get_anti_constraint(model, constraint.rhs)

    if lhs_anti_constraint === nothing || rhs_anti_constraint === nothing
        return nothing
    end

    T = parametric_type(model.inner)
    fct = MOIU.operate(vcat, T, lhs_anti_constraint.fct, rhs_anti_constraint.fct)
    return anti_bool_constraint(constraint, fct, lhs_anti_constraint, rhs_anti_constraint)
end

"""
    anti_bool_constraint(::AndConstraint, fct, lhs_constraint::Constraint, rhs_constraint::Constraint)

Return the OrConstraint with the already anti constraints `lhs_constraint` and `rhs_constraint`
::AndConstraint is just for dispatching ;)
"""
function anti_bool_constraint(::AndConstraint, fct, lhs_constraint::Constraint, rhs_constraint::Constraint)
    set = OrSet{typeof(lhs_constraint.fct), typeof(rhs_constraint.fct)}(lhs_constraint.set, rhs_constraint.set)

    internals = ConstraintInternals(0,fct,set,get_indices(fct))
    return OrConstraint(internals, lhs_constraint, rhs_constraint)
end

function anti_bool_constraint(::OrConstraint, fct, lhs_constraint::Constraint, rhs_constraint::Constraint)
    set = AndSet{typeof(lhs_constraint.fct), typeof(rhs_constraint.fct)}(lhs_constraint.set, rhs_constraint.set)

    internals = ConstraintInternals(0,fct,set,get_indices(fct))
    return AndConstraint(internals, lhs_constraint, rhs_constraint)
end

"""
    MOI.add_constraint(
        model::Optimizer,
        vars::MOI.VectorOfVariables,
        set::MOI.AbstractVectorSet,
    )

Add all kinds of vector of variables constraints like:
TableConstraint and AllDifferentConstraint
"""
function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::MOI.AbstractVectorSet,
)
    check_inbounds(model, vars)
    com = model.inner

    internals = create_interals(com, vars, set)

    constraint = init_constraint_struct(set, internals)

    add_constraint!(model, constraint)
    if set isa AllDifferentSetInternal
        com.info.n_constraint_types.alldifferent += 1
    elseif set isa TableSetInternal
        com.info.n_constraint_types.table += 1
    elseif set isa EqualSetInternal
        com.info.n_constraint_types.equality += 1
    end

    return MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::MOI.EqualTo{T},
) where {T<:Real}
    check_inbounds(model, func)

    if length(func.terms) == 1
        vidx = func.terms[1].variable_index.value
        val = convert(Int, set.value / func.terms[1].coefficient)
        push!(model.inner.init_fixes, (vidx, val))
        return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(0)
    elseif length(func.terms) == 2 && set.value == zero(T)
        if func.terms[1].coefficient == -func.terms[2].coefficient
            # we have the form a == b
            vecOfvar = MOI.VectorOfVariables([
                func.terms[1].variable_index,
                func.terms[2].variable_index,
            ])
            com = model.inner
            internals = ConstraintInternals(
                length(com.constraints) + 1,
                vecOfvar,
                CS.EqualSetInternal(2),
                Int[v.value for v in vecOfvar.variables],
            )
            constraint = EqualConstraint(internals, ones(Int, 2))

            add_constraint!(model, constraint)
            com.info.n_constraint_types.equality += 1
            return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(length(model.inner.constraints))
        end
    end

    lc = new_linear_constraint(model, func, set)

    add_constraint!(model, lc)
    model.inner.info.n_constraint_types.equality += 1

    return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(length(model.inner.constraints))
end

function add_variable_less_than_variable_constraint(
    model::Optimizer,
    func::SAF{T},
    set::MOI.LessThan{T},
) where {T<:Real}
    reverse_order = false
    if func.terms[1].coefficient != 1.0 || func.terms[2].coefficient != -1.0
        if func.terms[1].coefficient == -1.0 && func.terms[2].coefficient == 1.0
            # rhs is lhs and other way around
            reverse_order = true
        end
    end

    com = model.inner

    if reverse_order
        lhs = func.terms[2].variable_index.value
        rhs = func.terms[1].variable_index.value
    else
        lhs = func.terms[1].variable_index.value
        rhs = func.terms[2].variable_index.value
    end

    internals = ConstraintInternals(
        length(model.inner.constraints) + 1, # constraint idx
        func,
        set,
        [lhs, rhs],
    )
    svc = SingleVariableConstraint(internals, lhs, rhs)

    add_constraint!(model, svc)
    com.info.n_constraint_types.inequality += 1

    return MOI.ConstraintIndex{SAF{T},typeof(set)}(length(model.inner.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::MOI.LessThan{T},
) where {T<:Real}
    check_inbounds(model, func)

    # support for a <= b which is written as a-b <= 0
    # currently only supports if the coefficients are 1 and -1
    if set.upper == 0.0 &&
       length(func.terms) == 2 &&
       abs(func.terms[1].coefficient) == 1.0 &&
       abs(func.terms[2].coefficient) == 1.0 &&
       func.terms[1].coefficient == -func.terms[2].coefficient
        return add_variable_less_than_variable_constraint(model, func, set)
    end

    lc = new_linear_constraint(model, func, set)

    add_constraint!(model, lc)
    model.inner.info.n_constraint_types.inequality += 1

    return MOI.ConstraintIndex{SAF{T},typeof(set)}(length(model.inner.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::NotEqualTo{T},
) where {T<:Real}
    com = model.inner
    com.info.n_constraint_types.notequal += 1

    # support for cx != a where a and c are constants
    if length(func.terms) == 1
        # reformulate to x != a where x is variable and a a constant
        rhs = set.value / func.terms[1].coefficient
        if isapprox_discrete(com, rhs)
            rhs = get_approx_discrete(rhs)
            rm!(
                com,
                com.search_space[func.terms[1].variable_index.value],
                rhs;
                changes = false,
            )
        end
        return MOI.ConstraintIndex{SAF{T},NotEqualTo{T}}(0)
    end

    lc = new_linear_constraint(model, func, set)

    add_constraint!(model, lc)

    return MOI.ConstraintIndex{SAF{T},NotEqualTo{T}}(length(com.constraints))
end
function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::Strictly{T, MOI.LessThan{T}},
) where {T<:Real}
    check_inbounds(model, func)

    lc = new_linear_constraint(model, func, set)

    add_constraint!(model, lc)
    model.inner.info.n_constraint_types.inequality += 1

    return MOI.ConstraintIndex{SAF{T},typeof(set)}(length(model.inner.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::IS,
) where {A,T<:Real,ASS<:MOI.AbstractScalarSet,IS<:MOI.IndicatorSet{A,ASS}}
    com = model.inner
    com.info.n_constraint_types.indicator += 1

    indices = get_indices(func)

    internals = ConstraintInternals(
        length(com.constraints)+1,
        func,
        MOI.IndicatorSet{A}(set.set),
        indices,
    )

    activator_internals = get_activator_internals(A, indices)

    lc = get_inner_constraint(func, set, set.set)

    constraint = IndicatorConstraint(internals, activator_internals, lc)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},MOI.IndicatorSet{A,ASS}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::IS,
) where {A,S,IS<:CS.IndicatorSet{A,S}}
    com = model.inner
    com.info.n_constraint_types.indicator += 1

    internals = create_interals(com, vars, set)

    inner_constraint = get_inner_constraint(vars, set, set.set)

    indices = internals.indices
    activator_internals = get_activator_internals(A, indices)
    constraint =
        IndicatorConstraint(internals, activator_internals, inner_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{MOI.VectorOfVariables,CS.IndicatorSet{A,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::IS,
) where {T,A,S<:MOI.AbstractVectorSet,IS<:CS.IndicatorSet{A,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1

    internals = create_interals(com, func, set)

    inner_constraint = get_inner_constraint(func, set, set.set)
    indices = internals.indices

    activator_internals = get_activator_internals(A, indices)

    constraint =
        IndicatorConstraint(internals, activator_internals, inner_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},CS.IndicatorSet{A,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::RS,
) where {A,S<:MOI.AbstractScalarSet,T<:Real,RS<:ReifiedSet{A,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1

    indices = get_indices(func)

    internals = ConstraintInternals(
        length(com.constraints)+1,
        func,
        typeof(set)(set.set, set.dimension),
        indices,
    )

    # for normal linear constraints
    lc = get_inner_constraint(func, set, set.set)
    anti_lc = get_anti_constraint(model, lc)

    activator_internals = get_activator_internals(A, indices)
    constraint = ReifiedConstraint(internals, activator_internals, lc, anti_lc)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},CS.ReifiedSet{A,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::RS,
) where {A,S<:MOI.AbstractVectorSet,RS<:CS.ReifiedSet{A,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1
    internals = create_interals(com, vars, set)

    inner_constraint = get_inner_constraint(vars, set, set.set)
    anti_constraint = get_anti_constraint(model, inner_constraint)
    indices = internals.indices
    activator_internals = get_activator_internals(A, indices)
    constraint =
        ReifiedConstraint(internals, activator_internals, inner_constraint, anti_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{MOI.VectorOfVariables,CS.ReifiedSet{A,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::RS,
) where {T,A,S<:MOI.AbstractVectorSet,RS<:CS.ReifiedSet{A,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1

    internals = create_interals(com, func, set)

    inner_constraint = get_inner_constraint(func, set, set.set)
    anti_constraint = get_anti_constraint(model, inner_constraint)
    indices = internals.indices
    activator_internals = get_activator_internals(A, indices)
    constraint =
        ReifiedConstraint(internals, activator_internals, inner_constraint, anti_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},CS.ReifiedSet{A,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::BS,
) where {T,BS<:BoolSet}
    com = model.inner
    internals = create_interals(com, func, set)
    constraint = init_constraint_struct(set, internals)
    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},typeof(set)}(length(com.constraints))
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMMaterials.jl/blob/master/LICENSE

# # 3D Beam with ideal plastic material model from Materials.jl

using JuliaFEM, FEMMaterials, Materials, FEMBase, LinearAlgebra, Plots
import FEMMaterials: Continuum3D, MecaMatSo
pkg_dir = dirname(dirname(pathof(FEMMaterials)))

# ## Let's read the discretized geometry and create boundary conditions
#
# File `plactic_beam.inp` is created with 3rd party meshing tool.
# File contains surface sets `BC1`, `BC2` and `PRESSURE` as well
# element set `Body1`
mesh = abaqus_read_mesh(joinpath(pkg_dir,"examples","data_3dbeam","plastic_beam.inp"))
beam_elements = create_elements(mesh, "Body1")
bc_elements_1 = create_nodal_elements(mesh, "BC1")
bc_elements_2 = create_nodal_elements(mesh, "BC2")
trac_elements = create_surface_elements(mesh, "PRESSURE")

for j in 1:3
    update!(bc_elements_1, "displacement $j", 0.0)
end
update!(bc_elements_2, "displacement 1", 0.0)
update!(bc_elements_2, "displacement 2", 0.0)
update!(trac_elements, "surface pressure", 0.0 => 0.00)
update!(trac_elements, "surface pressure", 1.0 => 2.70)
trac = Problem(Continuum3D, "traction", 3)
bc = Problem(Dirichlet, "fix displacement", 3, "displacement")
add_elements!(trac, trac_elements)
add_elements!(bc, bc_elements_1)
add_elements!(bc, bc_elements_2)

# ## Next, we set the material properties for each element set
#
# In this example we only have the one element set `Body1` in variable `beam_elements`

update!(beam_elements, "youngs_modulus", 200.0e3)
update!(beam_elements, "poissons_ratio", 0.3)
update!(beam_elements, "yield_stress", 100.0)

beam = Problem(Continuum3D, "plastic beam", 3)
beam.properties.material_model = :IdealPlastic
add_elements!(beam, beam_elements)

# ## And next, we setup the analysis
#
# `t0`, `t1`, and `dt` are the start time, end time, and time step respectively.

analysis = Analysis(MecaMatSo, "solve problem")
analysis.properties.max_iterations = 50
analysis.properties.t0 = 0.0
analysis.properties.t1 = 1.0
analysis.properties.dt = 0.05

# ## Writing the results needs to be setup as well

xdmf = Xdmf("3dbeam_results_output"; overwrite=true)
add_results_writer!(analysis, xdmf)

# ## Finally, adding the problems together and running the analysis
#
# All earlier defined problems are added together. Also result file need to be
# closed to flush everything from the writing buffer to the file.
add_problems!(analysis, beam, trac, bc)
run!(analysis)
close(xdmf)

# ## The first post-processing step is to calculate maximum von Mises stresses
#
# `vmis` contains all integration points stresses and `vmis_` just the maximum.
# Let's plot the maximum von Mises stress as a function of time

tim = 0.0:0.05:1.0
vmis_ = []
for t in tim
    vmis = []
    for element in beam_elements
        for ip in get_integration_points(element)
            s11, s22, s33, s12, s23, s31 = ip("stress", t)
            push!(vmis, sqrt(1/2*((s11-s22)^2 + (s22-s33)^2 + (s33-s11)^2 + 6*(s12^2+s23^2+s31^2))))
        end
    end
    push!(vmis_, maximum(vmis))
end
plot(tim,vmis_)
#-
png("max_vonmises_stress_as_a_function_of_time")

# ![max_vonmises_stress_as_a_function_of_time][vonmises]
# [vonmises]: https://raw.githubusercontent.com/JuliaFEM/FEMMaterials.jl/master/notebooks/max_vonmises_stress_as_a_function_of_time.png

# ## The second post-processing step is to collect displacements
#
# Here as an example node number 96 second degree of freedom displacement is
# extracted.

u2_96 = []
for t in tim
    push!(u2_96, beam("displacement", t)[96][2])
end
plot(tim,u2_96)
#-
png("node_96_displacement_as_a_function_of_time")

# ![node_96_displacement_as_a_function_of_time][displacement]
# [displacement]: https://raw.githubusercontent.com/JuliaFEM/FEMMaterials.jl/master/notebooks/node_96_displacement_as_a_function_of_time.png


text:
@register Base.getindex(x,i::Integer) false
@register Base.getindex(x,i) # define one and only one promotion rule
@register Base.binomial(n,k)

Base.sign(x::Symbolic) = Term{Int}(sign, [x])
Base.sign(x::Num) = Num(sign(value(x)))
derivative(::typeof(sign), args::NTuple{1,Any}, ::Val{1}) = 0
@register Base.signbit(x)::Bool
derivative(::typeof(signbit), args::NTuple{1,Any}, ::Val{1}) = 0
derivative(::typeof(abs), args::NTuple{1,Any}, ::Val{1}) = IfElse.ifelse(signbit(args[1]),-one(args[1]),one(args[1]))
function derivative(::typeof(min), args::NTuple{2,Any}, ::Val{1})
    x, y = args
    IfElse.ifelse(x < y, one(x), zero(x))
end
function derivative(::typeof(min), args::NTuple{2,Any}, ::Val{2})
    x, y = args
    IfElse.ifelse(x < y, zero(y), one(y))
end
function derivative(::typeof(max), args::NTuple{2,Any}, ::Val{1})
    x, y = args
    IfElse.ifelse(x > y, one(x), zero(x))
end
function derivative(::typeof(max), args::NTuple{2,Any}, ::Val{2})
    x, y = args
    IfElse.ifelse(x > y, zero(y), one(y))
end

@register Base.rand(x)
@register Base.randn(x)

@register Distributions.pdf(dist,x)
@register Distributions.logpdf(dist,x)
@register Distributions.cdf(dist,x)
@register Distributions.logcdf(dist,x)
@register Distributions.quantile(dist,x)

@register Distributions.Uniform(mu,sigma) false
@register Distributions.Normal(mu,sigma) false

@register ∈(x::Num, y::AbstractArray)
@register ∪(x, y)
@register ∩(x, y)
@register ∨(x, y)
@register ∧(x, y)
@register ⊆(x, y)

LinearAlgebra.norm(x::Num, p::Real) = abs(x)


text:
using OmegaCore
using OmegaTest
using OmegaCore.Traits

using Test
using Distributions
using Random

function testtraits()
  ω = SimpleΩ(Dict())
  ω2 = OmegaCore.appendscope(ω, [1,])
  # traits(typeof(ω2.tags))
  g(traits(typeof(ω2)))
end

function testsimple1()
  ω = defΩ()()
  ω = OmegaCore.tagrng(ω, Random.GLOBAL_RNG)
  ω2 = defΩ()()
  ω2 = OmegaCore.tagrng(ω2, Random.GLOBAL_RNG)
  x = 1 ~ Normal(0, 1)
  map(x, [ω, ω2])
end

function testsimple2()
  x = 1 ~ Normal(0, 1)
  rng = Random.GLOBAL_RNG
  ΩT = defΩ()
  y = OmegaCore.condvar(x)
  ω =  OmegaCore.OmegaRejectionSample.condomegasample1(rng, ΩT, y, OmegaCore.RejectionSample)
  x(ω)
end

function testsimple3a()
  x = 1 ~ Normal(0, 1)
  rng = Random.GLOBAL_RNG
  ΩT = defΩ()
  y = OmegaCore.condvar(x)
  ω =  OmegaCore.OmegaRejectionSample.condomegasample1(rng, ΩT, y, OmegaCore.RejectionSample)
  map(x, typeof(ω)[ω, ω])
end

function testsimple3()
  x = 1 ~ Normal(0, 1)
  rng = Random.GLOBAL_RNG
  ΩT = defΩ()
  y = OmegaCore.condvar(x)
  ω =  OmegaCore.OmegaRejectionSample.condomegasample1(rng, ΩT, y, OmegaCore.RejectionSample)
  x.([ω, ω])
end

function testsimple4()
  x = 1 ~ Normal(0, 1)
  rng = Random.GLOBAL_RNG
  ΩT = defΩ()
  y = OmegaCore.condvar(x)
  ω =  OmegaCore.OmegaRejectionSample.condomegasample(rng, ΩT, y, 5, OmegaCore.RejectionSample)
  map(x, ω)
end

@testset "infer types" begin
  @test isinferred(testsimple1)
  @test isinferred(testsimple2)
  @test isinferred(testsimple3)
  @test isinferred(testsimple3a)
  @test isinferred(testsimple4)
end


text:
# # Getting started
#
# ## Installation
#
# To install the package, hit `]` in the Julia REPL to enter the package manager mode, then run
#
# ```julia
# pkg> add ProximalAlgorithms
# ```
#
# to get the latest stable version, or
#
# ```julia
# pkg> add ProximalAlgorithms#master
# ```
#
# the latest development version (`master` branch).
#
# ## Overview
#
# The methods implemented in ProximalAlgorithms are commonly referred to as (you've guessed it) *proximal algorithms*,
# in that they rely on the *proximal operator* (or *mapping*) to deal with non-differentiable terms in the objective.
# Loosely speaking, the algorithms in this package can be used to solve problems of the form
#
# ```math
# \operatorname*{minimize}_x\ \sum_{i=1}^N f_i(x)
# ```
#
# where ``N`` depends on the specific algorithm, together with specific assumptions on the terms ``f_i`` (like smoothness, convexity, strong convexity).
# The problem above is solved by iteratively accessing specific *first order information* on the terms ``f_i``,
# like their gradient ``\nabla f_i`` or their proximal mapping ``\operatorname{prox}_{f_i}``:
# ```math
# \mathrm{prox}_{\gamma f_i}(x) = \arg\min_z \left\{ f_i(z) + \tfrac{1}{2\gamma}\|z-x\|^2 \right\}
# ```
# The literature on proximal operators and algorithms is vast: for an overview, one can refer to [Parikh2014](@cite), [Beck2017](@cite).
#
# To evaluate these first-order primitives, ProximalAlgorithms does the following:
# * ``\nabla f_i`` falls back to using automatic differentiation (as provided by [Zygote](https://github.com/FluxML/Zygote.jl)).
# * ``\operatorname{prox}_{f_i}`` relies on the intereface of [ProximalOperators](https://github.com/JuliaFirstOrder/ProximalOperators.jl).
# Both of the above can be implemented for custom function types, as [documented here](@ref custom_terms).
#
#md # !!! note
#md #
#md #     Each of the [implemented algorithms](@ref problems_algorithms) assumes a different structure of the objective to be optimized
#md #     (e.g., a specific number ``N`` of terms), with specific assumptions (e.g., smoothness, convexity).
#md #     Furthermore, multiple algorithms can often be applied to the same problem (possibly through appropriate reformulation) and are expected
#md #     to perform differently; sometimes, even the *same algorithm* can be applied in multiple ways to the same problem, by grouping (*splitting*)
#md #     the terms in different ways.
#md #
#md #     Because of these reasons, ProximalAlgorithms **does not offer a modeling language** to automagically minimize any objective:
#md #     rather, the user is expected to formulate their problem by providing the right objective terms to the algorithm of choice.
#md #     Please refer to the [this section of the manual](@ref problems_algorithms) for information on what terms can be provided and under which assumptions.
#
# ## [Example: box constrained quadratic](@id box_qp)
#
# As a four-line example, consider the minimization of a 2D quadratic function subject to box constraints:

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms

quadratic_cost(x) = dot([3.4 1.2; 1.2 89.1] * x, x) / 2 + dot([2.3, 99.9], x)
box_indicator = ProximalOperators.IndBox(-1, 1)

ffb = ProximalAlgorithms.FastForwardBackward(maxit=1000, tol=1e-5, verbose=true)
solution, iterations = ffb(zeros(2), f=quadratic_cost, g=box_indicator)

# Here, we defined the cost function `quadratic_cost`, and the constraint indicator `box_indicator`.
# Then we set up the optimization algorithm as the fast forward-backward splitting method (fast proximal gradient),
# with options for the maximum number of iterations, termination tolerance, verbosity.
# Finally, running the algorithm requires an initial point and the objective terms.
#
# We can verify the correctness of the solution by checking that the negative gradient is orthogonal to the constraints, pointing outwards:

-ProximalAlgorithms.gradient(quadratic_cost, solution)[1]

# Or by plotting the solution against the cost function and constraint:

using Gadfly

contour = layer(
    z=(x,y) -> quadratic_cost([x, y]), xmin=[-2], xmax=[2], ymin=[-2], ymax=[2],
    Geom.contour(levels=[quadratic_cost(solution) + k * 12 for k in 0:100]),
)
point = layer(x=[solution[1]], y=[solution[2]])
box = layer(x=[-1, -1, 1, 1, -1], y=[1, -1, -1, 1, 1], Geom.path)
plot(contour, box, point, Guide.xlabel(nothing), Guide.ylabel(nothing))

# ## [Iterator interface](@id iterator_interface)
#
# Under the hood, all algorithms are implemented as standard Julia iterators. In the above example,
# the `FastForwardBackward` algorithm internally uses a [`FastForwardBackwardIteration`](@ref ProximalAlgorithms.FastForwardBackwardIteration)
# iterator object.
# Each iterator outputs the complete iteration state: interacting directly with this object allows for more
# fine-grained inspection and control over the algorithm, for example to report the algorithm's progress in a
# custom way (cost decrease, or whatever other metric), or to interrupt the optimization based on a
# custom stopping criterion.
#
#md # !!! note
#md #
#md #     Since each algorithm iterator type has its own logic, it will also have its own dedicated state structure.
#md #     Interacting with the state then requires being familiar with its structure, and with the nature of its attributes.
#md #
#
# For example, the optimization can be customized as follows: here we print the value of the cost function
# every ten iterations, and stop on a custom condition on the fixed-point residual norm.

ffbiter = ProximalAlgorithms.FastForwardBackwardIteration(x0=zeros(2), f=quadratic_cost, g=box_indicator)

for (k, state) in enumerate(ffbiter)
    if mod(k, 10) == 0
        println("iteration #$k: $(state.f_x)")
    end
    if norm(state.res) / state.gamma <= 1e-4
        break
    end
end

#md # !!! warning
#md #
#md #     To save on allocations, most (if not all) algorithms re-use state objects when iterating,
#md #     by updating the state *in place* instead of creating a new one. For this reason, one **should not**
#md #     - mutate the state object in any way, as this may corrupt the algorithm's logic,
#md #     - `collect` the sequence of states, since this will result in an array of identical objects.
#md #
# ## [Example: sparse linear regression](@id sparse_linreg)
#
# Let's look at a least squares regression problem with L1 regularization:
# we will use the "diabetes dataset" (see [here](https://www4.stat.ncsu.edu/~boos/var.select/)), so let's start by loading the data.

using HTTP

splitlines(s) = split(s, "\n")
splitfields(s) = split(s, "\t")
parsefloat64(s) = parse(Float64, s)

function load_diabetes_dataset()
    res = HTTP.request("GET", "https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt")
    lines = res.body |> String |> strip |> splitlines
    return hcat((line |> splitfields .|> parsefloat64 for line in lines[2:end])...)'
end

data = load_diabetes_dataset()

training_input = data[1:end-100, 1:end-1]
training_label = data[1:end-100, end]

test_input = data[end-99:end, 1:end-1]
test_label = data[end-99:end, end]

n_training, n_features = size(training_input)

# Now we can set up the optimization problem we want to solve: we will minimize the mean squared error
# for a linear model, appropriately scaled so that the features in the training data are normally distributed
# ("standardization", this is known to help the optimization process).
#
# After some simple manipulation, this standardized linear model can be implemented as follows:

using LinearAlgebra
using Statistics

input_loc = mean(training_input, dims=1) |> vec
input_scale = std(training_input, dims=1) |> vec

linear_model(wb, input) = input * wb[1:end-1] .+ wb[end]

function standardized_linear_model(wb, input)
    w_scaled = wb[1:end-1] ./ input_scale
    wb_scaled = vcat(w_scaled, wb[end] - dot(w_scaled, input_loc))
    return linear_model(wb_scaled, input)
end

# The loss term in the cost is then the following. Note that this is a regular Julia function:
# since the algorithm we will apply requires its gradient, automatic differentiation will
# do the work for us.

mean_squared_error(label, output) = mean((output .- label) .^ 2) / 2

training_loss(wb) = mean_squared_error(training_label, standardized_linear_model(wb, training_input))

# As regularization we will use the L1 norm, implemented in [ProximalOperators](https://github.com/JuliaFirstOrder/ProximalOperators.jl):

using ProximalOperators

reg = ProximalOperators.NormL1(1)

# We want to minimize the sum of `training_loss` and `reg`, and for this task we can use `FastForwardBackward`,
# which implements the fast proximal gradient method (also known as fast forward-backward splitting, or FISTA).
# Therefore we construct the algorithm, then apply it to our problem by providing a starting point,
# and the objective terms `f=training_loss` (smooth) and `g=reg` (non smooth).

using ProximalAlgorithms

ffb = ProximalAlgorithms.FastForwardBackward()
solution, iterations = ffb(zeros(n_features + 1), f=training_loss, g=reg)

# We can now check how well the trained model performs on the test portion of our data.

test_output = standardized_linear_model(solution, test_input)
mean_squared_error(test_label, test_output)


text:
# TODO: this is really a module for face-centered diamond-cubic
#       and over time it should become much more general of course

module Si

using JuLIP, JuLIP.Potentials, QuadGK, ForwardDiff, LinearAlgebra

import JuLIPMaterials

using JuLIPMaterials.CLE: elastic_moduli, voigt_moduli,
         fourth_order_basis, sextic_roots, A_coefficients, D_coefficients,
         little_a

const CLE = JuLIPMaterials.CLE
const FCC = JuLIPMaterials.FCC
# CauchyBorn = JuLIPMaterials.CauchyBorn

"""
`si110_plane(s::AbstractString) -> at::ASEAtoms, b, xcore, a `

Generates a unit cell for an FCC crystal with orthogonal cell vectors chosen
such that the F1 direction is the burgers vector and the F3 direction the normal
to the standard edge dislocation:
   b = F1 ~ a1;  ν = F3 ~ a2-a3
The third cell vector: F2 ~ a * e1. Here, ~ means they are rotated from the
a1, a2, a3 directions. This cell contains two atoms.

Returns
* `at`: Atoms object
* `b` : Burgers vector
* `xcore` : a core-offset (to add to any lattice position)
* `a` : lattice parameter
"""
function si110_plane(s::Symbol;
                     a = cell(bulk(s, cubic=true))[1,1] )
   # TODO: can si110_plane be combined with FCC.fcc_110_plane?
   @assert s == :Si
   # for Si, default a = 5.43
   # TODO: ensure only that s is an FCC species
   # ====================================================
   # construct the cell matrix
   F = a * JMat( [ sqrt(2)/2  0    0;
                         0    1    0;
                         0    0    sqrt(2)/2 ] )
   X = a * [ JVec([         0.0,   0.0,          0.0 ]),
             JVec([ 0.5/sqrt(2),   0.5,  0.5/sqrt(2) ]),
             JVec([         0.0, -0.25,  0.5/sqrt(2) ]),
             JVec([ 0.5/sqrt(2),  0.25,          0.0 ]) ]
   # construct ASEAtoms
   at = Atoms(s, X)
   set_cell!(at, F')
   # compute a burgers vector in these coordinates
   b = a * sqrt(2)/2 * JVec([1.0,0.0,0.0])
   # compute a core-offset (to add to any lattice position)
   xcore = [-.7, 1.0, 0] * a / 5.43
   # return the information
   return at, b, xcore, a
end

"a fully equilibrated SW potential"
function sw_eq()
    T(σ, at) = tr(stress(StillingerWeber(σ=σ), at))
    at = bulk(:Si, pbc=true)
    r0 = 2.09474
    r1 = r0 - 0.1
    s0, s1 = T(r0, at), T(r1, at)
    while (abs(s1) > 1e-8) && abs(r0 - r1) > 1e-8
        rnew = (r0 * s1 - r1 * s0) / (s1 - s0)
        r0, r1 = r1, rnew
        s0, s1 = s1, T(rnew, at)
    end
    return StillingerWeber(σ=r1)
end

function si110_cluster(sym, L; kwargs...)
    # @assert isodd(L)   # TODO: why?
    atu, b, _, a = si110_plane(sym; kwargs...)
    at = atu * (L, L, 1)
    set_pbc!(at, (false, false, true))
    b = b[1]
    X = positions(at)
    # This choice picks the lower left and upper right atom (not site) positions
    #   TODO I don't like the [-1, 0.1, 0] at all
    xcore = (1/2)*(X[length(X)-2]+X[3]) + [-1, 0.1, 0]
    return at, b, xcore
end

# """
# a function that identifies multi-lattice structure in 2 layers of bulk-Si
# (yes - very restrictive but will do for now!)
# """
# function si_multilattice(at)
#     TOL = 0.2
#     J0 = Int[]
#     J1 = Int[]
#     Jdel = Int[]
#     for (i, j, r, R, _) in sites(at, rnn("Si")+0.1)
#         foundneig = false
#         for (jj, RR) in zip(j, R)
#             if (RR[1] == 0.0) && (abs(RR[2] - 1.3575) < TOL)
#                 # neighbour above >> make (i, jj) a site
#                 push!(J0, i)
#                 push!(J1, jj)
#                 foundneig = true
#                 break
#             elseif (RR[1] == 0.0) && (abs(RR[2] + 1.3575) < TOL)
#                 # neighbour below >> (jj, i) is a site that will be pushed when i ↔ jj
#                 foundneig = true
#                 break
#             end
#         end
#         if !foundneig
#             # i has no neighbour above or below >> probably we just get rid of it
#             push!(Jdel, i)
#         end
#     end
#     return J0, J1, Jdel
# end
#
#
#
# function symml_displacement!(at, u)
#     I0, I1, Idel = si_multilattice(at)
#     @assert isempty(Idel)  # if Idel is not empty then (for now) we don't know what to do
#     X = positions(at)
#     W = CauchyBorn.WcbQuad(calculator(at))   # TODO: generalise this to general calculators
#     F0 = cell(W.at)'
#     p0 = W(F0)
#     # transformation matrices
#     Tp = [0 1/√2  -1/√2; 1 0 0; 0 1/√2 1/√2]
#     Tm = diagm([1,1,-1]) * Tp
#
#     for (i0, i1) in zip(I0, I1)   # each pair (i0, i1) corresponds to a ML lattice site
#         x0, x1 = X[i0], X[i1]
#         x1[3] > x0[3] ? T = Tp : T = Tm
#         x̄ = 0.5 * (x0 + x1)   # centre of mass of the bond
#         U, ∇U = u(x̄)          # displacement and displacement gradient
#         F = T' * (I + ∇U) * T
#         q = T * (W(F * F0) - p0)    # construct the shift corresponding to F = Id + ∇U
#         X[i0], X[i1] = x0 + U - 0.5 * q, x1 + U + 0.5 * q
#     end
#     set_positions!(at, X)
#     return at
# end
#
# function ml_displacement!(at, u)
#     I0, I1, Idel = si_multilattice(at)
#     @assert isempty(Idel)  # if Idel is not empty then (for now) we don't know what to do
#     X = positions(at)
#     W = CauchyBorn.WcbQuad(calculator(at))
#
#     # transformation matrices
#     Tp = [0 1/√2  -1/√2; 1 0 0; 0 1/√2 1/√2]
#     Tm = diagm([1,1,-1]) * Tp
#
#     F0 = cell(W.at)'  # get reference information
#     p0 = W(F0)
#
#     for (i0, i1) in zip(I0, I1)   # each pair (i0, i1) corresponds to a ML lattice site
#         x0, x1 = X[i0], X[i1]
#         x1[3] > x0[3] ? T = Tp : T = Tm
#         U, ∇U = u(x0)            # displacement and displacement gradient
#         F = T' * (I + ∇U) * T
#         q = T * (W(F * F0) - p0)   # construct the shift corresponding to F = Id + ∇U
#         X[i0], X[i1] = x0 + U, x1 + U + q
#     end
#     set_positions!(at, X)
#     return at
# end

"""
`edge110` : produces a high-quality CLE solution for an edge dislocation
in bulk silicon, 110 orientation.

## Keyword Arguments

* a : lattice parameter, default value is the default lattice parameter of the species, other allows values are `:equilibrate`
"""
function edge110(sym::Symbol, R::Real;
                  truncate=true, cle=:anisotropic, ν=0.25,
                  calc=sw_eq(), symmetrize = nothing,  # switch back to sym = true
                  TOL=1e-4, zDir=1,
                  eos_correction = true,
                  a = cell(bulk(sym, cubic=true))[1,1])

   @assert sym == :Si

   # compute the lattice parameter
   if a == :equilibrate
       atu = bulk(sym, cubic=true, pbc=true)
       set_calculator!(atu, calc)
       variablecell!(atu)
       minimise!(atu)
       a = cell(atu)[1,1]
   end

   # setup undeformed geometry
   at, b, x0 = si110_cluster(sym, ceil(Int, R/rnn(:Si)+3); a = a)
   set_calculator!(at, calc)

   # if cle == :anisotropic
   #    W = CauchyBorn.WcbQuad(calc)
   #    C = elastic_moduli(W)
   #    Cv = round.(voigt_moduli(C), 8)
   #    U = CLE.EdgeCubic(b, Cv, a, x0 = x0)
   # elseif cle == :isotropic
   #    @assert sym == nothing
   #    U = xx -> ([CLE.u_edge_isotropic(xx[1]-x0[1], xx[2]-x0[2], b, ν)...; 0.0], nothing)
   # else
   #    error("unknown `cle` option")
   # end

  @assert symmetrize == nothing
  U = xx -> ([CLE.u_edge_isotropic(xx[1]-x0[1], xx[2]-x0[2], b, ν)...; 0.0], nothing)

   if symmetrize == true
      @assert cle == :anisotropic
      symml_displacement!(at, U)
   elseif symmetrize == false
      @assert cle == :anisotropic
      ml_displacement!(at, U)
   elseif symmetrize == nothing
      X = positions(at)
      for (i,x) in enumerate(X)
         X[i] += U(x)[1]
      end
      set_positions!(at, X)
   else
      error("edge110: unknown parameters symmetrize = $(symmetrize)")
   end

   if truncate
      r = [ norm(x[1:2] - x0[1:2]) for x in positions(at) ]
      Ikeep = findall( r .<= R )
      X = positions(at)[Ikeep]
      C = cell(at)
      at = Atoms(:Si, X)
      set_cell!(at, C)
   end

   return at, x0
end

end


text:
# # Color palette swapping
# We can hide an image in another one by switching color palettes.
# Let's try to hide `img_secret` in `img`:
using Images
using DitherPunk
using Clustering
using IndirectArrays
using TestImages

img = testimage("peppers")
img_secret = testimage("airplaneF16");

# First we dither to eight colors:
ncolors = 8
d = dither(img, FloydSteinberg(), ncolors)
#
d_secret = dither(img_secret, FloydSteinberg(), ncolors)

# Both of these images are IndirectArrays, which means they contain a color scheme
d.values
# and a matrix of indices of type `UInt8` pointing to the color in the color scheme:
d.index

# Since `typemax(UInt8)` is 255 and therefore larger than `ncolor^2`,
# we can fit both images into a single index matrix by modifying the color schemes.
# To recover the secret image, it is then only necessary to swap the new color schemes.
combineindex(a, b) = UInt8((a - 1) * ncolors + b)
function uncombineindex(c)
    b = (c - 1) % ncolors + 1
    a = Int((c - b) / ncolors) + 1
    return a, b
end

function combine_images(ia1::IndirectArray, ia2::IndirectArray)
    index = map(t -> combineindex(t...), zip(ia1.index, ia2.index))
    cs1 = similar(ia1.values, ncolors^2)
    cs2 = similar(ia2.values, ncolors^2)

    for i in 1:(ncolors^2)
        a, b = uncombineindex(i)
        cs1[i] = ia1.values[a]
        cs2[i] = ia2.values[b]
    end
    return index, cs1, cs2
end;

# Let's test this on our dithered images
index, cs, cs_secret = combine_images(d, d_secret);

# When using the color scheme `cs`, we recover the peppers:
IndirectArray(index, cs)
# However, when using the secret color scheme `cs_secret`, the F-16 airplane appears:
IndirectArray(index, cs_secret)

# This certainly isn't a very secure way to hide an image in another one.
# However, it can be useful when trying to fit multiple images on low memory microcontrollers.
# As a matter of fact, this example was inspired by
# [Mark Ferrari's GDC 2016 talk](https://youtu.be/aMcJ1Jvtef0?t=717)
# where he used a similar technique for this very purpose.


text:
using Flux

model = Chain(
  Dense(10, 5, σ),
  Dense(5, 2),
  softmax)

@show(model(rand(10)))


text:
using ComplexValues
const Spherical = ComplexValues.Spherical
const Polar = ComplexValues.Polar

using Test

@testset "Constructions of Polar" begin
	@test Polar(-5+1im) isa Number
	@test Polar(1.0im) isa Number
	@test Polar(-2) isa Number
	@test Polar(10.0f0) isa Number
	@test Polar(2.0f0,pi) isa Number
	@test Polar(pi,pi) isa Number
end

@testset "Constructions of Spherical" begin
	@test Spherical(-5+1im) isa Number
	@test Spherical(1.0im) isa Number
	@test Spherical(-2) isa Number
	@test Spherical(10.0f0) isa Number
	@test Spherical(2.0f0,pi) isa Number
	@test Spherical(pi/2,pi) isa Number
	@test S2coord(Spherical(1)) ≈ [1,0,0]
end

@testset "Conversions between Polar, Spherical" begin
	for z in [-5+1im,1.0im,-2,10f0,Polar(Inf,-pi/5),Spherical(1,.5)]
		@test Polar(Spherical(z)) ≈ z
	end
	for z in [-5+1im,1.0im,-2,10f0]
		@test Polar{Float64}(Spherical{Float64}(z)) ≈ z
	end
end

@testset "Conversions in and out of Complex" begin
	for z in [-5+1im,1.0im,-2,10f0,Polar(Inf,-pi/5),Spherical(1,.5)]
		@test Polar(Complex(z)) ≈ Polar(z)
		@test Complex(Polar(z)) ≈ Complex(z)
		@test Spherical(Complex(z)) ≈ Spherical(z)
		@test Complex(Spherical(z)) ≈ Complex(z)
	end
end

@testset "Zero and infinity for $T" for T in [Polar,Spherical]
    z = T(Inf)
    @test isinf(z)
    @test isinf(abs(z))
    @test iszero(1/z)
    @test iszero(abs(inv(z)))
    @test iszero(Complex(4/z))
    z = T(0)
    @test iszero(z)
    @test isinf(3im/z)
    z = zero(T)
    @test iszero(z)
    @test isinf(3im/z)
end

@testset "Binary functions on $S,$T" for S in [Polar,Spherical], T in [Polar,Spherical]
	u = S(3.0+5.0im)
	v = T(-5+1im)
	@test Complex(u+v) ≈ Complex(u)+Complex(v)
	@test Complex(u*v) ≈ Complex(u)*Complex(v)
	@test Complex(u/v) ≈ Complex(u)/Complex(v)
	@test Complex(u-v) ≈ Complex(u)-Complex(v)
	@test Complex(u^v) ≈ Complex(u)^Complex(v)
end

@testset "Unary functions on $T" for T in [Polar,Spherical]
    u = T(3.0+5.0im)
    @test Complex(u-4) ≈ Complex(u)-4
    @test Complex(4-u) ≈ 4-Complex(u)
    @test Complex(cos(exp(u))) ≈ cos(exp(Complex(u)))
    @test Complex(sqrt(inv(u))) ≈ sqrt(inv(Complex(u)))
end


text:
"""
    sim_values(::Type{T}) where T<:SpatQSimSpec
    sim_values(::T) where T<:SpatQSimSpec

Return the values of the parameters that change for each simulation study.
"""
sim_values(::T) where T<:SpatQSimSpec = sim_values(T)
sim_values(::Type{<:QDevScalingSpec}) = 10.0 .^ (-3:0.5:-0.5)
sim_values(::Type{<:SharedQSpec}) = 0:0.2:1
sim_values(::Type{<:PrefIntensitySpec}) = [0, 1, 2, 4, 8, 16]
sim_values(::Type{<:DensityDependentQSpec}) = 0:0.25:1.25
sim_values(::Type{<:HabQSpec}) = 1.75 .^ (-2:3)

function sim_value_idx(spec::SpatQSimSpec)
    findfirst(isapprox(sim_value(spec)),
              sim_values(spec))
end


text:
function make_pullback(ℓ, autodiff::AD.AbstractBackend)
    function ∇ℓ(x)
        res = AD.value_and_gradient(autodiff, ℓ, x)
        return (first(res), first(last(res)))
    end
end

# This is custom since Zygote and AD have some weird performanc regression currently
function make_pullback(ℓ, ::AD.ReverseRuleConfigBackend)
    function ∇ℓ(x)
        f, b = AD.Zygote.pullback(ℓ, x)
        return (f, first(b(1.0)))
    end
end

function make_pullback(ℓ, autodiff::Function)
    function ∇ℓ(x)
        return (ℓ(x), autodiff(x))
    end
end


text:
using Plots

include("../src/conical_radon.jl")
include("../src/inv_conical_radon.jl")
include("../src/util.jl")
include("../src/phantom.jl")

function innerprod(f,g)
    a = 1
    return sum(f .* conj(g)) / ((2a)^2)
end

ϕ = π/4
μ = 1

# target = rand(100,100)^0
target = phantom(100,2)

Cμf = conical_radon(target, ϕ, μ)
plot(heatmap(target),heatmap(reverse(Cμf, dims=1)), size = (1000,400))

begin
X,Z = size(Cμf)
n = length(size(Cμf))
f = zeros(ComplexF64,X,Z)

# @time for κ in [100], m in [100]
@time for t in 1:10000
    κ = .1randn()
    m = .1randn()

    Φ_κm = zeros(ComplexF64,X,Z)

    for x in 1:X, z in 1:Z
        Φ_κm[x,z] = Φ(κ,m, x,z)
    end
    # cplx_heatmap(Φ_κm)
    f += (Φ_κm .* (innerprod(Cμf,Φ_κm)/λ(κ,m, μ,ϕ)))
end
cplx_heatmap(f)
end

# Φ_κm = zeros(ComplexF64,I,J)
# for x in 1:I, z in 1:J
#     Φ_κm[x,z] = Φ(10.3,25.5, x/200-1,z/200-1)
# end
# cplx_heatmap(Φ_κm)


text:
"""
InitialConditions(;
                    initial_thickness = fill!(Array{Float64}(undef,1,1),NaN),
                    initial_viscosity = fill!(Array{Float64}(undef,1,1),NaN),
                    initial_temperature = fill!(Array{Float64}(undef,1,1),NaN),
                    initial_damage = fill!(Array{Float64}(undef,1,1),NaN))

Construct a WAVI.jl initial conditions object.
Unpassed arguments default to 1x1 nan matrix; unspecified initial conditions are overwritten by default values specified in Params structure when model is assembled.

Keyword arguments
=================

- 'initial_thickness': (nx x ny) matrix defining ice thickness at t = 0
- 'initial_viscosity': (nx x ny x nz) matrix defining viscosity on sigma levels at t = 0
- 'initial_temperature': (nx x ny x nz) matrix defining temperature on sigma levels at t = 0
- 'initial_damage': (nx x ny x nz) matrix defining ice damage at t = 0
"""
@with_kw struct InitialConditions{T <: Real}
    initial_thickness::Array{T,2} = fill!(Array{Float64}(undef,1,1),NaN)
    initial_viscosity::Array{T,3} = fill!(Array{Float64}(undef,1,1,1),NaN)
    initial_temperature::Array{T,3} = fill!(Array{Float64}(undef,1,1,1),NaN)
    initial_damage::Array{T,3} = fill!(Array{Float64}(undef,1,1,1),NaN)
end


text:
using MonteCarloMeasurements

# Expected value of x log-weighted by ℓ
function expect(x,ℓ)
    w = exp(ℓ - maximum(ℓ))
    mean(x*w) / mean(w)
end

x = Particles(5000,TDist(3))
ℓ = logpdf(Normal(),x) - logpdf(TDist(3), x)

expect(x,ℓ)
expect(x^2,ℓ)
expect(x^3,ℓ)
expect(x^4,ℓ)

# Now sample z ~ Normal()

z = Particles(5000,Normal())

expect(z,ℓ)
expect(z^2,ℓ)
expect(z^3,ℓ)
expect(z^4,ℓ)

# And now a new weighted variable

y = Particles(5000,TDist(3))
ℓ += logpdf(Normal(),y) - logpdf(TDist(3), y)

expect(y,ℓ)
expect(y^2,ℓ)
expect(y^3,ℓ)
expect(y^4,ℓ)

expect(x+y,ℓ)
expect((x+y)^2,ℓ)

expect(x+z,ℓ)
expect((x+z)^2,ℓ)


text:
module Quadratures

import GaussQuadrature
import StaticArrays: SVector, SMatrix, MMatrix
import LinearAlgebra: Diagonal

export QuadratureStyle,
    GLL, GL, polynomial_degree, degrees_of_freedom, quadrature_points

"""
   QuadratureStyle

Quadrature style supertype. See sub-types:
 - [`GLL`](@ref)
 - [`GL`](@ref)
 - [`Uniform`](@ref)
"""
abstract type QuadratureStyle end

"""
    polynomial_degree(QuadratureStyle) -> Int

Returns the polynomial degree of the `QuadratureStyle` concrete type
"""
function polynomial_degree end

"""
    degrees_of_freedom(QuadratureStyle) -> Int

Returns the degrees_of_freedom of the `QuadratureStyle` concrete type
"""
function degrees_of_freedom end

"""
    points, weights = quadrature_points(::Type{FT}, quadrature_style)

The points and weights of the quadrature rule in floating point type `FT`.
"""
function quadrature_points end

"""
    GLL{Nq}()

Gauss-Legendre-Lobatto quadrature using `Nq` quadrature points.
"""
struct GLL{Nq} <: QuadratureStyle end

polynomial_degree(::GLL{Nq}) where {Nq} = Nq - 1
degrees_of_freedom(::GLL{Nq}) where {Nq} = Nq

@generated function quadrature_points(::Type{FT}, ::GLL{Nq}) where {FT, Nq}
    points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.both)
    :($(SVector{Nq}(points)), $(SVector{Nq}(weights)))
end

"""
    GL{Nq}()

Gauss-Legendre quadrature using `Nq` quadrature points.
"""
struct GL{Nq} <: QuadratureStyle end

polynomial_degree(::GL{Nq}) where {Nq} = Nq - 1
degrees_of_freedom(::GL{Nq}) where {Nq} = Nq

@generated function quadrature_points(::Type{FT}, ::GL{Nq}) where {FT, Nq}
    points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.neither)
    :($(SVector{Nq}(points)), $(SVector{Nq}(weights)))
end

"""
    Uniform{Nq}()

Uniformly-spaced quadrature.
"""
struct Uniform{Nq} <: QuadratureStyle end

@generated function quadrature_points(::Type{FT}, ::Uniform{Nq}) where {FT, Nq}
    points = SVector{Nq}(range(-1 + 1 / Nq, step = 2 / Nq, length = Nq))
    weights = SVector{Nq}(ntuple(i -> 2 / Nq, Nq))
    :($points, $weights)
end

"""
    barycentric_weights(x::AbstractVector)

The barycentric weights associated with the array of point locations `x`:

```math
w[i] = \\frac{1}{\\prod_{j \\ne i} (x[i] - x[j])
```
Reference:
  - [Berrut2004](@cite) equation 3.2
"""
function barycentric_weights(r::SVector{Nq, T}) where {Nq, T}
    SVector{Nq}(ntuple(Nq) do i
        w = one(T)
        for j in 1:Nq
            if j != i
                w *= (r[j] - r[i])
            end
        end
        inv(w)
    end)
end
@generated function barycentric_weights(
    ::Type{FT},
    quadstyle::QuadratureStyle,
) where {FT}
    barycentric_weights(quadrature_points(FT, quadstyle())[1])
end

function interpolation_matrix(
    points_to::SVector{Nto},
    points_from::SVector{Nfrom},
) where {Nto, Nfrom}
    T = eltype(points_to)
    bw = barycentric_weights(points_from)
    M = zeros(MMatrix{Nto, Nfrom, T, Nto * Nfrom})
    for i in 1:Nto
        x_to = points_to[i]
        skip_row = false
        for j in 1:Nfrom
            if x_to == points_from[j]
                # assign to one to avoid singularity condition
                M[i, j] = one(T)
                # skip over the equal boundry condition
                skip_row = true
            end
            skip_row && break
        end
        skip_row && continue
        w = bw ./ (x_to .- points_from)
        M[i, :] .= w ./ sum(w)
    end
    return SMatrix(M)
end

@generated function interpolation_matrix(
    ::Type{FT},
    quadto::QuadratureStyle,
    quadfrom::QuadratureStyle,
) where {FT}
    interpolation_matrix(
        quadrature_points(FT, quadto())[1],
        quadrature_points(FT, quadfrom())[1],
    )
end

"""
    V = orthonormal_poly(points, quad)

`V[i,j]` contains the `j-1`th Legendre polynomial evaluated at `points[i]`.
i.e. it is the mapping from the modal to the nodal representation.
"""
function orthonormal_poly(
    points::SVector{Np, FT},
    quad::GLL{Nq},
) where {FT, Np, Nq}
    N = Nq - 1
    a, b = GaussQuadrature.legendre_coefs(FT, N)
    if N == 0
        return SMatrix{Np, 1}(ntuple(x -> b[1], Np))
    end
    return SMatrix{Np, Nq}(GaussQuadrature.orthonormal_poly(points, a, b))
end

function spectral_filter_matrix(
    quad::GLL{Nq},
    Σ::SVector{Nq, FT},
) where {Nq, FT}
    points, _ = quadrature_points(FT, quad)
    V = orthonormal_poly(points, quad)
    return V * Diagonal(Σ) / V
end

function cutoff_filter_matrix(
    ::Type{FT},
    quad::GLL{Nq},
    Nc::Integer,
) where {FT, Nq}
    Σ = SVector(ntuple(i -> i <= Nc ? FT(1) : FT(0), Nq))
    return spectral_filter_matrix(quad, Σ)
end

"""
    differentiation_matrix(r::SVector{Nq, T},
                           wb=barycentric_weights(r)::SVector{Nq,T}) where {Nq, T}

The spectral differentiation matrix for a polynomial of degree `Nq` - 1 defined on the
points `r` with associated barycentric weights `wb`.

```math
D_{i,j} = \\begin{cases}
    \\frac{w_j / w_i}{x_i - x_j} \\text{ if } i \\ne j \\\\
    -\\sum_{k \\ne j} D_{k,j} \\text{ if } i = j
\\end{cases}
```

Reference:
 - [Berrut2004](@cite),#https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
"""
function differentiation_matrix(
    r::SVector{Nq, T},
    wb = barycentric_weights(r)::SVector{Nq, T},
) where {Nq, T}
    SMatrix{Nq, Nq, T, Nq * Nq}(
        begin
            if i == j
                D = zero(T)
                for l in 1:Nq
                    if l != i
                        D += one(T) / (r[i] - r[l])
                    end
                end
                D
            else
                (wb[i] / wb[j]) / (r[j] - r[i])
            end
        end for j in 1:Nq, i in 1:Nq
    )
end
@generated function differentiation_matrix(
    ::Type{FT},
    quadstyle::QuadratureStyle,
) where {FT}
    differentiation_matrix(quadrature_points(FT, quadstyle())[1])
end

end # module


text:
import Pkg
Pkg.activate(".")

using Revise
using SimpleSDMLayers
using Plots
using GBIF
using DataFrames

# Get world temperature data
temperature = worldclim(1)

kingfisher = GBIF.taxon("Megaceryle alcyon", strict=true)
kf_occurrences = occurrences(kingfisher)
# Get at least 200 occurrences
while length(kf_occurrences) < 200
    occurrences!(kf_occurrences)
    @info "$(length(kf_occurrences)) occurrences"
end
kf_occurrences

# Clip layer to occurrences
temperature_clip = clip(temperature, kf_occurrences)

presabs1 = mask(temperature_clip, kf_occurrences, Bool)

# Plot occurrences
contour(temperature_clip, fill = true)
scatter!(longitudes(kf_occurrences), latitudes(kf_occurrences))

## Test new mask feature ####

presabs = mask(temperature_clip, kf_occurrences, Bool)
abund = mask(temperature_clip, kf_occurrences, Float64)

presabs.grid
abund.grid

plot(presabs)
plot(abund, c = :BuPu)

presabs2 = convert(Float64, presabs)
presabs.grid
presabs2.grid
abund.grid

tmp = copy(abund)
replace!(tmp.grid, 0 => nothing)
tmp.grid |> unique
abund.grid |> unique
plot(convert(Float64, abund), c = :lightgrey)
plot!(convert(Float64, tmp))

tmp = copy(presabs)
replace!(tmp.grid, 0 => nothing)
tmp.grid |> unique
presabs.grid |> unique
plot(convert(Float64, presabs), c = :lightgrey)
plot!(convert(Float64, tmp))

tmp = copy(presabs)
tmp.grid
tmp = convert(Int, tmp)
tmp.grid
replace!(tmp.grid, 0 => nothing)
tmp.grid |> unique
presabs.grid |> unique
plot(convert(Float64, presabs), c = :lightgrey)
plot!(convert(Float64, tmp))

# Test removezeros argument
presabs1 = mask(temperature_clip, kf_occurrences, Bool)
presabs2 = mask(temperature_clip, kf_occurrences, Bool, removezeros = true)
presabs1.grid |> unique
presabs2.grid |> unique

## Others

grid = [true false false]
grid = [true false NaN]
grid = [true false nothing]
grid = [true false missing]
grid = [1 0 nothing]
grid = [1 0 NaN]

heatmap(grid)

layer = temperature_clip
layer = presabs
lg = copy(layer.grid)

eltype(layer) <: Number

lg = eltype(lg) <: Number ? lg : convert(Matrix{Union{eltype(lg), Float64}}, lg)
replace!(lg, nothing => NaN)
lg = convert(Matrix{Float64}, lg)
longitudes(layer), latitudes(layer), lg

# Sliding window
buffered = slidingwindow(presabs, maximum, 50.0)
plot(buffered)

## Test new replace/replace!

presabs
presabs.grid

tmp = copy(presabs)
replace!(tmp, 0 => nothing)
tmp.grid
tmp.grid |> unique
plot(convert(Float64, tmp), c = :viridis)

## Test DataFrames mask

kf_df = DataFrame(kf_occurrences)
select!(kf_df, [:key, :latitude, :longitude])

# Basic mask

layer = similar(temperature_clip, Bool)
layer = similar(temperature_clip, Float32)

presabs = mask(temperature_clip, kf_df, Bool)
abund   = mask(temperature_clip, kf_df, Float32)

plot(convert(Float32, presabs), c = :viridis)
plot(abund, c = :viridis)

# Latitude/longitude colnames

latitude = :latitude
longitude = :longitude

kf_df2 = rename(kf_df, :longitude => :lon, :latitude => :lat)
mask(temperature_clip, kf_df2, Bool) # should fail
mask(temperature_clip, kf_df2, Float32) # should fail

mask(temperature_clip, kf_df2, Bool; latitude = :lat, longitude = :lon) # should work
mask(temperature_clip, kf_df2, Float32; latitude = :lat, longitude = :lon) # should work

mask(temperature_clip, kf_df, Bool) # should still work
mask(temperature_clip, kf_df, Float32) # should still work

# Allocations

kf_df
kf_df2 = copy(kf_df)

kf_df2.key[1] = 1
kf_df2
kf_df2[!, :key] .= 1
kf_df2
kf_df

kf_df2 = copy(kf_df)
tmp = mask(temperature_clip, kf_df2, Float32)
tmp.grid
isequal(kf_df, kf_df2)

lat = kf_df[1, :latitude]
lon = kf_df[1, :longitude]
temperature_clip[lon, lat]

tmp = SimpleSDMResponse(kf_df, :key, temperature_clip)
tmp.grid

tempdf = DataFrame(temperature_clip)
templayer = SimpleSDMPredictor(tempdf, :values, temperature_clip)

isequal(templayer, temperature_clip)
isequal(templayer.grid, temperature_clip.grid)

df = tempdf
lats = df[!, latitude]
lons = df[!, longitude]
uniquelats = unique(lats)
uniquelons = unique(lons)
grid = Array{Any}(nothing, size(layer))
lats_idx = [SimpleSDMLayers._match_latitude(layer, lat) for lat in lats]
lons_idx = [SimpleSDMLayers._match_longitude(layer, lon) for lon in lons]
DataFrame(lats = lats_idx, lons = lons_idx)
unique(DataFrame, [:lats, :lons])

## Filter unique sites
tempdf = DataFrame(temperature_clip)
filter!(x -> !isnothing(x.values), tempdf)

# Benchmarks
smallmultidf = vcat(copy(kf_df), copy(kf_df))

mediummultidf = vcat(copy(kf_df), copy(kf_df), copy(kf_df), copy(kf_df),
                     copy(kf_df), copy(kf_df), copy(kf_df), copy(kf_df))

midmediummultidf = vcat(copy(mediummultidf), copy(mediummultidf), copy(mediummultidf),
                        copy(mediummultidf), copy(mediummultidf), copy(mediummultidf))

bigmultidf = vcat(copy(tempdf), copy(tempdf))

hugemultidf = vcat(copy(bigmultidf), copy(bigmultidf), copy(bigmultidf), copy(bigmultidf),
                   copy(bigmultidf), copy(bigmultidf), copy(bigmultidf), copy(bigmultidf))

# 400 obs (2 x 200)
@time mask(temperature_clip, smallmultidf, Bool) # without unique filtering
# 0.029148 seconds (35.14 k allocations: 5.966 MiB)
# 0.008258 seconds (22.03 k allocations: 5.249 MiB)
@time mask(temperature_clip, smallmultidf, Bool) # with unique filtering
# 0.120525 seconds (75.84 k allocations: 6.536 MiB)
# 0.003245 seconds (12.56 k allocations: 3.843 MiB)

# 1600 obs (8 x 200)
@time mask(temperature_clip, mediummultidf, Bool) # without unique filtering
# 0.038381 seconds (101.17 k allocations: 14.591 MiB)
# 0.030026 seconds (88.04 k allocations: 13.873 MiB, 49.52% gc time)
@time mask(temperature_clip, mediummultidf, Bool) # with unique filtering
# 0.112739 seconds (79.48 k allocations: 6.614 MiB)
# 0.003388 seconds (16.16 k allocations: 3.921 MiB)

# 10 000 obs
@time mask(temperature_clip, midmediummultidf, Bool) # with unique filtering
# 0.118171 seconds (541.20 k allocations: 72.088 MiB, 12.66% gc time)
# 0.092093 seconds (528.03 k allocations: 71.368 MiB, 15.12% gc time)
@time mask(temperature_clip, midmediummultidf, Bool) # with unique filtering
# 0.126579 seconds (103.51 k allocations: 7.160 MiB)
# 0.004906 seconds (40.16 k allocations: 4.466 MiB)

# 100 000 obs
@time mask(temperature_clip, bigmultidf, Bool)
# 0.928445 seconds (7.37 M allocations: 993.580 MiB, 11.33% gc time)
# 0.914822 seconds (7.37 M allocations: 993.580 MiB, 12.57% gc time)
@time mask(temperature_clip, bigmultidf, Bool)
# 0.591735 seconds (3.75 M allocations: 505.812 MiB, 9.40% gc time)
# 0.488477 seconds (3.69 M allocations: 503.119 MiB, 10.87% gc time)

# 1 000 000 obs
@time mask(temperature_clip, hugemultidf, Bool)
# 7.427266 seconds (58.99 M allocations: 7.747 GiB, 11.92% gc time)
# 6.763797 seconds (58.98 M allocations: 7.746 GiB, 10.16% gc time)
@time mask(temperature_clip, hugemultidf, Bool)
# 0.610696 seconds (3.75 M allocations: 528.286 MiB, 9.03% gc time)
# 0.464374 seconds (3.69 M allocations: 525.593 MiB, 8.72% gc time)

# Make sure result is the same
res1 = mask(temperature_clip, smallmultidf, Bool)
res2 = mask(temperature_clip, smallmultidf, Bool)
isequal(res1, res2) # nope, but whatever
isequal(res1.grid, res2.grid) # yup, that's what's important


text:
using ggplot

x = seq(-5,5, by = 0.2)
y = sin(x)
z = cos(y)

table = data.frame(x = x, sin =y, cos = z)

print("preview of the input data:")
print(table, max.print = 13)

bitmap(plt, file = `${@dir}/line_3colors_julia.png`):

	plt = ggplot(table, padding = "padding: 225px 600px 300px 250px;")
	plt += geom_line( aes(x = "x", y = "sin"), width = 8)
	plt += geom_line( aes(x = "x", y = "cos"), width = 8)
	plt += geom_line( aes(x = "x", y = "~sin + cos"), width = 8)
    plt += ggtitle("line overlaps")
    plt += xlab("time(ticks)")
    plt += ylab("signal data")

    plot(plt, width = 2600, height = 1600)

end

max_stars_repo_name: rsharp-lang/ggplot

max_issues_repo_name: rsharp-lang/ggplot

max_forks_repo_name: rsharp-lang/ggplot


text:
#------------------------------------------------------------------------------

# computes an invariant subspaces of the given matrices
#
# returns a common invariant subspace represented by basis vectors
#
# Gleb: I guess we can say that we guarantee that if there is a proper nonempty over Q,
#       it will be found
function invariant_subspace_1(matrices::AbstractArray{T}) where {T<:AbstractSparseMatrix}
    isempty(matrices) && error("empty system invariants are ill-defined")

    # generate a basis for the Algebra
    algebra = find_basis(deepcopy(matrices))

    # find the radical of the Algebra
    @info "computing the radical.."
    @time radical = find_radical_sup(algebra)

    # find an invariant subspace
    if length(radical) != 0 && length(radical) < dim(algebra)
        # Gleb: why `length(radical) < dim(algebra)` ?
        # Alex: we want radical to be nontrivial here
        @info "radical is nontrivial, computing the general kernel"
        invariant = general_kernel(map(to_dense, radical))
        invariant = [
            from_dense([invariant[:, i]...], Nemo.QQ)
            for i in 1:size(invariant, 2)
        ]
    else
        @info "radical is trivial, using randomized algorithm"
        if dim(algebra) != 0
            invariant = invariant_subspace_semisimple(algebra)
        else
            invariant = fullspace(size(first(matrices), 1), base_ring(first(matrices)))
        end
    end

    if isempty(invariant)
        @warn "no invariant subspaces"
    else
        @info "$(length(invariant)) dimensional subspace found"
    end

    for v in invariant
        scale!(v, lcm(v))
    end

    invariant
end

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

# computes invariant subspaces using the given hint,
# , i.e, an array of vectors generating a subspace
function __invariant_subspace_2(matrices::AbstractArray{T}, hint) where {T<:AbstractSparseMatrix}
    isempty(matrices) && error("empty system invariants are ill-defined")

    V = linear_span!(hint)

    new_pivots = collect(keys(V.echelon_form))
    i = 0

    while !isempty(new_pivots)
        pivots_to_process = sort(deepcopy(new_pivots))
        empty!(new_pivots)

        for pivot in pivots_to_process
            for (vectidx, vect) in enumerate(matrices)

                product = vect * V.echelon_form[pivot]

                i += 1

                if !iszero(product)
                    new_pivot = eat_sparsik!(V, product)
                    if new_pivot != reduced
                        push!(new_pivots, new_pivot)
                    end
                end

            end
        end

    end

    return basis(V)
end

function invariant_subspace_2(matrices::AbstractArray{T}) where {T<:AbstractSparseMatrix}
    n = size(first(matrices), 1)
    ground = base_ring(first(matrices))

    is = rand(1:n, 2)

    vs = [ unit_sparsik(n, i, ground) for i in is ]

    @info "hint" vs
    return __invariant_subspace_2(matrices, vs)
end

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

function __many_invariant_subspaces(
        As::AbstractArray{T},
        find_invariant)   where {T<:AbstractSparseMatrix}

    n = size(first(As), 1)
    ground = base_ring(first(As))

    # search for a subspace
    V = find_invariant(As)

    # no subspaces found
    if length(V) == 0
        return []
    end

    toreturn = [ V ]

    @info "found $(length(V))-dim subspace in ambient $(size(first(As), 1))-dim"

    # restrict
    if length(V) > 1 && length(V) < n
        As_V = restrict(As, V)
        As_V_sparse = map(x-> from_dense(x, ground), map(Array, As_V))
        subspaces = __many_invariant_subspaces(As_V_sparse, find_invariant)
        append!(toreturn, map(vs -> lift(vs, V), subspaces))
    end

    # factorize
    if length(V) + 1 < n # < n
        As_V = factorize(As, V)
        if !isempty(As_V)
            As_V_sparse = map(x-> from_dense(x, ground), map(Array, As_V))
            subspaces = __many_invariant_subspaces(As_V_sparse, find_invariant)
            lifted = map(
                vs -> lift(vs, augment_subspace(linear_span!(V))),
                subspaces)
            lifted = map(
                vs -> append!(vs, V),
                lifted)
            append!(toreturn, lifted)
        end
    end

    for vs in toreturn
        if !check_invariance!(deepcopy(As), deepcopy(vs))
            @assert false
        end
    end

    return toreturn
end

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

# Because matrices arising here are small,
# we mostly use dense Nemo algebra

# restricts the given matrices `As` to the subspace
# generated by vectors from `vs`
# Returns new, restriced representations of the matrices
function restrict(As::AbstractArray, vs)
    # Gleb: TODO for the future: the restriction is implicitly computed by check_invariance, so
    # it could be done there (and in the sparse way!)
    n = size(first(As), 1)
    ground = base_ring(first(As))

    MSpace = AbstractAlgebra.Generic.MatrixSpace(ground, n, length(vs))

    # Asvs = [A*f1 A*f2 .. Afd]
    # for fi in vs
    Asvs = [
        MSpace(hcat([ to_dense(A * v) for v in vs ]...))
        for A in As
    ]

    # matrix = [v1 v2 .. vn]
    # for vi in vs
    matrix = MSpace(hcat([ to_dense(v) for v in vs ]...))

    # solve for each A from As
    #   A*fi = a1*f1 + a2*f2 + ... + ad*fd
    [
        AbstractAlgebra.Generic.solve(matrix, Avs)
        for Avs in Asvs
    ]
end

# factorized the given matrices `As` by the subspace
# generated by vectors from `vs`
# Returns new, factorized representations of the matrices
function factorize(As::AbstractArray, vs)
    n = size(first(As), 1)
    ground = base_ring(first(As))

    MSpace = AbstractAlgebra.Generic.MatrixSpace(ground, n, n - length(vs))
    MSpaceplus = AbstractAlgebra.Generic.MatrixSpace(ground, n, n)

    # complement ⊕ vs = Rn
    complement = augment_subspace(linear_span!(vs))

    # Asvs = [A*f1 A*f2 .. Afd]
    # for fi in complement
    Asvs = [
        MSpace(hcat([to_dense(A * v) for v in complement ]...))
        for A in As
    ]

    # matrix = [f1 f2 .. fd e1 .. en]
    # for fi in complement
    # for ei in vs
    matrix = MSpaceplus(hcat(
        map(to_dense, complement)...,
        map(to_dense, vs)...)
    )

    # solve for each A from As
    #   A*fi = a1*f1 + a2*f2 + ... + ad*fd + e
    #   where e ∈ vs
    ans = [
        AbstractAlgebra.Generic.solve(matrix, Avs)[1:end-length(vs), :]
        for Avs in Asvs
    ]
    return ans
end

# lifts vectors `vs` to the subspace generated by vectors `Vs`
# Returns new, lifted representations of the vectors
function lift(vs, Vs)
    n = dim(first(Vs))
    ground = base_ring(first(Vs))

    lifted = []

    for v in vs
        lifted_vector = zero_sparsik(n, ground)
        for (i, t) in v
            lifted_vector = lifted_vector + scale(Vs[i], t)
        end
        push!(lifted, lifted_vector)
    end

    lifted
end

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