text:
Base.@propagate_inbounds function Base.setindex!(
    A::AbstractArray{T,N},
    val,
    i::Vararg{AbstractIndex,N}
   ) where {T,N}
    setindex!(A, to_indices(A, i))
end

@propagate_inbounds function Base.setindex!(a::IndicesArray, X, inds...)
    return unsafe_setindex!(parent(a), X, to_indices(axes(a), inds))
end

unsafe_setindex!(a, X, inds::Tuple) = @inbounds(setindex!(a, X, inds...))

#=
# Base.@propagate_inbounds function Base.setindex!(A::IndicesArray) end

for f in (:getindex, :view, :dotview)
    @eval begin
        @propagate_inbounds function Base.$f(A::NamedDimsArray; named_inds...)
            inds = order_named_inds(A; named_inds...)
            return Base.$f(A, inds...)
        end

        @propagate_inbounds function Base.$f(a::NamedDimsArray, inds::Vararg{<:Integer})
            # Easy scalar case, will just return the element
            return Base.$f(parent(a), inds...)
        end

        @propagate_inbounds function Base.$f(a::NamedDimsArray, ci::CartesianIndex)
            # Easy scalar case, will just return the element
            return Base.$f(parent(a), ci)
        end

        @propagate_inbounds function Base.$f(a::NamedDimsArray, inds...)
            # Some nonscalar case, will return an array, so need to give that names.
            data = Base.$f(parent(a), inds...)
            L = remaining_dimnames_from_indexing(dimnames(a), inds)
            return NamedDimsArray{L}(data)
        end
    end
end

############################################
# setindex!
@propagate_inbounds function Base.setindex!(a::NamedDimsArray, value; named_inds...)
    inds = order_named_inds(a; named_inds...)
    return setindex!(a, value, inds...)
end

@propagate_inbounds function Base.setindex!(a::NamedDimsArray, value, inds...)
    return setindex!(parent(a), value, inds...)
end
=#


text:
# Autogenerated wrapper script for eccodes_jll for armv7l-linux-gnueabihf
export eccodes

