text:
using HMatrices
using Test
using StaticArrays
using LinearAlgebra
using HMatrices: ACA, PartialACA, TSVD, RkMatrix

include(joinpath(HMatrices.PROJECT_ROOT,"test","testutils.jl"))

using Random
Random.seed!(1)

@testset "Scalar" begin
    T = ComplexF64
    m,n = 100,100
    X = rand(SVector{3,Float64},m)
    Y = map(i->SVector(10,0,0) + rand(SVector{3,Float64}),1:n)
    K = helmholtz_matrix(X,Y,1.0)
    M = Matrix(K)
    irange,jrange = 1:m,1:n
    @testset "aca_full" begin
        atol = 1e-5
        aca = ACA(atol=atol)
        R = aca(K,irange,jrange)
        @test norm(Matrix(R) - M) < atol
        rtol = 1e-5
        aca = ACA(rtol=rtol)
        R   = aca(M,irange,jrange)
        norm(Matrix(R) - M)
        @test norm(Matrix(R) - M) < rtol*norm(M)
        r = 10
        aca = ACA(rank=r)
        R = aca(M,irange,jrange)
        @test rank(R) == r
    end
    @testset "aca_partial" begin
        atol = 1e-5
        aca = PartialACA(atol=atol)
        R = aca(M,irange,jrange)
        @test norm(Matrix(R) - M) < atol
        rtol = 1e-5
        aca = PartialACA(rtol=rtol)
        R = aca(M,irange,jrange)
        @test norm(Matrix(R) - M) < rtol*norm(M)
        r = 10
        aca = PartialACA(rank=r)
        R = aca(M,irange,jrange)
        @test rank(R) == r

        # test fast update of frobenius norm
        m,n = 10000,1000
        r = 10
        T = ComplexF64
        A = [rand(T,m) for _ in 1:r]
        B = [rand(T,n) for _ in 1:r]
        R = RkMatrix(A,B)
        old_norm = norm(Matrix(R),2)
        push!(A,rand(T,m))
        push!(B,rand(T,n))
        Rnew = RkMatrix(A,B)
        new_norm = norm(Matrix(Rnew),2)
        @test new_norm ≈ HMatrices._update_frob_norm(old_norm,A,B)

        # test simple case where things are not compressible
        A  = rand(2,2)
        comp = PartialACA(rtol=1e-5)
        @test comp(A,1:2,1:2) ≈ A
    end
    @testset "truncated svd" begin
        atol = 1e-5
        tsvd = TSVD(atol=atol)
        R    = tsvd(M,irange,jrange)
        # the inequality below is guaranteed to be true  for  the spectral norm
        # i.e. (the `opnorm` with `p=2`).
        @test opnorm(Matrix(R) - M) < atol
        rtol = 1e-5
        tsvd = TSVD(rtol=rtol)
        R    = tsvd(M,irange,jrange)
        @test opnorm(Matrix(R) - M) < rtol*opnorm(M)
        r = 10
        tsvd = TSVD(rank=r)
        R  = tsvd(M,irange,jrange)
        @test rank(R) == r
    end
end

@testset "Tensorial" begin
    T = SMatrix{3,3,ComplexF64,9}
    # T = SMatrix{3,3,Float64,9}
    m,n = 100,100
    X = rand(SVector{3,Float64},m)
    Y = map(i->SVector(10,0,0) + rand(SVector{3,Float64}),1:n)
    K = elastosdynamic_matrix(X,Y,1.0,2.0,1.0,1.0)
    # K = ElastostaticMatrix(X,Y,1.0,2.0)
    M = Matrix(K)
    irange,jrange = 1:m,1:n
    @testset "aca_full" begin
        atol = 1e-5
        aca = ACA(atol=atol)
        R = aca(K,irange,jrange);
        @test norm(Matrix(R) - M) < atol
        rtol = 1e-5
        aca = ACA(rtol=rtol)
        R   = aca(M,irange,jrange)
        norm(Matrix(R) - M)
        @test norm(Matrix(R) - M) < rtol*norm(M)
        r = 10
        aca = ACA(rank=r)
        R = aca(M,irange,jrange);
        @test rank(R) == r
    end
    @testset "aca_partial" begin
        atol = 1e-5
        aca = PartialACA(atol=atol)
        R = aca(M,irange,jrange)
        @test norm(Matrix(R) - M) < atol
        rtol = 1e-5
        aca = PartialACA(rtol=rtol)
        R = aca(M,irange,jrange)
        @test norm(Matrix(R) - M) < rtol*norm(M)
        r = 10
        aca = PartialACA(rank=r)
        R = aca(M,irange,jrange);
        @test rank(R) == r

        # test fast update of frobenius norm
        m,n = 10000,1000
        r = 10
        T = ComplexF64
        A = [rand(T,m) for _ in 1:r]
        B = [rand(T,n) for _ in 1:r]
        R = RkMatrix(A,B)
        old_norm = norm(Matrix(R),2)
        push!(A,rand(T,m))
        push!(B,rand(T,n))
        Rnew = RkMatrix(A,B)
        new_norm = norm(Matrix(Rnew),2)
        @test new_norm ≈ HMatrices._update_frob_norm(old_norm,A,B)

        # test simple case where things are not compressible
        A  = rand(2,2)
        comp = PartialACA(rtol=1e-5)
        @test comp(A,1:2,1:2) ≈ A
    end
end

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
export converter

mutable struct FromQEOutput
    o::PyObject
end
FromQEOutput(inp_file_list, inp_static, inp_q_points) =
    FromQEOutput(qha_input_maker.FromQEOutput(inp_file_list, inp_static, inp_q_points))

@pyinterface FromQEOutput

function converter(inp_file_list, inp_static, inp_q_points)
    conv = FromQEOutput(inp_file_list, inp_static, inp_q_points)
    conv.read_file_list()
    conv.read_static()
    conv.read_q_points()
    conv.read_frequency_files()
    conv.write_to_file()
end

ext: jl


text:
@testset "Spin tests" begin
    x = (0.0, 0.0, 0.0)
    N = 1000000
    for i in 1:N
        x = x .+ SpinMC.uniformOnSphere()
    end
    @test all(isapprox.(x ./ N, (0.0, 0.0, 0.0), atol=1.0e-2))

    M = [1.1 2.2 3.3; 4.4 5.5 6.6; 7.7 8.8 9.9]
    s1 = (1.0, 2.0, 3.0)
    s2 = (4.0, 5.0, 6.0)
    @test SpinMC.exchangeEnergy(s1, SpinMC.InteractionMatrix(M), s2) ≈ 607.2

    a1 = (3/2, sqrt(3)/2)
    a2 = (3/2, -sqrt(3)/2)
    basis = [(0.0, 0.0), (1.0, 0.0)]
    uc = UnitCell(a1, a2)
    b1 = addBasisSite!(uc, basis[1])
    b2 = addBasisSite!(uc, basis[2])
    addInteraction!(uc, b1, b2, [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], (0, 0))
    addInteraction!(uc, b1, b2, [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], (0, -1))
    addInteraction!(uc, b1, b2, [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], (-1, 0))
    lattice = Lattice(uc, (4,4))
    for i in 1:length(lattice)
        setSpin!(lattice, i, (0.0,0.0,1.0))
    end
    @test getEnergy(lattice) ≈ -48.0
    @test SpinMC.getEnergyDifference(lattice, 1, (0.0, 0.0, -1.0)) ≈ 6.0
    setSpin!(lattice, 1, (0.0, 0.0, -1.0))
    @test getEnergy(lattice) ≈ -42.0

    @test getMagnetization(lattice) ≈ [0.0, 0.0, 30.0/32.0]
    @test all(getCorrelation(lattice, 1) .≈ [(i == 1 ? 1.0 : -1.0) for i in 1:32])
    @test all(getCorrelation(lattice, 2) .≈ [(i == 1 ? -1.0 : 1.0) for i in 1:32])
end


text:
using Intervals: Beginning, Ending, LeftEndpoint, RightEndpoint, contiguous, overlaps,
    isbounded, isunbounded

function unique_paired_permutation(v::Vector{T}) where T
    results = Tuple{T, T}[]
    for (i, j) in unique(sort([x, y]) for x in eachindex(v), y in eachindex(v))
        push!(results, (v[i], v[j]))
    end
    return results
end

const INTERVAL_TYPES = [Interval, AnchoredInterval{Ending}, AnchoredInterval{Beginning}]

# Defines what conversions are possible. Tests for invalid conversions will be skipped.
viable_convert(::Type{Interval}, interval::AbstractInterval) = true

function viable_convert(::Type{AnchoredInterval{Beginning}}, interval::AbstractInterval)
    return isbounded(interval) && isfinite(first(interval))
end

function viable_convert(::Type{AnchoredInterval{Ending}}, interval::AbstractInterval)
    return isbounded(interval) && isfinite(last(interval))
end

@testset "comparisons: $A vs. $B" for (A, B) in unique_paired_permutation(INTERVAL_TYPES)

    # Compare two intervals which are non-overlapping:
    # Visualization of the finite case:
    #
    # [12]
    #     [45]
    @testset "non-overlapping" begin
        test_intervals = product(
            [
                Interval{Closed, Closed}(1, 2),
                Interval{Closed, Closed}(-Inf, 2),
                Interval{Closed, Closed}(-∞, 2),
                Interval{Unbounded,Closed}(nothing, 2),
            ],
            [
                Interval{Closed, Closed}(4, 5),
                Interval{Closed, Closed}(4, Inf),
                Interval{Closed, Closed}(4, ∞),
                Interval{Closed,Unbounded}(4, nothing),
            ],
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            earlier = convert(A, a)
            later = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(b))
            expected_overlap = Interval{promote_type(eltype(a), eltype(b))}()

            @test earlier != later
            @test !isequal(earlier, later)
            @test hash(earlier) != hash(later)

            @test isless(earlier, later)
            @test !isless(later, earlier)

            @test earlier < later
            @test !(later < earlier)

            @test earlier ≪ later
            @test !(later ≪ earlier)

            @test !issubset(earlier, later)
            @test !issubset(later, earlier)

            @test intersect(earlier, later) == expected_overlap
            @test_throws ArgumentError merge(earlier, later)
            @test union([earlier, later]) == [earlier, later]
            @test !overlaps(earlier, later)
            @test !contiguous(earlier, later)
            @test superset([earlier, later]) == expected_superset
        end
    end

    # Compare two intervals which "touch" but both intervals do not include that point:
    # Visualization of the finite case:
    #
    # (123)
    #   (345)
    @testset "touching open/open" begin
        test_intervals = product(
            [
                Interval{Open, Open}(1, 3),
                Interval{Open, Open}(-Inf, 3),
                Interval{Open, Open}(-∞, 3),
                Interval{Unbounded,Open}(nothing, 3),
            ],
            [
                Interval{Open, Open}(3, 5),
                Interval{Open, Open}(3, Inf),
                Interval{Open, Open}(3, ∞),
                Interval{Open,Unbounded}(3, nothing),
            ],
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            earlier = convert(A, a)
            later = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(b))
            expected_overlap = Interval{promote_type(eltype(a), eltype(b))}()

            @test earlier != later
            @test !isequal(earlier, later)
            @test hash(earlier) != hash(later)

            @test isless(earlier, later)
            @test !isless(later, earlier)

            @test earlier < later
            @test !(later < earlier)

            @test earlier ≪ later
            @test !(later ≪ earlier)

            @test !issubset(earlier, later)
            @test !issubset(later, earlier)

            @test intersect(earlier, later) == expected_overlap
            @test_throws ArgumentError merge(earlier, later)
            @test union([earlier, later]) == [earlier, later]
            @test !overlaps(earlier, later)
            @test !contiguous(earlier, later)
            @test superset([earlier, later]) == expected_superset
        end
    end

    # Compare two intervals which "touch" and the later interval includes that point:
    # Visualization of the finite case:
    #
    # (123)
    #   [345]
    @testset "touching open/closed" begin
         test_intervals = product(
            [
                Interval{Open, Open}(1, 3),
                Interval{Open, Open}(-Inf, 3),
                Interval{Open, Open}(-∞, 3),
                Interval{Unbounded,Open}(nothing, 3),
            ],
            [
                Interval{Closed, Closed}(3, 5),
                Interval{Closed, Closed}(3, Inf),
                Interval{Closed, Closed}(3, ∞),
                Interval{Closed,Unbounded}(3, nothing),
            ],
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            earlier = convert(A, a)
            later = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(b))
            expected_overlap = Interval{promote_type(eltype(a), eltype(b))}()

            @test earlier != later
            @test !isequal(earlier, later)
            @test hash(earlier) != hash(later)

            @test isless(earlier, later)
            @test !isless(later, earlier)

            @test earlier < later
            @test !(later < earlier)

            @test earlier ≪ later
            @test !(later ≪ earlier)

            @test !issubset(earlier, later)
            @test !issubset(later, earlier)

            @test intersect(earlier, later) == expected_overlap
            @test merge(earlier, later) == expected_superset
            @test union([earlier, later]) == [expected_superset]
            @test !overlaps(earlier, later)
            @test contiguous(earlier, later)
            @test superset([earlier, later]) == expected_superset
        end
    end

    # Compare two intervals which "touch" and the earlier interval includes that point:
    # Visualization of the finite case:
    #
    # [123]
    #   (345)
    @testset "touching closed/open" begin
        test_intervals = product(
            [
                Interval{Closed, Closed}(1, 3),
                Interval{Closed, Closed}(-Inf, 3),
                Interval{Closed, Closed}(-∞, 3),
                Interval{Unbounded,Closed}(nothing, 3),
            ],
            [
                Interval{Open, Open}(3, 5),
                Interval{Open, Open}(3, Inf),
                Interval{Open, Open}(3, ∞),
                Interval{Open,Unbounded}(3, nothing),
            ],
        )

        @testset "$a vs $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            earlier = convert(A, a)
            later = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(b))
            expected_overlap = Interval{promote_type(eltype(a), eltype(b))}()

            @test earlier != later
            @test !isequal(earlier, later)
            @test hash(earlier) != hash(later)

            @test isless(earlier, later)
            @test !isless(later, earlier)

            @test earlier < later
            @test !(later < earlier)

            @test earlier ≪ later
            @test !(later ≪ earlier)

            @test !issubset(earlier, later)
            @test !issubset(later, earlier)

            @test intersect(earlier, later) == expected_overlap
            @test merge(earlier, later) == expected_superset
            @test union([earlier, later]) == [expected_superset]
            @test !overlaps(earlier, later)
            @test contiguous(earlier, later)
            @test superset([earlier, later]) == expected_superset
        end
    end

    # Compare two intervals which "touch" and both intervals include that point:
    # Visualization of the finite case:
    #
    # [123]
    #   [345]
    @testset "touching closed/closed" begin
        test_intervals = product(
            [
                Interval{Closed, Closed}(1, 3),
                Interval{Closed, Closed}(-Inf, 3),
                Interval{Closed, Closed}(-∞, 3),
                Interval{Unbounded,Closed}(nothing, 3),
            ],
            [
                Interval{Closed, Closed}(3, 5),
                Interval{Closed, Closed}(3, Inf),
                Interval{Closed, Closed}(3, ∞),
                Interval{Closed,Unbounded}(3, nothing),
            ],
        )

        @testset "$a vs $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            earlier = convert(A, a)
            later = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(b))
            expected_overlap = Interval{Closed, Closed}(last(a), first(b))

            @test earlier != later
            @test !isequal(earlier, later)
            @test hash(earlier) != hash(later)

            @test isless(earlier, later)
            @test !isless(later, earlier)

            @test earlier < later
            @test !(later < earlier)

            @test !(earlier ≪ later)
            @test !(later ≪ earlier)

            @test !issubset(earlier, later)
            @test !issubset(later, earlier)

            @test intersect(earlier, later) == expected_overlap
            @test merge(earlier, later) == expected_superset
            @test union([earlier, later]) == [expected_superset]
            @test overlaps(earlier, later)
            @test !contiguous(earlier, later)
            @test superset([earlier, later]) == expected_superset
        end
    end

    # Compare two intervals which overlap
    # Visualization of the finite case:
    #
    # [1234]
    #  [2345]
    @testset "overlapping" begin
        test_intervals = product(
            [
                Interval{Closed, Closed}(1, 4),
                Interval{Closed, Closed}(-Inf, 4),
                Interval{Closed, Closed}(-∞, 4),
                Interval{Unbounded,Closed}(nothing, 4),
            ],
            [
                Interval{Closed, Closed}(2, 5),
                Interval{Closed, Closed}(2, Inf),
                Interval{Closed, Closed}(2, ∞),
                Interval{Closed,Unbounded}(2, nothing),
            ],
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            earlier = convert(A, a)
            later = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(b))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(a))

            @test earlier != later
            @test !isequal(earlier, later)
            @test hash(earlier) != hash(later)

            @test isless(earlier, later)
            @test !isless(later, earlier)

            @test earlier < later
            @test !(later < earlier)

            @test !(earlier ≪ later)
            @test !(later ≪ earlier)

            @test !issubset(earlier, later)
            @test !issubset(later, earlier)

            @test intersect(earlier, later) == expected_overlap
            @test merge(earlier, later) == expected_superset
            @test union([earlier, later]) == [expected_superset]
            @test overlaps(earlier, later)
            @test !contiguous(earlier, later)
            @test superset([earlier, later]) == expected_superset
        end
    end

    @testset "equal ()/()" begin
        test_intervals = (
            [
                Interval{Open, Open}(l, u),
                Interval{Open, Open}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(a))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(b))

            @test a == b
            @test isequal(a, b)
            @test hash(a) == hash(b)

            @test !isless(a, b)
            @test !isless(a, b)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal [)/()" begin
        test_intervals = (
            [
                Interval{Closed, Open}(l, u),
                Interval{Open, Open}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(a))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(b))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test isless(a, b)
            @test !isless(b, a)

            @test a < b
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test !issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal (]/()" begin
        test_intervals = (
            [
                Interval{Open, Closed}(l, u),
                Interval{Open, Open}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(a))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(b))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test !isless(a, b)
            @test !isless(b, a)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test !issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal []/()" begin
        test_intervals = (
            [
                Interval{Closed, Closed}(l, u),
                Interval{Open, Open}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(a))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(b))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test isless(a, b)
            @test !isless(b, a)

            @test a < b
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test !issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal [)/[]" begin
        test_intervals = (
            [
                Interval{Closed, Open}(l, u),
                Interval{Closed, Closed}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(b), RightEndpoint(b))
            expected_overlap = Interval(LeftEndpoint(a), RightEndpoint(a))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test !isless(a, b)
            @test !isless(b, a)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test issubset(a, b)
            @test !issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal (]/[]" begin
        test_intervals = (
            [
                Interval{Open, Closed}(l, u),
                Interval{Closed, Closed}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(b), RightEndpoint(b))
            expected_overlap = Interval(LeftEndpoint(a), RightEndpoint(a))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test !isless(a, b)
            @test isless(b, a)

            @test !(a < b)
            @test b < a

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test issubset(a, b)
            @test !issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal []/[]" begin
        test_intervals = (
            [
                Interval{Closed, Closed}(l, u),
                Interval{Closed, Closed}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(b), RightEndpoint(b))
            expected_overlap = Interval(LeftEndpoint(a), RightEndpoint(a))

            @test a == b
            @test isequal(a, b)
            @test hash(a) == hash(b)

            @test !isless(a, b)
            @test !isless(b, a)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal [/(" begin
        test_intervals = (
            [
                Interval{Closed,Unbounded}(l, u),
                Interval{Open,Unbounded}(l, u),
            ]
            for (l, u) in product((1, -Inf, -∞), (nothing,))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(a))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(b))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test isless(a, b)
            @test !isless(b, a)

            @test a < b
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test !issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal ]/)" begin
        test_intervals = (
            [
                Interval{Unbounded,Closed}(l, u),
                Interval{Unbounded,Open}(l, u),
            ]
            for (l, u) in product((nothing,), (5, Inf, ∞))
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(a), RightEndpoint(a))
            expected_overlap = Interval(LeftEndpoint(b), RightEndpoint(b))

            @test a != b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            @test !isless(a, b)
            @test !isless(b, a)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test !issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal unbounded" begin
        test_intervals = (
            [
                Interval{Int,Unbounded,Unbounded}(nothing, nothing),
                Interval{Int,Unbounded,Unbounded}(nothing, nothing),
            ],
        )

        @testset "$a vs. $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            a = convert(A, a)
            b = convert(B, b)
            expected_superset = Interval(LeftEndpoint(b), RightEndpoint(b))
            expected_overlap = Interval(LeftEndpoint(a), RightEndpoint(a))

            @test a == b
            @test isequal(a, b)
            @test hash(a) == hash(b)

            @test !isless(a, b)
            @test !isless(b, a)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    @testset "equal -0.0/0.0" begin
        # Skip tests when we're comparing ending and beginning anchored intervals
        if Set((A, B)) != Set((AnchoredInterval{Ending}, AnchoredInterval{Beginning}))
            a = convert(A, Interval(0.0, -0.0))
            b = convert(B, Interval(-0.0, 0.0))
            expected_superset = Interval(0.0, 0.0)
            expected_overlap = Interval(0.0, 0.0)

            @test a == b
            @test !isequal(a, b)
            @test hash(a) != hash(b)

            # All other comparison should still work as expected
            @test !isless(a, b)
            @test !isless(b, a)

            @test !(a < b)
            @test !(b < a)

            @test !(a ≪ b)
            @test !(b ≪ a)

            @test issubset(a, b)
            @test issubset(b, a)

            @test intersect(a, b) == expected_overlap
            @test merge(a, b) == expected_superset
            @test union([a, b]) == [expected_superset]
            @test overlaps(a, b)
            @test !contiguous(a, b)
            @test superset([a, b]) == expected_superset
        end
    end

    # Compare two intervals where the first interval is contained by the second
    # Visualization of the finite case:
    #
    #  [234]
    # [12345]
    @testset "containing" begin
        test_intervals = product(
            [
                Interval{Closed, Closed}(2, 4),
            ],
            [
                Interval{Closed, Closed}(1, 5),
                Interval{Closed, Closed}(1, Inf),
                Interval{Closed, Closed}(1, ∞),
                Interval{Closed, Closed}(-Inf, 5),
                Interval{Closed, Closed}(-∞, 5),
                Interval{Closed, Closed}(-Inf, Inf),
                Interval{Closed, Closed}(-∞, ∞),
                Interval{Closed, Unbounded}(1, nothing),
                Interval{Closed, Unbounded}(-Inf, nothing),
                Interval{Closed, Unbounded}(-∞, nothing),
                Interval{Unbounded, Closed}(nothing, 5),
                Interval{Unbounded, Closed}(nothing, Inf),
                Interval{Unbounded, Closed}(nothing, ∞),
                Interval{Unbounded, Unbounded}(nothing, nothing),
            ],
        )

        @testset "$a vs $b" for (a, b) in test_intervals
            viable_convert(A, a) || continue
            viable_convert(B, b) || continue

            smaller = convert(A, a)
            larger = convert(B, b)
            expected_superset = Interval(larger)
            expected_overlap = Interval(smaller)

            @test smaller != larger
            @test !isequal(smaller, larger)
            @test hash(smaller) != hash(larger)

            @test !isless(smaller, larger)
            @test isless(larger, smaller)

            @test !(smaller < larger)
            @test larger < smaller

            @test !(smaller ≪ larger)
            @test !(larger ≪ smaller)

            @test issubset(smaller, larger)
            @test !issubset(larger, smaller)

            @test intersect(smaller, larger) == expected_overlap
            @test merge(smaller, larger) == expected_superset
            @test union([smaller, larger]) == [expected_superset]
            @test overlaps(smaller, larger)
            @test !contiguous(smaller, larger)
            @test superset([smaller, larger]) == expected_superset
        end
    end
end


text:
const CURAND_STATUS_SUCCESS = (UInt32)(0)
const CURAND_STATUS_VERSION_MISMATCH = (UInt32)(100)
const CURAND_STATUS_NOT_INITIALIZED = (UInt32)(101)
const CURAND_STATUS_ALLOCATION_FAILED = (UInt32)(102)
const CURAND_STATUS_TYPE_ERROR = (UInt32)(103)
const CURAND_STATUS_OUT_OF_RANGE = (UInt32)(104)
const CURAND_STATUS_LENGTH_NOT_MULTIPLE = (UInt32)(105)
const CURAND_STATUS_DOUBLE_PRECISION_REQUIRED = (UInt32)(106)
const CURAND_STATUS_LAUNCH_FAILURE = (UInt32)(201)
const CURAND_STATUS_PREEXISTING_FAILURE = (UInt32)(202)
const CURAND_STATUS_INITIALIZATION_FAILED = (UInt32)(203)
const CURAND_STATUS_ARCH_MISMATCH = (UInt32)(204)
const CURAND_STATUS_INTERNAL_ERROR = (UInt32)(999)


text:
# jupyter.jl

import Base: show, mimewritable
import Base.Markdown: plain, Code, List, tag, htmlesc
import Requests: get
import URIParser: escape
import JSON: print
using Compat

Base.mimewritable(::Type{MIME"text/markdown"}, result::QueryResult) = true

function path_of_image_from_param(param::Dict)
    lhs = param["lhs"]
    simpler(t::Tuple) = last(t)
    simpler(a::Any) = a
    simple = join(map(simpler, lhs), '.')
    (escape(sprint(print,lhs)), simple)
end

function look_like_view(str::AbstractString)
    m = match(r"(0x[\dabcdef]*); frame = \(([\d\.]*) ([\d\.]*); ([\d\.]*) ([\d\.]*)\);", str)
    if !isa(m, RegexMatch)
        m = match(r"(0x[\dabcdef]*); frame = {{([\d\.]*), ([\d\.]*)}, {([\d\.]*), ([\d\.]*)}};", str)
    end
    if isa(m, RegexMatch)
        (address, xywh) = (m.captures[1], m.captures[2:end])
        (true, address, map(n->parse(Float32,n), xywh))
    else
        (false, nothing, nothing)
    end
end

function image_scale(str::AbstractString)
    (isview, address, xywh) = look_like_view(str)
    if isview
        (x,y,w,h) = xywh
        scale = 1
        if w > 1000 && h > 500
            scale = 1/5
        elseif w > 100 && h > 50
            scale = 1/2
        end
        "width:$(w*scale)px; height: $(h*scale)px;"
    else
        ""
    end
end

function save_image_url_to_file(url::AbstractString, path::AbstractString, r::AbstractString)
    isempty(path) && return url
    isdir("images") || mkdir("images")
    image_path = "images/$(path)-$r.png"
    f = open(image_path, "w")
    write(f, get(url).data)
    close(f)
    image_path
end

function escape_html(str::AbstractString)
    s = replace(str, "&", "&")
    s = replace(s, ">", ">")
    replace(s, "<", "<")
end

function save_image(url::AbstractString, alt::AbstractString)
    randhash = randstring(10)
    src = "$url&r=$randhash"
    save_to_local = true
    if haskey(ENV, "SWIFTER_SAVE_IMAGE")
        save_to_local = !(ENV["SWIFTER_SAVE_IMAGE"] in [false, "false", "0"])
    end
    if save_to_local # && !isempty(el)
        image_path = save_image_url_to_file(src, alt, randhash)
        src = "$image_path?r=$randhash"
    end
    src
end

function show_view_vector(stream::IO, mime::MIME"text/markdown", appvec::Vector{App}, vec::Vector)
    for (idx,el) in enumerate(vec)
        isa(el, Void) && continue
        (isview, address, xywh) = look_like_view(el)
        if isview
            write(stream, escape_html(el))
            app = appvec[idx]
            url = "$(app.url)/image?address=$address"
            tag(stream, :img, :src=>save_image(url, address), :alt=>address, :style=>image_scale(el))
        else
            plain(stream, escape_html(el))
        end
    end
end

@compat function Base.show(stream::IO, mime::MIME"text/markdown", vec::Vector{QueryResult})
    show_view_vector(stream, mime, map(v->v.app, vec), map(v->v.info.value, vec))
end

@compat function Base.show(stream::IO, mime::MIME"text/markdown", result::QueryResult)
    if isa(result.info.value, AbstractArray)
        show_view_vector(stream, mime, [result.app for x in 1:length(result.info.value)], result.info.value)
    elseif isa(result.info.value, Void)
        # nothing
    else
        plain(stream, Code(string(result.info.value)))
    end
    if :view == result.info.typ
        (path,simple) = path_of_image_from_param(result.param)
        app = result.app
        url = "$(app.url)/image?path=$path"
        tag(stream, :img, :src=>save_image(url, simple), :alt=>simple, :style=>image_scale(result.info.value))
    end
end


text:
abstract type AbstractThunk <: AbstractDifferential end

Base.Broadcast.broadcastable(x::AbstractThunk) = broadcastable(unthunk(x))

@inline function Base.iterate(x::AbstractThunk)
    val = unthunk(x)
    element, state = iterate(val)
    return element, (val, state)
end

@inline function Base.iterate(::AbstractThunk, (val, state))
    element, new_state = iterate(val, state)
    return element, (val, new_state)
end

#####
##### `Thunk`
#####

"""
    Thunk(()->v)
A thunk is a deferred computation.
It wraps a zero argument closure that when invoked returns a differential.
`@thunk(v)` is a macro that expands into `Thunk(()->v)`.

Calling a thunk, calls the wrapped closure.
If you are unsure if you have a `Thunk`, call [`unthunk`](@ref) which is a no-op when the
argument is not a `Thunk`.
If you need to unthunk recursively, call [`extern`](@ref), which also externs the differial
that the closure returns.

```jldoctest
julia> t = @thunk(@thunk(3))
Thunk(var"#4#6"())

julia> extern(t)
3

julia> t()
Thunk(var"#5#7"())

julia> t()()
3
```

### When to `@thunk`?
When writing `rrule`s (and to a lesser exent `frule`s), it is important to `@thunk`
appropriately.
Propagation rules that return multiple derivatives may not have all deriviatives used.
 By `@thunk`ing the work required for each derivative, they then compute only what is needed.

#### How do thunks prevent work?
If we have `res = pullback(...) = @thunk(f(x)), @thunk(g(x))`
then if we did `dx + res[1]` then only `f(x)` would be evaluated, not `g(x)`.
Also if we did `Zero() * res[1]` then the result would be `Zero()` and `f(x)` would not be evaluated.

#### So why not thunk everything?
`@thunk` creates a closure over the expression, which (effectively) creates a `struct`
with a field for each variable used in the expression, and call overloaded.

Do not use `@thunk` if this would be equal or more work than actually evaluating the expression itself.

For more details see the manual section [on using thunks effectively](http://www.juliadiff.org/ChainRulesCore.jl/dev/writing_good_rules.html#Use-Thunks-appropriately-1)
"""
struct Thunk{F} <: AbstractThunk
    f::F
end

"""
    @thunk expr

Define a [`Thunk`](@ref) wrapping the `expr`, to lazily defer its evaluation.
"""
macro thunk(body)
    # Basically `:(Thunk(() -> $(esc(body))))` but use the location where it is defined.
    # so we get useful stack traces if it errors.
    func = Expr(:->, Expr(:tuple), Expr(:block, __source__, body))
    return :(Thunk($(esc(func))))
end

"""
    unthunk(x)

On `AbstractThunk`s this removes 1 layer of thunking.
On any other type, it is the identity operation.

In contrast to [`extern`](@ref) this is nonrecursive.
"""
@inline unthunk(x) = x

@inline extern(x::AbstractThunk) = extern(unthunk(x))

# have to define this here after `@thunk` and `Thunk` is defined
Base.conj(x::AbstractThunk) = @thunk(conj(unthunk(x)))

(x::Thunk)() = x.f()
@inline unthunk(x::Thunk) = x()

Base.show(io::IO, x::Thunk) = println(io, "Thunk($(repr(x.f)))")

"""
    InplaceableThunk(val::Thunk, add!::Function)

A wrapper for a `Thunk`, that allows it to define an inplace `add!` function.

`add!` should be defined such that: `ithunk.add!(Δ) = Δ .+= ithunk.val`
but it should do this more efficently than simply doing this directly.
(Otherwise one can just use a normal `Thunk`).

Most operations on an `InplaceableThunk` treat it just like a normal `Thunk`;
and destroy its inplacability.
"""
struct InplaceableThunk{T<:Thunk, F} <: AbstractThunk
    val::T
    add!::F
end

unthunk(x::InplaceableThunk) = unthunk(x.val)
(x::InplaceableThunk)() = unthunk(x)

function Base.show(io::IO, x::InplaceableThunk)
    println(io, "InplaceableThunk($(repr(x.val)), $(repr(x.add!)))")
end


text:
@with_kw_noshow struct TopPushK{S<:AbstractSurrogate, I<:Integer, T<:Real} <: AbstractTopPushK{S}
    K::I
    C::T = 1
    l::S = Hinge(ϑ = 1.0)

    @assert K >= 1 "The vaule of `K` must be greater or equal to 1."
end

TopPushK(K::Int; kwargs...) = TopPushK(; K = K, kwargs...)
show(io::IO, model::TopPushK) = print(io, "TopPushK($(model.K), $(model.C), $(model.l))")

@with_kw_noshow struct TopPush{S<:AbstractSurrogate, T<:Real} <: AbstractTopPushK{S}
    C::T = 1
    l::S = Hinge(ϑ = 1.0)
end

show(io::IO, model::TopPush) = print(io, "TopPush($(model.C), $(model.l))")

@with_kw_noshow struct τFPL{S<:AbstractSurrogate, R<:Real, T<:Real} <: AbstractTopPushK{S}
    τ::R
    C::T = 1
    l::S = Hinge(ϑ = 1.0)

    @assert 0 < τ < 1  "The vaule of `τ` must lay in the interval (0,1)."
end

τFPL(τ::Real; kwargs...) = τFPL(; τ = τ, kwargs...)
show(io::IO, model::τFPL) = print(io, "τFPL($(model.τ), $(model.C), $(model.l))")

getK(model::TopPushK, data) = model.K
getK(model::TopPush, data) = 1
getK(model::τFPL, data::Primal) = max(1, round(Int, model.τ * data.nneg))
getK(model::τFPL, data::Dual) = max(1, round(Int, model.τ * data.nβ))

# -------------------------------------------------------------------------------
# Primal problem - General solver
# -------------------------------------------------------------------------------
function optimize(solver::General, model::AbstractTopPushK, data::Primal)

    Xpos = @view data.X[data.ind_pos, :]
    Xneg = @view data.X[data.ind_neg, :]

    w = Convex.Variable(data.dim)
    t = Convex.Variable()
    y = Convex.Variable(data.npos)
    z = Convex.Variable(data.nneg)
    K = getK(model, data)

    objective = Convex.sumsquares(w)/2 + model.C * Convex.sum(model.l.value_exact.(y))
    constraints = [y == t + Convex.sum(z)/K - Xpos*w,
                   z >= Xneg*w - t,
                   z >= 0]

    problem = Convex.minimize(objective, constraints)
    Convex.solve!(problem, solver.solver)

    return vec(w.value), t.value + sum(z.value)/K
end

function optimize(solver::General, model::TopPush, data::Primal)

    Xpos = @view data.X[data.ind_pos, :]
    Xneg = @view data.X[data.ind_neg, :]

    w = Convex.Variable(data.dim)
    t = Convex.Variable()
    y = Convex.Variable(data.npos)

    objective = Convex.sumsquares(w)/2 + model.C * Convex.sum(model.l.value_exact.(y))
    constraints = [y == t - Xpos*w,
                   t >= Convex.maximum(Xneg*w)]

    problem = Convex.minimize(objective, constraints)
    Convex.solve!(problem, solver.solver)

    return vec(w.value), t.value
end

# -------------------------------------------------------------------------------
# Primal problem - Gradient descent solver
# -------------------------------------------------------------------------------
function initialization(model::AbstractTopPushK, data::Primal)
    w = zeros(eltype(data.X), data.dim)
    Δ = zero(w)
    s = scores(model, data, w)
    return w, s, Δ
end

function objective(model::AbstractTopPushK, data::Primal, w, t, s = scores(model, data, w))
    return w'*w/2 + model.C * sum(model.l.value.(t .- s[data.ind_pos]))
end

function threshold(model::AbstractTopPushK, data::Primal, s)
    K = getK(model, data)
   return mean(partialsort(s[data.ind_neg], 1:K, rev = true))
end

function threshold(::TopPush, data::Primal, s)
   return maximum(s[data.ind_neg])
end

function gradient!(model::AbstractTopPushK, data::Primal, w, s, Δ)
    K = getK(model, data)
    ind_t = partialsortperm(s[data.ind_neg], 1:K, rev = true)
    t     = mean(s[data.ind_neg[ind_t]])

    ∇l = model.l.gradient.(t .- s[data.ind_pos])
    ∇t = vec(mean(data.X[data.ind_neg[ind_t], :], dims = 1))

    Δ .= w .+ model.C .* (sum(∇l)*∇t .- data.X[data.ind_pos,:]'*∇l)
    return t
end

function gradient!(model::TopPush, data::Primal, w, s, Δ)
    t, ind_t = findmax(s[data.ind_neg])

    ∇l = model.l.gradient.(t .- s[data.ind_pos])
    ∇t = vec(data.X[data.ind_neg[ind_t], :])

    Δ .= w .+ model.C .* (sum(∇l)*∇t .- data.X[data.ind_pos,:]'*∇l)
    return t
end

# -------------------------------------------------------------------------------
# Dual problem - General solver
# -------------------------------------------------------------------------------
# Hinge loss
function optimize(solver::General, model::M, data::Dual{<:DTrain}; ε::Real = 1e-6) where {M<:AbstractTopPushK{<:Hinge}}

    α = Convex.Variable(data.nα)
    β = Convex.Variable(data.nβ)
    K = getK(model, data)
    KK = data.K + ε .* I

    objective   = - Convex.quadform(vcat(α, β), KK)/2 + Convex.sum(α)/model.l.ϑ
    constraints = [Convex.sum(α) == Convex.sum(β),
                   α <= model.l.ϑ*model.C,
                   α >= 0,
                   β >= 0]

    if M <: TopPushK || M <: τFPL
        push!(constraints, β <= Convex.sum(α)/K)
    end

    problem = Convex.maximize(objective, constraints)
    Convex.solve!(problem, solver.solver)

    return vec(α.value), vec(β.value)
end

# Truncated quadratic loss
function optimize(solver::General, model::M, data::Dual{<:DTrain}; ε::Real = 1e-6) where {M<:AbstractTopPushK{<:Quadratic}}

    α = Convex.Variable(data.nα)
    β = Convex.Variable(data.nβ)
    K = getK(model, data)
    KK = data.K + ε .* I

    objective   = - Convex.quadform(vcat(α, β), KK)/2 +
                    Convex.sum(α)/model.l.ϑ - Convex.sumsquares(α)/(4*model.C*model.l.ϑ^2)
    constraints = [Convex.sum(α) == Convex.sum(β),
                   α >= 0,
                   β >= 0]

    if M <: TopPushK || M <: τFPL
        push!(constraints, β <= Convex.sum(α)/K)
    end

    problem = Convex.maximize(objective, constraints)
    Convex.solve!(problem, solver.solver)

    return vec(α.value), vec(β.value)
end

# -------------------------------------------------------------------------------
# Dual problem - Graient descent solver
# -------------------------------------------------------------------------------
function initialization(model::AbstractTopPushK, data::Dual{<:DTrain})
    αβ   = rand(eltype(data.K), data.nα + data.nβ)
    α, β = @views αβ[data.ind_α], αβ[data.ind_β]
    projection!(model, data, α, β)

    s  = data.K * vcat(α, β)
    return α, β, αβ, s
end

# Hinge loss
function objective(model::AbstractTopPushK{<:Hinge}, data::Dual{<:DTrain}, α, β, s = data.K * vcat(α, β))
    - s'*vcat(α, β)/2 + sum(α)/model.l.ϑ
end

function gradient!(model::AbstractTopPushK{<:Hinge}, data::Dual{<:DTrain}, α, β, s, Δ)
    Δ             .= .- s
    Δ[data.ind_α] .+= 1/model.l.ϑ
end

function projection!(model::M, data::Dual{<:DTrain}, α, β) where {M<:AbstractTopPushK{<:Hinge}}
    K = getK(model, data)
    αs, βs = projection(α, β, model.l.ϑ*model.C, K)
    α .= αs
    β .= βs
    return α, β
end

# Truncated quadratic loss
function objective(model::AbstractTopPushK{<:Quadratic}, data::Dual{<:DTrain}, α, β, s = data.K * vcat(α, β))
    - s'*vcat(α, β)/2 + sum(α)/model.l.ϑ - sum(abs2, α)/(4*model.C*model.l.ϑ^2)
end

function gradient!(model::AbstractTopPushK{<:Quadratic}, data::Dual{<:DTrain}, α, β, s, Δ)
    Δ             .= .- s
    Δ[data.ind_α] .+= 1/model.l.ϑ .- α/(2*model.l.ϑ^2*model.C)
end

function projection!(model::M, data::Dual{<:DTrain}, α, β) where {M<:AbstractTopPushK{<:Quadratic}}
    K = getK(model, data)
    αs, βs = projection(α, β, K)
    α .= αs
    β .= βs
    return α, β
end

function threshold(model::AbstractTopPushK, data::Dual, s)
    K = getK(model, data)
   return mean(partialsort(.- s[data.ind_β], 1:K, rev = true))
end

function threshold(::TopPush, data::Dual, s)
   return maximum(.- s[data.ind_β])
end

function primal_objective(model::AbstractTopPushK, data::Dual{<:DTrain}, α, β, s = data.K * vcat(α, β))

    K = getK(model, data)
    t = threshold(model, data, s)
    z = max.(.- s[data.ind_β] .- t, 0)
    y = t + sum(z)/K  .- s[data.ind_α]
    w_norm = s'*vcat(α, β)/2
    return w_norm + model.C * sum(model.l.value.(y))
end

# -------------------------------------------------------------------------------
# Dual problem - Coordinate descent solver
# -------------------------------------------------------------------------------
function loss(::AbstractTopPushK, ::Dual{<:DTrain}, a::Real, b::Real, Δ::Real)
    a*Δ^2/2 + b*Δ
end

function select_k(::AbstractTopPushK, data::Dual{<:DTrain}, α, β)
    rand(1:(data.nα + data.nβ))
end

function apply!(::AbstractTopPushK, data::Dual{<:DTrain}, best::BestUpdate, α, β, αβ, s, αsum, βsort)
    βsorted!(data, best, β, βsort)
    if best.k <= data.nα && best.l > data.nα
        αsum .+= best.Δ
    end
    αβ[best.k] = best.vars[1]
    αβ[best.l] = best.vars[2]
    scores!(data, best, s)
end

function apply!(::TopPush, data::Dual{<:DTrain}, best::BestUpdate, α, β, αβ, s, αsum, βsort)
    αβ[best.k] = best.vars[1]
    αβ[best.l] = best.vars[2]
    scores!(data, best, s)
end

# Hinge loss
function rule_αα!(model::AbstractTopPushK{<:Hinge}, data::Dual{<:DTrain}, best::BestUpdate, k, l, α, β, s, αsum, βsort)

    αk, αl = α[k], α[l]
    C, ϑ   = model.C, model.l.ϑ

    a = - data.K[k,k] + 2*data.K[k,l] - data.K[l,l]
    b = - s[k] + s[l]
    Δ = solution(a, b, max(-αk, αl - ϑ*C), min(ϑ*C - αk, αl))

    vars = (αk = αk + Δ, αl = αl - Δ)
    L    = loss(model, data, a, b, Δ)
    update!(best, k, l, Δ, L, vars)
end

function rule_αβ!(model::M, data::Dual{<:DTrain}, best::BestUpdate, k, l, α, β, s, αsum, βsort) where {M<:AbstractTopPushK{<:Hinge}}

    αk, βl = α[k], β[l - data.nα]
    C, ϑ   = model.C, model.l.ϑ
    K = getK(model, data)

    a = - data.K[k,k] - 2*data.K[k,l] - data.K[l,l]
    b = - s[k] - s[l] + 1/ϑ

    if K == 1
        Δ = solution(a, b, max(-αk, -βl), C*ϑ - αk)
    else
        βmax = find_βmax(βsort, β, l - data.nα)
        Δ    = solution(a, b, max(-αk, -βl, K*βmax - αsum[1]), min(C*ϑ - αk, (αsum[1] - K*βl)/(K-1)))
    end

    vars = (αk = αk + Δ, βl = βl + Δ)
    L    = loss(model, data, a, b, Δ)
    update!(best, k, l, Δ, L, vars)
end

function rule_ββ!(model::M, data::Dual{<:DTrain}, best::BestUpdate, k, l, α, β, s, αsum, βsort) where {M<:AbstractTopPushK{<:Hinge}}

    βk, βl = β[k - data.nα], β[l - data.nα]
    C, ϑ   = model.C, model.l.ϑ
    K = getK(model, data)

    a = - data.K[k,k] + 2*data.K[k,l] - data.K[l,l]
    b = - s[k] + s[l]

    if K == 1
        Δ = solution(a, b, -βk, βl)
    else
        Δ = solution(a, b, max(-βk, βl - αsum[1]/K), min(αsum[1]/K - βk, βl))
    end

    vars = (βk = βk + Δ, βl = βl - Δ)
    L    = loss(model, data, a, b, Δ)
    update!(best, k, l, Δ, L, vars)
end

# Truncated quadratic loss
function rule_αα!(model::AbstractTopPushK{<:Quadratic}, data::Dual{<:DTrain}, best::BestUpdate, k, l, α, β, s, αsum, βsort)

    αk, αl = α[k], α[l]
    C, ϑ   = model.C, model.l.ϑ

    a = - data.K[k,k] + 2*data.K[k,l] - data.K[l,l] - 1/(C*ϑ^2)
    b = - s[k] + s[l] - (αk - αl)/(2*C*ϑ^2)
    Δ = solution(a, b, - αk, αl)

    vars = (αk = αk + Δ, αl = αl - Δ)
    L    = loss(model, data, a, b, Δ)
    update!(best, k, l, Δ, L, vars)
end

function rule_αβ!(model::M, data::Dual{<:DTrain}, best::BestUpdate, k, l, α, β, s, αsum, βsort) where {M<:AbstractTopPushK{<:Quadratic}}

    αk, βl = α[k], β[l - data.nα]
    C, ϑ   = model.C, model.l.ϑ
    K = getK(model, data)

    a = - data.K[k,k] - 2*data.K[k,l] - data.K[l,l] - 1/(2*C*ϑ^2)
    b = - s[k] - s[l] + 1/ϑ - αk/(2*C*ϑ^2)

    if K == 1
        Δ = solution(a, b, max(-αk, -βl), Inf)
    else
        βmax = find_βmax(βsort, β, l - data.nα)
        Δ    = solution(a, b, max(-αk, -βl, K*βmax - αsum[1]), (αsum[1] - K*βl)/(K-1))
    end

    vars = (αk = αk + Δ, βl = βl + Δ)
    L    = loss(model, data, a, b, Δ)
    update!(best, k, l, Δ, L, vars)
end

function rule_ββ!(model::M, data::Dual{<:DTrain}, best::BestUpdate, k, l, α, β, s, αsum, βsort) where {M<:AbstractTopPushK{<:Quadratic}}

    βk, βl = β[k - data.nα], β[l - data.nα]
    C, ϑ   = data.n, model.C, model.l.ϑ
    K = getK(model, data)

    a = - data.K[k,k] + 2*data.K[k,l] - data.K[l,l]
    b = - s[k] + s[l]

    if K == 1
        Δ = solution(a, b, -βk, βl)
    else
        Δ = solution(a, b, max(-βk, βl - αsum[1]/K), min(αsum[1]/K - βk, βl))
    end

    vars = (βk = βk + Δ, βl = βl - Δ)
    L    = loss(model, data, a, b, Δ)
    update!(best, k, l, Δ, L, vars)
end


text:
function geom_point(cc::Cairo.CairoContext, ctxs::Vector{<:AbstractRefFrame})

    x, y = user2dev(ctxs,
                    getparam(ctxs, :x; fail = true),
                    getparam(ctxs, :y; fail = true))

    shape = getparam(ctxs, :shape, "circle")
    color = getparam(ctxs, :color, "black")
    fill = getparam(ctxs, :fill, nothing)
    alpha = getparam(ctxs, :alpha, 1.0)
    size = getparam(ctxs, :size, 2.0)

    points(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray;
           alpha = alpha, color = color, fill = fill, shape = shape, size = size)
end

function geom_line(cc::Cairo.CairoContext, ctxs::Vector{<:AbstractRefFrame})

    x, y = user2dev(ctxs,
                    getparam(ctxs, :x; fail = true),
                    getparam(ctxs, :y; fail = true))

    color = getparam(ctxs, :color, "black")
    alpha = getparam(ctxs, :alpha, 1.0)
    linetype = getparam(ctxs, :linetype, "solid")
    size = getparam(ctxs, :size, 1.0)

    lines(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray;
           alpha = alpha, color = color, size = size, linetype = linetype)
end

function geom_col(cc::Cairo.CairoContext, ctxs::Vector{<:AbstractRefFrame})

    x, y = user2dev(ctxs,
                    getparam(ctxs, :x; fail = true),
                    getparam(ctxs, :y; fail = true))

    basex, basey = user2dev(ctxs, 0, 0)

    color = getparam(ctxs, :color, "black")
    alpha = getparam(ctxs, :alpha, 1.0)
    linetype = getparam(ctxs, :linetype, "solid")
    size = getparam(ctxs, :size, 1.0)

    cols(cc::Cairo.CairoContext, x::AbstractArray, y::AbstractArray, basey;
           alpha = alpha, color = color, size = size, linetype = linetype)
end

const geometry_funcs = Dict(:point => geom_point,
                            :line => geom_line,
                            :col => geom_col)

function geometry(cc::Cairo.CairoContext, ctxs::Vector{<:AbstractRefFrame})
    geom = getparam(ctxs, :geom, :point)
    @assert haskey(geometry_funcs, geom) "Unknown geometry: $geom"
    geometry_funcs[geom](cc, ctxs)
end


text:
export mis1

"""
Solving MIS problem with simple branching algorithm.
"""
function mis1(eg::EliminateGraph)
    N = nv(eg)
    if N == 0
        return 0
    else
        vmin, dmin = mindegree_vertex(eg)
        return 1 + neighborcover_mapreduce(y->eliminate(mis1, eg, NeighborCover(y)), max, eg, vmin)
    end
end


text:
const _arg_desc = KW(

# series args
:label 			=> "String type. The label for a series, which appears in a legend.  If empty, no legend entry is added.",
:seriescolor 		=> "Color Type. The base color for this series.  `:auto` (the default) will select a color from the subplot's `color_palette`, based on the order it was added to the subplot",
:seriesalpha 		=> "Number in [0,1]. The alpha/opacity override for the series.  `nothing` (the default) means it will take the alpha value of the color.",
:seriestype 		=> "Symbol. This is the identifier of the type of visualization for this series. Choose from $(_allTypes) or any series recipes which are defined.",
:linestyle 			=> "Symbol. Style of the line (for path and bar stroke).  Choose from $(_allStyles)",
:linewidth         	=> "Number. Width of the line (in pixels)",
:linecolor         	=> "Color Type. Color of the line (for path and bar stroke).  `:match` will take the value from `:seriescolor`, (though histogram/bar types use `:black` as a default).",
:linealpha         	=> "Number in [0,1]. The alpha/opacity override for the line.  `nothing` (the default) means it will take the alpha value of linecolor.",
:fillrange         	=> "Number or AbstractVector.  Fills area between fillrange and y for line-types, sets the base for bar/stick types, and similar for other types.",
:fillcolor         	=> "Color Type. Color of the filled area of path or bar types.  `:match` will take the value from `:seriescolor`.",
:fillalpha         	=> "Number in [0,1]. The alpha/opacity override for the fill area.  `nothing` (the default) means it will take the alpha value of fillcolor.",
:markershape       	=> "Symbol, Shape, or AbstractVector. Choose from $(_allMarkers).",
:markercolor       	=> "Color Type. Color of the interior of the marker or shape. `:match` will take the value from `:seriescolor`.",
:markeralpha       	=> "Number in [0,1]. The alpha/opacity override for the marker interior.  `nothing` (the default) means it will take the alpha value of markercolor.",
:markersize        	=> "Number or AbstractVector. Size (radius pixels) of the markers.",
:markerstrokestyle 	=> "Symbol. Style of the marker stroke (border).  Choose from $(_allStyles)",
:markerstrokewidth 	=> "Number. Width of the marker stroke (border. in pixels)",
:markerstrokecolor 	=> "Color Type. Color of the marker stroke (border).  `:match` will take the value from `:foreground_color_subplot`.",
:markerstrokealpha 	=> "Number in [0,1]. The alpha/opacity override for the marker stroke (border).  `nothing` (the default) means it will take the alpha value of markerstrokecolor.",
:bins              	=> "Integer, NTuple{2,Integer}, AbstractVector or Symbol. Default is :auto (the Freedman-Diaconis rule). For histogram-types, defines the approximate number of bins to aim for, or the auto-binning algorithm to use (:sturges, :sqrt, :rice, :scott or :fd). For fine-grained control pass a Vector of break values, e.g. `range(minimum(x), stop = maximum(x), length = 25)`",
:smooth            	=> "Bool.  Add a regression line?",
:group             	=> "AbstractVector. Data is split into a separate series, one for each unique value in `group`.",
:x                 	=> "Various. Input data. First Dimension",
:y                 	=> "Various. Input data. Second Dimension",
:z                 	=> "Various. Input data. Third Dimension. May be wrapped by a `Surface` for surface and heatmap types.",
:marker_z          	=> "AbstractVector, Function `f(x,y,z) -> z_value`, or Function `f(x,y) -> z_value`, or nothing. z-values for each series data point, which correspond to the color to be used from a markercolor gradient.",
:line_z          	=> "AbstractVector, Function `f(x,y,z) -> z_value`, or Function `f(x,y) -> z_value`, or nothing. z-values for each series line segment, which correspond to the color to be used from a linecolor gradient.  Note that for N points, only the first N-1 values are used (one per line-segment).",
:fill_z                 => "Matrix{Float64} of the same size as z matrix, which specifies the color of the 3D surface; the default value is `nothing`.",
:levels            	=> "Integer, NTuple{2,Integer}, or AbstractVector. Levels or number of levels (or x-levels/y-levels) for a contour type.",
:orientation       	=> "Symbol.  Horizontal or vertical orientation for bar types.  Values `:h`, `:hor`, `:horizontal` correspond to horizontal (sideways, anchored to y-axis), and `:v`, `:vert`, and `:vertical` correspond to vertical (the default).",
:bar_position      	=> "Symbol.  Choose from `:overlay` (default), `:stack`.  (warning: May not be implemented fully)",
:bar_width              => "nothing or Number. Width of bars in data coordinates. When nothing, chooses based on x (or y when `orientation = :h`).",
:bar_edges              => "Bool.  Align bars to edges (true), or centers (the default)?",
:xerror                 => "AbstractVector or 2-Tuple of Vectors. x (horizontal) error relative to x-value.  If 2-tuple of vectors, the first vector corresponds to the left error (and the second to the right)",
:yerror            	=> "AbstractVector or 2-Tuple of Vectors. y (vertical) error relative to y-value.  If 2-tuple of vectors, the first vector corresponds to the bottom error (and the second to the top)",
:ribbon            	=> "Number or AbstractVector.  Creates a fillrange around the data points.",
:quiver            	=> "AbstractVector or 2-Tuple of vectors.  The directional vectors U,V which specify velocity/gradient vectors for a quiver plot.",
:arrow             	=> "nothing (no arrows), Bool (if true, default arrows), Arrow object, or arg(s) that could be style or head length/widths.  Defines arrowheads that should be displayed at the end of path line segments (just before a NaN and the last non-NaN point).  Used in quiverplot, streamplot, or similar.",
:normalize         	=> "Bool or Symbol. Histogram normalization mode. Possible values are: false/:none (no normalization, default), true/:pdf (normalize to a discrete Probability Density Function, where the total area of the bins is 1), :probability (bin heights sum to 1) and :density (the area of each bin, rather than the height, is equal to the counts - useful for uneven bin sizes).",
:weights           	=> "AbstractVector. Used in histogram types for weighted counts.",
:show_empty_bins    => "Bool. Whether empty bins in a 2D histogram are colored as 0 (true), or transparent (the default).",
:contours          	=> "Bool. Add contours to the side-grids of 3D plots?  Used in surface/wireframe.",
:contour_labels         => "Bool. Show labels at the contour lines?",
:match_dimensions  	=> "Bool. For heatmap types... should the first dimension of a matrix (rows) correspond to the first dimension of the plot (x-axis)?  The default is false, which matches the behavior of Matplotlib, Plotly, and others.  Note: when passing a function for z, the function should still map `(x,y) -> z`.",
:subplot           	=> "Integer (subplot index) or Subplot object.  The subplot that this series belongs to.",
:series_annotations     => "AbstractVector of String or PlotText.  These are annotations which are mapped to data points/positions.",
:primary                => "Bool.  Does this count as a 'real series'?  For example, you could have a path (primary), and a scatter (secondary) as 2 separate series, maybe with different data (see sticks recipe for an example).  The secondary series will get the same color, etc as the primary.",
:hover 	        	=> "nothing or vector of strings. Text to display when hovering over each data point.",
:colorbar_entry   => "Bool.  Include this series in the color bar?  Set to `false` to exclude.",

# plot args
:plot_title               => "String. Title for the whole plot (not the subplots) (Note: Not currently implemented)",
:background_color         => "Color Type.  Base color for all backgrounds.",
:background_color_outside => "Color Type or `:match` (matches `:background_color`).  Color outside the plot area(s)",
:foreground_color         => "Color Type.  Base color for all foregrounds.",
:size                     => "NTuple{2,Int}. (width_px, height_px) of the whole Plot",
:pos                      => "NTuple{2,Int}. (left_px, top_px) position of the GUI window (note: currently unimplemented)",
:window_title             => "String. Title of the standalone gui-window.",
:show                     => "Bool.  Should this command open/refresh a GUI/display?  This allows displaying in scripts or functions without explicitly calling `display`",
:layout                   => "Integer (number of subplots), NTuple{2,Integer} (grid dimensions), AbstractLayout (for example `grid(2,2)`), or the return from the `@layout` macro.  This builds the layout of subplots.",
:link                     => "Symbol.  How/whether to link axis limits between subplots. Values: `:none`, `:x` (x axes are linked by columns), `:y` (y axes are linked by rows), `:both` (x and y are linked), `:all` (every subplot is linked together regardless of layout position).",
:overwrite_figure         => "Bool.  Should we reuse the same GUI window/figure when plotting (true) or open a new one (false).",
:html_output_format       => "Symbol.  When writing html output, what is the format?  `:png` and `:svg` are currently supported.",
:tex_output_standalone    => "Bool. When writing tex output, should the source include a preamble for a standalone document class.",
:inset_subplots           => "nothing or vector of 2-tuple (parent,bbox).  optionally pass a vector of (parent,bbox) tuples which are the parent layout and the relative bounding box of inset subplots",
:dpi 			  => "Number.  Dots Per Inch of output figures",
:thickness_scaling        => "Number. Scale for the thickness of all line elements like lines, borders, axes, grid lines, ... defaults to 1.",
:display_type             => "Symbol (`:auto`, `:gui`, or `:inline`).  When supported, `display` will either open a GUI window or plot inline.",
:extra_kwargs 	          => "Either one of (`:plot`, `:subplot`, `:series`) to specify for which element extra keyword args are collected or a KW (Dict{Symbol,Any}) to pass a map of extra keyword args which may be specific to a backend. Default: `:series`.\n Example: `pgfplotsx(); scatter(1:5, extra_kwargs=Dict(:subplot=>Dict(\"axis line shift\" => \"10pt\"))`",
:fontfamily               => "String or Symbol.  Default font family for title, legend entries, tick labels and guides",
:warn_on_unsupported      => "Bool. Warn on unsupported attributes, series types and marker shapes",

# subplot args
:title                    => "String. Subplot title.",
:titlelocation           => "Symbol. Position of subplot title. Values: `:left`, `:center`, `:right`",
:titlefontfamily          => "String or Symbol. Font family of subplot title.",
:titlefontsize            => "Integer. Font pointsize of subplot title.",
:titlefonthalign          => "Symbol. Font horizontal alignment of subplot title: :hcenter, :left, :right or :center",
:titlefontvalign          => "Symbol. Font vertical alignment of subplot title: :vcenter, :top, :bottom or :center",
:titlefontrotation        => "Real. Font rotation of subplot title",
:titlefontcolor           => "Color Type. Font color of subplot title",
:background_color_subplot => "Color Type or `:match` (matches `:background_color`).  Base background color of the subplot.",
:background_color_legend  => "Color Type or `:match` (matches `:background_color_subplot`).  Background color of the legend.",
:background_color_inside  => "Color Type or `:match` (matches `:background_color_subplot`).  Background color inside the plot area (under the grid).",
:foreground_color_subplot => "Color Type or `:match` (matches `:foreground_color`).  Base foreground color of the subplot.",
:foreground_color_legend  => "Color Type or `:match` (matches `:foreground_color_subplot`).  Foreground color of the legend.",
:foreground_color_title   => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of subplot title.",
:color_palette            => "Vector of colors (cycle through) or color gradient (generate list from gradient) or `:auto` (generate a color list using `Colors.distiguishable_colors` and custom seed colors chosen to contrast with the background).  The color palette is a color list from which series colors are automatically chosen.",
:legend                   => "Bool (show the legend?) or (x,y) tuple or Symbol (legend position) or angle or (angle,inout) tuple. Bottom left corner of legend is placed at (x,y).  Symbol values: `:none`; `:best`; `:inline`; `:inside`; `:legend`; any valid combination of `:(outer ?)(top/bottom ?)(right/left ?)`, i.e.: `:top`, `:topright`, `:outerleft`, `:outerbottomright` ... (note: only some may be supported in each backend). Legend is positioned at (angle degrees) (so (90,:outer) is roughly equivalent to :outertop), close to the inside of the axes or the outside if inout=:outer.",
:legendfontfamily         => "String or Symbol. Font family of legend entries.",
:legendfontsize           => "Integer. Font pointsize of legend entries.",
:legendfonthalign         => "Symbol. Font horizontal alignment of legend entries: :hcenter, :left, :right or :center",
:legendfontvalign         => "Symbol. Font vertical alignment of legend entries: :vcenter, :top, :bottom or :center",
:legendfontrotation       => "Real. Font rotation of legend entries",
:legendfontcolor          => "Color Type. Font color of legend entries",
:legendtitle              => "String. Legend title.",
:legendtitlefontfamily    => "String or Symbol. Font family of the legend title.",
:legendtitlefontsize           => "Integer. Font pointsize the legend title.",
:legendtitlefonthalign         => "Symbol. Font horizontal alignment of the legend title: :hcenter, :left, :right or :center",
:legendtitlefontvalign         => "Symbol. Font vertical alignment of the legend title: :vcenter, :top, :bottom or :center",
:legendtitlefontrotation       => "Real. Font rotation of the legend title",
:legendtitlefontcolor          => "Color Type. Font color of the legend title",
:colorbar                 => "Bool (show the colorbar?) or Symbol (colorbar position).  Symbol values: `:none`, `:best`, `:right`, `:left`, `:top`, `:bottom`, `:legend` (matches legend value) (note: only some may be supported in each backend)",
:clims 					  => "`:auto`, NTuple{2,Number}, or a function that takes series data in and returns NTuple{2,Number}.  Fixes the limits of the colorbar.",
:legendfont               => "Font. Font of legend items.",
:legendtitlefont          => "Font. Font of the legend title.",
:annotations              => "(x,y,text) tuple(s).  Can be a single tuple or a list of them.  Text can be String or PlotText (created with `text(args...)`)  Add one-off text annotations at the x,y coordinates.",
:projection               => "Symbol or String.  '3d' or 'polar'",
:aspect_ratio             => "Symbol (:equal or :none) or Number. Plot area is resized so that 1 y-unit is the same size as `aspect_ratio` x-units. With `:none`, images inherit aspect ratio of the plot area.",
:margin                   => "Measure (multiply by `mm`, `px`, etc). Base for individual margins... not directly used.  Specifies the extra padding around subplots.",
:left_margin              => "Measure (multiply by `mm`, `px`, etc) or `:match` (matches `:margin`).  Specifies the extra padding to the left of the subplot.",
:top_margin               => "Measure (multiply by `mm`, `px`, etc) or `:match` (matches `:margin`).  Specifies the extra padding on the top of the subplot.",
:right_margin             => "Measure (multiply by `mm`, `px`, etc) or `:match` (matches `:margin`).  Specifies the extra padding to the right of the subplot.",
:bottom_margin            => "Measure (multiply by `mm`, `px`, etc) or `:match` (matches `:margin`).  Specifies the extra padding on the bottom of the subplot.",
:subplot_index            => "Integer.  Internal (not set by user).  Specifies the index of this subplot in the Plot's `plt.subplot` list.",
:colorbar_title           => "String.  Title of colorbar.",
:framestyle               => "Symbol.  Style of the axes frame. Choose from $(_allFramestyles)",
:camera                   => "NTuple{2, Real}. Sets the view angle (azimuthal, elevation) for 3D plots",

# axis args
:guide     		  => "String. Axis guide (label).",
:guide_position		  => "Symbol. Position of axis guides: :top, :bottom, :left or :right",
:lims      		  => "NTuple{2,Number} or Symbol. Force axis limits.  Only finite values are used (you can set only the right limit with `xlims = (-Inf, 2)` for example). `:round` widens the limit to the nearest round number ie. [0.1,3.6]=>[0.0,4.0]",
:ticks     		  => "Vector of numbers (set the tick values), Tuple of (tickvalues, ticklabels), or `:auto`",
:scale     		  => "Symbol. Scale of the axis: `:none`, `:ln`, `:log2`, `:log10`",
:rotation  		  => "Number. Degrees rotation of tick labels.",
:flip      		  => "Bool.  Should we flip (reverse) the axis?",
:formatter                => "Function, :scientific, :plain or :auto. A method which converts a number to a string for tick labeling.",
:tickfontfamily           => "String or Symbol. Font family of tick labels.",
:tickfontsize             => "Integer. Font pointsize of tick labels.",
:tickfonthalign           => "Symbol. Font horizontal alignment of tick labels: :hcenter, :left, :right or :center",
:tickfontvalign           => "Symbol. Font vertical alignment of tick labels: :vcenter, :top, :bottom or :center",
:tickfontrotation         => "Real. Font rotation of tick labels",
:tickfontcolor            => "Color Type. Font color of tick labels",
:guidefontfamily          => "String or Symbol. Font family of axes guides.",
:guidefontsize            => "Integer. Font pointsize of axes guides.",
:guidefonthalign          => "Symbol. Font horizontal alignment of axes guides: :hcenter, :left, :right or :center",
:guidefontvalign          => "Symbol. Font vertical alignment of axes guides: :vcenter, :top, :bottom or :center",
:guidefontrotation        => "Real. Font rotation of axes guides",
:guidefontcolor           => "Color Type. Font color of axes guides",
:foreground_color_axis    => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of axis ticks.",
:foreground_color_border  => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of plot area border (spines).",
:foreground_color_text    => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of tick labels.",
:foreground_color_guide   => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of axis guides (axis labels).",
:mirror                   => "Bool.  Switch the side of the tick labels (right or top).",
:grid                     => "Bool, Symbol, String or `nothing`. Show the grid lines? `true`, `false`, `:show`, `:hide`, `:yes`, `:no`, `:x`, `:y`, `:z`, `:xy`, ..., `:all`, `:none`, `:off`",
:foreground_color_grid    => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of grid lines.",
:gridalpha                => "Number in [0,1]. The alpha/opacity override for the grid lines.",
:gridstyle                => "Symbol. Style of the grid lines. Choose from $(_allStyles)",
:gridlinewidth            => "Number. Width of the grid lines (in pixels)",
:foreground_color_minor_grid => "Color Type or `:match` (matches `:foreground_color_subplot`). Color of minor grid lines.",
:minorgrid                => "Bool. Adds minor grid lines and ticks to the plot. Set minorticks to change number of gridlines",
:minorticks               => "Integer. Intervals to divide the gap between major ticks into",
:minorgridalpha           => "Number in [0,1]. The alpha/opacity override for the minorgrid lines.",
:minorgridstyle           => "Symbol. Style of the minor grid lines. Choose from $(_allStyles)",
:minorgridlinewidth       => "Number. Width of the minor grid lines (in pixels)",
:tick_direction           => "Symbol.  Direction of the ticks. `:in` or `:out`",
:showaxis                 => "Bool, Symbol or String.  Show the axis. `true`, `false`, `:show`, `:hide`, `:yes`, `:no`, `:x`, `:y`, `:z`, `:xy`, ..., `:all`, `:off`",
:widen                    => "Bool. Widen the axis limits by a small factor to avoid cut-off markers and lines at the borders. Defaults to `true`.",
:draw_arrow               => "Bool. Draw arrow at the end of the axis.",
)


text:
@testset "ImagePreprocessing" begin
    encfns = [
        () -> ImagePreprocessing(),
        () -> ImagePreprocessing(buffered=false),
        () -> ImagePreprocessing(T=Float64),
        () -> ImagePreprocessing(augmentations=FastAI.augs_lighting()),
        () -> ImagePreprocessing(C=Gray{N0f8}, means=SVector(0.), stds=SVector(1.)),
    ]
    for encfn in encfns
        enc = encfn()
        block = Image{2}()
        img = rand(RGB{N0f8}, 10, 10)
        testencoding(enc, block, img)

        ctx = Validation()
        outblock = encodedblock(enc, block)
        a = encode(enc, ctx, block, img)
        rimg = decode(enc, ctx, outblock, a)
        if eltype(rimg) <: RGB
            @test img ≈ rimg
        end
    end

    @testset "3D" begin
        enc = ImagePreprocessing()
        block = Image{3}()
        img = rand(RGB{N0f8}, 10, 10, 10)
        testencoding(enc, block, img)

        enc = ImagePreprocessing(buffered = false)
        block = Image{3}()
        img = rand(RGB{N0f8}, 10, 10, 10)
        testencoding(enc, block, img)
    end
end


text:
using BayesianDataFusion

if nprocs() <= 1
  addprocs(8)
  @everywhere using BayesianDataFusion
end

W = sprand(1_500_000, 1000, 0.01)
#W = sprand(150_000, 1000, 0.01)
rd = RelationData(W, class_cut = 0.5)
assignToTest!(rd.relations[1], 50)

############ num_latent = 10 ############
## single threaded
blas_set_num_threads(1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10);

blas_set_num_threads(8)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10);

blas_set_num_threads(16)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10);

## multi-threaded (1)
blas_set_num_threads(1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2], latent_blas_threads=1);

## multi-threaded (2)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=1);
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=2);