using libpng_jll
using OpenJpeg_jll
JLLWrappers.@generate_wrapper_header("eccodes")
JLLWrappers.@declare_library_product(eccodes, "libeccodes.so")
function __init__()
    JLLWrappers.@generate_init_header(libpng_jll, OpenJpeg_jll)
    JLLWrappers.@init_library_product(
        eccodes,
        "lib/libeccodes.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using BenchmarkTools
using RiemannTheta
using Distances

######### some testing values

Ω1 = [ 1.690983006 + 0.9510565162im 1.5+0.363271264im ;
      1.5+0.363271264im 1.309016994+0.9510565162im ]
T1 = Matrix(cholesky(imag.(Ω1)))

Ω2 = -1/(2π * im) * [ 111.207 96.616 ; 96.616 83.943 ]
T2 = Matrix(cholesky(imag.(Ω2)))

srand(0)
tmp = 5*rand(10,10) - 2.5
Ω3 = 5*rand(10,10) - 2.5 + (tmp * tmp') * im
T3 = Matrix(cholesky(imag.(Ω3)))

derivs1 = [ rand(ComplexF64, 2) for i in 1:1 ]
derivs2 = [ rand(ComplexF64, 10) for i in 1:5 ]

ϵ1, ϵ2 = 1e-3, 1e-8

############## lll_reduce.jl  ###############

############## radius.jl  ###############

@btime radius($ϵ1, $T1) # 3.2 μs
@btime radius($ϵ1, $T1, $derivs1) # 95 μs
@btime radius($ϵ2, $T3) # 540 μs
@btime radius($ϵ2, $T3, $derivs2) # 770 μs

############## innerpoints.jl  ###############

R = radius(ϵ1, T1)
@btime innerpoints(T1, R) # 4.6 μs

R = radius(ϵ2, T3)
@btime innerpoints(T3, R) # 770 μs

##########   riemanntheta function     ###########

zs1 = [ComplexF64[0.5, 0.]]
@btime riemanntheta(zs1, Ω1, eps=ϵ1) # 50.512 μs (269 allocations: 24.03 KiB)
@btime riemanntheta(zs1, Ω1, eps=ϵ2) # 58.536 μs (390 allocations: 35.91 KiB)

R = RiemannTheta.radius(ϵ1, Matrix(cholesky(imag.(Ω1))), Vector{ComplexF64}[], 5.)
RiemannTheta.innerpoints(Matrix(cholesky(imag.(Ω1))), R)
R = RiemannTheta.radius(ϵ2, Matrix(cholesky(imag.(Ω1))), Vector{ComplexF64}[], 5.)
RiemannTheta.innerpoints(Matrix(cholesky(imag.(Ω1))), R)

zs2 = [ ComplexF64[x, 2x] for x in -1:0.01:1 ]
@btime riemanntheta(zs2, Ω1, eps=ϵ1) # 608.379 μs (1097 allocations: 487.64 KiB)
@btime riemanntheta(zs2, Ω1, eps=ϵ2) # 930.824 μs (1122 allocations: 793.72 KiB)

@btime riemanntheta(zs2, Ω1, eps=ϵ1, derivs=derivs1) # 3.149 ms (61344 allocations: 3.45 MiB)
@btime riemanntheta(zs2, Ω1, eps=ϵ2, derivs=derivs1) # 4.760 ms (96737 allocations: 5.42 MiB)

srand(0)
zs3 = [ rand(ComplexF64, 10) for i in 1:20 ]
@btime riemanntheta(zs3, Ω3, eps=ϵ1) # 1.095 ms (8566 allocations: 852.06 KiB)
@btime riemanntheta(zs3, Ω3, eps=ϵ2) # 1.921 ms (12017 allocations: 1.08 MiB)

R = RiemannTheta.radius(ϵ1, Matrix(cholesky(imag.(Ω3))), Vector{ComplexF64}[], 5.)
RiemannTheta.innerpoints(Matrix(cholesky(imag.(Ω3))), R)
R = RiemannTheta.radius(ϵ2, Matrix(cholesky(imag.(Ω3))), Vector{ComplexF64}[], 5.)
RiemannTheta.innerpoints(Matrix(cholesky(imag.(Ω3))), R)

@btime riemanntheta(zs3, Ω3, eps=ϵ1, derivs=derivs2) # 97.773 ms (1004788 allocations: 112.40 MiB)
@btime riemanntheta(zs3, Ω3, eps=ϵ2, derivs=derivs2) # 200.701 ms (2147103 allocations: 240.10 MiB)

# comparable to openRT timings ?
# openRT at 90 ms
srand(0)
zs4 = [ ComplexF64[Complex(rand(),0.), Complex(rand(),0.)] for i in 1:10000 ]
Ω4 = [ 3im 0. ; 0. 2im ]
@btime riemanntheta(zs4, Ω4, eps=1e-8) #  22.784 ms (20724 allocations: 20.02 MiB)

Profile.clear()
@profile for i in 1:100 ; riemanntheta(zs3, Ω3, eps=ϵ1) ; end
Profile.print(format=:flat)

Count File                        Line Function
36 ...Theta\src\finite_sum.jl   124 finite_sum(::Array{Float64,2}, ::Ar...
 1 ...Theta\src\finite_sum.jl   141 finite_sum_small(::Array{Float64,2}...
 1 ...Theta\src\finite_sum.jl   142 finite_sum_small(::Array{Float64,2}...
25 ...Theta\src\finite_sum.jl   149 finite_sum_small(::Array{Float64,2}...
 2 ...Theta\src\finite_sum.jl   150 finite_sum_small(::Array{Float64,2}...
 2 ...Theta\src\finite_sum.jl   152 finite_sum_small(::Array{Float64,2}...
 2 ...Theta\src\finite_sum.jl   153 finite_sum_small(::Array{Float64,2}...
 2 ...Theta\src\finite_sum.jl   154 finite_sum_small(::Array{Float64,2}...
 1 ...Theta\src\finite_sum.jl   156 finite_sum_small(::Array{Float64,2}...

 1 ...heta\src\innerpoints.jl    25 (::RiemannTheta.#_innerpoints#8{Arr...
12 ...heta\src\innerpoints.jl    26 (::RiemannTheta.#_innerpoints#8{Arr...
 9 ...heta\src\innerpoints.jl    32 (::RiemannTheta.#_innerpoints#8{Arr...
102 ...heta\src\innerpoints.jl    34 (::RiemannTheta.#_innerpoints#8{Arr...
 8 ...heta\src\innerpoints.jl    19 innerpoints(::Array{Float64,2}, ::F...
22 ...heta\src\innerpoints.jl    41 innerpoints(::Array{Float64,2}, ::F...

 1 ...RiemannTheta\src\lll.jl    52 lll(::Array{Array{Float64,1},1})
18 ...RiemannTheta\src\lll.jl    53 lll(::Array{Array{Float64,1},1})
 2 ...RiemannTheta\src\lll.jl    54 lll(::Array{Array{Float64,1},1})
33 ...RiemannTheta\src\lll.jl    55 lll(::Array{Array{Float64,1},1})
 1 ...RiemannTheta\src\lll.jl    59 lll(::Array{Array{Float64,1},1})
 1 ...RiemannTheta\src\lll.jl    60 lll(::Array{Array{Float64,1},1})
 1 ...RiemannTheta\src\lll.jl    65 lll(::Array{Array{Float64,1},1})
 1 ...RiemannTheta\src\lll.jl    68 lll(::Array{Array{Float64,1},1})
 1 ...RiemannTheta\src\lll.jl    69 lll(::Array{Array{Float64,1},1})
 1 ...RiemannTheta\src\lll.jl    70 lll(::Array{Array{Float64,1},1})
 3 ...RiemannTheta\src\lll.jl    34 (::RiemannTheta.#☼#1{Array{Array{...
 8 ...RiemannTheta\src\lll.jl    36 (::RiemannTheta.#☼#1{Array{Array{...
36 ...RiemannTheta\src\lll.jl    37 (::RiemannTheta.#☼#1{Array{Array{...
 2 ...RiemannTheta\src\lll.jl    38 (::RiemannTheta.#☼#1{Array{Array{...

 1 ...iemannTheta\src\main.jl   113 #25
 1 ...iemannTheta\src\main.jl    76 #oscillatory_part#23(::Float64, ::A...
60 ...iemannTheta\src\main.jl    80 #oscillatory_part#23(::Float64, ::A...
30 ...iemannTheta\src\main.jl    81 #oscillatory_part#23(::Float64, ::A...
36 ...iemannTheta\src\main.jl    83 #oscillatory_part#23(::Float64, ::A...
 3 ...iemannTheta\src\main.jl   148 #riemanntheta#28(::Float64, ::Array...
229 ...iemannTheta\src\main.jl   149 #riemanntheta#28(::Float64, ::Array...
 4 ...iemannTheta\src\main.jl   152 #riemanntheta#28(::Float64, ::Array...
 2 ...iemannTheta\src\main.jl   110 exponential_part(::Array{Array{Comp...
 1 ...iemannTheta\src\main.jl   113 exponential_part(::Array{Array{Comp...

60 ...mannTheta\src\radius.jl    69 radius(::Float64, ::Array{Float64,2...

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

######### optim   ###########


text:
# Test constraint deletion
@testset "JuMP.delete (Constraints)" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, inf, Infinite(par))
    @variable(m, rv, SemiInfinite(inf, 0.5))
    @variable(m, pt, Point(inf, 0.5))
    @variable(m, x)
    data = TestData(par, 0, 1)
    meas = measure(inf + par - x, data)
    d1 = @deriv(inf, par)
    @constraint(m, c1, par - inf + pt + 2x - rv + meas - d1 <= 1)
    @constraint(m, c2, inf <= 9, DomainRestrictions(par => 0))
    @constraint(m, c3, [inf, x] in MOI.Zeros(2))
    # test normal deletion
    @test isa(delete(m, c1), Nothing)
    @test !is_valid(m, c1)
    @test !used_by_constraint(par)
    @test !(index(c1) in InfiniteOpt._constraint_dependencies(inf))
    @test !used_by_constraint(pt)
    @test !(index(c1) in InfiniteOpt._constraint_dependencies(x))
    @test !used_by_constraint(rv)
    @test !used_by_constraint(meas)
    @test !used_by_constraint(d1)
    @test length(InfiniteOpt._data_dictionary(c1)) == 2
    # test assertion error
    @test_throws AssertionError delete(m, c1)
    # test restricted deletion
    @test isa(delete(m, c2), Nothing)
    @test !is_valid(m, c2)
    @test !(index(c2) in InfiniteOpt._constraint_dependencies(inf))
    @test length(InfiniteOpt._data_dictionary(c2)) == 1
    @test isempty(m.constraint_restrictions)
    # test vector constraint deletion
    @test isa(delete(m, c3), Nothing)
    @test !is_valid(m, c3)
    @test !used_by_constraint(inf)
    @test !used_by_constraint(x)
    @test isempty(InfiniteOpt._data_dictionary(c3))
end

# Test deleting variable information
@testset "Variable Information" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, 0 <= inf1 <= 1, Infinite(par), Bin)
    @variable(m, inf2 == 1, Infinite(par), Int)
    @variable(m, pt1, Point(inf1, 0.5))
    @variable(m, pt2, Point(inf2, 0.5))
    @variable(m, 0 <= gb1 <= 1, Bin)
    @variable(m, gb2 == 1, Int)
    @variable(m, 0 <= d1 <= 1, Deriv(inf1, par))
    @variable(m, d2 == 1, Deriv(inf2, par))
    # test delete_lower_bound
    @testset "JuMP.delete_lower_bound" begin
        # test with infinite variable
        @test has_lower_bound(inf1)
        @test isa(delete_lower_bound(inf1), Nothing)
        @test !has_lower_bound(inf1)
        @test_throws ErrorException delete_lower_bound(inf2)
        # test with point variable
        @test has_lower_bound(pt1)
        @test isa(delete_lower_bound(pt1), Nothing)
        @test !has_lower_bound(pt1)
        @test_throws ErrorException delete_lower_bound(pt2)
        # test with finite variable
        @test has_lower_bound(gb1)
        @test isa(delete_lower_bound(gb1), Nothing)
        @test !has_lower_bound(gb1)
        @test_throws ErrorException delete_lower_bound(gb2)
        # test with derivative
        @test has_lower_bound(d1)
        @test isa(delete_lower_bound(d1), Nothing)
        @test !has_lower_bound(d1)
        @test_throws ErrorException delete_lower_bound(d2)
    end
    # test delete_upper_bound
    @testset "JuMP.delete_upper_bound" begin
        # test with infinite variable
        @test has_upper_bound(inf1)
        @test isa(delete_upper_bound(inf1), Nothing)
        @test !has_upper_bound(inf1)
        @test_throws ErrorException delete_upper_bound(inf2)
        # test with point variable
        @test has_upper_bound(pt1)
        @test isa(delete_upper_bound(pt1), Nothing)
        @test !has_upper_bound(pt1)
        @test_throws ErrorException delete_upper_bound(pt2)
        # test with finite variable
        @test has_upper_bound(gb1)
        @test isa(delete_upper_bound(gb1), Nothing)
        @test !has_upper_bound(gb1)
        @test_throws ErrorException delete_upper_bound(gb2)
        # test with derivative
        @test has_upper_bound(d1)
        @test isa(delete_upper_bound(d1), Nothing)
        @test !has_upper_bound(d1)
        @test_throws ErrorException delete_upper_bound(d2)
    end
    # test unfix
    @testset "JuMP.unfix" begin
        # test with infinite variable
        @test is_fixed(inf2)
        @test isa(unfix(inf2), Nothing)
        @test !is_fixed(inf2)
        @test_throws ErrorException unfix(inf1)
        # test with point variable
        @test is_fixed(pt2)
        @test isa(unfix(pt2), Nothing)
        @test !is_fixed(pt2)
        @test_throws ErrorException unfix(pt1)
        # test with finite variable
        @test is_fixed(gb2)
        @test isa(unfix(gb2), Nothing)
        @test !is_fixed(gb2)
        @test_throws ErrorException unfix(gb1)
        # test with derivative
        @test is_fixed(d2)
        @test isa(unfix(d2), Nothing)
        @test !is_fixed(d2)
        @test_throws ErrorException unfix(d1)
    end
    # test unset_binary
    @testset "JuMP.unset_binary" begin
        # test with infinite variable
        @test is_binary(inf1)
        @test isa(unset_binary(inf1), Nothing)
        @test !is_binary(inf1)
        @test_throws ErrorException unset_binary(inf2)
        # test with point variable
        @test is_binary(pt1)
        @test isa(unset_binary(pt1), Nothing)
        @test !is_binary(pt1)
        @test_throws ErrorException unset_binary(pt2)
        # test with finite variable
        @test is_binary(gb1)
        @test isa(unset_binary(gb1), Nothing)
        @test !is_binary(gb1)
        @test_throws ErrorException unset_binary(gb2)
    end
    # test unset_integer
    @testset "JuMP.unset_integer" begin
        # test with infinite variable
        @test is_integer(inf2)
        @test isa(unset_integer(inf2), Nothing)
        @test !is_integer(inf2)
        @test_throws ErrorException unset_integer(inf1)
        # test with point variable
        @test is_integer(pt2)
        @test isa(unset_integer(pt2), Nothing)
        @test !is_integer(pt2)
        @test_throws ErrorException unset_integer(pt1)
        # test with finite variable
        @test is_integer(gb2)
        @test isa(unset_integer(gb2), Nothing)
        @test !is_integer(gb2)
        @test_throws ErrorException unset_integer(gb1)
    end
end

# Test scalar parameter deletion
@testset "Scalar Parameters" begin
    # initialize model
    m = InfiniteModel()
    # setup the parameters
    @infinite_parameter(m, par in [0, 1])
    @infinite_parameter(m, par2 in [0, 1])
    @infinite_parameter(m, pars[1:2] in [0, 1])
    @infinite_parameter(m, par3 in [0, 1])
    @finite_parameter(m, fin == 42)
    @finite_parameter(m, fin2 == 42)
    dpar = dispatch_variable_ref(par)
    dpars = dispatch_variable_ref.(pars)
    dpar2 = dispatch_variable_ref(par2)
    dpar3 = dispatch_variable_ref(par3)
    dfin = dispatch_variable_ref(fin)
    dfin2 = dispatch_variable_ref(fin2)
    # setup the variables
    @variable(m, inf, Infinite(par))
    @variable(m, inf2, Infinite(par, par2))
    @variable(m, inf3, Infinite(par, pars))
    @variable(m, inf4, Infinite(par, pars))
    @variable(m, pt, Point(inf, 0.5))
    pt2 = @variable(m, variable_type = Point(inf2, 0.5, 0.5))
    @variable(m, pt3, Point(inf3, 0, [0, 0]))
    @variable(m, x)
    var = build_variable(error, inf4, Dict{Int, Float64}(1 => 0.5), check = false)
    rv = add_variable(m, var)
    dinf = dispatch_variable_ref(inf)
    dinf2 = dispatch_variable_ref(inf2)
    dinf3 = dispatch_variable_ref(inf3)
    dinf4 = dispatch_variable_ref(inf4)
    dpt = dispatch_variable_ref(pt)
    dpt2 = dispatch_variable_ref(pt2)
    dpt3 = dispatch_variable_ref(pt3)
    dx = dispatch_variable_ref(x)
    drv = dispatch_variable_ref(rv)
    # setup the measure
    data = TestData(par, 0, 1)
    mref = measure(inf + par - x + rv + par2 + par3 + fin, data)
    mref2 = measure(par2, data)
    dmref = dispatch_variable_ref(mref)
    dmref2 = dispatch_variable_ref(mref2)
    # setup a derivative
    d1 = @deriv(mref, par2)
    # setup the constraints
    @constraint(m, con, inf2 + inf4 - par2 + par3 + fin <= 0)
    constr = ScalarConstraint(par2, MOI.GreaterThan(0.))
    con2 = add_constraint(m, constr)
    @constraint(m, con3, [par2, par2] in MOI.Zeros(2))
    # setup objective
    set_objective(m, MOI.MIN_SENSE, fin)
    # test _check_param_in_data
    @testset "_check_param_in_data" begin
        @test_throws ErrorException InfiniteOpt._check_param_in_data(par, data)
        @test isa(InfiniteOpt._check_param_in_data(par2, data), Nothing)
        @test isa(InfiniteOpt._check_param_in_data(pars[1], data), Nothing)
        @test_throws ErrorException InfiniteOpt._check_param_in_data(par, BadData())
    end
    # test _update_measures
    @testset "_update_measures" begin
        @test InfiniteOpt._update_measures(m, par2) isa Nothing
        @test isequal_canonical(measure_function(dmref), inf + par - x + rv + par3 + fin)
        @test isequal(measure_function(dmref2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
        @test InfiniteOpt._object_numbers(dmref) == [2, 3, 4]
        @test InfiniteOpt._object_numbers(dmref2) == []
        @test InfiniteOpt._parameter_numbers(dmref) == [2, 3, 4, 5]
        @test InfiniteOpt._parameter_numbers(dmref2) == []
        # undo changes
        meas = Measure(inf + par - x + rv + par2 + par3 + fin, data, [2, 4], [2, 5], false)
        InfiniteOpt._set_core_variable_object(dmref, meas)
        meas = Measure(par2, data, [2], [2], false)
        InfiniteOpt._set_core_variable_object(dmref2, meas)
    end
    # test _update_constraints
    @testset "_update_constraints" begin
        @test InfiniteOpt._update_constraints(m, par2) isa Nothing
        @test isequal_canonical(jump_function(constraint_object(con)), inf2 + inf4 + par3 + fin)
        @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
        @test isempty(setdiff(InfiniteOpt._object_numbers(con), [1, 2, 3, 4]))
        @test InfiniteOpt._object_numbers(con2) == []
        @test !is_valid(m, con3)
        # undo changes
        constr = ScalarConstraint(inf2 + inf4 - par2 + par3 + fin, MOI.LessThan(0.))
        InfiniteOpt._set_core_constraint_object(con, constr)
        constr = ScalarConstraint(par2, MOI.GreaterThan(0.))
        InfiniteOpt._set_core_constraint_object(con2, constr)
        unregister(m, :con3)
        @constraint(m, con3, [par2, par2] in MOI.Zeros(2))
    end
    # test _update_number_list
    @testset "_update_number_list" begin
        list = [1, 2, 3, 4, 5]
        @test InfiniteOpt._update_number_list([2, 3], list) isa Nothing
        @test list == [1, 2, 3]
        list = [1, 3, 4, 5]
        @test InfiniteOpt._update_number_list([2], list) isa Nothing
        @test list == [1, 2, 3, 4]
    end
    # test JuMP.delete for IndependentParameters
    @testset "JuMP.delete (IndependentParameterRef)" begin
        # add parameter function for object number updating
        pfunc = parameter_function(sin, par)
        # test used by infinite variable
        @test_throws ErrorException delete(m, par2)
        @test delete(m, inf2) isa Nothing
        @test !is_valid(m, inf2)
        @test !is_valid(m, pt2)
        # test parameter function dependency
        push!(InfiniteOpt._parameter_function_dependencies(par2), ParameterFunctionIndex(1))
        @test_throws ErrorException delete(m, par2)
        popfirst!(InfiniteOpt._parameter_function_dependencies(par2))
        # test normal usage
        @test isa(delete(m, par2), Nothing)
        @test !is_valid(m, par2)
        @test !is_valid(m, d1)
        @test isequal_canonical(measure_function(dmref), inf + par - x + rv + par3 + fin)
        @test isequal_canonical(measure_function(dmref2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
        @test isequal(parameter_refs(inf4), (par, pars))
        @test isequal(parameter_refs(rv), (pars,))
        @test isequal_canonical(jump_function(constraint_object(con)), inf4 + par3 + fin)
        @test isequal(jump_function(constraint_object(con2)),zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
        expected = [IndependentParameterIndex(1), DependentParametersIndex(1),
                    IndependentParameterIndex(3)]
        @test InfiniteOpt._param_object_indices(m) == expected
        @test InfiniteOpt._last_param_num(m) == 4
        @test InfiniteOpt._object_number(dpar) == 1
        @test InfiniteOpt._object_number(dpars[1]) == 2
        @test InfiniteOpt._object_number(dpar3) == 3
        @test InfiniteOpt._parameter_number(dpar3) == 4
        @test InfiniteOpt._parameter_number(dpars[2]) == 3
        @test InfiniteOpt._parameter_numbers(dinf4) == [1, 2, 3]
        @test isempty(setdiff(InfiniteOpt._object_numbers(inf3), [1, 2]))
        @test InfiniteOpt._parameter_numbers(dinf) == [1]
        @test isempty(setdiff(InfiniteOpt._object_numbers(drv), [1, 2]))
        @test InfiniteOpt._object_numbers(dmref) == [3]
        @test InfiniteOpt._parameter_numbers(dmref) == [4]
        @test isempty(setdiff(InfiniteOpt._object_numbers(con), [1, 2, 3]))
        @test isempty(InfiniteOpt._object_numbers(con2))
        @test !is_valid(m, con3)
        # test special measure case with single parameter (possible through
        # multiple deletions) and with single parameter in constraint
        @test isa(delete(m, par3), Nothing)
        @test isequal_canonical(measure_function(dmref), inf + par - x + rv + fin)
        @test isequal_canonical(jump_function(constraint_object(con)), inf4 + fin)
        @test isempty(setdiff(InfiniteOpt._object_numbers(con), [1, 2]))
        expected = [IndependentParameterIndex(1), DependentParametersIndex(1)]
        @test InfiniteOpt._param_object_indices(m) == expected
        @test InfiniteOpt._last_param_num(m) == 3
        @test InfiniteOpt._parameter_number(dpars[2]) == 3
        @test InfiniteOpt._parameter_numbers(dinf4) == [1, 2, 3]
        # test invalid parameter
        @test_throws AssertionError delete(m, par2)
        @test_throws AssertionError delete(m, par3)
        # test measure data check
        @test_throws ErrorException delete(m, par)
    end
    # test delete for finite parameters
    @testset "JuMP.delete (FiniteParameterRef)" begin
        # test first case
        @test delete(m, fin) isa Nothing
        @test isequal_canonical(measure_function(dmref), inf + par - x + rv)
        @test isequal_canonical(jump_function(constraint_object(con)), inf4 + 0.)
        @test isequal(objective_function(m), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
        @test objective_sense(m) == MOI.FEASIBILITY_SENSE
        # test error
        @test_throws AssertionError delete(m, fin)
        # test different objective
        @objective(m, Min, fin2 + x + 1)
        @test delete(m, fin2) isa Nothing
        @test isequal_canonical(objective_function(m), x + 1)
        @test objective_sense(m) == MOI.MIN_SENSE
    end
end

# test deletion of dependent parameters
@testset "Dependent Parameters" begin
    # initialize model
    m = InfiniteModel()
    # setup the parameters
    @infinite_parameter(m, par in [0, 1])
    @infinite_parameter(m, pars[1:2] in [0, 1])
    @infinite_parameter(m, par2 in [0, 1])
    dpar = dispatch_variable_ref(par)
    dpars = dispatch_variable_ref.(pars)
    dpar2 = dispatch_variable_ref(par2)
    # setup the variables
    @variable(m, inf, Infinite(par))
    @variable(m, inf2, Infinite(par, par2))
    @variable(m, inf3, Infinite(par, pars))
    @variable(m, pt, Point(inf, 0.5))
    pt2 = @variable(m, variable_type = Point(inf2, 0.5, 0.5))
    @variable(m, pt3, Point(inf3, 0, [0, 0]))
    @variable(m, x)
    dinf = dispatch_variable_ref(inf)
    dinf2 = dispatch_variable_ref(inf2)
    dinf3 = dispatch_variable_ref(inf3)
    dpt = dispatch_variable_ref(pt)
    dpt2 = dispatch_variable_ref(pt2)
    dpt3 = dispatch_variable_ref(pt3)
    dx = dispatch_variable_ref(x)
    # setup the measure
    data = TestData(par, 0, 1)
    mref = measure(inf + par - x + pars[1], data)
    data = TestData(pars[2], 0, 1)
    mref2 = measure(inf2, data)
    dmref = dispatch_variable_ref(mref)
    dmref2 = dispatch_variable_ref(mref2)
    # setup the derivatives
    d1 = @deriv(mref, pars[1])
    # setup the constraints
    @constraint(m, con, inf2 - par2 + pars[1] + pars[2] <= 0)
    @constraint(m, con2, pars in MOI.Zeros(2))
    # test delete for dependent parameters
    @testset "JuMP.delete (DependentParameterRefs)" begin
        # test measure error
        @test_throws ErrorException delete(m, pars)
        mindex = index(dmref2)
        filter!(e -> e != mindex, InfiniteOpt._measure_dependencies(inf2))
        filter!(e -> e != mindex, InfiniteOpt._measure_dependencies(pars[2]))
        InfiniteOpt._delete_data_object(dmref2)
        # test infinite variable error
        @test_throws ErrorException delete(m, pars)
        @test delete(m, inf3) isa Nothing
        @test !is_valid(m, inf3)
        @test !is_valid(m, pt3)
        # test parameter function dependency
        data = InfiniteOpt._data_object(first(pars))
        push!(data.parameter_func_indices, ParameterFunctionIndex(1))
        @test_throws ErrorException delete(m, pars)
        empty!(data.parameter_func_indices)
        # test regular
        @test delete(m, pars) isa Nothing
        @test !is_valid(m, d1)
        @test isequal(parameter_refs(dinf), (par,))
        @test string(dinf) == "inf(par)"
        @test parameter_values(dpt2) == (0.5, 0.5)
        @test isequal_canonical(measure_function(dmref), inf + par - x)
        @test isequal_canonical(jump_function(constraint_object(con)), inf2 - par2)
        expected = [IndependentParameterIndex(1), IndependentParameterIndex(2)]
        @test InfiniteOpt._param_object_indices(m) == expected
        @test InfiniteOpt._last_param_num(m) == 2
        @test InfiniteOpt._object_number(dpar) == 1
        @test InfiniteOpt._object_number(dpar2) == 2
        @test InfiniteOpt._parameter_number(dpar2) == 2
        @test isempty(setdiff(InfiniteOpt._object_numbers(dinf2), [1, 2]))
        @test InfiniteOpt._parameter_numbers(dinf) == [1]
        @test InfiniteOpt._object_numbers(dmref) == []
        @test InfiniteOpt._parameter_numbers(dmref) == []
        @test isempty(setdiff(InfiniteOpt._object_numbers(con), [1, 2]))
        @test !is_valid(m, con2)
        # test assertion error
        @test_throws AssertionError delete(m, pars)
    end
end

# Test semi_infinite variable deletion
@testset "JuMP.delete (SemiInfinite Variables)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @infinite_parameter(m, par2 in [0, 1])
    @variable(m, inf, Infinite(par, par2))
    @variable(m, pt, Point(inf, 0.5, 0.5))
    @variable(m, x)
    var = build_variable(error, inf, Dict{Int, Float64}(2 => 0.5), check = false)
    rv = add_variable(m, var)
    var = build_variable(error, inf, Dict{Int, Float64}(2 => 0), check = false)
    rv2 = add_variable(m, var)
    data = TestData(par, 0, 1)
    meas = measure(inf + par - x + rv, data)
    meas2 = measure(rv2, data)
    d1 = @deriv(rv, par)
    @constraint(m, con, x + rv <= 0)
    constr = ScalarConstraint(rv2, MOI.GreaterThan(0.))
    con2 = add_constraint(m, constr)
    @constraint(m, con3, [rv, rv] in MOI.Zeros(2))
    # test normal deletion
    @test isa(delete(m, rv), Nothing)
    @test !is_valid(m, d1)
    @test isequal_canonical(measure_function(meas), inf + par - x)
    @test InfiniteOpt._object_numbers(meas) == [2]
    @test isequal_canonical(jump_function(constraint_object(con)), x + 0)
    @test InfiniteOpt._object_numbers(con) == []
    @test InfiniteOpt._semi_infinite_variable_dependencies(inf) == [JuMP.index(rv2)]
    @test !haskey(InfiniteOpt._data_dictionary(m, SemiInfiniteVariable), JuMP.index(rv))
    @test !is_valid(m, rv)
    # test deletion of special cases
    @test isa(delete(m, rv2), Nothing)
    @test isequal(measure_function(meas2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(meas2) == []
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(con2) == []
    @test InfiniteOpt._semi_infinite_variable_dependencies(inf) == []
    @test !haskey(InfiniteOpt._data_dictionary(m, SemiInfiniteVariable), JuMP.index(rv2))
    @test isempty(m.semi_lookup)
    # test error
    @test_throws AssertionError delete(m, rv)
    @test_throws AssertionError delete(m, rv2)
end

# Test variable deletion
@testset "JuMP.delete (Finite Variables)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, 0 <= x <= 1, Bin)
    @variable(m, y == 1, Int)
    data = TestData(par, 0, 0)
    meas1 = measure(x + y + par, data)
    meas2 = measure(y, data)
    @constraint(m, con1, x + y + par <= 0)
    con2 = add_constraint(m, ScalarConstraint(y, MOI.LessThan(0.)))
    @constraint(m, con3, [x, y] in MOI.Zeros(2))
    @constraint(m, con4, max(par, x) <= 0)
    @objective(m, Min, x + y)
    # test deletion of x
    @test isa(delete(m, x), Nothing)
    @test num_constraints(m) == 5
    @test isequal_canonical(measure_function(meas1), y + par)
    @test isequal_canonical(jump_function(constraint_object(con1)), y + par)
    @test isequal_canonical(objective_function(m), y + 0)
    @test isequal_canonical(jump_function(constraint_object(con4)), max(par, 0))
    @test !is_valid(m, con3)
    @test !haskey(InfiniteOpt._data_dictionary(m, FiniteVariable), JuMP.index(x))
    # test deletion of y
    set_objective_function(m, y)
    @test isa(delete(m, y), Nothing)
    @test num_constraints(m) == 3
    @test isequal_canonical(measure_function(meas1), par + 0)
    @test isequal(measure_function(meas2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test isequal_canonical(jump_function(constraint_object(con1)), par + 0)
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test isequal(objective_function(m), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test !haskey(InfiniteOpt._data_dictionary(m, FiniteVariable), JuMP.index(y))
    # test errors
    @test_throws AssertionError delete(m, x)
    @test_throws AssertionError delete(m, y)
end

# Test variable deletion
@testset "JuMP.delete (Point Variables)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, inf, Infinite(par))
    @variable(m, 0 <= x <= 1, Point(inf, 0), Bin)
    @variable(m, y == 1, Point(inf, 1), Int)
    data = TestData(par, 0, 1)
    meas1 = measure(x + y + par, data)
    meas2 = measure(y, data)
    @constraint(m, con1, x + y + par <= 0)
    con2 = add_constraint(m, ScalarConstraint(y, MOI.LessThan(0.)))
    @constraint(m, con3, [x, y] in MOI.Zeros(2))
    @objective(m, Min, x + y)
    # test deletion of x
    @test isa(delete(m, x), Nothing)
    @test num_constraints(m) == 4
    @test isequal_canonical(measure_function(meas1), y + par)
    @test InfiniteOpt._object_numbers(meas1) == []
    @test isequal_canonical(jump_function(constraint_object(con1)), y + par)
    @test InfiniteOpt._object_numbers(con1) == [1]
    @test isequal_canonical(objective_function(m), y + 0)
    @test !haskey(InfiniteOpt._data_dictionary(m, PointVariable), JuMP.index(x))
    @test !is_valid(m, con3)
    # test deletion of y
    set_objective_function(m, y)
    @test isa(delete(m, y), Nothing)
    @test num_constraints(m) == 2
    @test isequal_canonical(measure_function(meas1), par + 0)
    @test isequal(measure_function(meas2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test isequal_canonical(jump_function(constraint_object(con1)), par + 0)
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test isequal(objective_function(m), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test !haskey(InfiniteOpt._data_dictionary(m, PointVariable), JuMP.index(y))
    @test isempty(m.point_lookup)
    # test errors
    @test_throws AssertionError delete(m, x)
    @test_throws AssertionError delete(m, y)
end

# Test variable deletion
@testset "JuMP.delete (Infinite Variables)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, 0 <= x <= 1, Infinite(par), Bin)
    @variable(m, y == 1, Infinite(par), Int)
    @variable(m, x0, Point(x, 0))
    var = build_variable(error, x, Dict{Int, Float64}(1 => 0.5), check = false)
    rv = add_variable(m, var)
    data = TestData(par, 0, 1)
    meas1 = measure(x + y + par, data)
    meas2 = measure(y, data)
    d1 = @deriv(x, par)
    @constraint(m, con1, x + y + par <= 0)
    con2 = add_constraint(m, ScalarConstraint(y, MOI.LessThan(0.)))
    @constraint(m, con3, [x, y] in MOI.Zeros(2))
    # test deletion of x
    @test isa(delete(m, x), Nothing)
    @test num_constraints(m) == 4
    @test isequal_canonical(measure_function(meas1), y + par)
    @test InfiniteOpt._object_numbers(meas1) == []
    @test isequal_canonical(jump_function(constraint_object(con1)), y + par)
    @test InfiniteOpt._object_numbers(con1) == [1]
    @test InfiniteOpt._infinite_variable_dependencies(par) == [index(y)]
    @test !is_valid(m, rv)
    @test !is_valid(m, x0)
    @test !is_valid(m, d1)
    @test !haskey(InfiniteOpt._data_dictionary(m, InfiniteVariable), JuMP.index(x))
    @test !is_valid(m, con3)
    # test deletion of y
    @test isa(delete(m, y), Nothing)
    @test num_constraints(m) == 2
    @test isequal_canonical(measure_function(meas1), par + 0)
    @test isequal(measure_function(meas2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(meas1) == []
    @test isequal_canonical(jump_function(constraint_object(con1)), par + 0)
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(con1) == [1]
    @test InfiniteOpt._object_numbers(con2) == []
    @test InfiniteOpt._infinite_variable_dependencies(par) == []
    @test !haskey(InfiniteOpt._data_dictionary(m, InfiniteVariable), JuMP.index(y))
    # test errors
    @test_throws AssertionError delete(m, x)
    @test_throws AssertionError delete(m, y)
end

# Test derivative deletion
@testset "JuMP.delete (Derivatives)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @variable(m, 0 <= x <= 1, Infinite(par), Bin)
    @variable(m, y == 1, Infinite(par), Int)
    d1 = @deriv(x, par)
    d2 = @deriv(y, par)
    d3 = @deriv(d1, par)
    @variable(m, dx0, Point(d1, 0))
    var = build_variable(error, d1, Dict{Int, Float64}(1 => 0.5), check = false)
    rv = add_variable(m, var)
    data = TestData(par, 0, 1)
    meas1 = measure(d1 + y + par, data)
    meas2 = measure(d2, data)
    @constraint(m, con1, d1 + d2 + par <= 0)
    con2 = add_constraint(m, ScalarConstraint(d2, MOI.LessThan(0.)))
    cref = @constraint(m, d1 == 0)
    push!(InfiniteOpt._derivative_constraint_dependencies(d1), index(cref))
    InfiniteOpt._set_has_derivative_constraints(par, true)
    @constraint(m, con3, [d1, d2] in MOI.Zeros(2))
    # test deletion of d1
    @test isa(delete(m, d1), Nothing)
    @test num_constraints(m) == 7
    @test isequal_canonical(measure_function(meas1), y + par)
    @test InfiniteOpt._object_numbers(meas1) == []
    @test isequal_canonical(jump_function(constraint_object(con1)), d2 + par)
    @test InfiniteOpt._object_numbers(con1) == [1]
    @test InfiniteOpt._derivative_dependencies(par) == [index(d2)]
    @test !is_valid(m, rv)
    @test !is_valid(m, dx0)
    @test !is_valid(m, d3)
    @test !is_valid(m, cref)
    @test !is_valid(m, con3)
    @test !haskey(InfiniteOpt._data_dictionary(m, Derivative), JuMP.index(d1))
    @test !has_derivative_constraints(par)
    @test !haskey(m.deriv_lookup, (x, par))
    # test deletion of d2
    @test isa(delete(m, d2), Nothing)
    @test num_constraints(m) == 7
    @test isequal_canonical(measure_function(meas1), y + par)
    @test isequal(measure_function(meas2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(meas1) == []
    @test isequal_canonical(jump_function(constraint_object(con1)), par + 0)
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(con1) == [1]
    @test InfiniteOpt._object_numbers(con2) == []
    @test InfiniteOpt._derivative_dependencies(par) == []
    @test !haskey(InfiniteOpt._data_dictionary(m, Derivative), JuMP.index(d2))
    @test !haskey(m.deriv_lookup, (y, par))
    # test errors
    @test_throws AssertionError delete(m, d1)
    @test_throws AssertionError delete(m, d2)
end

# Test infinite parameter function deletion
@testset "JuMP.delete (Parameter Function)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    f = parameter_function(sin, par)
    d = deriv(f, par)
    d = deriv(f, par)
    rv = add_variable(m, build_variable(error, f, Dict(1 => 0.)))
    data = TestData(par, 0, 0)
    meas1 = measure(f + par, data)
    meas2 = measure(f, data)
    @constraint(m, con1, f + par <= 0)
    con2 = add_constraint(m, ScalarConstraint(f, MOI.LessThan(0.)))
    @constraint(m, con3, [f, f] in MOI.Zeros(2))
    # test deletion of x
    @test isa(delete(m, f), Nothing)
    @test num_constraints(m) == 2
    @test isequal_canonical(measure_function(meas1), par + 0)
    @test isequal(measure_function(meas2), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test isequal_canonical(jump_function(constraint_object(con1)), par + 0)
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test !haskey(InfiniteOpt._data_dictionary(f), JuMP.index(f))
    @test !is_valid(m, f)
    @test !is_valid(m, d)
    @test !is_valid(m, rv)
    @test !is_valid(m, con3)
    # test errors
    @test_throws AssertionError delete(m, f)
end

# Test variable deletion
@testset "JuMP.delete (Measures)" begin
    # intialize the model
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 1])
    @infinite_parameter(m, par2 in [0, 1])
    @infinite_parameter(m, pars[1:2] in [0, 1])
    @variable(m, w, Infinite(pars))
    @variable(m, x, Infinite(par))
    @variable(m, y, Infinite(par2))
    @variable(m, x0, Point(x, 0))
    var = build_variable(error, x, Dict{Int, Float64}(1 => 0.5), check = false)
    rv = add_variable(m, var)
    data = TestData(par, 0, 1)
    data2 = TestData(pars, [0, 0], [1, 1])
    meas = measure(x, data)
    meas1 = measure(x + x0 + rv + par + meas, data)
    meas2 = measure(x + x0 + rv + par, data)
    meas3 = measure(meas1 + x0, data)
    meas4 = measure(meas2, data)
    meas5 = measure(w, data2)
    meas6 = measure(x + y, data)
    d1 = @deriv(meas6, par2)
    @constraint(m, con1, x0 + meas1 <= 0)
    con2 = add_constraint(m, ScalarConstraint(meas2, MOI.LessThan(0.)))
    @constraint(m, con3, [meas1, meas1] in MOI.Zeros(2))
    @objective(m, Min, meas1 + x0)
    # test deletion of meas1
    @test isa(delete(m, meas1), Nothing)
    @test isequal_canonical(measure_function(meas3), x0 + 0)
    @test InfiniteOpt._object_numbers(meas3) == []
    @test isequal_canonical(jump_function(constraint_object(con1)), x0 + 0)
    @test InfiniteOpt._object_numbers(con1) == []
    @test isequal_canonical(objective_function(m), x0 + 0)
    @test InfiniteOpt._measure_dependencies(x) == [JuMP.index(meas), JuMP.index(meas2), JuMP.index(meas6)]
    @test InfiniteOpt._measure_dependencies(y) == [JuMP.index(meas6)]
    @test length(InfiniteOpt._measure_dependencies(par)) == 5
    @test InfiniteOpt._measure_dependencies(rv) == [JuMP.index(meas2)]
    @test !haskey(InfiniteOpt._data_dictionary(m, Measure), JuMP.index(meas1))
    @test !is_valid(m, con3)
    # test deletion of meas2
    set_objective_function(m, meas2)
    @test isa(delete(m, meas2), Nothing)
    @test isequal(measure_function(meas4), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(meas4) == []
    @test isequal(jump_function(constraint_object(con2)), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._object_numbers(con2) == []
    @test isequal(objective_function(m), zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}))
    @test InfiniteOpt._measure_dependencies(x) == [JuMP.index(meas), JuMP.index(meas6)]
    @test InfiniteOpt._measure_dependencies(y) == [JuMP.index(meas6)]
    @test length(InfiniteOpt._measure_dependencies(par)) == 4
    @test InfiniteOpt._measure_dependencies(rv) == []
    @test !haskey(InfiniteOpt._data_dictionary(m, Measure), JuMP.index(meas2))
    # test deletion of meas5
    @test isa(delete(m, meas5), Nothing)
    @test InfiniteOpt._measure_dependencies(pars[1]) == []
    @test InfiniteOpt._measure_dependencies(pars[2]) == []
    # test deletion of meas6
    @test isa(delete(m, meas6), Nothing)
    @test !is_valid(m, d1)
    # test errors
    @test_throws AssertionError delete(m, meas1)
    @test_throws AssertionError delete(m, meas2)
end


text:
using LinearMaps
using Base.Test

import Base: *

A=2*rand(Complex128,(20,10)).-1
v=rand(Complex128,10)
w=rand(Complex128,20)

# test wrapped map for matrix
M=LinearMap(A)
@test M*v==A*v

# test transposition and full
@test M'*w==A'*w
@test M.'*w==A.'*w

@test full(M)==A
@test full(M')==A'
@test full(M.')==A.'

# test function map
F=LinearMap(cumsum,2)
@test full(F)==[1. 0.;1. 1.]

N=100
F=LinearMap(fft,N;isreal=false)/sqrt(N)
U=full(F) # will be a unitary matrix
@test_approx_eq U'*U eye(N)

F=LinearMap(cumsum,10)
@test F*v==cumsum(v)
@test_throws ErrorException F'*v

# test linear combinations
A=2*rand(Complex128,(10,10)).-1
M=LinearMap(A)
v=rand(Complex128,10)

@test full(3*M)==3*A
@test full(M+A)==2*A
@test full(-M)==-A
@test full(3*M'-F)==3*A'-full(F)
@test (3*M-1im*F)'==3*M'+1im*F'

@test_approx_eq (2*M'+3*I)*v (2*A'+3*I)*v

# test composition
@test (F*F)*v==F*(F*v)
@test (F*A)*v==F*(A*v)
@test_approx_eq full(M*M.') A*A.'
@test !isposdef(M*M.')
@test isposdef(M*M')
@test isposdef(F.'*F)
@test isposdef((M*F)'*M*F)
@test (M*F).'==F.'*M.'

L=3*F+1im*A+F*M'*F
LF=3*full(F)+1im*A+full(F)*full(M)'*full(F)
@test_approx_eq full(L) LF

# test inplace operations
w=similar(v)
Base.A_mul_B!(w,L,v)
@test_approx_eq w LF*v

# test new type
type SimpleFunctionMap <: AbstractLinearMap{Float64}
    f::Function
    N::Int
end

Base.size(A::SimpleFunctionMap)=(A.N,A.N)
Base.issym(A::SimpleFunctionMap)=false
*(A::SimpleFunctionMap,v::Vector)=A.f(v)

F=SimpleFunctionMap(cumsum,10)
@test size(F,1)==10
w=similar(v)
Base.A_mul_B!(w,F,v)
@test w==F*v
@test_throws MethodError F'*v
@test_throws MethodError F.'*v
@test_throws MethodError Ac_mul_B!(w,F,v)
@test_throws MethodError At_mul_B!(w,F,v)


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
abstract Instance

include("pool.jl")
include("graph.jl")
include("source.jl")
include("instance.jl")
include("evaluate.jl")
include("operators.jl")
include("matmul.jl")
include("convolutions.jl")
include("activation.jl")
include("reduce.jl")
include("limits.jl")
include("mutators.jl")
include("index.jl")
include("update.jl")
include("autodiff.jl")
include("devec.jl")
include("snapshots.jl")


text:
using Pkg
Pkg.activate("./")
# load pkgs
using Ipopt, HarmonicPowerModels, PowerModels
using JuMP #adding avoids problems with Revise

# pkg const
const _PMs = PowerModels
const _HPM = HarmonicPowerModels

# set the solver
solver = Ipopt.Optimizer

# path to the data
path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_3bus2gen.m")

# load data
data  = _PMs.parse_file(path)
# sanity check
pm_fundamental = _PMs.instantiate_model(data, _PMs.IVRPowerModel, _PMs.build_opf_iv);
result_fundamental = optimize_model!(pm_fundamental, optimizer=solver)
result_fundamental["termination_status"]

# set up current ratings consistently
#assume minimum voltage
vmmin = 0.8
for (b,branch) in data["branch"]
    branch["c_rating"] = branch["rate_a"]/vmmin
end
for (d,load) in data["load"]
    load["c_rating"] = abs(load["pd"] + im* load["qd"])/vmmin
end
for (g,gen) in data["gen"]
    gen["c_rating"] = abs(gen["pmax"] + im* gen["qmax"])/vmmin
end
hdata = _HPM.replicate(data)

##
# solve the hopf
# result = run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
pm = _PMs.instantiate_model(hdata, _PMs.IVRPowerModel, _HPM.build_hopf_iv; ref_extensions=[_HPM.ref_add_xfmr!]);
result = optimize_model!(pm, optimizer=solver)

##
##
for (n,nw) in result["solution"]["nw"]
    for (i,bus) in nw["bus"]
            bus["vm"] =  abs(bus["vr"] +im*  bus["vi"])
            bus["va"] =  angle(bus["vr"] +im*  bus["vi"])*180/pi
    end
end
println("Harmonic 3")
_PMs.print_summary(result["solution"]["nw"]["2"])
println("Harmonic 1")
_PMs.print_summary(result["solution"]["nw"]["1"])
result["objective"]
result["termination_status"]


text:
using Random
Random.seed!(5)
using Statistics, ZigZagBoomerang, LinearAlgebra, Test, SparseArrays
parallel = Threads.nthreads() > 1
const D = 2
#n = 1_000_000
n = (2<<19) - D
d = D + n
if parallel
    κ = Threads.nthreads() # number of threads
    Δ = 0.01 # how much samples per epoch
    d2 = d÷κ
    partition = ZigZagBoomerang.Partition(κ, d)
end

K = 20
T = 10.0

const σ2 = 1.0
const γ0 = 0.1

β = randn(D)
X = [randn(n) for i in 1:D]
t = range(0.0, 1.0, length=n)

x = sin.(2pi*t)
y = sum(X .* β) + x + sqrt(σ2)*randn(n)

dia = -0.5ones(n-1)
Γ = 100*sparse(SymTridiagonal(1.0ones(n), dia))^2
GX0 = [X[1] X[2] ones(n)]
GX = [X[1] X[2]]

# rough point estimate, trend as constant or piecewise constant
xpost = 0*x
βpost = (GX0\(y))[1:D]
for i in 1:3 # 3 em steps
    global βpost, xpost
    xpost = (I'*I + σ2*Γ)\(I'*y - sum(X .* βpost))
    βpost = (GX\(y - xpost))
end

# rough estimate of posterior precision matrix using order of standard errors
Γpost = [I(2)*(γ0 + n) zeros(2,n); zeros(2,n)' (Γ + I/σ2)]

if parallel # if we want to do it parallel, we have to break some dependencies in our posterior proxy
    for k1 in 0:κ-1
        for k2 in 0:κ-1
            k2 == k1 && continue
            abs(k2-k1) > 1 && continue
            Γpost[(k1*d2 + 1):(k1+1)*d2,(k2*d2 + 1):(k2+1)*d2] *= 0
        end
    end
end
dropzeros!(Γpost)
μpost = [βpost; xpost]

bias = [0.0, 0.0]
for i in 1:D
    for j = 1:n # use an unbiased estimate of the log-likelihood
        global bias
        bias[i] += X[i][j]*(xpost[j] + X[1][j]*βpost[1] + X[2][j]*βpost[2] - y[j])/σ2
    end
end
bias

# Compute an unbiased estimate of the i'th partial derivative of the negative loglikelihood in a smart way
function ∇ϕ(u, i_, Γ, X, y, (xhat, βhat, bias), K = 20)
    β = u[1:D]
    x = @view u[D+1:end]
    if i_ > D
        i = i_ - D
        a = ZigZagBoomerang.idot(Γ, i, x) + (x[i] - y[i])/σ2
        for di in 1:D
            a += β[di]*X[di][i]/σ2
        end
        return a
    else #x'*C*x +  (a+B*x)'*A*(B*x+a) on http://www.matrixcalculus.org
        i = i_
        a = γ0*β[i]
        s = n/K/σ2
        for k in 1:K # use an unbiased estimate of the log-likelihood with control variate
            j = rand(1:n)
            a += s*X[i][j]*(x[j] + X[1][j]*β[1] + X[2][j]*β[2] - y[j])
            a -= s*X[i][j]*(xhat[j] + X[1][j]*βhat[1] + X[2][j]*βhat[2] - y[j])
        end
        return a + bias[i]
    end
end

t0 = 0.0
u = μpost
x0 = copy(u) + 0.01*randn(d) # jiggle the starting point to see convergence
θ0 = rand([-1.0, 1.0], d)
θ0[1] = θ0[2] = 5*sqrt(1/n)

# Graphical structure of posterior
G2 = [i+D => D .+ rowvals(Γ)[nzrange(Γ, i)] for i in eachindex(y)]
G1 = [1 => [1:d;], 2 => [1:d;]]
G = [G1; G2]

# precision bounds
c = 1e-9*ones(d)
dt = T/100

Z = ZigZag(Γpost, μpost)
if parallel
    trc, _ = @time ZigZagBoomerang.parallel_spdmp(partition, ∇ϕ, t0, x0, θ0, T, c, G, Z, Γ, X, y, (xpost, βpost, bias), K;  adapt=true, progress=true, Δ=Δ)
else
    trc, _ = @time ZigZagBoomerang.spdmp(∇ϕ, t0, x0, θ0, T, c, G, Z, Γ, X, y, (xpost, βpost, bias), K; structured=true, adapt=true, progress=true)
end
J = [1,2, n÷2]
subtrc = subtrace(trc, J)

#ts, xs = ZigZagBoomerang.sep(collect(discretize(subtrc, dt)))
ts, xs = ZigZagBoomerang.sep(subtrc)

# posterior mean
u = mean(trc)
β̂ = u[1:D]
x̂ = u[D+1:end]

β - β̂
x - x̂

fig = Figure()
ax = fig[1,1:3] = Axis(fig, title="x")
fig[2,1] = Axis(fig, title="β1")
fig[2,2] = Axis(fig, title="β2")
fig[2,3] = Axis(fig, title="x[$(J[end]-2)]")
jj = 50
scatter!(ax, t[1:jj:end], y[1:jj:end], markersize=0.1)
lines!(ax, t[1:jj:end], x̂[1:jj:end], color=:lightblue)
lines!(ax, t[1:jj:end], x[1:jj:end], color=:green, linewidth=3.0)
lines!(fig[2,1], ts, getindex.(xs, 1))
lines!(fig[2,1], ts, fill(β[1], length(ts)), color=:green)
lines!(fig[2,2], ts, getindex.(xs, 2))
lines!(fig[2,2], ts, fill(β[2], length(ts)), color=:green)
lines!(fig[2,3], ts, getindex.(xs, 3))
lines!(fig[2,3], ts, fill(x[J[end]-2], length(ts)), color=:green)
display(fig)

using FileIO
save("million.png", fig)


text:
######################################################################
#
#   Lasso Solution
#
######################################################################

struct LassoSolution{T, S}
  x::SparseIterate{T}
  residuals::Vector{T}
  penalty::S
  σ::Union{T, Nothing}

  LassoSolution{T, S}(x::SparseIterate{T}, residuals::AbstractVector{T}, penalty::ProxL1{T}, σ::T) where {T, S} =
    new(x, residuals, penalty, σ)
  LassoSolution{T, S}(x::SparseIterate{T}, residuals::AbstractVector{T}, penalty::ProxL1{T}) where {T, S} =
    new(x, residuals, penalty, nothing)
end

######################################################################
#
#   Lasso Interface
#
######################################################################

function lasso(
  X::StridedMatrix{T},
  y::StridedVector{T},
  λ::T,
  options::CDOptions=CDOptions()) where {T<:AbstractFloat}

  x = SparseIterate( size(X, 2) )
  f = CDLeastSquaresLoss(y, X)
  g = ProxL1(λ)
  coordinateDescent!(x, f, g, options)

  LassoSolution{T, typeof(g)}(x, f.r, g, std(f.r))
end

function lasso(
  X::StridedMatrix{T},
  y::StridedVector{T},
  λ::T,
  ω::Array{T},
  options::CDOptions=CDOptions()) where {T<:AbstractFloat}

  x = SparseIterate( size(X, 2) )
  f = CDLeastSquaresLoss(y, X)
  g = ProxL1(λ, ω)
  coordinateDescent!(x, f, g, options)

  LassoSolution{T, typeof(g)}(x, f.r, g, std(f.r))
end

######################################################################
#
#   Sqrt-Lasso Interface
#
######################################################################

function sqrtLasso(
  X::StridedMatrix{T},
  y::StridedVector{T},
  λ::T,
  options::CDOptions=CDOptions();
  standardizeX=true) where {T<:AbstractFloat}

  p = size(X, 2)
  x = SparseIterate( p )
  f = CDSqrtLassoLoss(y, X)
  if standardizeX
    stdX = Array{T}(p)
    _stdX!(stdX, X)
    g = ProxL1(λ, stdX)
  else
    g = ProxL1(λ)
  end
  coordinateDescent!(x, f, g, options)

  LassoSolution{T, typeof(g)}(x, f.r, g, std(f.r))
end

function sqrtLasso(
  X::StridedMatrix{T},
  y::StridedVector{T},
  λ::T,
  ω::Array{T},
  options::CDOptions=CDOptions()) where {T<:AbstractFloat}

  p = size(X, 2)
  x = SparseIterate( p )
  f = CDSqrtLassoLoss(y, X)
  g = ProxL1(λ, ω)
  coordinateDescent!(x, f, g, options)

  LassoSolution{T, typeof(g)}(x, f.r, g, std(f.r))
end

######################################################################
#
#   Scaled Lasso Interface
#
######################################################################

function scaledLasso!(
  x::SparseIterate{T},
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  λ::T,
  ω::AbstractVector{T},
  options::IterLassoOptions=IterLassoOptions()
  )  where {T<:AbstractFloat}

  n, p = size(X)
  f = CDLeastSquaresLoss(y,X)

  # initializa σ
  if options.initProcedure == :Screening
    σ = _findInitSigma!(X, y, options.sinit, f.r)
  elseif options.initProcedure == :InitStd
    σ = options.σinit
  elseif options.initProcedure == :WarmStart
    initialize!(f, x)
    σ = std(f.r)
  else
    throw(ArgumentError("Incorrect initialization Symbol"))
  end

  g = ProxL1(λ * σ, ω)
  for iter=1:options.maxIter
    coordinateDescent!(x, f, g, options.optionsCD)
    σnew = sqrt( sum(abs2, f.r) / n )

    if abs(σnew - σ) / σ < options.optTol
      break
    end
    σ = σnew
    g = ProxL1(λ * σ, ω)
  end

  LassoSolution{T, typeof(g)}(x, f.r, g, std(f.r))
end

######################################################################
#
#   Lasso Path Interface
#
######################################################################

function feasibleLasso!(
  x::SparseIterate{T},
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  λ0::T,
  options::IterLassoOptions=IterLassoOptions()
  )  where {T<:AbstractFloat}

  n, p = size(X)
  f = CDLeastSquaresLoss(y,X)
  Γ = Array{T}(p)              # stores loadings
  Γold = Array{T}(p)

  # initializa residuals
  if options.initProcedure == :Screening
    _findInitResiduals!(X, y, options.sinit, f.r)
  elseif options.initProcedure == :InitStd
    σ = options.σinit
    _stdX!(Γ, X)
    coordinateDescent!(x, f, ProxL1(λ0*σ, Γ), options.optionsCD)
  elseif options.initProcedure == :WarmStart
    initialize!(f, x)
  else
    throw(ArgumentError("Incorrect initialization Symbol"))
  end
  _getLoadings!(Γ, X, f.r)

  g = ProxL1(λ0, Γ)
  for iter=1:options.maxIter
    copy!(Γold, Γ)

    coordinateDescent!(x, f, g, options.optionsCD)
    _getLoadings!(Γ, X, f.r)

    if maximum(abs.(Γold  - Γ)) / maximum(Γ) < options.optTol
      break
    end
  end

  LassoSolution(x, f.r, g, std(f.r))
end

######################################################################
#
#   Lasso Path Interface
#
######################################################################

struct LassoPath{T<:AbstractFloat}
  λpath::Vector{T}
  βpath::Vector{SparseIterate{T,1}}
end

function refitLassoPath(
  path::LassoPath{T},
  X::StridedMatrix{T},
  Y::StridedVector{T}) where {T<:AbstractFloat}

  λpath = path.λpath
  βpath = path.βpath

  out = Dict{Vector{Int64},Vector{Float64}}()
  for i=1:length(λpath)
    S = findall(!iszero, βpath[i])
    if haskey(out, S)
      continue
    end
    out[S] = X[:, S] \ Y
  end
  out
end

# λArr is in decreasing order
function LassoPath(
  X::StridedMatrix{T},
  Y::StridedVector{T},
  λpath::Vector{T},
  options=CDOptions();
  max_hat_s=Inf, standardizeX::Bool=true) where {T<:AbstractFloat}

  n, p = size(X)
  stdX = Array{T}(undef, p)
  if standardizeX
    _stdX!(stdX, X)
  else
    fill!(stdX, one(T))
  end

  x = SparseIterate(T, p)
  f = CDLeastSquaresLoss(Y, X)

  numλ  = length(λpath)
  βpath = Vector{SparseIterate{T}}(undef, numλ)

  for indλ=1:numλ
    coordinateDescent!(x, f, ProxL1(λpath[indλ], stdX), options)
    βpath[indλ] = copy(x)
    if nnz(x) > max_hat_s
      resize!(λpath, indλ)
      break
    end
  end

  LassoPath{T}(copy(λpath), βpath)
end


text:
function print_documentation(io::IO, node::ExprNode, indent, options, members::Bool=false; kwargs...)
    print_documentation(io, node.cursor, indent, options, members; kwargs...)
end

const ESCAPE_PATTERN = r"""(\$|\\|"(?=""))"""

function print_documentation(io::IO, cursor::CLCursor, indent, options, members::Bool=false;
        prologue::Vector{String}=String[], epilogue::Vector{String}=String[])

    fold_single_line_comment = get(options, "fold_single_line_comment", false)
    extract_c_comment_style = get(options, "extract_c_comment_style", "disable")
    show_c_function_prototype = get(options, "show_c_function_prototype", false)
    ids = get(options, "DAG_ids", Dict())

    if extract_c_comment_style == "doxygen"
        doc = format_doxygen(cursor, options, members)
    elseif extract_c_comment_style == "raw"
        doc = format_raw(cursor, members)
    else
        return
    end

    if show_c_function_prototype && cursor isa Clang.CLFunctionDecl
        prototype = ["### Prototype", "```c", get_prototype(cursor), "```"]
    else
        prototype = String[]
    end

    # Do not print """ if no doc
    all(isempty, doc) && isempty(prototype) && return
    doc = [prologue; doc; prototype; epilogue]
    # _ is already escaped with one slash, here we add one more
    doc = replace.(doc, ESCAPE_PATTERN=>s"\\\1")
    last(doc) == "" && pop!(doc)
    if length(doc) == 1 && fold_single_line_comment
        line = only(doc)
        println(io, indent, '"'^3, line, '"'^3)
    else
        println(io, indent * '"'^3)
        for line in doc
            println(io, indent, line)
        end
        println(io, indent * '"'^3)
    end
end

function get_prototype(node)
    code = Clang.getSourceCode(node.cursor)
    code = replace(code, r"/\*.*?\*/|//.*?$"m=>"")
    code = replace(code, r"\s+"m=>" ")
    "$code;"
end

escape_underscore(word) = replace(word, "_"=>"\\_")

"Replace identifiers with @ref"
function gen_automatic_links(text, ids)
    result = []
    # ...`id`(...)
    ESCAPE_PATTERN = r"^(.*?)(%)?(\*{1,2}|`?)([\w_][\w\d_]*)\3((?:\(.*?\))?.*)$"
    while true
        m = match(ESCAPE_PATTERN, text)

        isnothing(m) && break
        prev, escape, marker, word, next = m[1], m[2], m[3], m[4], m[5]
        text = next

        if isnothing(escape) && Symbol(word) in keys(ids)
            ref = "[`$word`](@ref)"
        else
            word = marker == "`" ? word : escape_underscore(word)
            ref = marker * word * marker
        end

        push!(result, prev, ref)
    end
    push!(result, text)
    return join(result)
end

function format_raw(cursor, members=false)
    comment = Clang.getRawCommentText(cursor)
    doc = strip_comment_markers(comment)
    member_docs = members ? get_member_doc(format_raw, cursor) : String[]
    isempty(member_docs) ? doc : [doc; ""; member_docs]
end

function strip_comment_markers(s::AbstractString)::Vector
    # Assume the comments are either /**/ or //,
    # although Clang may associate multiple consecutive comments as a whole

    # // /*
    first_line_pattern = r"^\s*(/{2,}!?|/\*+!?)\s*(.*)$"
    # // *
    middle_line_pattern = r"^\s*(/{2,}!?\s?|\*?\s?)(.*)$"
    # // */
    last_line_pattern = r"^\s*(/{2,}!?\s*(.*)|(.*?)\*+/)$"
    # // /**/
    single_line_pattern = r"^\s*(/\*+<?\s*(.*)\s*\*/|/{2,}!?<?\s*(.*))$"

    lines = split(s, '\n')
    if length(lines) == 0
        return String[]
    elseif length(lines) == 1
        m = match(single_line_pattern, only(lines))
        if isnothing(m)
            return [s]
        else
            return [isnothing(m[2]) ? m[3] : m[2]]
        end
    else
        first = lines[1]
        last = lines[end]
        middle = lines[2:end-1]
        stripped = SubString{String}[]
        m_first = match(first_line_pattern, first)
        m_last = match(last_line_pattern, last)
        if isnothing(m_first)
            push!(stripped, first)
        elseif !isempty(m_first[2])
            push!(stripped, m_first[2])
        end
        for line in middle
            m = match(middle_line_pattern, line)
            if isnothing(m)
                push!(stripped, line)
            else
                push!(stripped, m[2])
            end
        end
        if isnothing(m_last)
            push!(stripped, last)
        else
            m = isnothing(m_last[2]) ? m_last[3] : m_last[2]
            isempty(m) || push!(stripped, m)
        end
        return stripped
    end
end

function join_row(row)
    parts = ["|"]
    for i in row
        push!(parts, " ", i, " |")
    end
    join(parts)
end

function render_table(table::AbstractMatrix{<:AbstractString})
    size(table, 2) <= 1 && return String[]
    slash_width = count.(ESCAPE_PATTERN, table) .* textwidth('\\')
    raw_width = textwidth.(table)
    widths = maximum(raw_width .+ slash_width, dims=2)
    cell_widths = max.(widths, 4)
    lines = String[]
    titles = rpad.(table[:, 1], cell_widths  .- slash_width[:, 1])

    push!(lines, join_row(titles))
    length_row = rpad.(":", cell_widths, '-')
    push!(lines, join_row(length_row))
    last_col = fill("", size(table, 1))
    for (col,slash_col) in Iterators.drop(zip(eachcol(table), eachcol(slash_width)), 1)
        for i in eachindex(col)
            if col[i] == last_col[i]
                col[i] = ""
                slash_col[i] = 0
            else
                last_col[i] = col[i]
            end
        end
        row = rpad.(col, cell_widths .- slash_col)
        push!(lines, join_row(row))
    end
    lines
end

get_member_doc(f, cursor::CLCursor) = String[]

function get_member_doc(f, cursor::CLEnumDecl)
    table = ["Enumerator", "Note"]
    for c in children(cursor)
        c isa CLEnumConstantDecl || continue
        name = escape_underscore(spelling(c))
        doc = join(f(c), ' ')
        doc = replace(doc, '\n'=>' ')
        isempty(doc) || push!(table, name, doc)
    end
    return render_table(reshape(table, 2, :))
end

function get_member_doc(f, cursor::CLStructDecl)
    table = ["Field", "Note"]
    for c in children(cursor)
        c isa CLFieldDecl || continue
        name = escape_underscore(spelling(c))
        doc = join(f(c), ' ')
        doc = replace(doc, '\n'=>' ')
        isempty(doc) || push!(table, name, doc)
    end
    return render_table(reshape(table, 2, :))
end

function format_doxygen(cursor, options, members=false)
    comment = Clang.getParsedComment(cursor)
    child_nodes = children(comment)
    # Collect parameters
    parameters, paragraphs = [], []
    returns = missing
    seealso = String[]
    for c in child_nodes
        if c isa Clang.ParamCommand
            push!(parameters, c)
        elseif c isa Clang.BlockCommand
            name = Clang.getCommandName(c)
            if name in ["returns", "return"]
                returns = c
            elseif name in ["sa", "see"]
                append!(seealso, format_block(Clang.getParagraph(c), options))
            else
                push!(paragraphs, c)
            end
        elseif c isa Union{Clang.Paragraph,Clang.VerbatimLine,Clang.VerbatimBlockCommand}
            push!(paragraphs, c)
        else
            @info "Unhandled toplevel element: $c $(children(c))"
        end
    end

    member_docs = members ? get_member_doc(x->format_doxygen(x, options, false), cursor) : String[]

    lines = String[]
    for p in paragraphs
        Clang.isWhiteSpace(p) && continue
        t = format_block(p, options)
        append!(lines, t)
        push!(lines, "")
    end
    if !isempty(parameters)
        push!(lines, "### Parameters")
        for p in parameters
            append!(lines, format_parameter(p, options))
            push!(lines, "")
        end
    end
    if !ismissing(returns)
        push!(lines, "### Returns")
        append!(lines, format_block(Clang.getParagraph(returns), options))
        push!(lines, "")
    end
    if !isempty(member_docs)
        append!(lines, member_docs)
    end
    if !isempty(seealso)
        push!(lines, "### See also")
        push!(lines, join(strip.(seealso), ", "), "")
    end
    lines
end

function format_parameter(p, options)
    name = Clang.getParamName(p)
    dir = parameter_pass_direction_name(Clang.getDirection(p))
    content = format_inline.(children(p), Ref(options))
    ["* `$name`:$dir$(content[1])"; @view content[2:end]]
end

"""
    format_block(comment) -> Vector{String}
Format block elements, return a collection of lines.
"""
function format_block end

function format_block(t::Clang.VerbatimLine, options)
    # this command is usually a result of parse error
    ["`" * Clang.getText(t) * "`"]
end

function format_block(x::Clang.Paragraph, options)
    t = format_inline(x, options)
    # Remove leading space
    t = replace(t, r"^\s+"=>"")
    [t]
end

function format_block(t::Clang.VerbatimBlockCommand, options)
    lines = Clang.getText.(children(t))
    # Clang does not parse its own documentation right
    # This is a quick fix for it
    lines = filter(x->!occursin('\n', x), lines)
    ["```c++"; lines; "```"]
end

function format_block(x::Clang.BlockCommand, options)
    name = Clang.getCommandName(x)
    args = join(Clang.getArguments(x), ' ')
    content = only(format_block(Clang.getParagraph(x), options))
    name in ["li", "arg"] && return ["* $content"]
    name in ["brief", "details"] && return [content]
    name in ["note", "warning"] && return ["!!! $name", "", "    $content"]
    ["\\$name$args $content"]
end

"""
    format_inline(comment) -> String
Format inline elements, return a string.
"""
function format_inline end

function format_inline(t::Clang.Text, options)
    text = Clang.getText(t)
    # Fold consecutive space
    text = replace(text, r"(\s)\s+"=>s"\1")
    gen_automatic_links(text, get(options, "DAG_ids", Dict()))
end

function format_inline(t::Clang.VerbatimLine, options)
    "`" * Clang.getText(t) * "`"
end

function format_inline(t::Clang.VerbatimBlockCommand, options)
    lines = Clang.getText.(children(t))
    # Clang does not parse its own documentation right
    # This is a quick fix for it
    lines = filter(x->!occursin('\n', x), lines)
    ["```c++"; lines; "```"]
end

"""Render a "word" argument respecting C's identifier definition."""
function render_argument(raw, options, marker_begin, marker_end=marker_begin)
    # libclang thinks "(\c code)" is (`code)` instead of (`code`)
    m = match(r"^([\w\d_ ]+)(.*)$", raw)
    ids = get(options, "DAG_ids", Dict())
    if isnothing(m)
        "$marker_begin$raw$marker_end"
    else
        code = m[1]
        text = m[2]
        # Replace code with a link
        if Symbol(code) in keys(ids)
            "[`$code`](@ref)$text"
        else
            "$marker_begin$code$marker_end$text"
        end
    end
end

function format_inline(c::Clang.InlineCommand, options)
    k = Clang.getRenderKind(c)
    args = join(Clang.getArguments(c), ' ')
    if k == Clang.CXCommentInlineCommandRenderKind_Normal
        return args
    elseif k == Clang.CXCommentInlineCommandRenderKind_Bold
        return render_argument(args, options, "**")
    elseif k == Clang.CXCommentInlineCommandRenderKind_Monospaced
        return render_argument(args, options, "`")
    elseif k == Clang.CXCommentInlineCommandRenderKind_Emphasized
        return render_argument(args, options, "*")
    elseif k == Clang.CXCommentInlineCommandRenderKind_Anchor
        return ""
    end
end

function format_inline(p::Clang.Paragraph, options)
    join(format_inline.(children(p), Ref(options)))
end

function parameter_pass_direction_name(d)
    ismissing(d) && return ""
    d == Clang.CXCommentParamPassDirection_In && return raw"\[in\]"
    d == Clang.CXCommentParamPassDirection_Out && return raw"\[out\]"
    d == Clang.CXCommentParamPassDirection_InOut && return raw"\[in,out\]"
end

function format_inline(p::Clang.ParamCommand, options)
    name = Clang.getParamName(p)
    name = escape_underscore(name)
    dir = parameter_pass_direction_name(Clang.getDirection(p))
    content = join(format_inline.(children(p), Ref(options)))
    "\\param$dir $name$content"
end

function format_inline(c::Clang.BlockCommand, options)
    name = Clang.getCommandName(c)
    args = join(Clang.getArguments(c), ' ')
    content = format_inline(Clang.getParagraph(c), options)
    "\\$name$args$content"
end

function format_inline(t::Union{Clang.HTMLStartTag, Clang.HTMLEndTag}, options)
    Clang.getAsString(t)
end


text:
module MySQLDatabaseAdapter

import Revise
import MySQL, DataFrames, DataStreams, Logging
using SearchLight, SearchLight.Database, SearchLight.Exceptions

export DatabaseHandle, ResultHandle

#
# Setup
#

const DB_ADAPTER = MySQL
const DEFAULT_PORT = 3306

const COLUMN_NAME_FIELD_NAME = :Field

const DatabaseHandle = DB_ADAPTER.MySQLHandle
const ResultHandle   = DataStreams.Data.Rows

const TYPE_MAPPINGS = Dict{Symbol,Symbol}( # Julia => MySQL
  :char       => :CHARACTER,
  :string     => :VARCHAR,
  :text       => :TEXT,
  :integer    => :INTEGER,
  :int        => :INTEGER,
  :float      => :FLOAT,
  :decimal    => :DECIMAL,
  :datetime   => :DATETIME,
  :timestamp  => :TIMESTAMP,
  :time       => :TIME,
  :date       => :DATE,
  :binary     => :BLOB,
  :boolean    => :BOOLEAN,
  :bool       => :BOOLEAN
)

function db_adapter()::Symbol
  Symbol(DB_ADAPTER)
end

#
# Connection
#

"""
    connect(conn_data::Dict)::DatabaseHandle

Connects to the database and returns a handle.
"""
function connect(conn_data::Dict) :: DatabaseHandle
    MySQL.connect(conn_data["host"], conn_data["username"], get(conn_data, "password", "");
                  db    = conn_data["database"],
                  port  = get(conn_data, "port", 3306),
                  opts  = Dict(getfield(MySQL.API, Symbol(k))=>v for (k,v) in conn_data["options"]))
end

"""
    disconnect(conn::DatabaseHandle)::Nothing

Disconnects from database.
"""
function disconnect(conn::DatabaseHandle) :: Nothing
  MySQL.disconnect(conn)
end

#
# Utility
#

"""
    table_columns_sql(table_name::String)::String

Returns the adapter specific query for SELECTing table columns information corresponding to `table_name`.
"""
function table_columns_sql(table_name::String) :: String
  "SHOW COLUMNS FROM `$table_name`"
end

"""
    create_migrations_table(table_name::String)::Bool

Runs a SQL DB query that creates the table `table_name` with the structure needed to be used as the DB migrations table.
The table should contain one column, `version`, unique, as a string of maximum 30 chars long.
Returns `true` on success.
"""
function create_migrations_table(table_name::String)::Bool
  "CREATE TABLE `$table_name` (
    `version` varchar(30) NOT NULL DEFAULT '',
    PRIMARY KEY (`version`)
  ) ENGINE=InnoDB DEFAULT CHARSET=utf8" |> SearchLight.Database.query

  @info "Created table $table_name"

  true
end

#
# Data sanitization
#

"""
    escape_column_name(c::String, conn::DatabaseHandle)::String

Escapes the column name.

# Examples
```julia
julia>
```
"""
function escape_column_name(c::String, conn::DatabaseHandle) :: String
  """`$(replace(c, "`"=>"-"))`"""
end

"""
    escape_value{T}(v::T, conn::DatabaseHandle)::T

Escapes the value `v` using native features provided by the database backend.

# Examples
```julia
julia>
```
"""
function escape_value(v::T, conn::DatabaseHandle)::T where {T}
  isa(v, Number) ? v : "'$(MySQL.escape(conn, string(v)))'"
end

function mysql_escape(s::String, conn::DatabaseHandle)
  MySQL.escape(conn, s)
end

#
# Query execution
#

"""
    query(sql::String, suppress_output::Bool, conn::DatabaseHandle)::DataFrames.DataFrame

Executes the `sql` query against the database backend and returns a DataFrame result.

# Examples:
```julia
julia> query(SearchLight.to_fetch_sql(Article, SQLQuery(limit = 5)), false, Database.connection)

2017-01-16T21:36:21.566 - info: SQL QUERY: SELECT \"articles\".\"id\" AS \"articles_id\", \"articles\".\"title\" AS \"articles_title\", \"articles\".\"summary\" AS \"articles_summary\", \"articles\".\"content\" AS \"articles_content\", \"articles\".\"updated_at\" AS \"articles_updated_at\", \"articles\".\"published_at\" AS \"articles_published_at\", \"articles\".\"slug\" AS \"articles_slug\" FROM \"articles\" LIMIT 5

  0.000985 seconds (16 allocations: 576 bytes)

5×7 DataFrames.DataFrame
...
```
"""
function query(sql::String, suppress_output::Bool, conn::DatabaseHandle) :: DataFrames.DataFrame
  try
    _result = if suppress_output || ( ! SearchLight.config.log_db && ! SearchLight.config.log_queries )
      DB_ADAPTER.Query(conn, sql)
    else
      @info sql
      @time DB_ADAPTER.Query(conn, sql)
    end

    result = if startswith(sql, "INSERT ")
      DataFrames.DataFrame(SearchLight.LAST_INSERT_ID_LABEL => last_insert_id(conn))
    elseif startswith(sql, "ALTER ") || startswith(sql, "CREATE ") || startswith(sql, "DROP ") || startswith(sql, "DELETE ") || startswith(sql, "UPDATE ")
      DataFrames.DataFrame(:result => "OK")
    else
      DataFrames.DataFrame(_result)
    end
  catch ex
    @error """MySQL error when running "$(escape_string(sql))" """

    if isa(ex, MySQL.MySQLInternalError) && (ex.errno == 2013 || ex.errno == 2006)
      @warn ex, " ...Attempting reconnection"
      query(sql, suppress_output, SearchLight.Database.connect())
    else
      rethrow(ex)
    end
  end
end

"""
"""
function to_find_sql(m::Type{T}, q::SearchLight.SQLQuery, joins::Union{Nothing,Vector{SearchLight.SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  sql::String = ( "$(to_select_part(m, q.columns, joins)) $(to_from_part(m)) $(to_join_part(m, joins)) $(to_where_part(q.where)) " *
                      "$(to_group_part(q.group)) $(to_having_part(q.having)) $(to_order_part(m, q.order)) " *
                      "$(to_limit_part(q.limit)) $(to_offset_part(q.offset))") |> strip
  replace(sql, r"\s+"=>" ")
end

const to_fetch_sql = to_find_sql

"""
"""
function to_store_sql(m::T; conflict_strategy = :error)::String where {T<:SearchLight.AbstractModel}
  uf = SearchLight.persistable_fields(m)

  sql = if ! SearchLight.ispersisted(m) || (SearchLight.ispersisted(m) && conflict_strategy == :update)
    pos = findfirst(x -> x == SearchLight.primary_key_name(m), uf)
    pos > 0 && splice!(uf, pos)

    fields = SearchLight.SQLColumn(uf)
    vals = join( map(x -> string(SearchLight.to_sqlinput(m, Symbol(x), getfield(m, Symbol(x)))), uf), ", ")

    "INSERT INTO $(m._table_name) ( $fields ) VALUES ( $vals )" *
        if ( conflict_strategy == :error ) ""
        elseif ( conflict_strategy == :ignore ) " ON DUPLICATE KEY UPDATE `$(SearchLight.primary_key_name(m))` = `$(SearchLight.primary_key_name(m))`"
        elseif ( conflict_strategy == :update && getfield(m, Symbol(SearchLight.primary_key_name(m))).value !== nothing )
          " ON DUPLICATE KEY UPDATE $(update_query_part(m))"
        else ""
        end
  else
    "UPDATE $(m._table_name) SET $(update_query_part(m))"
  end

  sql
end

"""
"""
function delete_all(m::Type{T}; truncate::Bool = true, reset_sequence::Bool = true, cascade::Bool = false)::Nothing where {T<:SearchLight.AbstractModel}
  _m::T = m()
  (truncate ? "TRUNCATE $(_m._table_name)" : "DELETE FROM $(_m._table_name)") |> SearchLight.query

  nothing
end

"""
"""
function delete(m::T)::T where {T<:SearchLight.AbstractModel}
  SearchLight.ispersisted(m) || throw(SearchLight.Exceptions.NotPersistedException(m))

  "DELETE FROM $(m._table_name) WHERE $(SearchLight.primary_key_name(m)) = '$(m.id.value)'" |> SearchLight.query

  m.id = SearchLight.DbId()

  m
end

"""
"""
function count(m::Type{T}, q::SearchLight.SQLQuery = SearchLight.SQLQuery())::Int where {T<:SearchLight.AbstractModel}
  count_column = SearchLight.SQLColumn("COUNT(*) AS __cid", raw = true)
  q = SearchLight.clone(q, :columns, push!(q.columns, count_column))

  SearchLight.finddf(m, q)[1, Symbol("__cid")]
end

"""
"""
function update_query_part(m::T)::String where {T<:SearchLight.AbstractModel}
  update_values = join(map(x -> "$(string(SearchLight.SQLColumn(x))) = $(string(SearchLight.to_sqlinput(m, Symbol(x), getfield(m, Symbol(x)))))", SearchLight.persistable_fields(m)), ", ")

  " $update_values WHERE $(m._table_name).$(SearchLight.primary_key_name(m)) = '$(m.id.value)'"
end

"""
"""
function column_data_to_column_name(column::SearchLight.SQLColumn, column_data::Dict{Symbol,Any}) :: String
  "$(SearchLight.to_fully_qualified(column_data[:column_name], column_data[:table_name])) AS $(isempty(column_data[:alias]) ? SearchLight.to_sql_column_name(column_data[:column_name], column_data[:table_name]) : column_data[:alias])"
end

"""
"""
function to_select_part(m::Type{T}, cols::Vector{SearchLight.SQLColumn}, joins::Union{Nothing,Vector{SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  "SELECT " * SearchLight.Database._to_select_part(m, cols, joins)
end

"""
"""
function to_from_part(m::Type{T})::String where {T<:SearchLight.AbstractModel}
  "FROM " * SearchLight.Database.escape_column_name(m()._table_name)
end

function to_where_part(w::Vector{SearchLight.SQLWhereEntity}) :: String
  where = isempty(w) ?
          "" :
          "WHERE " * (string(first(w).condition) == "AND" ? "TRUE " : "FALSE ") * join(map(wx -> string(wx), w), " ")

  replace(where, r"WHERE TRUE AND "i => "WHERE ")
end

"""
"""
function to_order_part(m::Type{T}, o::Vector{SearchLight.SQLOrder})::String where {T<:SearchLight.AbstractModel}
  isempty(o) ?
    "" :
    "ORDER BY " * join(map(x -> (! SearchLight.is_fully_qualified(x.column.value) ? SearchLight.to_fully_qualified(m, x.column) : x.column.value) * " " * x.direction, o), ", ")
end

"""
"""
function to_group_part(g::Vector{SearchLight.SQLColumn}) :: String
  isempty(g) ?
    "" :
    " GROUP BY " * join(map(x -> string(x), g), ", ")
end

"""
"""
function to_limit_part(l::SearchLight.SQLLimit) :: String
  l.value != "ALL" ? "LIMIT " * (l |> string) : ""
end

"""
"""
function to_offset_part(o::Int) :: String
  o != 0 ? "OFFSET " * (o |> string) : ""
end

"""
"""
function to_having_part(h::Vector{SearchLight.SQLWhereEntity}) :: String
  having =  isempty(h) ?
            "" :
            "HAVING " * (string(first(h).condition) == "AND" ? "TRUE " : "FALSE ") * join(map(w -> string(w), h), " ")

  replace(having, r"HAVING TRUE AND "i => "HAVING ")
end

"""
"""
function to_join_part(m::Type{T}, joins::Union{Nothing,Vector{SQLJoin{N}}} = nothing)::String where {T<:SearchLight.AbstractModel, N<:Union{Nothing,SearchLight.AbstractModel}}
  joins === nothing && return ""

  _m::T = m()
  join_part = ""

  join_part * join( map(x -> string(x), joins), " " )
end

"""
    cast_type(v::Bool)::Union{Bool,Int,Char,String}

Converts the Julia type to the corresponding type in the database.
"""
function cast_type(v::Bool)::Int
  v ? 1 : 0
end

"""
"""
function create_table_sql(f::Function, name::String, options::String = "") :: String
  "CREATE TABLE `$name` (" * join(f()::Vector{String}, ", ") * ") $options" |> strip
end

"""
"""
function column_sql(name::String, column_type::Symbol, options::String = ""; default::Any = nothing, limit::Union{Int,Nothing,String} = nothing, not_null::Bool = false)::String
  "`$name` $(TYPE_MAPPINGS[column_type] |> string) " *
    (isa(limit, Int) || isa(limit, String) ? "($limit)" : "") *
    (default === nothing ? "" : " DEFAULT $default ") *
    (not_null ? " NOT NULL " : "") *
    " " * options
end

"""
"""
function column_id_sql(name::String = "id", options::String = ""; constraint::String = "", nextval::String = "") :: String
  "`$name` INT NOT NULL AUTO_INCREMENT PRIMARY KEY $options"
end

"""
"""
function add_index_sql(table_name::String, column_name::String; name::String = "", unique::Bool = false, order::Symbol = :none) :: String
  name = isempty(name) ? SearchLight.Database.index_name(table_name, column_name) : name
  "CREATE $(unique ? "UNIQUE" : "") INDEX `$(name)` ON `$table_name` (`$column_name`)"
end

"""
"""
function add_column_sql(table_name::String, name::String, column_type::Symbol; default::Any = nothing, limit::Union{Int,Nothing} = nothing, not_null::Bool = false) :: String
  "ALTER TABLE `$table_name` ADD $(column_sql(name, column_type, default = default, limit = limit, not_null = not_null))"
end

"""
"""
function drop_table_sql(name::String) :: String
  "DROP TABLE `$name`"
end

"""
"""
function remove_column_sql(table_name::String, name::String, options::String = "") :: String
  "ALTER TABLE `$table_name` DROP COLUMN `$name` $options"
end

"""
"""
function remove_index_sql(table_name::String, name::String, options::String = "") :: String
  "DROP INDEX `$name` ON `$table_name` $options"
end

"""
"""
function rand(m::Type{T}; limit = 1)::Vector{T} where {T<:SearchLight.AbstractModel}
  SearchLight.find(m, SearchLight.SQLQuery(limit = SearchLight.SQLLimit(limit), order = [SearchLight.SQLOrder("rand()", raw = true)]))
end

"""
"""
function last_insert_id(conn)
  MySQL.insertid(conn)
end

end


text:
function sparse_regression(X, Theta, lambda; maxiters=20)

    N = size(X)[2]

    # initial guess
    C = Theta \ X

    # find coefficients and apply threshold
    i = abs.(C) .< lambda
    j = Vector{Bool}(1 .- i)
    C[i] .= 0
    sparsity = sum(i[:])

    for _=1:maxiters

        for k=1:N
            if all(i[:,k])
                C[:,k] .= 0
            else
                C[j[:,k],k] = Theta[:,j[:,k]] \ X[:,k]
            end
        end

        # find coefficients and apply threshold
        i = abs.(C) .< lambda
        j = Vector{Bool}(1 .- i)
        C[i] .= 0

        # if the sparsity hasn't changed, stop iteration
        new_sparsity = sum(i[:])
        if new_sparsity == sparsity
            break
        end

        sparisty = new_sparsity
    end

    return C
end


text:
using CuArrays, CUDAnative, CUDAdrv

function add!(a,b,c)
    local tid = 1
    while (tid <= min(length(a), length(b), length(c)))
        c[tid] = a[tid] + b[tid]
        tid += 1
    end
    return nothing
end

function main()

    a = CuArrays.CuArray(fill(0, 10))
    b = CuArrays.CuArray(fill(0, 10))
    c = CuArrays.CuArray(fill(0, 10))

    for i in 1:10
        a[i] = -i
        b[i] = i * i
    end

    # IMPORTANT NOTE TO SELF
    # you can pass tuples to represent a grid to blocks, just like in CUDA C <3
    @cuda blocks=10 add!(a,b,c)

    a=Array(a)
    b=Array(b)
    c=Array(c)

    for i in 1:length(a)
        println(a[i], " + ", b[i], " = ", c[i])
    end
end

main()


text:
using JuLIP, ForwardDiff, StaticArrays, BenchmarkTools

# simple EAM potential
f(R) = sqrt( 1.0 + sum( exp(-norm(r)) for r in R ) )

nneigs = 30
R0 = [ @SVector rand(3) for n = 1:nneigs ]

f_fd(R) =  ForwardDiff.gradient( T -> f(vecs(T)),  mat(R)[:] ) |> vecs

function f_d(R::AbstractVector{<:JVec})
   ∇f = zeros(JVecF, length(R))
   ρ̄ = sum( exp(-norm(r)) for r in R )
   dF = 0.5 * (1.0 + ρ̄)^(-0.5)
   return [ dF * (- exp(-norm(r))) * r/norm(r)  for r in R ]
end

@assert f_d(R0) ≈ f_fd(R0)

@btime f($R0);
@btime f_d($R0);
@btime f_fd($R0);

# f1(R::AbstractMatrix) =
#       sqrt( 1.0 + sum( exp.(-sqrt.(sum(abs2, R, 1))) ) )
# f1_rd(R::Matrix) =
#       reshape(
#             ReverseDiff.gradient( S -> f1(reshape(S, 3, length(S) ÷ 3)), R[:] ),
#             size(R) )
# R1 = mat(R0)
# f1_rd(R1)
#
# @btime f1($R1);
# @btime f1_rd($R1);


text:
using BenchmarkTools: @btime
import Random
using Distributions: AbstractMixtureModel, MixtureModel, LogNormal, Normal, pdf, ncomponents, probs, component, components, ContinuousUnivariateDistribution
using Test

# v0.22.1
function current_master(d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    v = 0.0
    @inbounds for i in eachindex(p)
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            v += pdf(c, x) * pi
        end
    end
    return v
end

# compute the overhead of having a mixture
function evaluate_manual_pdf(distributions, priors, x)
    return sum(pdf(d, x) * p for (d, p) in zip(distributions, priors))
end

function improved_version(d, x)
    p = probs(d)
    return sum(enumerate(p)) do (i, pi)
        if pi > 0
            @inbounds c = component(d, i)
            pdf(c, x) * pi
        else
            zero(eltype(p))
        end
    end
end

function improved_no_inbound(d, x)
    p = probs(d)
    return sum(enumerate(p)) do (i, pi)
        if pi > 0
           c = component(d, i)
           pdf(c, x) * pi
        else
            zero(eltype(p))
        end
    end
end

function ifelse_version(d, x)
    ps = probs(d)
    cs = components(d)
    P = eltype(ps)
    return sum(ifelse(ps[i] > 0, ps[i] * pdf(cs[i], x), zero(P)) for i in eachindex(ps))
end

function forloop(d, x)
    ps = probs(d)
    cs = components(d)
    s = zero(eltype(ps))
    @inbounds for i in eachindex(ps)
        if ps[i] > 0
            s += ps[i] * pdf(cs[i], x)
        end
    end
    return s
end

function indexed_sum_comp(d, x)
    ps = probs(d)
    cs = components(d)
    @inbounds sum(ps[i] * pdf(cs[i], x) for i in eachindex(ps) if ps[i] > 0)
end

function indexed_boolprod(d, x)
    ps = probs(d)
    cs = components(d)
    @inbounds sum((ps[i] > 0) * (ps[i] * pdf(cs[i], x)) for i in eachindex(ps))
end

function indexed_boolprod_noinbound(d, x)
    ps = probs(d)
    cs = components(d)
    return sum((ps[i] > 0) * (ps[i] * pdf(cs[i], x)) for i in eachindex(ps))
end

function sumcomp_cond(d, x)
    ps = probs(d)
    cs = components(d)
    s = zero(eltype(ps))
    @inbounds sum(ps[i] * pdf(cs[i], x) for i in eachindex(ps) if ps[i] > 0)
end

distributions = [
    Normal(-1.0, 0.3),
    Normal(0.0, 0.5),
    Normal(3.0, 1.0),
]

priors = [0.25, 0.25, 0.5]

gmm_normal = MixtureModel(distributions, priors)

@info "Small Gaussian mixture"

Random.seed!(42)
for x in rand(5)
    @info "sampling $x"
    @info "evaluate_manual_pdf"
    vman = @btime evaluate_manual_pdf($distributions, $priors, $x)
    @info "current_master"
    vmaster =  @btime current_master($gmm_normal, $x)
    @info "improved_version"
    v1 = @btime improved_version($gmm_normal, $x)
    @info "improved_no_inbound"
    v2 = @btime improved_no_inbound($gmm_normal, $x)
    @info "ifelse_version"
    v3 = @btime ifelse_version($gmm_normal, $x)
    @info "forloop"
    v4 = @btime forloop($gmm_normal, $x)
    @info "indexed_sum_comp"
    v5 = @btime indexed_sum_comp($gmm_normal, $x)
    @info "indexed_boolprod"
    v6 = @btime indexed_boolprod($gmm_normal, $x)
    @info "indexed_boolprod_noinbound"
    v6 = @btime indexed_boolprod_noinbound($gmm_normal, $x)
    @info "sumcomp_cond"
    v7 = @btime sumcomp_cond($gmm_normal, $x)
    @assert vman ≈ vmaster ≈ v1 ≈ v2 ≈ v3 ≈ v4 ≈ v5 ≈ v6 ≈ v7
    @info "==================="
end

large_normals = [Normal(rand(), rand()) for _ in 1:1000]
large_probs = [rand() for _ in 1:1000]
large_probs .= large_probs ./ sum(large_probs)

gmm_large = MixtureModel(large_normals, large_probs)

@info "Large Gausian mixture"

Random.seed!(42)
for x in rand(5)
    @info "sampling $x"
    @info "evaluate_manual_pdf"
    vman = @btime evaluate_manual_pdf($large_normals, $large_probs, $x)
    @info "current_master"
    vmaster =  @btime current_master($gmm_large, $x)
    @info "improved_version"
    v1 = @btime improved_version($gmm_large, $x)
    @info "improved_no_inbound"
    v2 = @btime improved_no_inbound($gmm_large, $x)
    @info "ifelse_version"
    v3 = @btime ifelse_version($gmm_large, $x)
    @info "forloop"
    v4 = @btime forloop($gmm_large, $x)
    @info "indexed_sum_comp"
    v5 = @btime indexed_sum_comp($gmm_large, $x)
    @info "indexed_boolprod"
    v6 = @btime indexed_boolprod($gmm_large, $x)
    @info "indexed_boolprod_noinbound"
    v6 = @btime indexed_boolprod_noinbound($gmm_large, $x)
    @info "sumcomp_cond"
    v7 = @btime sumcomp_cond($gmm_large, $x)
    @assert vman ≈ vmaster ≈ v1 ≈ v2 ≈ v3 ≈ v4 ≈ v5 ≈ v6 ≈ v7
    @info "==================="
end

large_het = append!(
    ContinuousUnivariateDistribution[Normal(rand(), rand()) for _ in 1:1000],
    ContinuousUnivariateDistribution[LogNormal(rand(), rand()) for _ in 1:1000],
)

large_het_probs = [rand() for _ in 1:2000]
large_het_probs .= large_het_probs ./ sum(large_het_probs)

gmm_het = MixtureModel(large_het, large_het_probs)

@info "Heterogeneous distributions"

Random.seed!(42)
for x in rand(5)
    @info "sampling $x"
    @info "evaluate_manual_pdf"
    vman = @btime evaluate_manual_pdf($large_het, $large_het_probs, $x)
    @info "current_master"
    vmaster =  @btime current_master($gmm_het, $x)
    @info "improved_version"
    v1 = @btime improved_version($gmm_het, $x)
    @info "improved_no_inbound"
    v2 = @btime improved_no_inbound($gmm_het, $x)
    @info "ifelse_version"
    v3 = @btime ifelse_version($gmm_het, $x)
    @info "forloop"
    v4 = @btime forloop($gmm_het, $x)
    @info "indexed_sum_comp"
    v5 = @btime indexed_sum_comp($gmm_het, $x)
    @info "indexed_boolprod"
    v6 = @btime indexed_boolprod($gmm_het, $x)
    @info "indexed_boolprod_noinbound"
    v6 = @btime indexed_boolprod_noinbound($gmm_het, $x)
    @info "sumcomp_cond"
    v7 = @btime sumcomp_cond($gmm_het, $x)
    @assert vman ≈ vmaster ≈ v1 ≈ v2 ≈ v3 ≈ v4 ≈ v5 ≈ v6 ≈ v7
    @info "==================="
end

@info "Test with one NaN"

distributions = [
    Normal(-1.0, 0.3),
    Normal(0.0, 0.5),
    Normal(3.0, 1.0),
    Normal(NaN, 1.0),
]

priors = [0.25, 0.25, 0.5, 0.0]

gmm_normal = MixtureModel(distributions, priors)

Random.seed!(42)
for x in rand(1)
    @info "sampling $x"
    @info "current_master"
    vmaster =  @btime current_master($gmm_normal, $x)
    @info "improved_version"
    v1 = @btime improved_version($gmm_normal, $x)
    @info "improved_no_inbound"
    v2 = @btime improved_no_inbound($gmm_normal, $x)
    @info "ifelse_version"
    v3 = @btime ifelse_version($gmm_normal, $x)
    @info "forloop"
    v4 = @btime forloop($gmm_normal, $x)
    @info "indexed_sum_comp"
    v5 = @btime indexed_sum_comp($gmm_normal, $x)
    @info "indexed_boolprod"
    v6 = @btime indexed_boolprod($gmm_normal, $x)
    @info "indexed_boolprod_noinbound"
    v7 = @btime indexed_boolprod_noinbound($gmm_normal, $x)
    @info "sumcomp_cond"
    v8 = @btime sumcomp_cond($gmm_normal, $x)
    @test vmaster ≈ v1 ≈ v2 ≈ v3 ≈ v4 ≈ v5 ≈ v6 ≈ v7 ≈ v8
end


text:
include("../imports.jl")

@testset ExtendedTestSet "Affine" begin
    P = CoordinateTransformations.Translation(10, 10)
    tfm = Project(P)

    keypoints = Keypoints([SVector(1., 1)], (50, 50))
    image = Image(rand(RGB, 50, 50))

    indices = boundsranges(getbounds(keypoints))

    @testset ExtendedTestSet "Affine Keypoints" begin
        @test getprojection(tfm, getbounds(keypoints)) == P
        @test_nowarn project(P, keypoints, indices)
        tkeypoints = apply(tfm, keypoints)
        @test itemdata(tkeypoints) == [SVector(11, 11)]
        @test getbounds(tkeypoints)[1] == SVector(10, 10)
    end

    @testset ExtendedTestSet "Affine Image" begin
        @test_nowarn project(P, image, indices)
        timage = apply(tfm, image)
        @test itemdata(timage)[11, 11] ≈ itemdata(image)[1, 1]
        @test getbounds(timage)[1] == SVector(10, 10)
    end

    tfmcropped = tfm |> CenterCrop((50, 50))

    @testset ExtendedTestSet "Cropped, Image" begin
        tfmcropped = tfm |> CenterCrop((50, 50))
        @test_nowarn apply(tfmcropped, keypoints)
        tkeypoints = apply(tfmcropped, keypoints)
        @test getbounds(tkeypoints)[1] == SVector(10, 10)
        @test itemdata(tkeypoints)[1] == SVector(11, 11)
    end

    @testset ExtendedTestSet "Cropped with PinOrigin, Image" begin
        tfmcropped = tfm |> PinOrigin() |> CenterCrop((50, 50))
        @test_nowarn apply(tfmcropped, keypoints)
        tkeypoints = apply(tfmcropped, keypoints)
        @test getbounds(tkeypoints)[1] == SVector(0, 0)
    end

    @testset ExtendedTestSet "Composed" begin
        tfms = Project(Translation(10, 10)) |> Project(Translation(-10, -10))
        randstate = getrandstate(tfms)
        @test getprojection(tfms, getbounds(image); randstate) == Translation(0, 0)
        @test itemdata(apply(tfms, image)) == itemdata(image)
        @test itemdata(apply(tfms, keypoints)) == itemdata(keypoints)
    end

    @testset ExtendedTestSet "ScaleKeepAspect" begin
        image = Image(zeros(RGB, 100, 50))
        keypoints = Keypoints([SVector(0., 0)], (50, 50))

        imbounds = getbounds(image)
        kpbounds = getbounds(keypoints)
        tfm1 = ScaleKeepAspect((50, 50))
        @test getprojection(tfm1, imbounds) ≈ getprojection(ScaleRatio((1, 1)), imbounds)
        @test getprojection(tfm1, kpbounds) ≈ getprojection(ScaleRatio((1, 1)), kpbounds,)

        tfm2 = ScaleKeepAspect((25, 25))
        @test getprojection(tfm2, imbounds) ≈ getprojection(ScaleRatio((1 / 2, 1 / 2)), imbounds)
        @test getprojection(tfm2, kpbounds) ≈ getprojection(ScaleRatio((1 / 2, 1 / 2)), kpbounds)
    end

    @testset ExtendedTestSet "Scale" begin
        keypoints = Keypoints([SVector(0., 0)], (50, 50))
        image = Image(rand(RGB, 50, 50))

        @testset ExtendedTestSet "ScaleFixed" begin
            tfm = ScaleFixed((25, 25))
            @test_nowarn apply(tfm, image)
            @test_nowarn apply(tfm, keypoints)
            timage = apply(tfm, image)
            tkeypoints = apply(tfm, keypoints)
            @test boundsranges(getbounds(timage)) == (0:25, 0:25)
            @test getbounds(timage) == getbounds(tkeypoints)

        end
        @testset ExtendedTestSet "ScaleRatio" begin
            tfm = ScaleRatio((1/2, 1/2))
            @test_nowarn apply(tfm, image)
            @test_nowarn apply(tfm, keypoints)
            timage = apply(tfm, image)
            tkeypoints = apply(tfm, keypoints)
            @test boundsranges(getbounds(timage)) == (0:25, 0:25)
            @test getbounds(timage) == getbounds(tkeypoints)
        end
    end

    @testset ExtendedTestSet "Rotate" begin
        tfm = Rotate(10)
        image = Image(rand(RGB, 50, 50))
        @test_nowarn apply(tfm, image)

    end

    @testset ExtendedTestSet "Reflect" begin
        tfm = Reflect(10)
        image = Image(rand(RGB, 50, 50))
        @test_nowarn apply(tfm, image)

    end

    @testset ExtendedTestSet "CenterResizeCrop" begin
        tfm = CenterResizeCrop((25, 40))
        keypoints = Keypoints([SVector(0., 0)], (50, 50))
        image = Image(rand(RGB, 50, 50))
        timage, tkeypoints = apply(tfm, (image, keypoints))

        @test getbounds(timage) ≈ getbounds(tkeypoints)
        @test boundsranges(getbounds(timage)) == (1:25, 1:40)
    end

    @testset ExtendedTestSet "CroppedAffine inplace" begin
        image1 = Image(rand(RGB, 100, 100))
        image2 = Image(rand(RGB, 100, 100))
        tfm = ScaleFixed((50, 50))
        buffer = makebuffer(tfm, image1)

        @test_nowarn apply!(buffer, tfm, image2)
    end
end

@testset ExtendedTestSet "Big pipeline" begin
    @testset ExtendedTestSet "2D" begin
        sz = (100, 100)
        items = (
            Image(rand(RGB, sz)),
            Keypoints(rand(SVector{2, Float32}, 50), sz),
            MaskBinary(rand(Bool, sz)),
            MaskMulti(rand(1:8, sz)),
        )

        tfms = compose(
            Rotate(10),
            FlipX(), FlipY(),
            ScaleRatio((.8, .8)),
            RandomResizeCrop((50, 50)),
        )
        @test_nowarn apply(tfms, items)
    end

    @testset ExtendedTestSet "3D" begin
        sz = (50, 50, 50)
        items = (
            Image(rand(RGB, sz)),
            Keypoints(rand(SVector{3, Float32}, 50), sz),
            MaskBinary(rand(Bool, sz)),
            MaskMulti(rand(1:8, sz)),
        )

        tfms = compose(
            ScaleRatio((.8, .8, .8)),
            RandomResizeCrop((25, 25, 25)),
        )
        @test_nowarn apply(tfms, items)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4655bdc2-1ed7-11ec-022d-7df0b8fdb907
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 661cc8fe-44f3-422a-91f5-96956b748847
begin
	using SqState
	using Plots
	using QuantumStateBase
	using Flux
	using LinearAlgebra
end

# ╔═╡ c55d7af1-0cd3-4b0a-bbd5-0048f15ea851
m = SqState.get_model("model_q2ρ1_70")

# ╔═╡ 8f82e42f-8c2c-4bd0-ab20-95b0930670dc
# args = SqState.rand_arg((0, 1), (0, 2π), (0, 1))
args = (0.9,  4.17646, 0.5, 0.3, 0.852858, 0.147142);

# ╔═╡ f0cfd6f3-1e55-46cf-bd58-937e52cb7daa
begin
	dim = 70
	state = SqState.construct_state_sqth_th(args..., 1000)
	d = Float32.(rand(state, 4096, IsGaussian))
	ρ = m(reshape(d[1, :], :, 1, 1))
	ρ = reshape(ρ[:, 1, 1] + im * ρ[:, 2, 1], dim, dim)
end;

# ╔═╡ 08ddc8c4-57c7-4b2f-baeb-c424c3227065
show_dim = 35

# ╔═╡ e9c31f3c-4f06-43b4-9982-4f69e997cde8
begin
	m1 = maximum(real.(state.𝛒))
	heatmap(real.(state.𝛒)[1:show_dim, 1:show_dim], color=:coolwarm, clim=(-m1, m1))
end

# ╔═╡ 6c0a3a23-335c-4661-8a2f-3729b3f0a1de
begin
	m2 = maximum(real.(ρ))
	heatmap(real.(ρ)[1:show_dim, 1:show_dim], color=:coolwarm, clim=(-m2, m2))
end

# ╔═╡ 6a68f473-cb16-4635-a3b0-9aae0e18d86d
function fidelity(ρ1, ρ2)
	sqrt_ρ1 = sqrt(ρ1)
	return tr(sqrt(sqrt_ρ1*ρ2*sqrt_ρ1))^2
end

# ╔═╡ 00f15b97-52e0-43aa-a2ee-136c7a8a6e95
fidelity(ρ, state.𝛒[1:dim, 1:dim])

# ╔═╡ 35a5ee23-a2df-47c6-a096-a5e1c60328b6
# begin
# 	f_dim = 70

# 	f = 0
# 	for _ in 1:50
# 		f_state = SqState.construct_state(
# 			SqState.rand_arg((0, 1), (0, 2π), (0, 1))...,
# 			1000
# 		)
# 		f_ρ = f_state.𝛒[1:f_dim, 1:f_dim]

# 		f_ρ̂ = m(reshape(Float32.(rand(state, 4096, IsGaussian))[1, :], :, 1, 1))
# 		f_ρ̂ = reshape(f_ρ̂[:, 1, 1] + im * f_ρ̂[:, 2, 1], f_dim, f_dim)
# 		f += fidelity(f_ρ̂, f_ρ)
# 	end
# end

# ╔═╡ 0171f5c1-2b14-4694-b8b0-14e50591db4b
# f/50

# ╔═╡ Cell order:
# ╠═4655bdc2-1ed7-11ec-022d-7df0b8fdb907
# ╠═661cc8fe-44f3-422a-91f5-96956b748847
# ╠═c55d7af1-0cd3-4b0a-bbd5-0048f15ea851
# ╠═8f82e42f-8c2c-4bd0-ab20-95b0930670dc
# ╠═f0cfd6f3-1e55-46cf-bd58-937e52cb7daa
# ╠═08ddc8c4-57c7-4b2f-baeb-c424c3227065
# ╟─e9c31f3c-4f06-43b4-9982-4f69e997cde8
# ╟─6c0a3a23-335c-4661-8a2f-3729b3f0a1de
# ╟─6a68f473-cb16-4635-a3b0-9aae0e18d86d
# ╠═00f15b97-52e0-43aa-a2ee-136c7a8a6e95
# ╟─35a5ee23-a2df-47c6-a096-a5e1c60328b6
# ╟─0171f5c1-2b14-4694-b8b0-14e50591db4b


text:
"""
    Entry{ET}

    value: data
    isinverted: flag indicating if matrix value has been inverted
"""
mutable struct Entry{ET}
    value::ET
    isinverted::Bool

    function Entry{T}(dims...;
        static = true) where T
        static ? value = szeros(T,dims...) : value = zeros(T,dims...)
        new{typeof(value)}(value, false)
    end

    Entry(dims...; static = true) = Entry{Float64}(dims...; static = static)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, entry::Entry)
    summary(io, entry)
    println(io, "")
    println(io, "value:      "*string(entry.value))
    println(io, "isinverted: "*string(entry.isinverted))
end

function Base.zero(::Entry{ET}) where ET
    dims = [ET.parameters[1].parameters...]
    return Entry{ET.parameters[2]}(dims...)
end

Base.zero(::Type{Dojo.Entry}) = nothing