############ num_latent = 30 ############
blas_set_num_threads(1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=30);

blas_set_num_threads(2)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=30);

blas_set_num_threads(4)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=30);

## multi-threaded (1)
r1 = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2], latent_blas_threads=1);

## multi-threaded (2)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=1);
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3], latent_blas_threads=2);

## multi-threaded (4)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2,3,4,5], latent_blas_threads=1);

## multi-threaded (8)
r = macau(rd, burnin = 2, psamples = 2, num_latent=10, latent_pids=[2:9;], latent_blas_threads=1);


text:
# Contains the rocket set up and structs

#=
Now we want to set up the rocket dynamics. We have a linear system:
xk+1 = Axk + Buk + G

where xk = [sk; vk] (the position and velocity, respectively). If sk is n×1,
A is 2n×2n, B is 2n×n, and G is a 2n×1 vector.

Thus, we get the system:
[A B -I] * [xk; uk; xk+1] = G
which we can stack as

[A B -I 0 0 0 0 0 0]
[0 0 A B -I 0 0 0 0]
[0 0 0 0 A B -I 0 0]
[0 0 0 0 0 0 A B -I]
=#

using SparseArrays, LinearAlgebra

# include("trajectory-setup.jl")

"""
    rocket_simple(mass, isp, grav, deltaTime)

Barebones "spherical" rocket. This is to say, this struct only holds the mass
and the isp (specific impulse) of the rocket.
"""
struct rocket_simple
    mass
    isp
    grav
    deltaTime
end

@doc raw"""
    rocketDynamics(r::rocket_simple, nDim::Int64)

We have a linear dynamics system, that when continuous is
```math
\frac{dx}{dt} = Ax + Bu + g
```
If we have `y = [x; u; g]`, then,
```math
\frac{dy}{dt} = [A B I; 0 0 0; 0 0 0] y = W y
```
which has solution
```math
y(t) = e^{W t} y0
```
where W is square. In a discrete system,
```math
yk+1 = e^{W Δt} yk = E yk
```

Note that due to all the zeros in W, we get:
```math
xk+1 = E11 xk + E12 uk + E13 g
```
```math
uk+1 = E21 xk + E22 uk + E23 g = uk
```
```math
g = E31 xk + E32 uk + E33 g = g
```

Thus, once we calculate E, we extract E11, E12, and E13 to a system of the
form:
```math
x_{k+1} = Aₑ x_k + Bₑ u_k + Gₑ
```
We compute and return Aₑ, Bₑ, and Gₑ for a [`rocket_simple`](@ref).

For a simple rocket, we have

```math
A = [0 \ I; 0 \ 0] \quad B = [0; \frac{1}{m} I] \quad G = [0; -g]
```

We then take the exponent and return a resulting sparse matrices

"""
function rocketDynamics(r::rocket_simple, nDim::Int64)

    # Start off with a dense representation bc it is needed for the
    # exponentiation.
    z0 = zeros(nDim, nDim)
    aMat = [z0 I; z0 z0]
    bMat = [z0; -(1/r.mass) * I]
    gMat = [z0; I]
    # gVec = [spzeros(nDim, 1); r.grav]

    # Top Row
    abg = [aMat bMat gMat]
    abg_size = size(abg)

    # We need to pad the remaining rows to make a square matrix
    remainingRows = abg_size[2] - abg_size[1]
    abgTot = [abg; zeros(remainingRows, abg_size[2])]

    # Take the exponent e^(W Δt)
    eMat = exp(abgTot * r.deltaTime)

    # Extract the effective matrices
    rowsA = size(aMat, 1)
    colsA = size(aMat, 2)
    colsB = size(bMat, 2)

    aEffMat = sparse(eMat[1:rowsA, 1:colsA])
    bEffMat = sparse(eMat[1:rowsA, colsA + 1:colsA + colsB])
    gEffMat = sparse(eMat[1:rowsA, colsA + colsB + 1:end])

    return aEffMat, bEffMat, gEffMat
end

@doc raw"""
    rocketDynamicsStack(r::rocket_simple, nDim::Int64, NSteps::Int64)

We want to take the per discretization dynamics and write it as a full matrix.

Specifically, we get the system:
`[A B -I] * [xk; uk; xk+1] = G`
which we can stack as

```math
[A \ B \ -I \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]
```

```math
[0 \ 0 \ A \ B \ -I \ 0 \ 0 \ 0 \ 0]
```

```math
[0 \ 0 \ 0 \ 0 \ A \ B \ -I \ 0 \ 0]
```

```math
[0 \ 0 \ 0 \ 0 \ 0 \ 0 \ A \ B \ -I]
```

We return the above matrix and the stacked `G` matrix.

See: [`rocketDynamicsFull`](@ref) and [`rocketDynamics`](@ref)
"""
function rocketDynamicsStack(r::rocket_simple, nDim::Int64, NSteps::Int64)
    # First calculate the A and B matrices. The C vector is not needed in this
    # case.
    Ak, Bk, Gk = rocketDynamics(r::rocket_simple, nDim::Int64)
    ABI_unit = [Ak Bk -I]
    colLen = size(ABI_unit, 2)

    # At each step we go down 2 * n rows and move right 3 * n columns.
    sizeCols = 3 * NSteps * nDim + 2 * nDim
    sizeRows = 2 * nDim * NSteps

    # Initialize and empty sparse A matrix and G vector
    AStacked = spzeros(sizeRows, sizeCols)
    GStacked = spzeros(sizeRows, 1)

    rStart = 1
    cStart = 1
    for k in 1:NSteps
        rEnd = rStart + (2 * nDim - 1)
        cEnd = cStart + colLen - 1

        # First, we calculate that stacked A Matrix
        # println("Accessing $rStart:$rEnd and $cStart:$cEnd")
        AStacked[rStart:rEnd, cStart:cEnd] = ABI_unit

        # Second, we calculate the stacked G vector
        GStacked[rStart:rEnd] = -Gk * r.grav

        # Prepare for next iteration.
        rStart += 2 * nDim
        cStart += 3 * nDim
    end

    return AStacked, GStacked
end

"""
    rocketDynamicsFull(r::rocket_simple, x0::Array{Float64, 1},
                             xN::Array{Float64, 1}, NSteps::Int64,
                             noEnd::Bool = false)

We combine the stack at each `xk` with the initial and final conditions. Namely,

`[I0; AStack; IN] * XU = [x0; GStack; xN]`

where `XU` is the stacked trajectory vector (x states and u controls). We
return `[I0; AStack; IN]` and `[x0; GStack; xN]`

See [`rocketDynamicsStack`](@ref), [`rocketDynamics`](@ref),
[`rocket_simple`](@ref)
"""
function rocketDynamicsFull(r::rocket_simple, x0::Array{Float64, 1},
                                 xN::Array{Float64, 1}, NSteps::Int64,
                                 noEnd::Bool = false)
     if size(x0) != size(xN)
         error("x0 and xN must have the same size. " *
                 "Currently: x0 = $x0 and xN = $xN")
     end
     if mod(size(x0, 1), 2) != 0
         error("x0 and xN must have an even number of variables. " *
               "Example, x0 = [sx0; sy0; vx0; vy0]. Currently x0 = $x0.")
     end
    # The number of columns must be the same for all three. (Stacked vertically)
    nDim = Int64(size(x0, 1) / 2)
    sizeCols = 3 * NSteps * nDim + 2 * nDim
    sizeRowsI = 2 * nDim

    # Use the above function to determine the stacked matrix in the middle.
    AStacked, GStacked = rocketDynamicsStack(r, nDim, NSteps)
    I0 = [I spzeros(sizeRowsI, sizeCols - sizeRowsI)]

    if noEnd
        AFull = [I0; AStacked]
        BFull = [x0; GStacked]
    else
        IN = [spzeros(sizeRowsI, sizeCols - sizeRowsI) I]

        AFull = [I0; AStacked; IN]
        BFull = [x0; GStacked; xN]
    end

    return AFull, BFull
end

function rocketDynamicsFull(r::rocket_simple, x0::Array{Int64, 1},
                                 xN::Array{Int64, 1}, NSteps::Int64)
    return rocketDynamicsFull(r, float.(x0), float.(xN), NSteps)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
module Causal

using Spec
using DocStringExtensions
using Cassette
using Callbacks: donothing

using ..Space
import ..Space: combinetags
import ..Omega
using ..IDS: ID
using ..NonDet: RandVar, constant, proj
import ..NonDet: apl, ppapl, id, params, ciid, uid

include("replace.jl")
include("counterfactual.jl")
include("causes.jl")

export replace,
       cf,
       iscausebf
end


text:
module TestMapFunctions

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running map function tests")

##################################
## map for columnwise iterators ##
##################################

tm = TimeData.testcase(TimeData.Timematr, 1)

## eachcol
##--------

## multiplication
td3 = map(x -> x.*3.0, eachcol(tm)) # function acts on DataArray
tm3 = TimeData.Timematr(td3.vals, td3.idx)
expOut = tm.*3
@test TimeData.isequal(tm3, expOut)

## isna for each column
tdisna = map(x -> isna(x), eachcol(tm)) # function acts on DataArray
@test tdisna.vals[:, 1] == convert(DataArray, falses(4))

## cumulated sum for each column
tdaggr = map(x -> cumsum(x), eachcol(tm))
@test tdaggr.vals[:, 1] == @data([100, 220, 330, 500])

## eachvar
##--------

## multiplication
td3 = map(x -> x.*3.0, TimeData.eachvar(tm))
tm3 = TimeData.Timematr(td3.vals, td3.idx)
expOut = tm.*3
@test TimeData.isequal(tm3, expOut)

## isna for each column
tdisna = map(x -> TimeData.isnaElw(x), TimeData.eachvar(tm))
@test tdisna.vals[:, 1] == convert(DataArray, falses(4))

## cumulated sum for each column
tdaggr = map(x -> cumsum(x, 1), TimeData.eachvar(tm))
@test tdaggr.vals[:, 1] == @data([100, 220, 330, 500])

###############################
## map for rowwise iterators ##
###############################

tm = TimeData.testcase(TimeData.Timematr, 1)

## eachrow
##--------

## identity function
tdIden = map(x -> x, TimeData.eachrow(tm))
tmIden = TimeData.Timematr(tdIden.vals, tdIden.idx)
@test TimeData.isequal(tm, tmIden)

## value of second column everywhere
tdDup = map(x -> [x[1, 2], x[1, 2]], TimeData.eachrow(tm))
tmDup = TimeData.Timematr(tdDup.vals, tdDup.idx)
expOut = TimeData.Timematr(DataFrame(prices1 = tm.vals[:, 2],
                                     prices2 = tm.vals[:, 2]),
                           TimeData.idx(tm))
@test TimeData.isequal(tmDup, expOut)

## eachdate
##---------

## multiplication
td3 = map(x -> (x.*3.5).vals, TimeData.eachdate(tm))
tm3 = TimeData.Timematr(td3.vals, td3.idx)
expOut = tm.*3.5
@test TimeData.isequal(tm3, expOut)

td3 = map(x -> x.*3.5, TimeData.eachdate(tm))
tm3 = TimeData.Timematr(td3.vals, td3.idx)
@test TimeData.isequal(tm3, expOut)

## cumulated sum for each column
tdaggr = map(x -> cumsum(x, 2), TimeData.eachdate(tm))
@test tdaggr.vals[:, 1] == tm.vals[:, 1]
@test tdaggr.vals[:, 2] == @data([210, 240, 210, 300])

tdaggr = map(x -> cumsum(x, 2).vals, TimeData.eachdate(tm))
@test tdaggr.vals[:, 1] == tm.vals[:, 1]
@test tdaggr.vals[:, 2] == @data([210, 240, 210, 300])

###################################
## map for elementwise iterators ##
###################################

## multiplication
tn = TimeData.testcase(TimeData.Timenum, 2)
td2 = map(x -> 2*x, TimeData.eachentry(tn))
tn2 = convert(TimeData.Timenum, td2)
expOut = tn.*2
@test isequal(tn2, expOut)

tn = TimeData.testcase(TimeData.Timenum, 2)
td2 = map(x -> 2.*get(x, 1, 1), TimeData.eachobs(tn))
## mathematical operations on DataArray{NAtype, 1} with single row are
## currently not working

tn2 = convert(TimeData.Timenum, td2)
expOut = tn.*2
@test isequal(tn2, expOut)

##################
## test chkVars ##
##################

## column sum greater 490
td = TimeData.testcase(TimeData.Timedata, 1)
actOut = TimeData.chkVars(x -> sum(x) > 490, eachcol(td))
expOut = DataFrame(prices1 = true, prices2 = false)
@test isequal(expOut, actOut)

## NA in column
td = TimeData.testcase(TimeData.Timedata, 3)
actOut = TimeData.chkVars(x -> any(isna(x)), eachcol(td))
expOut = DataFrame(color = true, value = true)
@test isequal(expOut, actOut)

## column mean greater than 120 using eachvar
tm = TimeData.testcase(TimeData.Timematr, 1)
actOut = TimeData.chkVars(x -> mean(x, 1)[1, 1] .> 120,
                          TimeData.eachvar(tm))
expOut = DataFrame(prices1 = true, prices2 = false)
@test isequal(expOut, actOut)

## NA at day 2010-01-01
td = TimeData.testcase(TimeData.Timedata, 3)
actOut = TimeData.chkVars(x -> TimeData.isnaElw(x[Date(2010,1,1)]),
                          TimeData.eachvar(td))
expOut = DataFrame(color = true, value = false)
@test isequal(expOut, actOut)

###################
## test chkDates ##
###################

## is NA in row?
td = TimeData.testcase(TimeData.Timedata, 3)
function anyNA(dfr::DataFrame)
    containsNA = false
    for ii=1:size(dfr, 2)
        if isna(dfr[1, ii])
            containsNA = true
            break
        end
    end
    return containsNA
end

tdHolds = TimeData.chkDates(x -> anyNA(x), TimeData.eachrow(td))
@test isequal(tdHolds.vals[:, 1], @data([true, false, true]))

## is row equal to given date
tdHolds = TimeData.chkDates(x -> TimeData.idx(x) == [Date(2010,1,1)],
                            TimeData.eachdate(td))

@test isequal(tdHolds.vals[:, 1], @data([true, false, false]))

#################
## test chkElw ##
#################

tn = TimeData.testcase(TimeData.Timenum, 2)
tnHolds = TimeData.chkElw(x -> x.>120, TimeData.eachentry(tn))
expOut = tn .> 120
@test isequal(tnHolds, expOut)

########################
## test collapseDates ##
########################

## eachcol
##--------

## column sum with NAs
tn = TimeData.testcase(TimeData.Timenum, 2)
df = TimeData.collapseDates(x -> sum(x), eachcol(tn))
expOut = DataFrame(prices1 = [NA], prices2 = [NA])
@test isequal(df, expOut)

## column sum skipping NAs
function sumSkipNA(da::DataArray)
    sumVal = 0
    for ii=1:length(da)
        if !isna(da[ii])
            sumVal += da[ii]
        end
    end
    return sumVal
end

tn = TimeData.testcase(TimeData.Timenum, 2)
df = TimeData.collapseDates(x -> sumSkipNA(x), eachcol(tn))
expOut = DataFrame(prices1 = 630, prices2 = 360)
@test isequal(df, expOut)

## eachvar
##--------

## column sum
tm = TimeData.testcase(TimeData.Timematr, 1)
df = TimeData.collapseDates(x -> sum(x, 1), TimeData.eachvar(tm))
expOut = DataFrame(prices1 = 500, prices2 = 460)
@test isequal(df, expOut)

## count NAs
function countNAs(td::TimeData.AbstractTimedata)
    return sum(TimeData.asArr(TimeData.isnaElw(td), Bool, false))
end

tn = TimeData.testcase(TimeData.Timenum, 2)
df = TimeData.collapseDates(x -> countNAs(x),
                            TimeData.eachvar(tn))
expOut = DataFrame(prices1 = 1, prices2 = 2)
@test isequal(df, expOut)

#######################
## test collapseVars ##
#######################

## eachrow
##--------

## count NAs
function countNAs(df::DataFrame)
    sumNAs = 0
    for ii=1:size(df, 2)
        if isna(df[1, ii])
            sumNAs += 1
        end
    end
    return sumNAs
end

tn = TimeData.testcase(TimeData.Timenum, 2)
td = TimeData.collapseVars(x -> countNAs(x),
                           eachrow(tn))

expOut = TimeData.Timedata(DataFrame(funcVal = [1, 0, 1, 0, 1]),
                           TimeData.idx(tn))

@test isequal(expOut, td)

## eachdate
##---------

## calculate mean
function meanDfRow(df::DataFrame)
    sumVal = 0
    sumNAs = 0
    for ii=1:size(df, 2)
        if !isna(df[1, ii])
            sumVal += df[1, ii]
        else
            sumNAs += 1
        end
    end
    return sumVal ./ (size(df, 2) - sumNAs)
end

tn = TimeData.testcase(TimeData.Timenum, 2)
td = TimeData.collapseVars(x -> meanDfRow(x.vals),
                           TimeData.eachdate(tn))

mVals = [110, 120, 140, 150, 200]
expOut = TimeData.Timedata(DataFrame(funcVal = mVals),
                           TimeData.idx(tn))

@test isequal(expOut, td)

## using already existing rowmeans function for Timenum
td = TimeData.collapseVars(x -> TimeData.rowmeans(x),
                           TimeData.eachdate(tn))
@test isequal(expOut, td)
end


text:
##  SCRIPT NAME:   objectiveFunction.jl
##  DATE:          January 2015
#
#   This script returns the objective function in each time period, that is:
#   μH uH(cH,K) + μW uW(cW,K) + β EV(a')
#   where the way it is written I assume the choice variables are a', cH, cW,
#   and K (i.e. next period's assets, private consumption, and public consumption).
#   It is the felicity the family enjoys if current state variables are: a
#   (assets), wH, and wW (wages), and the family chooses next period assets a',
#   private consumption cH and cW, and public consumption. The continuation value
#   is inputted in the function.
#
#
##  ALEXANDROS THELOUDIS, UCL

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

function objectiveFunction(choice_var::Vector, grad::Vector, A0, wH0, wW0, interpol_A!_EVal!)

  # Declare choice variables:
  A!        = choice_var[1]
  consH     = choice_var[2]
  consW     = choice_var[3]
  kons      = choice_var[4]

  # Calculate gradient:
  if length(grad) > 0
    (value,derivative) = valgrad(interpol_A!_EVal!, A!)

    grad[1] = beta * derivative
    grad[2] = paretoH * mUtilH_cH(consH, kons)
    grad[3] = paretoW * mUtilW_cW(consW, kons)
    grad[4] = paretoH * mUtilH_k(consH, kons) + paretoW * mUtilW_k(consW, kons)
  end

  # Calculate objective (with interpolated future value function):
  objective = paretoH * utilityH(consH, kons) + paretoW * utilityW(consW, kons) + beta * interpol_A!_EVal![A!]

  return objective

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using DifferentiableStateSpaceModels, SparseArrays, LinearAlgebra, Parameters, Test,
      TimerOutputs, BenchmarkTools

@testset "Dense SW First Order" begin
    m = @include_example_module(Examples.SW07)

    p = [10, 0.51, 10, 0, 0.7, 0.742, 0, 0, 0.24, 0.2696, 6.0144, 0.025, 1.5, 0.6361, 1.5,
         0.3243, 0.8087, 0.47, 0.6, 1.9423, 1.5, 1.488, 0.2347, 0.0593, 0.8762, 0.9977,
         0.5799, 0.9957, 0.7165, 0, 0, 0, 0.3982, 0.18]
    generate_perturbation(m, p)
    reset_timer!()
    sol = generate_perturbation(m, p)
    print_timer()
    @test sol.retcode == :Success
end