text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries

using Test
using LinearAlgebra

@testset "Tests for HomogeneousPolynomial and TaylorN" begin
    eeuler = Base.MathConstants.e

    @test HomogeneousPolynomial <: AbstractSeries
    @test HomogeneousPolynomial{Int} <: AbstractSeries{Int}
    @test TaylorN{Float64} <: AbstractSeries{Float64}

    set_variables("x", numvars=2, order=6)
    _taylorNparams = TaylorSeries.ParamsTaylorN(6, 2, String["x₁", "x₂"])
    @test _taylorNparams.order == get_order()
    @test _taylorNparams.num_vars == get_numvars()
    @test _taylorNparams.variable_names == get_variable_names()
    @test _taylorNparams.variable_symbols == get_variable_symbols()

    @test eltype(set_variables(Int, "x", numvars=2, order=6)) == TaylorN{Int}
    @test eltype(set_variables("x", numvars=2, order=6)) == TaylorN{Float64}
    @test eltype(set_variables(BigInt, "x y", order=6)) == TaylorN{BigInt}
    @test eltype(set_variables("x y", order=6)) == TaylorN{Float64}
    @test eltype(set_variables(Int, :x, numvars=2, order=6)) == TaylorN{Int}
    @test eltype(set_variables(:x, numvars=2, order=6)) == TaylorN{Float64}
    @test eltype(set_variables(BigInt, [:x,:y], order=6)) == TaylorN{BigInt}
    @test eltype(set_variables([:x,:y], order=6)) == TaylorN{Float64}
    @test typeof(show_params_TaylorN()) == Nothing
    @test typeof(show_monomials(2)) == Nothing

    @test TaylorSeries.coeff_table[2][1] == [1,0]
    @test TaylorSeries.index_table[2][1] == 7
    @test TaylorSeries.in_base(get_order(),[2,1]) == 15
    @test TaylorSeries.pos_table[4][15] == 2

    @test get_order() == 6
    @test get_numvars() == 2

    @test get_variables()[1].order == get_order()
    @test get_variables(2)[1].order == 2
    @test get_variables(3)[1] == TaylorN(1, order=3)
    @test get_variables(Int, 3)[1] == TaylorN(Int, 1, order=3)
    @test length(get_variables()) == get_numvars()

    x, y = set_variables("x y", order=6)
    @test axes(x) == axes(y) == ()
    @test axes(x[1]) == axes(y[2]) == ()
    @test size(x) == (7,)
    @test size(x[1]) == (2,)
    @test size(x[2]) == (3,)
    @test firstindex(x) == 0
    @test firstindex(x[end]) == 1
    @test lastindex(y) == get_order()
    @test eachindex(x) == 0:6
    @test iterate(x) == (HomogeneousPolynomial([0.0], 0), 1)
    @test iterate(y, 1) == (HomogeneousPolynomial([0.0, 1.0], 1), 2)
    @test iterate(x, 7) == nothing

    @test x.order == 6
    @test TaylorSeries.name_taylorNvar(1) == " x"
    @test TaylorSeries._params_TaylorN_.variable_names == ["x","y"]
    @test TaylorSeries._params_TaylorN_.variable_symbols == [:x, :y]
    @test get_variable_symbols() == [:x, :y]
    @test TaylorSeries.lookupvar(:x) == 1
    @test TaylorSeries.lookupvar(:α) == 0
    @test TaylorSeries.get_variable_names() == ["x", "y"]
    @test x == HomogeneousPolynomial(Float64, 1)
    @test x == HomogeneousPolynomial(1)
    @test y == HomogeneousPolynomial(Float64, 2)
    @test y == HomogeneousPolynomial(2)
    @test !isnan(x)

    set_variables("x", numvars=2, order=17)
    v = [1,2]
    @test typeof(TaylorSeries.resize_coeffsHP!(v,2)) == Nothing
    @test v == [1,2,0]
    @test_throws AssertionError TaylorSeries.resize_coeffsHP!(v,1)
    hpol_v = HomogeneousPolynomial(v)
    hpol_v[3] = 3
    @test v == [1,2,3]
    hpol_v[1:3] = 3
    @test v == [3,3,3]
    hpol_v[1:2:2] = 0
    @test v == [0,3,3]
    hpol_v[1:1:2] = [1,2]
    @test all(hpol_v[1:1:2] .== [1,2])
    @test v == [1,2,3]
    hpol_v[:] = zeros(Int, 3)
    @test hpol_v == 0

    tn_v = TaylorN(HomogeneousPolynomial(zeros(Int, 3)))
    tn_v[0] = 1
    @test tn_v == 1
    tn_v[0:1] = [0, 1]
    @test tn_v[0] == 0 && tn_v[1] == HomogeneousPolynomial(1, 1)
    tn_v[0:1] = [HomogeneousPolynomial(0, 0), HomogeneousPolynomial([0,1])]
    @test tn_v[0] == 0 && tn_v[1] == HomogeneousPolynomial([0,1], 1)
    tn_v[:] = [HomogeneousPolynomial(1, 0), HomogeneousPolynomial(0, 1), hpol_v]
    @test tn_v == 1
    tn_v[:] = 0
    @test tn_v == 0
    tn_v[:] = [3,1,0]
    @test tn_v == TaylorN([HomogeneousPolynomial(3, 0), HomogeneousPolynomial(1, 1)], 2)
    tn_v[0:2] = [HomogeneousPolynomial(3, 0), HomogeneousPolynomial(1, 1), HomogeneousPolynomial(0, 2)]
    @test tn_v == TaylorN([HomogeneousPolynomial(3, 0), HomogeneousPolynomial(1, 1)], 2)
    tn_v[0:2:2] = [0,0]
    @test tn_v == TaylorN(HomogeneousPolynomial(1, 1), 2)

    xH = HomogeneousPolynomial([1,0])
    yH = HomogeneousPolynomial([0,1],1)
    @test xH == convert(HomogeneousPolynomial{Float64},xH)
    @test HomogeneousPolynomial(0,0)  == 0
    xT = TaylorN(xH, 17)
    yT = TaylorN(Int, 2, order=17)
    zeroT = zero( TaylorN([xH],1) )
    @test zeroT.coeffs == zeros(HomogeneousPolynomial{Int}, 1)
    @test size(xH) == (2,)
    @test firstindex(xH) == 1
    @test lastindex(yH) == 2
    @test length(zeros(HomogeneousPolynomial{Int}, 1)) == 2
    @test one(HomogeneousPolynomial(1,1)) == HomogeneousPolynomial([1,1])
    uT = one(convert(TaylorN{Float64},yT))
    @test uT == one(HomogeneousPolynomial)
    @test uT == convert(TaylorN{Float64},uT)
    @test zeroT[0] == HomogeneousPolynomial(0, 0)
    @test uT[0] == HomogeneousPolynomial(1, 0)
    @test ones(xH,1) == [1, xH+yH]
    @test typeof(ones(xH,2)) == Array{HomogeneousPolynomial{Int},1}
    @test length(ones(xH,2)) == 3
    @test ones(HomogeneousPolynomial{Complex{Int}},0) ==
        [HomogeneousPolynomial([complex(1,0)], 0)]
    @test !isnan(uT)
    @test TaylorSeries.fixorder(xH,yH) == (xH,yH)
    @test_throws AssertionError TaylorSeries.fixorder(zeros(xH,0)[1],yH)

    @test constant_term(xT) == 0
    @test constant_term(uT) == 1.0
    @test constant_term(xT) == constant_term(yT)
    @test constant_term(xH) == xH
    @test linear_polynomial(xT) == xT
    @test linear_polynomial(xT+yT) == xT+yT
    @test linear_polynomial(uT) == zero(yT)

    @test get_order(zeroT) == 1
    @test xT[1][1] == 1
    @test yH[2] == 1
    @test getcoeff(xT,(1,0)) == getcoeff(xT,[1,0]) == 1
    @test getcoeff(yH,(1,0)) == getcoeff(yH,[1,0]) == 0
    @test typeof(convert(HomogeneousPolynomial,1im)) ==
        HomogeneousPolynomial{Complex{Int}}
    @test convert(HomogeneousPolynomial,1im) ==
        HomogeneousPolynomial([complex(0,1)], 0)
    @test convert(HomogeneousPolynomial{Int},[1,1]) == xH+yH
    @test convert(HomogeneousPolynomial{Float64},[2,-1]) == 2.0xH-yH
    @test typeof(convert(TaylorN,1im)) == TaylorN{Complex{Int}}
    @test convert(TaylorN, 1im) ==
        TaylorN([HomogeneousPolynomial([complex(0,1)], 0)], 0)
    @test convert(TaylorN{Float64}, yH) == 1.0*yT
    @test convert(TaylorN{Float64}, [xH,yH]) == xT+1.0*yT
    @test convert(TaylorN{Int}, [xH,yH]) == xT+yT
    @test promote(xH, [1,1])[2] == xH+yH
    @test promote(xH, yT)[1] == xT
    @test promote(xT, [xH,yH])[2] == xT+yT
    @test typeof(promote(im*xT,[xH,yH])[2]) == TaylorN{Complex{Int}}

    @test iszero(zeroT.coeffs)
    @test iszero(zero(xH))
    @test !iszero(uT)
    @test iszero(zeroT)

    @test eltype(xH) == Int
    @test length(xH) == 2
    @test zero(xH) == 0*xH
    @test one(yH) == xH+yH
    @test xH * true == xH
    @test false * yH == zero(yH)
    @test get_order(yH) == 1
    @test get_order(xT) == 17
    @test xT * true == xT
    @test false * yT == zero(yT)
    @test HomogeneousPolynomial([1.0])*xH == xH

    @test xT == TaylorN([xH])
    @test one(xT) == TaylorN(1,5)
    @test TaylorN(uT) == convert(TaylorN{Complex},1)
    @test get_numvars() == 2
    @test length(uT) == get_order()+1
    @test eltype(convert(TaylorN{Complex{Float64}},1)) == Complex{Float64}

    @test 1+xT+yT == TaylorN(1,1) + TaylorN([xH,yH],1)
    @test xT-yT-1 == TaylorN([-1,xH-yH])
    @test xT*yT == TaylorN([HomogeneousPolynomial([0,1,0],2)])
    @test (1/(1-xT))[3] == HomogeneousPolynomial([1.0],3)
    @test xH^20 == HomogeneousPolynomial([0], get_order())
    @test (yT/(1-xT))[4] == xH^3 * yH
    @test mod(1+xT,1) == +xT
    @test (rem(1+xT,1))[0] == 0
    @test mod(1+xT,1.0) == +xT
    @test (rem(1+xT,1.0))[0] == 0
    @test abs(1-xT)  == 1-xT
    @test abs(-1-xT)  == 1+xT
    @test differentiate(yH,1) == differentiate(xH, :x₂)
    @test derivative(mod2pi(2pi+yT^3),2) == derivative(yT^3, :x₂)
    @test derivative(yT^3, :x₂) == derivative(yT^3, (0,1))
    @test derivative(yT) == zeroT == derivative(yT, (1,0))
    @test derivative((0,1), yT) == 1
    @test -xT/3im == im*xT/3
    @test (xH/3im)' == im*xH/3
    @test xT/BigInt(3) == TaylorN(BigFloat,1)/3
    @test xT/complex(0,BigInt(3)) == -im*xT/BigInt(3)
    @test (xH/complex(0,BigInt(3)))' ==
        im*HomogeneousPolynomial([BigInt(1),0])/3
    @test evaluate(xH) == zero(eltype(xH))
    @test xH() == zero(eltype(xH))
    @test xH([1,1]) == evaluate(xH, [1,1])
    @test xH((1,1)) == xH(1, 1.0) == evaluate(xH, (1, 1.0)) == 1
    hp = -5.4xH+6.89yH
    @test hp([1,1]) == evaluate(hp, [1,1])
    vr = rand(2)
    @test hp(vr) == evaluate(hp, vr)

    @test integrate(yH,1) == integrate(xH, :x₂)
    p = (xT-yT)^6
    @test integrate(derivative(p, 1), 1, yT^6) == p
    @test integrate(derivative(p, :x₁), :x₁, yT^6) == p
    @test derivative(integrate(p, 2), 2) == p
    @test derivative(integrate(p, :x₂), :x₂) == p
    @test derivative(TaylorN(1.0, get_order())) == TaylorN(0.0, get_order())
    @test integrate(TaylorN(6.0, get_order()), 1) == 6xT
    @test integrate(TaylorN(0.0, get_order()), 2) == TaylorN(0.0, get_order())
    @test integrate(TaylorN(0.0, get_order()), 2, xT) == xT
    @test integrate(TaylorN(0.0, get_order()), :x₂, xT) == xT
    @test integrate(xT^17, 2) == TaylorN(0.0, get_order())
    @test integrate(xT^17, 1, yT) == yT
    @test integrate(xT^17, 1, 2.0) == TaylorN(2.0, get_order())
    @test integrate(xT^17, :x₁, 2.0) == TaylorN(2.0, get_order())
    @test_throws AssertionError integrate(xT, 1, xT)
    @test_throws AssertionError integrate(xT, :x₁, xT)
    @test_throws AssertionError derivative(xT, (1,))
    @test_throws AssertionError derivative(xT, (1,2,3))
    @test_throws AssertionError derivative(xT, (-1,2))
    @test_throws AssertionError derivative((1,), xT)
    @test_throws AssertionError derivative((1,2,3), xT)
    @test_throws AssertionError derivative((-1,2), xT)

    @test derivative(2xT*yT^2, (8,8)) == 0
    @test derivative((8,8), 2xT*yT^2) == 0
    @test derivative(2xT*yT^2, 1) == 2yT^2
    @test derivative((1,0), 2xT*yT^2) == 0
    @test derivative(2xT*yT^2, (1,2)) == 4*one(yT)
    @test derivative((1,2), 2xT*yT^2) == 4
    @test xT*xT^3 == xT^4
    txy = 1.0 + xT*yT - 0.5*xT^2*yT + (1/3)*xT^3*yT + 0.5*xT^2*yT^2
    @test getindex((1+TaylorN(1))^TaylorN(2),0:4) == txy.coeffs[1:5]
    @test ( (1+TaylorN(1))^TaylorN(2) )[:] == ( (1+TaylorN(1))^TaylorN(2) ).coeffs[:]
    @test txy.coeffs[:] == txy[:]
    @test txy.coeffs[:] == txy[0:end]
    txy[:] .= ( -1.0 + 3xT*yT - xT^2*yT + (4/3)*xT^3*yT + (1/3)*xT*yT^3 + 0.5*xT^2*yT^2 + 0.5*xT*yT^3 )[:]
    @test txy[:] == ( -1.0 + 3xT*yT - xT^2*yT + (4/3)*xT^3*yT + (1/3)*xT*yT^3 + 0.5*xT^2*yT^2 + 0.5*xT*yT^3 )[:]
    txy[2:end-1] .= ( 1.0 - xT*yT + 0.5*xT^2*yT - (2/3)*xT*yT^3 - 0.5*xT^2*yT^2  + 7*xT^3*yT )[2:end-1]
    @test txy[2:end-1] == ( 1.0 - xT*yT + 0.5*xT^2*yT - (2/3)*xT*yT^3 - 0.5*xT^2*yT^2  + 7*xT^3*yT )[2:end-1]

    a = -5.0 + sin(xT+yT^2)
    b = deepcopy(a)
    @test a[:] == a[0:end]
    @test a[:] == b[:]
    @test a[1:end] == b[1:end]
    @test a[end][:] == b[end][:]
    @test a[end][1:end] == b[end][1:end]
    a[end][:] .= rand.()
    rv = a[end][:]
    @test a[end][:] == rv
    @test a[end][:] != b[end][:]
    a[end][1:end] .= rand.()
    rv = a[end][1:end]
    @test a[end][1:end] == rv
    @test a[end][1:end] != b[end][1:end]
    @test a[0:2:end] == a.coeffs[1:2:end]
    a[0:1:end] .= 0.0
    @test a == zero(a)

    hp = HomogeneousPolynomial(1)^8
    rv1 = rand( length(hp) )
    hp[:] = rv1
    @test rv1 == hp[:]
    rv2 = rand( length(hp)-2 )
    hp[1:end-2] .= rv2
    @test hp[1:end-2] == rv2
    @test hp[end-1:end] == rv1[end-1:end]
    hp[3:4] .= 0.0
    @test hp[1:2] == rv2[1:2]
    @test hp[3:4] == zeros(2)
    @test hp[5:end-2] == rv2[5:end]
    @test hp[end-1:end] == rv1[end-1:end]
    hp[:] = 0.0
    @test hp[:] == zero(rv1)
    @test all(hp[end-1:1:end] .== 0.0)

    pol = sin(xT+yT*xT)+yT^2-(1-xT)^3
    q = deepcopy(pol)
    q[:] = 0.0
    @test get_order.(q[:]) == collect(0:q.order)
    @test q[:] == zero(q[:])
    q[:] .= pol.coeffs
    @test q == pol
    @test q[:] == pol[:]
    q[2:end-1] .= 0.0
    @test q[2:end-1] == zero.(q[2:end-1])
    @test q[1] == pol[1]
    @test q[end] == pol[end]
    # q[:] = pol.coeffs
    # zH0 = zero(HomogeneousPolynomial{Float64})
    q[:] = 1.0
    @test q[1] == HomogeneousPolynomial([1,0])
    @test q[2] == HomogeneousPolynomial([1,0,0])
    q[:] .= pol.coeffs
    q[2:end-1] = one.(q[2:end-1])
    @test q[2:end-1] == one.(q[2:end-1])
    @test q[2] == HomogeneousPolynomial([1,1,1])
    @test q[1] == pol[1]
    @test q[end] == pol[end]
    q[:] .= pol.coeffs
    zHall = zeros(HomogeneousPolynomial{Float64}, q.order)
    q[:] .= zHall
    @test q[:] == zHall
    q[:] .= pol.coeffs
    q[1:end-1] .= zHall[2:end-1]
    @test q[1:end-1] == zHall[2:end-1]
    q[:] .= pol.coeffs
    @test q[:] != zeros(q.order+1)
    q[:] .= zeros(q.order+1)
    @test q[:] == zeros(q.order+1)
    q[:] .= pol.coeffs
    q[1:end-1] .= zeros(q.order+1)[2:end-1]
    @test q != pol
    @test all(q[1:1:end-1] .== 0.0)
    @test q[1:end-1] == zeros(q.order+1)[2:end-1]
    @test q[0] == pol[0]
    @test q[end] == pol[end]
    q[:] .= pol.coeffs
    pol2 = cos(sin(xT)-yT^3*xT)-3yT^2+sqrt(1-xT)
    q[2:end-2] .= pol2.coeffs[3:end-2]
    @test q[0:1] == pol[0:1]
    @test q[2:end-2] == pol2[2:end-2]
    @test q[end-1:end] == pol[end-1:end]
    @test q[2:2:end-2] == pol2[2:2:end-2]
    @test q[end-1:1:end] == pol[end-1:1:end]
    q[end-2:2:end] .= [0.0, 0.0]
    @test q[end-2] == 0.0
    @test_throws AssertionError q[end-2:2:end] = [0.0, 0.0, 0.0]
    q[end-2:2:end] .= pol.coeffs[end-2:2:end]
    @test q[end-2] == pol[end-2]
    q[end-2:2:end] .= pol.coeffs[end-2:2:end]
    @test_throws AssertionError q[end-2:2:end] = pol.coeffs[end-1:2:end]

    @test_throws AssertionError yT^(-2)
    @test_throws AssertionError yT^(-2.0)
    @test (1+xT)^(3//2) == ((1+xT)^0.5)^3
    @test real(xH) == xH
    @test imag(xH) == zero(xH)
    @test conj(im*yH) == (im*yH)'
    @test conj(im*yT) == (im*yT)'
    @test real( exp(1im * xT)) == cos(xT)
    @test getcoeff(convert(TaylorN{Rational{Int}},cos(xT)),(4,0)) ==
        1//factorial(4)
    cr = convert(TaylorN{Rational{Int}},cos(xT))
    @test getcoeff(cr,(4,0)) == 1//factorial(4)
    @test imag((exp(yT))^(-1im)') == sin(yT)
    exy = exp( xT+yT )
    @test evaluate(exy) == 1
    @test evaluate(exy, 0.1im, 0.01im) == exp(0.11im)
    @test exy(0.1im, 0.01im) == exp(0.11im)
    @test evaluate(exy,(0.1im, 0.01im)) == exp(0.11im)
    @test exy((0.1im, 0.01im)) == exp(0.11im)
    @test exy(true, (0.1im, 0.01im)) == exp(0.11im)
    @test evaluate(exy, (0.1im, 0.01im), sorting=false) == exy(false, (0.1im, 0.01im))
    @test evaluate(exy, (0.1im, 0.01im), sorting=false) == exy(false, 0.1im, 0.01im)
    @test evaluate(exy,[0.1im, 0.01im]) == exp(0.11im)
    @test exy([0.1im, 0.01im]) == exp(0.11im)
    @test isapprox(evaluate(exy, (1,1)), eeuler^2)
    @test exy(:x₁, 0.0) == exp(yT)
    txy = tan(xT+yT)
    @test getcoeff(txy,(8,7)) == 929569/99225
    ptxy = xT + yT + (1/3)*( xT^3 + yT^3 ) + xT^2*yT + xT*yT^2
    @test getindex(tan(TaylorN(1)+TaylorN(2)),0:4) == ptxy.coeffs[1:5]
    @test evaluate(xH*yH, 1.0, 2.0) == (xH*yH)(1.0, 2.0) == 2.0
    @test evaluate(xH*yH, (1.0, 2.0)) == 2.0
    @test evaluate(xH*yH, [1.0, 2.0]) == 2.0
    @test ptxy(:x₁, -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2
    @test ptxy(:x₁ => -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2
    @test evaluate(ptxy, :x₁ => -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2
    @test evaluate(ptxy, :x₁, -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2
    v = zeros(Int, 2)
    @test evaluate!([xT, yT], ones(Int, 2), v) == nothing
    @test v == ones(2)
    @test evaluate!([xT, yT][1:2], ones(Int, 2), v) == nothing
    @test v == ones(2)
    A_TN = [xT 2xT 3xT; yT 2yT 3yT]
    @test evaluate(A_TN, ones(2)) == [1.0 2.0 3.0; 1.0 2.0 3.0]
    @test evaluate(A_TN) == [0.0 0.0 0.0; 0.0 0.0 0.0]
    @test A_TN() == [0.0  0.0  0.0; 0.0  0.0  0.0]
    @test (view(A_TN,:,:))() == [0.0 0.0 0.0; 0.0 0.0 0.0]
    t = Taylor1(10)
    @test A_TN([t,t^2]) == [t 2t 3t; t^2 2t^2 3t^2]
    @test view(A_TN, :, :)(ones(2)) == A_TN(ones(2))
    @test view(A_TN, :, 1)(ones(2)) == A_TN[:,1](ones(2))

    @test evaluate(sin(asin(xT+yT)), [1.0,0.5]) == 1.5
    @test evaluate(asin(sin(xT+yT)), [1.0,0.5]) == 1.5
    @test tan(atan(xT+yT)) == xT+yT
    @test atan(tan(xT+yT)) == xT+yT
    @test atan(sin(1+xT+yT), cos(1+xT+yT)) == atan(sin(1+xT+yT)/cos(1+xT+yT))
    @test constant_term(atan(sin(3pi/4+xT+yT), cos(3pi/4+xT+yT))) == 3pi/4
    @test asin(xT+yT) + acos(xT+yT) == pi/2

    @test -sinh(xT+yT) + cosh(xT+yT) == exp(-(xT+yT))
    @test  sinh(xT+yT) + cosh(xT+yT) == exp(xT+yT)
    @test evaluate(- sinh(xT+yT)^2 + cosh(xT+yT)^2 , rand(2)) == 1
    @test evaluate(- sinh(xT+yT)^2 + cosh(xT+yT)^2 , zeros(2)) == 1
    @test tanh(xT + yT + 0im) == -1im * tan((xT+yT)*1im)
    @test cosh(xT+yT) == real(cos(im*(xT+yT)))
    @test sinh(xT+yT) == imag(sin(im*(xT+yT)))

    xx = 1.0*zeroT
    TaylorSeries.add!(xx, 1.0*xT, 2yT, 1)
    @test xx[1] == HomogeneousPolynomial([1,2])
    TaylorSeries.add!(xx, 5.0, 0)
    @test xx[0] == HomogeneousPolynomial([5.0])
    TaylorSeries.add!(xx, -5.0, 1)
    @test xx[1] == zero(xx[1])
    TaylorSeries.subst!(xx, 1.0*xT, yT, 1)
    @test xx[1] == HomogeneousPolynomial([1,-1])
    TaylorSeries.subst!(xx, 5.0, 0)
    @test xx[0] == HomogeneousPolynomial([-5.0])
    TaylorSeries.subst!(xx, -5.0, 1)
    @test xx[1] == zero(xx[end])
    TaylorSeries.div!(xx, 1.0+xT, 1.0+xT, 0)
    @test xx[0] == HomogeneousPolynomial([1.0])
    TaylorSeries.pow!(xx, 1.0+xT, 0.5, 1)
    @test xx[1] == HomogeneousPolynomial([0.5,0.0])
    xx = 1.0*zeroT
    TaylorSeries.pow!(xx, 1.0+xT, 1.5, 0)
    @test xx[0] == HomogeneousPolynomial([1.0])
    TaylorSeries.pow!(xx, 1.0+xT, 1.5, 1)
    @test xx[1] == HomogeneousPolynomial([1.5,0.0])
    xx = 1.0*zeroT
    TaylorSeries.pow!(xx, 1.0+xT, 0, 0)
    @test xx[0] == HomogeneousPolynomial([1.0])
    TaylorSeries.pow!(xx, 1.0+xT, 1, 1)
    @test xx[1] == HomogeneousPolynomial([1.0,0.0])
    xx = 1.0*zeroT
    TaylorSeries.pow!(xx, 1.0+xT, 2, 0)
    @test xx[0] == HomogeneousPolynomial([1.0])
    TaylorSeries.pow!(xx, 1.0+xT, 2, 1)
    @test xx[1] == HomogeneousPolynomial([2.0,0.0])
    xx = 1.0*zeroT
    TaylorSeries.sqrt!(xx, 1.0+xT, 0)
    TaylorSeries.sqrt!(xx, 1.0+xT, 1)
    @test xx[0] == 1.0
    @test xx[1] == HomogeneousPolynomial([0.5,0.0])
    xx = 1.0*zeroT
    TaylorSeries.exp!(xx, 1.0*xT, 0)
    TaylorSeries.exp!(xx, 1.0*xT, 1)
    @test xx[0] == 1.0
    @test xx[1] == HomogeneousPolynomial([1.0,0.0])
    xx = 1.0*zeroT
    TaylorSeries.log!(xx, 1.0+xT, 0)
    TaylorSeries.log!(xx, 1.0+xT, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.sincos!(xx, cxx, 1.0*xT, 0)
    TaylorSeries.sincos!(xx, cxx, 1.0*xT, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    @test cxx[0] == 1.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.tan!(xx, 1.0*xT, cxx, 0)
    TaylorSeries.tan!(xx, 1.0*xT, cxx, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    @test cxx[0] == 0.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.asin!(xx, 1.0*xT, cxx, 0)
    TaylorSeries.asin!(xx, 1.0*xT, cxx, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    @test cxx[0] == 1.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.acos!(xx, 1.0*xT, cxx, 0)
    TaylorSeries.acos!(xx, 1.0*xT, cxx, 1)
    @test xx[0] == acos(0.0)
    @test xx[1] == HomogeneousPolynomial(-1.0,1)
    @test cxx[0] == 1.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.atan!(xx, 1.0*xT, cxx, 0)
    TaylorSeries.atan!(xx, 1.0*xT, cxx, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    @test cxx[0] == 1.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.sinhcosh!(xx, cxx, 1.0*xT, 0)
    TaylorSeries.sinhcosh!(xx, cxx, 1.0*xT, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    @test cxx[0] == 1.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)
    xx = 1.0*zeroT
    cxx = zero(xx)
    TaylorSeries.tanh!(xx, 1.0*xT, cxx, 0)
    TaylorSeries.tanh!(xx, 1.0*xT, cxx, 1)
    @test xx[0] == 0.0
    @test xx[1] == HomogeneousPolynomial(1.0,1)
    @test cxx[0] == 0.0
    @test cxx[1] == HomogeneousPolynomial(0.0,1)

    g1(x, y) = x^3 + 3y^2 - 2x^2 * y - 7x + 2
    g2(x, y) = y + x^2 - x^4
    f1 = g1(xT, yT)
    f2 = g2(xT, yT)
    @test TaylorSeries.gradient(f1) == [ 3*xT^2-4*xT*yT-TaylorN(7,0), 6*yT-2*xT^2 ]
    @test ∇(f2) == [2*xT - 4*xT^3, TaylorN(1,0)]
    @test TaylorSeries.jacobian([f1,f2], [2,1]) == TaylorSeries.jacobian( [g1(xT+2,yT+1), g2(xT+2,yT+1)] )
    jac = Array{Int}(undef, 2, 2)
    TaylorSeries.jacobian!(jac, [g1(xT+2,yT+1), g2(xT+2,yT+1)])
    @test jac == TaylorSeries.jacobian( [g1(xT+2,yT+1), g2(xT+2,yT+1)] )
    TaylorSeries.jacobian!(jac, [f1,f2], [2,1])
    @test jac == TaylorSeries.jacobian([f1,f2], [2,1])
    @test TaylorSeries.hessian( f1*f2 ) ==
        [derivative((2,0), f1*f2) derivative((1,1), (f1*f2));
         derivative((1,1), f1*f2) derivative((0,2), (f1*f2))] == [4 -7; -7 0]
    @test TaylorSeries.hessian( f1*f2, [xT, yT] ) ==
        [derivative(f1*f2, (2,0)) derivative((f1*f2), (1,1));
         derivative(f1*f2, (1,1)) derivative((f1*f2), (0,2))]
    @test [xT yT]*TaylorSeries.hessian(f1*f2)*[xT, yT] == [ 2*TaylorN((f1*f2)[2]) ]
    @test TaylorSeries.hessian(f1^2)/2 == [ [49,0] [0,12] ]
    @test TaylorSeries.hessian(f1-f2-2*f1*f2) == (TaylorSeries.hessian(f1-f2-2*f1*f2))'
    @test TaylorSeries.hessian(f1-f2,[1,-1]) == TaylorSeries.hessian(g1(xT+1,yT-1)-g2(xT+1,yT-1))
    hes = Array{Int}(undef, 2, 2)
    TaylorSeries.hessian!(hes, f1*f2)
    @test hes == TaylorSeries.hessian(f1*f2)
    @test [xT yT]*hes*[xT, yT] == [ 2*TaylorN((f1*f2)[2]) ]
    TaylorSeries.hessian!(hes, f1^2)
    @test hes/2 == [ [49,0] [0,12] ]
    TaylorSeries.hessian!(hes, f1-f2-2*f1*f2)
    @test hes == hes'
    hes1 = Array{Int}(undef, 2, 2)
    TaylorSeries.hessian!(hes1, f1-f2,[1,-1])
    TaylorSeries.hessian!(hes, g1(xT+1,yT-1)-g2(xT+1,yT-1))
    @test hes1 == hes

    use_show_default(true)
    aa = sqrt(2) * xH
    ab = sqrt(2) * TaylorN(2, order=1)
    @test string(aa) ==
        "HomogeneousPolynomial{Float64}([1.4142135623730951, 0.0], 1)"
    @test string(ab) ==
        "TaylorN{Float64}(HomogeneousPolynomial{Float64}" *
        "[HomogeneousPolynomial{Float64}([0.0], 0), " *
        "HomogeneousPolynomial{Float64}([0.0, 1.4142135623730951], 1)], 1)"
    @test string([aa, aa]) ==
        "HomogeneousPolynomial{Float64}[HomogeneousPolynomial{Float64}" *
        "([1.4142135623730951, 0.0], 1), HomogeneousPolynomial{Float64}" *
        "([1.4142135623730951, 0.0], 1)]"
    @test string([ab, ab]) == "TaylorN{Float64}[TaylorN{Float64}" *
        "(HomogeneousPolynomial{Float64}[HomogeneousPolynomial{Float64}([0.0], 0), " *
        "HomogeneousPolynomial{Float64}([0.0, 1.4142135623730951], 1)], 1), " *
        "TaylorN{Float64}(HomogeneousPolynomial{Float64}[HomogeneousPolynomial{Float64}" *
        "([0.0], 0), HomogeneousPolynomial{Float64}([0.0, 1.4142135623730951], 1)], 1)]"
    use_show_default(false)
    @test string(aa) == " 1.4142135623730951 x₁"
    @test string(ab) == " 1.4142135623730951 x₂ + 𝒪(‖x‖²)"
    displayBigO(false)
    @test string(-xH) == " - 1 x₁"
    @test string(xT^2) == " 1 x₁²"
    @test string(1im*yT) == " ( 1 im ) x₂"
    @test string(xT-im*yT) == "  ( 1 ) x₁ - ( 1 im ) x₂"
    @test string([ab, ab]) ==
        "TaylorN{Float64}[ 1.4142135623730951 x₂,  1.4142135623730951 x₂]"
    displayBigO(true)
    @test string(-xH) == " - 1 x₁"
    @test string(xT^2) == " 1 x₁² + 𝒪(‖x‖¹⁸)"
    @test string(1im*yT) == " ( 1 im ) x₂ + 𝒪(‖x‖¹⁸)"
    @test string(xT-im*yT) == "  ( 1 ) x₁ - ( 1 im ) x₂ + 𝒪(‖x‖¹⁸)"

    @test_throws ArgumentError abs(xT)
    @test_throws AssertionError 1/x
    @test_throws AssertionError zero(x)/zero(x)
    @test_throws ArgumentError sqrt(x)
    @test_throws AssertionError x^(-2)
    @test_throws ArgumentError log(x)
    @test_throws AssertionError cos(x)/sin(y)
    @test_throws BoundsError xH[20]
    @test_throws BoundsError xT[20]

    a = 3x + 4y +6x^2 + 8x*y
    @test typeof( norm(x) ) == Float64
    @test norm(x) > 0
    @test norm(a) == norm([3,4,6,8.0])
    @test norm(a, 4) == sum([3,4,6,8.0].^4)^(1/4.)
    @test norm(a, Inf) == 8.0
    @test norm((3.0 + 4im)*x) == abs(3.0 + 4im)

    @test TaylorSeries.rtoldefault(TaylorN{Int}) == 0
    @test TaylorSeries.rtoldefault(TaylorN{Float64}) == sqrt(eps(Float64))
    @test TaylorSeries.rtoldefault(TaylorN{BigFloat}) == sqrt(eps(BigFloat))
    @test TaylorSeries.real(TaylorN{Float64}) == TaylorN{Float64}
    @test TaylorSeries.real(TaylorN{Complex{Float64}}) == TaylorN{Float64}
    @test isfinite(a)
    @test a[0] ≈ a[0]
    @test a[1] ≈ a[1]
    @test a[2] ≈ a[2]
    @test a[3] ≈ a[3]
    @test a ≈ a
    @test a .≈ a
    b = deepcopy(a)
    b[2][3] = Inf
    @test !isfinite(b)
    b[2][3] = NaN
    @test !isfinite(b)
    b[2][3] = a[2][3]+eps()
    @test isapprox(a[2], b[2], rtol=eps())
    @test a ≈ b
    b[2][2] = a[2][2]+sqrt(eps())
    @test a[2] ≈ b[2]
    @test a ≈ b

    f11(a,b) = (a+b)^a - cos(a*b)*b
    f22(a) = (a[1] + a[2])^a[1] - cos(a[1]*a[2])*a[2]
    @test taylor_expand(f11, 1.0,2.0) == taylor_expand(f22, [1,2.0])
    @test evaluate(taylor_expand(x->x[1] + x[2], [1,2])) == 3.0
    f33(x,y) = 3x+y
    @test eltype(taylor_expand(f33,1,1)) == eltype(1)
    x,y = get_variables()
    xysq = x^2 + y^2
    update!(xysq,[1.0,-2.0])
    @test xysq == (x+1.0)^2 + (y-2.0)^2
    update!(xysq,[-1.0,2.0])
    @test xysq == x^2 + y^2

    #test function-like behavior for TaylorN
    @test exy() == 1
    @test exy([0.1im,0.01im]) == exp(0.11im)
    @test isapprox(exy([1,1]), eeuler^2)
    @test sin(asin(xT+yT))([1.0,0.5]) == 1.5
    @test asin(sin(xT+yT))([1.0,0.5]) == 1.5
    @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(rand(2)) == 1
    @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(zeros(2)) == 1
    #number of variables changed to 4...
    dx = set_variables("x", numvars=4, order=10)
    P = sin.(dx)
    v = [1.0,2,3,4]
    for i in 1:4
        @test P[i](v) == evaluate(P[i], v)
    end
    @test P.(fill(v, 4)) == fill(P(v), 4)
    F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))]
    Q = F(v+dx)
    @test Q.( fill(v, 4) ) == fill(Q(v), 4)
    vr = map(x->rand(4), 1:4)
    @test Q.(vr) == map(x->Q(x), vr)
    for i in 1:4
        @test P[i]() == evaluate(P[i])
        @test Q[i]() == evaluate(Q[i])
    end
    @test P() == evaluate.(P)
    @test P() == evaluate(P)
    @test Q() == evaluate.(Q)
    @test Q() == evaluate(Q)
    @test Q[1:3]() == evaluate(Q[1:3])

    dx = set_variables("x", numvars=4, order=10)
    for i in 1:4
        @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i]
        @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i]
    end
    p = sin(exp(dx[1]*dx[2])+dx[3]*dx[2])/(1.0+dx[4]^2)
    q = zero(p)
    TaylorSeries.deg2rad!(q, p, 0)
    @test q[0] == p[0]*(pi/180)
    # TaylorSeries.deg2rad!.(q, p, [1,3,5])
    # for i in [0,1,3,5]
    #     @test q[i] == p[i]*(pi/180)
    # end
    TaylorSeries.rad2deg!(q, p, 0)
    @test q[0] == p[0]*(180/pi)
    # TaylorSeries.rad2deg!.(q, p, [1,3,5])
    # for i in [0,1,3,5]
    #     @test q[i] == p[i]*(180/pi)
    # end
    xT = 5+TaylorN(Int, 1, order=10)
    yT = TaylorN(2, order=10)
    TaylorSeries.deg2rad!(yT, xT, 0)
    @test yT[0] == xT[0]*(pi/180)
    TaylorSeries.rad2deg!(yT, xT, 0)
    @test yT[0] == xT[0]*(180/pi)
end

@testset "Integrate for several variables" begin

    t, x, y = set_variables("t x y")

    @test integrate(t, 1) == 0.5*t^2
    @test integrate(t, 2) == t * x
    @test integrate(t, 3) == t * y
    @test integrate(x, 1) == t * x
    @test integrate(x, 2) == 0.5*x^2
    @test integrate(y, 2) == x * y

end


text:
cd(@__DIR__)

# Setup
using Pkg
!isdir("ImageMagick-Test") && mkdir("ImageMagick-Test")
Pkg.activate(joinpath(@__DIR__,"ImageMagick-Test"))
Pkg.develop(PackageSpec(path="/Users/ian/Documents/GitHub/FileIO.jl"))
Pkg.add("ImageMagick")
Pkg.instantiate()

!isdir("QuartzImageIO-Test") && mkdir("QuartzImageIO-Test")
Pkg.activate(joinpath(@__DIR__,"QuartzImageIO-Test"))
Pkg.develop(PackageSpec(path="/Users/ian/Documents/GitHub/FileIO.jl"))
Pkg.add("QuartzImageIO")
Pkg.instantiate()

!isdir("PNG-Test") && mkdir("PNG-Test")
Pkg.activate(joinpath(@__DIR__,"PNG-Test"))
Pkg.develop(PackageSpec(path="/Users/ian/Documents/GitHub/FileIO.jl"))
Pkg.develop(PackageSpec(path="/Users/ian/Documents/GitHub/PNG.jl"))
Pkg.instantiate()

Pkg.activate()

corebench = raw"""
    using BenchmarkTools
    tfirst = @elapsed begin
        using FileIO
        FileIO.save("tst.png", x)
    end
    rm("tst.png")
    println("Time to first save (including save): $(round(tfirst, digits=3)) seconds")

    print("Save @btime:")
    @btime FileIO.save("tst.png", x) teardown=(rm("tst.png"))
    FileIO.save("tst.png", x)

    print("Load @btime:")
    @btime FileIO.load("tst.png")
    rm("tst.png")
    println("")
    """

for imgstr in ["rand(RGB{N0f8}, 244, 244)", "rand(RGB{N0f8}, 20000, 4000)"]
    setup = """
    using ImageCore
    x = $imgstr
    """

    @info "Testing with $imgstr image -----------------------------"

    for pkgname  in ["ImageMagick", "QuartzImageIO", "PNG"]
        @info "Testing FileIO with $pkgname (dedicated environment, new instance)"
        cd(joinpath(@__DIR__,"$pkgname-Test"))
        run(`julia --project=Project.toml --color=yes -e $(string(setup, corebench))`)
    end
end

ext: jl


text:
#==============================================================================#
# Budgets.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/budgets-2016-10-20.normal.json
#==============================================================================#

__precompile__()

module Budgets

using AWSCore

"""
    using AWSSDK.Budgets.create_budget
    create_budget([::AWSConfig], arguments::Dict)
    create_budget([::AWSConfig]; AccountId=, Budget=, <keyword arguments>)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "CreateBudget", arguments::Dict)
    budgets([::AWSConfig], "CreateBudget", AccountId=, Budget=, <keyword arguments>)

# CreateBudget Operation

Creates a budget and, if included, notifications and subscribers.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget.

## `Budget = [ ... ]` -- *Required*
The budget object that you want to create.
```
 Budget = [
        "BudgetName" => <required> ::String,
        "BudgetLimit" =>  [
            "Amount" => <required> ::String,
            "Unit" => <required> ::String
        ],
        "CostFilters" =>  ::Dict{String,String},
        "CostTypes" =>  [
            "IncludeTax" =>  ::Bool,
            "IncludeSubscription" =>  ::Bool,
            "UseBlended" =>  ::Bool,
            "IncludeRefund" =>  ::Bool,
            "IncludeCredit" =>  ::Bool,
            "IncludeUpfront" =>  ::Bool,
            "IncludeRecurring" =>  ::Bool,
            "IncludeOtherSubscription" =>  ::Bool,
            "IncludeSupport" =>  ::Bool,
            "IncludeDiscount" =>  ::Bool,
            "UseAmortized" =>  ::Bool
        ],
        "TimeUnit" => <required> "DAILY", "MONTHLY", "QUARTERLY" or "ANNUALLY",
        "TimePeriod" =>  [
            "Start" =>  timestamp,
            "End" =>  timestamp
        ],
        "CalculatedSpend" =>  [
            "ActualSpend" => <required> [
                "Amount" => <required> ::String,
                "Unit" => <required> ::String
            ],
            "ForecastedSpend" =>  [
                "Amount" => <required> ::String,
                "Unit" => <required> ::String
            ]
        ],
        "BudgetType" => <required> "USAGE", "COST", "RI_UTILIZATION" or "RI_COVERAGE"
    ]
```

## `NotificationsWithSubscribers = [[ ... ], ...]`
A notification that you want to associate with a budget. A budget can have up to five notifications, and each notification can have one SNS subscriber and up to ten email subscribers. If you include notifications and subscribers in your `CreateBudget` call, AWS creates the notifications and subscribers for you.
```
 NotificationsWithSubscribers = [[
        "Notification" => <required> [
            "NotificationType" => <required> "ACTUAL" or "FORECASTED",
            "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
            "Threshold" => <required> double,
            "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
        ],
        "Subscribers" => <required> [[
            "SubscriptionType" => <required> "SNS" or "EMAIL",
            "Address" => <required> ::String
        ], ...]
    ], ...]
```

# Returns

`CreateBudgetResponse`

# Exceptions

`InvalidParameterException`, `InternalErrorException`, `CreationLimitExceededException` or `DuplicateRecordException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/CreateBudget)
"""
@inline create_budget(aws::AWSConfig=default_aws_config(); args...) = create_budget(aws, args)

@inline create_budget(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "CreateBudget", args)

@inline create_budget(args) = create_budget(default_aws_config(), args)

"""
    using AWSSDK.Budgets.create_notification
    create_notification([::AWSConfig], arguments::Dict)
    create_notification([::AWSConfig]; AccountId=, BudgetName=, Notification=, Subscribers=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "CreateNotification", arguments::Dict)
    budgets([::AWSConfig], "CreateNotification", AccountId=, BudgetName=, Notification=, Subscribers=)

# CreateNotification Operation

Creates a notification. You must create the budget before you create the associated notification.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget that you want to create a notification for.

## `BudgetName = ::String` -- *Required*
The name of the budget that you want AWS to notified you about. Budget names must be unique within an account.

## `Notification = [ ... ]` -- *Required*
The notification that you want to create.
```
 Notification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

## `Subscribers = [[ ... ], ...]` -- *Required*
A list of subscribers that you want to associate with the notification. Each notification can have one SNS subscriber and up to ten email subscribers.
```
 Subscribers = [[
        "SubscriptionType" => <required> "SNS" or "EMAIL",
        "Address" => <required> ::String
    ], ...]
```

# Returns

`CreateNotificationResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException`, `CreationLimitExceededException` or `DuplicateRecordException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/CreateNotification)
"""
@inline create_notification(aws::AWSConfig=default_aws_config(); args...) = create_notification(aws, args)

@inline create_notification(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "CreateNotification", args)

@inline create_notification(args) = create_notification(default_aws_config(), args)

"""
    using AWSSDK.Budgets.create_subscriber
    create_subscriber([::AWSConfig], arguments::Dict)
    create_subscriber([::AWSConfig]; AccountId=, BudgetName=, Notification=, Subscriber=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "CreateSubscriber", arguments::Dict)
    budgets([::AWSConfig], "CreateSubscriber", AccountId=, BudgetName=, Notification=, Subscriber=)

# CreateSubscriber Operation

Creates a subscriber. You must create the associated budget and notification before you create the subscriber.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` associated with the budget that you want to create a subscriber for.

## `BudgetName = ::String` -- *Required*
The name of the budget that you want to subscribe to. Budget names must be unique within an account.

## `Notification = [ ... ]` -- *Required*
The notification that you want to create a subscriber for.
```
 Notification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

## `Subscriber = [ ... ]` -- *Required*
The subscriber that you want to associate with a budget notification.
```
 Subscriber = [
        "SubscriptionType" => <required> "SNS" or "EMAIL",
        "Address" => <required> ::String
    ]
```

# Returns

`CreateSubscriberResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `CreationLimitExceededException`, `DuplicateRecordException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/CreateSubscriber)
"""
@inline create_subscriber(aws::AWSConfig=default_aws_config(); args...) = create_subscriber(aws, args)

@inline create_subscriber(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "CreateSubscriber", args)

@inline create_subscriber(args) = create_subscriber(default_aws_config(), args)

"""
    using AWSSDK.Budgets.delete_budget
    delete_budget([::AWSConfig], arguments::Dict)
    delete_budget([::AWSConfig]; AccountId=, BudgetName=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DeleteBudget", arguments::Dict)
    budgets([::AWSConfig], "DeleteBudget", AccountId=, BudgetName=)

# DeleteBudget Operation

Deletes a budget. You can delete your budget at any time.

**Deleting a budget also deletes the notifications and subscribers associated with that budget.**

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget that you want to delete.

## `BudgetName = ::String` -- *Required*
The name of the budget that you want to delete.

# Returns

`DeleteBudgetResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DeleteBudget)
"""
@inline delete_budget(aws::AWSConfig=default_aws_config(); args...) = delete_budget(aws, args)

@inline delete_budget(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DeleteBudget", args)

@inline delete_budget(args) = delete_budget(default_aws_config(), args)

"""
    using AWSSDK.Budgets.delete_notification
    delete_notification([::AWSConfig], arguments::Dict)
    delete_notification([::AWSConfig]; AccountId=, BudgetName=, Notification=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DeleteNotification", arguments::Dict)
    budgets([::AWSConfig], "DeleteNotification", AccountId=, BudgetName=, Notification=)

# DeleteNotification Operation

Deletes a notification.

**Deleting a notification also deletes the subscribers associated with the notification.**

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget whose notification you want to delete.

## `BudgetName = ::String` -- *Required*
The name of the budget whose notification you want to delete.

## `Notification = [ ... ]` -- *Required*
The notification that you want to delete.
```
 Notification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

# Returns

`DeleteNotificationResponse`

# Exceptions

`InvalidParameterException`, `InternalErrorException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DeleteNotification)
"""
@inline delete_notification(aws::AWSConfig=default_aws_config(); args...) = delete_notification(aws, args)

@inline delete_notification(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DeleteNotification", args)

@inline delete_notification(args) = delete_notification(default_aws_config(), args)

"""
    using AWSSDK.Budgets.delete_subscriber
    delete_subscriber([::AWSConfig], arguments::Dict)
    delete_subscriber([::AWSConfig]; AccountId=, BudgetName=, Notification=, Subscriber=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DeleteSubscriber", arguments::Dict)
    budgets([::AWSConfig], "DeleteSubscriber", AccountId=, BudgetName=, Notification=, Subscriber=)

# DeleteSubscriber Operation

Deletes a subscriber.

**Deleting the last subscriber to a notification also deletes the notification.**

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget whose subscriber you want to delete.

## `BudgetName = ::String` -- *Required*
The name of the budget whose subscriber you want to delete.

## `Notification = [ ... ]` -- *Required*
The notification whose subscriber you want to delete.
```
 Notification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

## `Subscriber = [ ... ]` -- *Required*
The subscriber that you want to delete.
```
 Subscriber = [
        "SubscriptionType" => <required> "SNS" or "EMAIL",
        "Address" => <required> ::String
    ]
```

# Returns

`DeleteSubscriberResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DeleteSubscriber)
"""
@inline delete_subscriber(aws::AWSConfig=default_aws_config(); args...) = delete_subscriber(aws, args)

@inline delete_subscriber(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DeleteSubscriber", args)

@inline delete_subscriber(args) = delete_subscriber(default_aws_config(), args)

"""
    using AWSSDK.Budgets.describe_budget
    describe_budget([::AWSConfig], arguments::Dict)
    describe_budget([::AWSConfig]; AccountId=, BudgetName=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DescribeBudget", arguments::Dict)
    budgets([::AWSConfig], "DescribeBudget", AccountId=, BudgetName=)

# DescribeBudget Operation

Describes a budget.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget that you want a description of.

## `BudgetName = ::String` -- *Required*
The name of the budget that you want a description of.

# Returns

`DescribeBudgetResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DescribeBudget)
"""
@inline describe_budget(aws::AWSConfig=default_aws_config(); args...) = describe_budget(aws, args)

@inline describe_budget(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DescribeBudget", args)

@inline describe_budget(args) = describe_budget(default_aws_config(), args)

"""
    using AWSSDK.Budgets.describe_budgets
    describe_budgets([::AWSConfig], arguments::Dict)
    describe_budgets([::AWSConfig]; AccountId=, <keyword arguments>)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DescribeBudgets", arguments::Dict)
    budgets([::AWSConfig], "DescribeBudgets", AccountId=, <keyword arguments>)

# DescribeBudgets Operation

Lists the budgets associated with an account.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budgets that you want descriptions of.

## `MaxResults = ::Int`
Optional integer. Specifies the maximum number of results to return in response.

## `NextToken = ::String`
The pagination token that indicates the next set of results to retrieve.

# Returns

`DescribeBudgetsResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException`, `InvalidNextTokenException` or `ExpiredNextTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DescribeBudgets)
"""
@inline describe_budgets(aws::AWSConfig=default_aws_config(); args...) = describe_budgets(aws, args)

@inline describe_budgets(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DescribeBudgets", args)

@inline describe_budgets(args) = describe_budgets(default_aws_config(), args)

"""
    using AWSSDK.Budgets.describe_notifications_for_budget
    describe_notifications_for_budget([::AWSConfig], arguments::Dict)
    describe_notifications_for_budget([::AWSConfig]; AccountId=, BudgetName=, <keyword arguments>)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DescribeNotificationsForBudget", arguments::Dict)
    budgets([::AWSConfig], "DescribeNotificationsForBudget", AccountId=, BudgetName=, <keyword arguments>)

# DescribeNotificationsForBudget Operation

Lists the notifications associated with a budget.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget whose notifications you want descriptions of.

## `BudgetName = ::String` -- *Required*
The name of the budget whose notifications you want descriptions of.

## `MaxResults = ::Int`
Optional integer. Specifies the maximum number of results to return in response.

## `NextToken = ::String`
The pagination token that indicates the next set of results to retrieve.

# Returns

`DescribeNotificationsForBudgetResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException`, `InvalidNextTokenException` or `ExpiredNextTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DescribeNotificationsForBudget)
"""
@inline describe_notifications_for_budget(aws::AWSConfig=default_aws_config(); args...) = describe_notifications_for_budget(aws, args)

@inline describe_notifications_for_budget(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DescribeNotificationsForBudget", args)

@inline describe_notifications_for_budget(args) = describe_notifications_for_budget(default_aws_config(), args)

"""
    using AWSSDK.Budgets.describe_subscribers_for_notification
    describe_subscribers_for_notification([::AWSConfig], arguments::Dict)
    describe_subscribers_for_notification([::AWSConfig]; AccountId=, BudgetName=, Notification=, <keyword arguments>)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "DescribeSubscribersForNotification", arguments::Dict)
    budgets([::AWSConfig], "DescribeSubscribersForNotification", AccountId=, BudgetName=, Notification=, <keyword arguments>)

# DescribeSubscribersForNotification Operation

Lists the subscribers associated with a notification.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget whose subscribers you want descriptions of.

## `BudgetName = ::String` -- *Required*
The name of the budget whose subscribers you want descriptions of.

## `Notification = [ ... ]` -- *Required*
The notification whose subscribers you want to list.
```
 Notification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

## `MaxResults = ::Int`
Optional integer. Specifies the maximum number of results to return in response.

## `NextToken = ::String`
The pagination token that indicates the next set of results to retrieve.

# Returns

`DescribeSubscribersForNotificationResponse`

# Exceptions

`InternalErrorException`, `NotFoundException`, `InvalidParameterException`, `InvalidNextTokenException` or `ExpiredNextTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/DescribeSubscribersForNotification)
"""
@inline describe_subscribers_for_notification(aws::AWSConfig=default_aws_config(); args...) = describe_subscribers_for_notification(aws, args)

@inline describe_subscribers_for_notification(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "DescribeSubscribersForNotification", args)

@inline describe_subscribers_for_notification(args) = describe_subscribers_for_notification(default_aws_config(), args)

"""
    using AWSSDK.Budgets.update_budget
    update_budget([::AWSConfig], arguments::Dict)
    update_budget([::AWSConfig]; AccountId=, NewBudget=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "UpdateBudget", arguments::Dict)
    budgets([::AWSConfig], "UpdateBudget", AccountId=, NewBudget=)

# UpdateBudget Operation

Updates a budget. You can change every part of a budget except for the `budgetName` and the `calculatedSpend`. When a budget is modified, the `calculatedSpend` drops to zero until AWS has new usage data to use for forecasting.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget that you want to update.

## `NewBudget = [ ... ]` -- *Required*
The budget that you want to update your budget to.
```
 NewBudget = [
        "BudgetName" => <required> ::String,
        "BudgetLimit" =>  [
            "Amount" => <required> ::String,
            "Unit" => <required> ::String
        ],
        "CostFilters" =>  ::Dict{String,String},
        "CostTypes" =>  [
            "IncludeTax" =>  ::Bool,
            "IncludeSubscription" =>  ::Bool,
            "UseBlended" =>  ::Bool,
            "IncludeRefund" =>  ::Bool,
            "IncludeCredit" =>  ::Bool,
            "IncludeUpfront" =>  ::Bool,
            "IncludeRecurring" =>  ::Bool,
            "IncludeOtherSubscription" =>  ::Bool,
            "IncludeSupport" =>  ::Bool,
            "IncludeDiscount" =>  ::Bool,
            "UseAmortized" =>  ::Bool
        ],
        "TimeUnit" => <required> "DAILY", "MONTHLY", "QUARTERLY" or "ANNUALLY",
        "TimePeriod" =>  [
            "Start" =>  timestamp,
            "End" =>  timestamp
        ],
        "CalculatedSpend" =>  [
            "ActualSpend" => <required> [
                "Amount" => <required> ::String,
                "Unit" => <required> ::String
            ],
            "ForecastedSpend" =>  [
                "Amount" => <required> ::String,
                "Unit" => <required> ::String
            ]
        ],
        "BudgetType" => <required> "USAGE", "COST", "RI_UTILIZATION" or "RI_COVERAGE"
    ]
```

# Returns

`UpdateBudgetResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException` or `NotFoundException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/UpdateBudget)
"""
@inline update_budget(aws::AWSConfig=default_aws_config(); args...) = update_budget(aws, args)

@inline update_budget(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "UpdateBudget", args)

@inline update_budget(args) = update_budget(default_aws_config(), args)

"""
    using AWSSDK.Budgets.update_notification
    update_notification([::AWSConfig], arguments::Dict)
    update_notification([::AWSConfig]; AccountId=, BudgetName=, OldNotification=, NewNotification=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "UpdateNotification", arguments::Dict)
    budgets([::AWSConfig], "UpdateNotification", AccountId=, BudgetName=, OldNotification=, NewNotification=)

# UpdateNotification Operation

Updates a notification.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget whose notification you want to update.

## `BudgetName = ::String` -- *Required*
The name of the budget whose notification you want to update.

## `OldNotification = [ ... ]` -- *Required*
The previous notification associated with a budget.
```
 OldNotification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

## `NewNotification = [ ... ]` -- *Required*
The updated notification to be associated with a budget.
```
 NewNotification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

# Returns

`UpdateNotificationResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException` or `DuplicateRecordException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/UpdateNotification)
"""
@inline update_notification(aws::AWSConfig=default_aws_config(); args...) = update_notification(aws, args)

@inline update_notification(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "UpdateNotification", args)

@inline update_notification(args) = update_notification(default_aws_config(), args)

"""
    using AWSSDK.Budgets.update_subscriber
    update_subscriber([::AWSConfig], arguments::Dict)
    update_subscriber([::AWSConfig]; AccountId=, BudgetName=, Notification=, OldSubscriber=, NewSubscriber=)

    using AWSCore.Services.budgets
    budgets([::AWSConfig], "UpdateSubscriber", arguments::Dict)
    budgets([::AWSConfig], "UpdateSubscriber", AccountId=, BudgetName=, Notification=, OldSubscriber=, NewSubscriber=)

# UpdateSubscriber Operation

Updates a subscriber.

# Arguments

## `AccountId = ::String` -- *Required*
The `accountId` that is associated with the budget whose subscriber you want to update.

## `BudgetName = ::String` -- *Required*
The name of the budget whose subscriber you want to update.

## `Notification = [ ... ]` -- *Required*
The notification whose subscriber you want to update.
```
 Notification = [
        "NotificationType" => <required> "ACTUAL" or "FORECASTED",
        "ComparisonOperator" => <required> "GREATER_THAN", "LESS_THAN" or "EQUAL_TO",
        "Threshold" => <required> double,
        "ThresholdType" =>  "PERCENTAGE" or "ABSOLUTE_VALUE"
    ]
```

## `OldSubscriber = [ ... ]` -- *Required*
The previous subscriber associated with a budget notification.
```
 OldSubscriber = [
        "SubscriptionType" => <required> "SNS" or "EMAIL",
        "Address" => <required> ::String
    ]
```

## `NewSubscriber = [ ... ]` -- *Required*
The updated subscriber associated with a budget notification.
```
 NewSubscriber = [
        "SubscriptionType" => <required> "SNS" or "EMAIL",
        "Address" => <required> ::String
    ]
```

# Returns

`UpdateSubscriberResponse`

# Exceptions

`InternalErrorException`, `InvalidParameterException`, `NotFoundException` or `DuplicateRecordException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/budgets-2016-10-20/UpdateSubscriber)
"""
@inline update_subscriber(aws::AWSConfig=default_aws_config(); args...) = update_subscriber(aws, args)

@inline update_subscriber(aws::AWSConfig, args) = AWSCore.Services.budgets(aws, "UpdateSubscriber", args)

@inline update_subscriber(args) = update_subscriber(default_aws_config(), args)

end # module Budgets

#==============================================================================#
# End of file
#==============================================================================#


text:
# Dantzig-Wolfe Decomposition
# Edward J. Xu
# 2019.5.28
# include("$(homedir())/Documents/Github/DantzigWolfeDecomposition/Case/SSCFLP/main.jl")
########################################################################################################################
push!(LOAD_PATH, "$(homedir())//Documents/Github/DantzigWolfeDecomposition/Case/SSCFLP")
cd("$(homedir())//Documents/Github/DantzigWolfeDecomposition/Case/SSCFLP")
using JuMP
# using CPLEX
using Gurobi
using LinearAlgebra
using Random
using DantzigWolfeOptim
include("FuncSub.jl")
include("FuncMas.jl")
include("FuncOptim.jl")
include("FuncData.jl")
# include("DantzigWolfeOptim.jl")
########################################################################################################################

function main(whi::Float64)
    ## 1,  Data
    (vec_c, mat_a, vec_b, matIndexSub, numSub, numXInSub, m_mat_h, vecRowMatD) = getData(whi)
    ## 2,  Optim
    vecWhiBranch = [1, 1]
    if whi == 1
        wheBranch = true
        vecWhiSub = [1, 1]
        vecWhiVar = [12, 6]
    elseif whi == 2
        wheBranch = true
        vecWhiSub = [2, 1]
        vecWhiVar = [6, 6]
    else
        wheBranch = false
        vecWhiSub = [1]
        vecWhiVar = [1]
        vecWhiBranch = [1]
    end
    doDantzigWolfeOptim(vec_c, mat_a, vec_b, matIndexSub, numSub, numXInSub,
        m_mat_h, vecRowMatD, wheBranch, vecWhiSub, vecWhiVar, vecWhiBranch
        )
end

main(7.2)


text:
@testset "Cycles" begin
    completedg = CompleteDiGraph(4)
    pathdg = PathDiGraph(5)
    triangle = random_regular_graph(3, 2)
    quadrangle = random_regular_graph(4, 2)
    pentagon = random_regular_graph(5, 2)

    for g in testdigraphs(pathdg)
        @test maxsimplecycles(g) == 0
        @test maxsimplecycles(g, false) == 84
        @test simplecyclescount(g) == 0
        @test length(simplecycles(g)) == 0
        @test isempty(simplecycles(g)) == true
        @test isempty(simplecycles_iter(g)) == true
        @test simplecycleslength(g) == (zeros(5), 0)
        @test simplecyclescount(g, 10) == 0
        @test isempty(simplecycles_iter(g, 10)) == true
        @test simplecycleslength(g, 10) == (zeros(5), 0)
    end

    @test maxsimplecycles(4) == 20

    for g in testdigraphs(completedg)
        @test maxsimplecycles(g) == 20
        @test length(simplecycles(g)) == 20
        @test simplecycles(g) == @inferred(simplecycles_iter(g))
        @test simplecyclescount(g) == 20
        @test simplecycleslength(g) == ([0, 6, 8, 6], 20)
        @test simplecyclescount(g, 10) == 10
        @test simplecycleslength(g, 10)[2] == 10
    end

    for g in testgraphs(triangle)
        trianglelengths, triangletotal = simplecycleslength(DiGraph(g))
        @test sum(trianglelengths) == triangletotal
    end

    for g in testgraphs(quadrangle)
        quadranglelengths, quadrangletotal = simplecycleslength(DiGraph(g))
        @test sum(quadranglelengths) == quadrangletotal
        @test simplecycles(DiGraph(g)) == @inferred(simplecycles_iter(DiGraph(g)))
    end

    for g in testgraphs(pentagon)
        pentagonlengths, pentagontotal = simplecycleslength(DiGraph(g))
        @test sum(pentagonlengths) == pentagontotal
    end

    selfloopg = DiGraph([
        0 1 0 0;
        0 0 1 0;
        1 0 1 0;
        0 0 0 1;
    ])

    for g in testdigraphs(selfloopg)
        cycles = simplecycles(g)
        @test [3] in cycles
        @test [4] in cycles
        @test [1, 2, 3] in cycles || [2, 3, 1] in cycles || [3, 1, 2] in cycles
        @test length(cycles) == 3

        cycles2 = simplecycles_iter(g)
        @test [3] in cycles2
        @test [4] in cycles2
        @test [1, 2, 3] in cycles2 || [2, 3, 1] in cycles2 || [3, 1, 2] in cycles2
        @test length(cycles2) == 3
    end
end


text:
function state_likelihood_model(Eₘ_list)
    # Generate the model matrix used to calculate likelihoods
    # The model matrix depends on the measurement operators
    sum(abs,data(sum(Eₘ_list))-I)<1E-14 ||
        throw(ArgumentError("Eₘ operators do not form a valid POVM!"))
    return reduce(vcat,transpose.(vec.(Array.(Eₘ_list))))
end

function mle_state_tomo(M,A)
    (nm,d²) = size(A); d = isqrt(d²)
    length(M) == nm || throw(ArgumentError("number of measurements does not match model matrix size!"))
    d^2 == d² || throw(ArgumentError("invalid number of columns in model matrix!"))
    # create function with autodiff
    t0 = [(k=(i-2d)÷d+1; k*(2d-k)+1 == i ? 1/√d : 0.0) for i=1:d^2]
    f = OnceDifferentiable(t -> loglikelihood_statetomo(t,M,A), t0; autodiff = :forward)
    # optimize parameters such that norm(t) == 1 to keep trace == 1
    return optimize(f, t0, ConjugateGradient(manifold=Optim.Sphere()), Optim.Options())
end

function loglikelihood_statetomo(t,M,A)
    # Binomial statistics for the measurement count probability, up to some irrelevant constant
    X = t_expectation_values(t,A)
    return sum(x -> ((p,m)=x; -m*log(real(p))), zip(X,M))
end

function t_expectation_values(t::AbstractVector,A::Matrix)
    # assemble ρ from the t-vector
    ρ = build_density_matrix(t)
    # calculate probabilities by multiplying the vectorized density operator with the model matrix
    return A*vec(ρ)
end

function build_density_matrix(t::AbstractVector)
    # reconstruct the density matrix from the t-parameters
    T = build_T_matrix(t)
    # ρ = T'*T
    return T'*T
end

function build_T_matrix(t::AbstractVector{S}) where S<:Real
    # build the T matrix from the t-parameters
    d = isqrt(length(t))
    d^2 == length(t) || throw(ArgumentError("invalid t-vector length!"))
    T = zeros(Complex{S},d,d)
    r,c,k = d,0,1
    @inbounds while k <= d^2
        r += 1
        if r > d
            c = c + 1
            r = c
            T[r,c] = complex(t[k])
            k += 1
        else
            T[r,c] = complex(t[k],t[k+1])
            k += 2
        end
    end
    return T
end


text:
using Images, MLDatasets, MXNet, ImageView

train_x, train_y = CIFAR10.traindata();
test_x,  test_y  = CIFAR10.testdata();

# preview the images from train_x dataset
preview_img = zeros((0..., size(train_x, 1, 3)...));

for i = 1:10
    preview_img = vcat(preview_img, train_x[:, :, :, i])
end

imshow(colorview(RGB, permutedims(preview_img, (3, 2, 1))))

### PREPARING THE DATASET

train_x = train_x ./ 1;
test_x = test_x ./ 1;

train_length = 40000
validation_length = 10000

train_data_array  = mx.zeros((size(train_x, 1, 2, 3)..., train_length...));
train_label_array = mx.zeros(train_length);

validation_data_array  = mx.zeros((size(train_x, 1, 2, 3)..., validation_length...));
validation_label_array = mx.zeros(validation_length);

test_data_array  = mx.zeros((size(train_x, 1, 2, 3)..., size(test_x, 4)...));
test_label_array = mx.zeros(size(test_x, 4));

for idx = 1:train_length
    train_data_array[idx:idx] = reshape(train_x[:, :, :, idx], (size(train_x, 1, 2,3 )..., 1...))
    train_label_array[idx:idx] = train_y[idx]
end

for idx = 1:validation_length
    validation_data_array[idx:idx] = reshape(train_x[:, :, :, train_length + idx], (size(train_x, 1, 2,3 )..., 1...))
    validation_label_array[idx:idx] = train_y[train_length + idx]
end

for idx = 1:size(test_x, 4)
    test_data_array[idx:idx] = reshape(test_x[:, :, :, idx], (size(test_x, 1, 2,3 )..., 1...))
    test_label_array[idx:idx] = test_y[idx]
end

train_data_provider = mx.ArrayDataProvider(:data => train_data_array, :label => train_label_array, batch_size = 100, shuffle = true);
validation_data_provider = mx.ArrayDataProvider(:data => validation_data_array, :label => validation_label_array, batch_size = 100, shuffle = true);
test_data_provider = mx.ArrayDataProvider(:data => test_data_array, :label => test_label_array, batch_size = 100);

### Simple 1-layer NN

arch = @mx.chain mx.Variable(:data) =>
  mx.Flatten() =>
  mx.FullyConnected(num_hidden=128) =>
  mx.Activation(name=:relu1, act_type=:relu) =>
  mx.FullyConnected(num_hidden=10) =>
  mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 50, initializer = mx.XavierInitializer());

### Simple 2-layer NN

arch = @mx.chain mx.Variable(:data) =>
  mx.FullyConnected(num_hidden=128) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=64) =>
  mx.Activation(act_type=:relu) =>
  mx.FullyConnected(num_hidden=10) =>
  mx.Activation(act_type=:relu) =>
  mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = validation_data_provider, n_epoch = 50, initializer = mx.XavierInitializer());

## Convolution Neural Network

arch = @mx.chain mx.Variable(:data) =>
        mx.Convolution(kernel=(3, 3), num_filter=32) =>
        mx.Activation(act_type=:relu) =>
        mx.Dropout(p = 0.25) =>
        mx.Pooling( kernel=(2, 2), pool_type=:max) =>
        mx.Flatten() =>
        mx.FullyConnected(num_hidden=256) =>
        mx.Activation(act_type=:relu) =>
        mx.FullyConnected(num_hidden=10) =>
        mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = test_data_provider, n_epoch = 50, initializer = mx.XavierInitializer());


text:
function case_count_obs(trace::Trace)
    T = get_args(trace)[1]
    addrs = Gen.select([:step => t => :obs => :infected for t in 1:T]...,
                       [:step => t => :obs => :recovered for t in 1:T]...)
    observations = Gen.get_selected(get_choices(trace), addrs)
    return observations
end

function location_obs(trace::Trace; frac::Float64=0.5)
    model, T = get_retval(trace), get_args(trace)[1]
    tracked = 1:Int(round(frac * nagents(model)))
    addrs = Gen.select([:step => t => :obs => :location => agt
                        for t in 1:T for agt in tracked]...)
    observations = Gen.get_selected(get_choices(trace), addrs)
    return observations
end

function migration_obs(trace::Trace)
    model, T = get_retval(trace), get_args(trace)[1]
    m_addrs = Gen.select(
        [:step => t => :agents => :step! => :agent_step! => 1 =>
         agt => :migrate! => :m for t in 1:T for agt in 1:nagents(model)]...
    )
    observations = Gen.get_selected(get_choices(trace), m_addrs)
    return observations
end

function obs_dataframe(tr::Trace)
    model, T = get_retval(tr), get_args(tr)[1]
    data = reduce(vcat, tr[:step => t => :obs]' for t in 1:T)
    cnames = reduce(vcat, [[Symbol(:I, i), Symbol(:R, i)] for i in 1:model.C])
    return DataFrame(data, cnames)
end

function plot_obs!(tr::Trace, cities=nothing, legend=false; kwargs...)
    model = get_retval(tr)
    cities = isnothing(cities) ? collect(1:model.C) : cities
    cnames = reduce(vcat, [[Symbol(:I, i), Symbol(:R, i)] for i in cities])
    df = obs_dataframe(tr)[:, cnames]
    if !haskey(kwargs, :ls)
        kwargs = Dict{Symbol,Any}(kwargs...)
        kwargs[:ls] = reduce(hcat, [[:solid :dash] for i in cities])
    end
    labels = legend == false ? fill("", length(cnames)) : names(df)
    plot!(Matrix(df); lw=2, labels=permutedims(labels),
          legend=legend, ylims=[0, maximum(model.Ns)*1.5], kwargs...)
end

function plot_obs(tr::Trace, cities=nothing, legend=:topleft; kwargs...)
    plot(); plot_obs!(tr, cities, legend; kwargs...)
end

function plot_obs!(trs::AbstractVector{<:Trace}, cities=nothing; kwargs...)
    for tr in trs plot_obs!(tr, cities, false; kwargs...) end
    return plot!()
end

function plot_obs!(trs::AbstractVector{<:Trace}, ws::AbstractVector{Float64},
                  cities=nothing; kwargs...)
    max_w = maximum(ws)
    for (tr, w) in zip(trs, ws)
        plot_obs!(tr, cities, false; alpha=w/max_w, kwargs...)
    end
    return plot!()
end

function plot_obs_3d!(tr::Trace, cities=nothing; kwargs...)
    model = get_retval(tr)
    cities = isnothing(cities) ? collect(1:model.C) : cities
    cnames = reduce(vcat, [[Symbol(:I, i), Symbol(:R, i)] for i in cities])
    df = obs_dataframe(tr)[:, cnames]
    T = size(df)[1]
    for (z, (iname, rname)) in enumerate(partition(cnames, 2))
        plot!(1:T, fill(z, T), df[!, iname]; label=(z==1 ? "Infected" : ""),
              lw=2, ls=:solid, color=:salmon, alpha=1-(z-1)*0.25)
        plot!(1:T, fill(z, T), df[!, rname]; label=(z==1 ? "Recovered" : ""),
              lw=2, ls=:dash, color=:royalblue, alpha=1-(z-1)*0.25)
    end
    ylims!(0.5, model.C+0.5)
    yticks!(1:model.C, ["City $i" for i in 1:model.C])
    return plot!(camera=(15,15), size=(500,300), legend=:topleft, kwargs...)
end


text:
using MacroTools: @forward, postwalk

"""
    TopoChain(topo::FuncTopo, layers...)

Similar to a Flux.Chain, with the addition of the use of an FuncTopo to define the order/structure of the functions called.

# Example
```jldoctest
julia> topo = @functopo x:x => a:x => b:(a, b) => c => o

julia> model = TopoChain(topo,
                Dense(32, 64),
                Dense(32, 64),
                (x, y) -> x .* y,
                Dropout(0.1))

TopoChain(Dense(32, 64), Dense(32, 64), #5, Dropout(0.1)) representing the following function composition:
function(x)
    a = Dense(32, 64)(x)
    b = Dense(32, 64)(x)
    c = #5(a, b)
    o = Dropout(0.1)(c)
    o
end
```
"""
struct TopoChain{T<:Tuple, FS}
    models::T
    topo::FuncTopo{FS}
    TopoChain(topo::FuncTopo{FS}, xs...) where FS = new{typeof(xs), FS}(xs, topo)
end

@generated function (s::TopoChain{TP, FS})(xs...) where {TP, FS}
    _code = functopo_impl(FS)
    n = fieldcount(TP)
    ms = [Symbol(:__model, i, :__) for i = 1:n]
    head = Expr(:(=), Expr(:tuple, ms...), :(s.models))
    pushfirst!(_code.args, head)
    code = postwalk(_code) do x
        if x isa Expr && x.head === :ref && x.args[1] === :model
            i = x.args[2]
            y = :($(ms[i]))
            return y
        else
            x
        end
    end
    return code
end

@forward TopoChain.models Base.getindex, Base.length

function Base.show(io::IO, s::TopoChain)
    print(io, "TopoChain(")
    join(io, s.models, ", ")
    print(io, ") representing the following function composition: \n")
    print_topo(io, s.topo; models=s.models)
end

"show the structure of the TopoChain function"
show_topochainfunc(s::TopoChain) = print_topo(s.topo; models=s.models)


text:
function +(a::String, b::String)
	a * b
end

function Echo(port)
	# port = 8080
	server = listen(port)
	println("listen @ " + string(port))

	while true
		conn = accept(server)
		@async begin
			try
				while true
					incoming = readline(conn)
					print("> ", incoming)
					write(conn, incoming)
					# incoming = replace(input, "\n", "")
					# if incoming == ":exit" || input == ":quit"
				end
			catch err
				print("error: $err")
			end

		end
	end
end

Echo(8080)

max_stars_repo_name: openwonk/echo

max_issues_repo_name: openwonk/echo

max_forks_repo_name: openwonk/echo


text:
# Run forward model

println("Time to load BM.jl")
@time using BITEModel
const BM=BITEModel
using Parameters
import PyPlot

# Parameters to adjust:
#######################

update_cache = true
verbose = true

glaciertype = [:svalbard, :itmix, :synth, :peninsula, :rgi][1]
plotyes = true

# a few test glaciers:
region_svalbard = 7
all_svalbard = readdlm("data/rgi-6.0-huss/outlines/xyzn/svalbard/svalbard.dat")
cutoff = (0,10^10) # km^2
westfonna_ll = [18.5, 79.7]
austfonna_ll = [20, 79.2]
svalbard_nrs = Int[all_svalbard[i,1] for i=1:size(all_svalbard,1) if
                   cutoff[1]<=all_svalbard[i,end]<=cutoff[2] && # sizes
                   # remove Austfonna and Westfonna:
                   !(all_svalbard[i,2]>westfonna_ll[1] && all_svalbard[i,3]>westfonna_ll[2]) &&
                   !(all_svalbard[i,2]>austfonna_ll[1] && all_svalbard[i,3]>austfonna_ll[2])]
with_radar = ["00027", "00408", "00409", "00427", "00428", "00429", "00434"][2:end] # 27 is on west or austfonna
with_good_radar = [409, 428] # ELFENBEINBREEN and SVEIGBREEN

brucebreen = BM.RGIGlacier("RGI60-07.01455", "Brucebreen") # area 6km^2, land, long lat: 17.291 78.485
ll = (17.291, 78.485)
rabotbreen = BM.RGIGlacier("RGI60-07.01478", "Rabotbreen") # area 70km^2, land
tunabreen = BM.RGIGlacier("RGI60-07.01458", "Tunabreen") # area 160km^2, tide-water
whatbreen = BM.RGIGlacier("RGI60-07.00434") # area 8km^2
gid = [whatbreen, rabotbreen, brucebreen, tunabreen][1]

# Run it
########
probs = [118]
errs = []
# nr =222
#         gid = BM.RGIGlacier(nr,region_svalbard)
#         gl,gb,pp,pm,pn,pl = BM.init_forward(gid, verbose;
#                                             update_cache=update_cache)
#         pm = BM.MPara(pm, iv_nye=true, iv_dist_exp2=0.01)
#         print("Time to run forward model:"); gc()
#         @time hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb =
#             BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

#         qtot1d, qd1d = postproc1d(gb, pp, pm)

faulty = [319,440,]
strange = [849,857,892,944] # dem wrong 3x, dem boring 944

for nr in svalbard_nrs
    #nr in merge(faulty,strange) && continue
    println("\n")
    try
        gid = BM.RGIGlacier(nr,region_svalbard)
        @show gid
        gl,gb,pp,pm,pn,pl = BM.init_forward(gid, verbose;
                                            update_cache=update_cache)
        pm = BM.MPara(pm, iv_nye=true, iv_dist_exp2=0.01)
        print("Time to run forward model:"); gc()
        @time fwdsol = BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

        BM.postproc1d!(fwdsol)

        if plotyes
            print("Time to plot:")
            p = BM.plot_run(fwdsol, reuse=true, plot2d=[:h,:iv,:ivlog][1])
            Plots.savefig(p, "figs/$(gid.id).png")
        end
    catch e
        e==InterruptException() && throw(e)
        @show e
        push!(errs, (nr, e))
    end
end


text:
@testset "singletons" begin
    @testset "alias" begin
        @test unconditional() == Unconditional()
        @test exact() == Exact()
    end

    @testset "show" begin
        @test sprint(show, Unconditional()) == "Unconditional"
        @test sprintcompact(Unconditional()) == "U"

        @test sprint(show, Exact()) == "Parallel"
        @test sprintcompact(Exact()) == "P"
    end
end

@testset "unspecifiedpr" begin
    us = UnspecifiedParallel()
    @test unspecifiedpr() === us
    @test unspecifiedpr(Unconditional()) === us
    @test unspecifiedpr(Unconditional(), Exact()) === us

    @test sprint(show, us) == "Unspecified{U,P}"
    @test sprint(show, MIME("text/plain"), us) ==
        "Parallel among unspecified treatment groups"
end

@testset "nevertreated" begin
    nt0 = NeverTreatedParallel([0], Unconditional(), Exact())
    nt0y = NeverTreatedParallel(Date(0), Unconditional(), Exact())
    nt1 = NeverTreatedParallel([0,1], Unconditional(), Exact())
    nt1y = NeverTreatedParallel([Date(0), Date(1)], Unconditional(), Exact())

    @testset "inner constructor" begin
        @test NeverTreatedParallel([1,1,0], Unconditional(), Exact()) == nt1
        @test_throws ErrorException NeverTreatedParallel([], Unconditional(), Exact())
        @test_throws MethodError NeverTreatedParallel([-0.5], Unconditional(), Exact())
    end

    @testset "without @formula" begin
        @test nevertreated(0) == nt0
        @test nevertreated([0,1]) == nt1
        @test nevertreated(0:1) == nt1
        @test nevertreated(Set([0,1])) == nt1
        @test nevertreated([0,1,1]) == nt1
        @test nevertreated(Date(0)) == nt0y
        @test nevertreated((Date(0), Date(1))) == nt1y

        @test nevertreated(0, Unconditional(), Exact()) == nt0
        @test nevertreated(0, c=Unconditional(), s=Exact()) == nt0
        @test nevertreated([0,1], Unconditional(), Exact()) == nt1
        @test nevertreated([0,1], c=Unconditional(), s=Exact()) == nt1
    end

    @testset "with @formula" begin
        @test nevertreated(term(0)) == nt0

        f = @formula(y ~ nevertreated(cb(0,1)))
        @test nevertreated(f.rhs.args_parsed[1]) == nt1

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(0)))
        t = parse_treat(f)[1]
        @test t.pr == nt0

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(cb(0,1))))
        t = parse_treat(f)[1]
        @test t.pr == nt1

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(0, Unconditional, Exact)))
        t = parse_treat(f)[1]
        @test t.pr == nt0

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(0, unconditional, exact)))
        t = parse_treat(f)[1]
        @test t.pr == nt0

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(0, unconditional)))
        t = parse_treat(f)[1]
        @test t.pr == nt0

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(0, exact)))
        t = parse_treat(f)[1]
        @test t.pr == nt0

        f = @formula(y ~ treat(g, ttreat(t, 0), nevertreated(exact, unconditional, 0)))
        t = parse_treat(f)[1]
        @test t.pr == nt0
    end

    @testset "show" begin
        @test sprint(show, nt0) == "NeverTreated{U,P}(0)"
        @test sprint(show, MIME("text/plain"), nt0) == """
            Parallel trends with any never-treated group:
              Never-treated groups: 0"""

        @test sprint(show, nt1) == "NeverTreated{U,P}(0, 1)"
        @test sprint(show, MIME("text/plain"), nt1) == """
            Parallel trends with any never-treated group:
              Never-treated groups: 0, 1"""
    end
end

@testset "notyettreated" begin
    ny0 = NotYetTreatedParallel([0], [0], Unconditional(), Exact())
    ny0y = NotYetTreatedParallel(Date(0), Date(0), Unconditional(), Exact())
    ny1 = NotYetTreatedParallel([0,1], [0], Unconditional(), Exact())
    ny1y = NotYetTreatedParallel([Date(0), Date(1)], [Date(0)], Unconditional(), Exact())
    ny2 = NotYetTreatedParallel([0,1], [0], Unconditional(), Exact())
    ny3 = NotYetTreatedParallel([0,1], [0,1], Unconditional(), Exact())

    @testset "inner constructor" begin
        @test NotYetTreatedParallel([1,1,0], [1,0,0], Unconditional(), Exact()) == ny3
        @test_throws ArgumentError NotYetTreatedParallel([0], [], Unconditional(), Exact())
        @test_throws ArgumentError NotYetTreatedParallel([], [0], Unconditional(), Exact())
        @test_throws ArgumentError NotYetTreatedParallel([-0.5], [-1], Unconditional(), Exact())
    end

    @testset "without @formula" begin
        @test notyettreated(0) == ny0
        @test notyettreated([0,1], 0) == ny1
        @test notyettreated(0:1, (0,)) == ny2
        @test notyettreated(Set([0,1]), 0:1) == ny3
        @test notyettreated([0,1,1], [0,1,1]) == ny3
        @test notyettreated(Date(0)) == ny0y
        @test notyettreated((Date(0), Date(1)), Date(0)) == ny1y

        @test notyettreated(0, 0, Unconditional(), Exact()) == ny0
        @test notyettreated(0, c=Unconditional(), s=Exact()) == ny0
        @test notyettreated([0,1], [0], Unconditional(), Exact()) == ny1
        @test notyettreated([0,1], 0, c=Unconditional(), s=Exact()) == ny1
    end

    @testset "with @formula" begin
        @test notyettreated(term(0)) == ny0

        f = @formula(y ~ notyettreated(cb(0,1), 0))
        @test notyettreated(f.rhs.args_parsed...) == ny1

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(0)))
        t = parse_treat(f)[1]
        @test t.pr == ny0

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(cb(0,1), 0)))
        t = parse_treat(f)[1]
        @test t.pr == ny2

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(cb(0,1), cb(0,1))))
        t = parse_treat(f)[1]
        @test t.pr == ny3

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(0, Unconditional, Exact)))
        t = parse_treat(f)[1]
        @test t.pr == ny0

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(0, unconditional, exact)))
        t = parse_treat(f)[1]
        @test t.pr == ny0

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(0, unconditional)))
        t = parse_treat(f)[1]
        @test t.pr == ny0

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(0, exact)))
        t = parse_treat(f)[1]
        @test t.pr == ny0

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(exact, unconditional, 0)))
        t = parse_treat(f)[1]
        @test t.pr == ny0

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(cb(0,1), 0, Unconditional, Exact)))
        t = parse_treat(f)[1]
        @test t.pr == ny2

        f = @formula(y ~ treat(g, ttreat(t, 0), notyettreated(exact, unconditional, cb(0,1), cb(0,1))))
        t = parse_treat(f)[1]
        @test t.pr == ny3
    end

    @testset "show" begin
        @test sprint(show, ny0) == "NotYetTreated{U,P}(0)"
        @test sprint(show, MIME("text/plain"), ny0) == """
            Parallel trends with any not-yet-treated group:
              Not-yet-treated groups: 0
              Treated since: 0"""

        @test sprint(show, ny1) == "NotYetTreated{U,P}(0, 1)"
        @test sprint(show, MIME("text/plain"), ny1) == """
            Parallel trends with any not-yet-treated group:
              Not-yet-treated groups: 0, 1
              Treated since: 0"""

        @test sprint(show, ny2) == "NotYetTreated{U,P}(0, 1)"
        @test sprint(show, MIME("text/plain"), ny2) == """
            Parallel trends with any not-yet-treated group:
              Not-yet-treated groups: 0, 1
              Treated since: 0"""

        @test sprint(show, ny3) == "NotYetTreated{U,P}(0, 1)"
        @test sprint(show, MIME("text/plain"), ny3) == """
            Parallel trends with any not-yet-treated group:
              Not-yet-treated groups: 0, 1
              Treated since: 0, 1"""
    end
end

@testset "istreated istreated!" begin
    nt = nevertreated(-1)
    @test istreated(nt, -1) == false
    @test istreated(nt, 0)

    out = BitVector(undef, 3)
    x = [1, -1, missing]
    @test istreated!(out, nt, x) == [true, false, true]
    x = ScaledArray([1, -1, missing], 1)
    @test istreated!(out, nt, x) == [true, false, true]
    x = ScaledArray([1, 0, missing], 1)
    @test istreated!(out, nt, x) == [true, true, true]

    ny = notyettreated(5)
    @test istreated(ny, 5) == false
    @test istreated(ny, 4)
    ny = notyettreated(5, 4)
    @test istreated(ny, 4)

    out = BitVector(undef, 3)
    x = [1, 5, missing]
    @test istreated!(out, ny, x) == [true, false, true]
    x = ScaledArray([1, 5, missing], 1)
    @test istreated!(out, ny, x) == [true, false, true]
    x = ScaledArray([1, 0, missing], 1)
    @test istreated!(out, ny, x) == [true, true, true]
end

@testset "termvars" begin
    @test termvars(Unconditional()) == Symbol[]
    @test termvars(Exact()) == Symbol[]
    @test termvars(unspecifiedpr()) == Symbol[]
    @test termvars(nevertreated(-1)) == Symbol[]
    @test termvars(notyettreated(5)) == Symbol[]

    @test_throws ErrorException termvars(TestParaCondition())
    @test_throws ErrorException termvars(TestParaStrength())
    @test_throws ErrorException termvars(PR)
end


text:
function instantiate(model_type::DataType, data, pools, param::Dict{Symbol, <:Any})
    debug(LOGGER, "Instantiating model of type $(model_type) with params $(param)")
    m = model_type(data, pools)
    set_param!(m, param)
    return m
end

function set_param!(model, param::Dict{Symbol, <:Any})
    has_changed = false
    for (k, v) in param
        k in keys(get_model_params(model)) || throw(ArgumentError("Model $(typeof(model)) has no parameter $k. Valid parameters are $(get_model_params(model))."))
        fieldtype = typeof(getfield(model, k))
        is_valid_param_value(model, Val{k}, v) || throw(ArgumentError("Parameter $k has value $v which is outside the valid range for model $(typeof(model))."))
        if getfield(model,k) != v
            has_changed = true
            setfield!(model, k, v)
        end
    end
    has_changed && invalidate_solution!(model)
    return nothing
end

function set_pools!(model, pools::Vector{Symbol})
    islabelenc(pools, learning_pool_enc) || throw(ArgumentError("Pools must be a LabelEncoding of $(learning_pool_enc.label)."))
    set_pools!(model, labelmap(pools))
    return nothing
end

function set_pools!(model, pools::Dict{Symbol, Vector{Int}})
    if !isdefined(model, :pools) ||model.pools != pools
        model.pools = pools
        invalidate_solution!(model)
    end
    return nothing
end

function set_adjust_K!(model, adjust_K::Bool)
    if adjust_K && !model.adjust_K
        model.adjust_K = adjust_K
        update_K!(model)
        invalidate_solution!(model)
    elseif !adjust_K && model.adjust_K
        model.adjust_K = adjust_K
        invalidate_solution!(model)
    end
    return nothing
end

is_K_adjusted(model)::Bool = model.adjust_K

calculate_kernel_matrix(model) = MLKernels.kernelmatrix(Val(:col), model.kernel_fct, model.data)

function calculate_kernel_matrix(model::SubSVDD)
     map(k -> MLKernels.kernelmatrix(Val(:col), model.kernel_fct[k], model.data[model.subspaces[k],:]), eachindex(model.subspaces))
 end

function update_K!(model)
    updated_K = calculate_kernel_matrix(model)
    if !isdefined(model, :K) || updated_K != model.K || (!isdefined(model, :K_adjusted) && model.adjust_K)
        debug(LOGGER, "[UPDATE_K] Updating Kernel matrix.")
        invalidate_solution!(model)
        model.K = updated_K
        if model.adjust_K
            model.K_adjusted = adjust_kernel_matrix(updated_K)
        end
    else
        debug(LOGGER, "[UPDATE_K] Nothing to update, kernel matrix has not changed.")
    end
    return nothing
end

function set_data!(model, data::Array{T, 2}) where T <: Real
    if model.data != data
        model.data = data
        update_K!(model)
        invalidate_solution!(model)
    end
    return nothing
end

function set_kernel!(model::OCClassifier, kernel_fct)
    if !isdefined(model, :kernel_fct) || model.kernel_fct != kernel_fct
        model.kernel_fct = kernel_fct
        update_K!(model)
        invalidate_solution!(model)
    end
    return nothing
end

function initialize!(model::OCClassifier, strategy::InitializationStrategy)
    debug(LOGGER, "[INITIALIZE] Initializing $(typeof(model)) with strategy: $(typeof(strategy))")
    C, kernel_fct = get_parameters(model, strategy)
    debug(LOGGER, "[INITIALIZE] Fitted parameters are: C = $(C) and kernel_fct = $(kernel_fct)")
    set_C!(model, C)
    set_kernel!(model, kernel_fct)
    model.state = model_initialized
    return nothing
end

get_kernel(model::OCClassifier) = model.kernel_fct

apply_update_strategy!(model, new_pools, query_ids, old_idx_remaining, new_idx_remaining) = nothing

function update_with_feedback!(model, new_data, new_pools::Vector{Symbol}, query_ids::Vector{Int},
    old_idx_remaining::Vector{Int},
    new_idx_remaining::Vector{Int})
    @assert size(new_data, 2) == length(new_pools)
    @assert islabelenc(new_pools, SVDD.learning_pool_enc)

    apply_update_strategy!(model, new_pools, query_ids, old_idx_remaining, new_idx_remaining)

    set_pools!(model, new_pools)
    set_data!(model, new_data)
    invalidate_solution!(model)
    return nothing
end


text:
module BlockAverage

export block_average

"""

```
block_average(x::AbstractVector{T}, by = mean) where T<:Real
```

Function that computes the block average and standard deviation as a function of block size for a time-dependent, possibly correlated, set of data.

Returns three vectors `avg`, `err`, and `sizes`.  The output `avg` and `err` vectors will contain the estimate of the mean value of the property and the *standard deviation of the mean,* for each possible block size (stored in `sizes`). Thus, `avg[i]` for `sizes[i] = N` will contain the mean of `by(x)` computed for blocks of `x` of length `N` (thus, with `length(x)÷N` samples).

The function defined by `by=` must be able to operate on a `view` of a subvector of `x`.  By default, `x` is assumed to contain the value of the property to be considered, and then `by = mean`. Nevertheless, one could compute other property from `x`.

Call with complete set of parameters:

```julia
avg, err, sizes = block_average(x::AbstractVector{T}; by = mean,
                                max_block_size::Int = length(x)÷10,
                                min_block_size::Int = 1 ) where T<:Real
```

## Examples

With `x = rand(1000)`:

1. `x` contains directly the property to be estimated:
```julia
avg, err, sizes = block_average(x)
```

2. The property to be estimated is the average of the square of each `x`:
```julia
func(x) = sum(x.^2)/length(x)
avg, err, sizes = block_average(x,by=func)
```

3. The property is the probabilty of finding an element of `x` greater than `0.5`:
```julia
func(x) = count( x .> 0.5 ) / length(x)
avg, err, sizes = block_average(x,by=func)
```

"""
function block_average(x::AbstractVector{T}; by = x -> sum(x)/length(x),
                       min_block_size::Int = 1,
                       max_block_size::Int = length(x)÷10) where T<:Real

  n = length(x)

  xmean = Vector{T}(undef,0)
  xerr = Vector{T}(undef,0)
  sizes = Vector{Int}(undef,0)

  ib = 0
  for block_size in min_block_size:max_block_size
    ib += 1
    nblocks = n÷block_size
    remaining = n%block_size

    # If the remaining is greater than the number of blocks, add
    # some points to each block to maximize the use of the data
    block_size_final = block_size
    if remaining >= nblocks
      block_size_final += remaining÷nblocks
    end

    # If the block size turned out to be the same as the last one
    # continue to next size
    if ib > 1 && block_size_final == sizes[end]
      continue
    end

    # Add new point to vectors
    push!(sizes,block_size_final)
    push!(xmean,zero(T))
    push!(xerr,zero(T))

    # Compute the mean of each block
    for i in 1:nblocks
      xblock = @view x[brange(i,block_size_final)]
      val = by(xblock)
      xmean[end] += val
    end

    # Compute the standard deviation of the mean (σ²/√N)
    if nblocks > 1
      xmean[end] /= nblocks
      for i in 1:nblocks
        xblock = @view x[brange(i,block_size_final)]
        val = by(xblock)
        xerr[end] += (val - xmean[end])^2
      end
      xerr[end] = sqrt(xerr[end]/((nblocks-1)*nblocks))
    end
  end

  return xmean, xerr, sizes
end

# Range of indices for block i of size block_size
function brange(i,block_size)
  ifirst = (i-1)*block_size+1
  ilast = ifirst + block_size - 1
  return ifirst:ilast
end

# Generate correlated data to test
function test_data(n)
  temperature = 1.
  x = Vector{Float64}(undef,n)
  x[1] = 0.
  u = 0.
  i = 1
  while i < n
    x_trial = x[i] - 0.1 + 0.2*rand()
    u_trial = x_trial^2
    if u_trial < u || exp((u-u_trial)/temperature) > 0.5
      i += 1
      x[i] = x_trial
      u = u_trial
    end
  end
  x
end

end


text:
function levenberg_marquardt(f::Function, g::Function, x0::Vector{Float64}; tolX=1e-8, tolG=1e-12,
                             maxIter=100, lambda=100.0, show_trace=false, quiet=false)
	# finds argmin sum(f(x).^2) using the Levenberg-Marquardt algorithm
	#          x
	# The function f should take an input vector of length n and return an output vector of length m
	# The function g is the Jacobian of f, and should be an m x n matrix
	# x0 is an initial guess for the solution
	# fargs is a tuple of additional arguments to pass to f
	# available options:
	#   tolX - search tolerance in x
	#   tolG - search tolerance in gradient
	#   maxIter - maximum number of iterations
	#   lambda - (inverse of) initial trust region radius
	#   show_trace - print a status summary on each iteration if true
	# returns: x, J
	#   x - least squares solution for x
	#   J - estimate of the Jacobian of f at x

	const MAX_LAMBDA = 1e16 # minimum trust region radius
	const MIN_LAMBDA = 1e-16 # maximum trust region radius
	const MIN_STEP_QUALITY = 1e-3
	const MIN_DIAGONAL = 1e-6 # lower bound on values of diagonal matrix used to regularize the trust region step
	converged = false
	iterCt = 0
	x = x0
	delta_x = copy(x0)
	f_calls = 0
	g_calls = 0
	fcur = f(x)
	f_calls += 1
	residual = norm(fcur)^2

	while ~converged && iterCt < maxIter
		J = g(x)
		g_calls += 1
		# we want to solve:
		#    argmin 0.5*||J(x)*delta_x + f(x)||^2 + lambda*||diagm(J'*J)*delta_x||^2
		# Solving for the minimum gives:
		#    (J'*J + lambda*DtD) * delta_x == -J^T * f(x), where DtD = diagm(sum(J.^2,1))
		# Where we have used the equivalence: diagm(J'*J) = diagm(sum(J.^2, 1))
		# It is additionally useful to bound the elements of DtD below to help
		# prevent "parameter evaporation".
		DtD = diagm(Float64[max(x, MIN_DIAGONAL) for x in sum(J.^2,1)])
		delta_x = ( J'*J + sqrt(lambda)*DtD ) \ -J'*fcur

		# if the linear assumption is valid, our new residual should be:
		predicted_residual = norm(J*delta_x + fcur)^2
		# check for numerical problems in solving for delta_x by ensuring that the predicted residual is smaller
		# than the current residual
		if predicted_residual > residual + 2max(eps(predicted_residual),eps(residual))
			warn("""Problem solving for delta_x: predicted residual increase.
                             $predicted_residual (predicted_residual) >
                             $residual (residual) + $(eps(predicted_residual)) (eps)""")
		end
		# try the step and compute its quality
		trial_f = f(x + delta_x)
		f_calls += 1
		trial_residual = norm(trial_f)^2
		# step quality = residual change / predicted residual change
		rho = (trial_residual - residual) / (predicted_residual - residual)
		if rho > MIN_STEP_QUALITY
			x += delta_x
			fcur = trial_f
			residual = trial_residual
			# increase trust region radius
			lambda = max(0.1*lambda, MIN_LAMBDA)
		else
			# decrease trust region radius
			lambda = min(10*lambda, MAX_LAMBDA)
		end
		iterCt += 1
		# check convergence criteria:
		# 1. Small gradient: norm(J^T * fcur, Inf) < tolG
		# 2. Small step size: norm(delta_x) < tolX
		converged = (norm(J' * fcur, Inf) < tolG) || (norm(delta_x) < tolX*(tolX + norm(x)))
	end
	x
end

function nlsfitworker(model::Function, y::Matrix{Float64},
                          x::Vector{Float64}, p0::Vector{Float64})
  # performs Levenberg-Marquardt fitting
  nparams = length(p0)
  n = size(y,2)
  params = zeros(nparams, n)
  resids = similar(y)
  for j in 1:n
    f(p) = model(x, p) - y[:,j]
    g = Calculus.jacobian(f, :forward)
    results = levenberg_marquardt(f, g, p0; tolX=1e-3, tolG=1e-4,
                                        lambda=10.0)
    resids[:,j] = f(results) / sqrt(norm(y[:,j]))
    params[:,j] = results
  end
  (params, resids)
end

function nlsfit(f::Function, y::Matrix{Float64}, idxs::Vector{Int},
                x::Vector{Float64}, p0::Vector{Float64})
  tic()
  nt, n = size(y)
  nw = nworkers()
  nidxs = length(idxs)
  params = zeros(length(p0), n)
  resids = zeros(nt, n)
  dof = nt - length(p0)
  if nw > 1
    reflist = Any[]
    idxperworker = round(Int, nidxs / nw)
    workeridxs = Any[idxs[(w-1)*idxperworker+1:min(w*idxperworker,nidxs)] for w in 1:nw]
    @dprint "fitting $nt x $idxperworker points on each of $nw workers"
    workerids = workers()
    for w in 1:nw
      r = @spawnat workerids[w] nlsfitworker(f, y[:,workeridxs[w]], x, p0)
      push!(reflist, r)
    end
    for w in 1:nw
      params[:,workeridxs[w]], resids[:,workeridxs[w]] = fetch(reflist[w])
    end
  else   # run on one CPU core
    @dprint "running $nt x $nidxs points on one CPU core"
    params[:,idxs], resids[:,idxs] = nlsfitworker(f, y[:,idxs], x, p0)
  end
  t = toq()
  vps = nidxs / t
  @dprint @sprintf "processed %d voxels in %.1f s (%.1f vox/s)\n" nidxs t vps
  (params, resids, dof)
end


text:
"""
    minimum_weight_perfect_matching(g, w::Dict{Edge,Real})
    minimum_weight_perfect_matching(g, w::Dict{Edge,Real}, cutoff)

Given a graph `g` and an edgemap `w` containing weights associated to edges,
returns a matching with the mimimum total weight among the ones containing
exactly `nv(g)/2` edges.

Edges in `g` not present in `w` will not be considered for the matching.

This function relies on the BlossomV.jl package, a julia wrapper
around Kolmogorov's BlossomV algorithm.

Eventually a `cutoff` argument can be given, to the reduce computational time
excluding edges with weights higher than the cutoff.

The returned object is of type `MatchingResult`.

In case of error try to change the optional argument `tmaxscale` (default is `tmaxscale=10`).
"""
function minimum_weight_perfect_matching end

function minimum_weight_perfect_matching(g::Graph, w::Dict{E,U}, cutoff, kws...) where {U<:Real, E<:Edge}
    wnew = Dict{E, U}()
    for (e, c) in w
        if c <= cutoff
            wnew[e] = c
        end
    end
    return minimum_weight_perfect_matching(g, wnew; kws...)
end

function minimum_weight_perfect_matching(g::Graph, w::Dict{E,U}; tmaxscale=10.) where {U<:AbstractFloat, E<:Edge}
    wnew = Dict{E, Int32}()
    cmax = maximum(values(w))
    cmin = minimum(values(w))
    tmax = typemax(Int32)  / tmaxscale # /10 is kinda arbitrary,
                                # hopefully high enough to not occur in overflow problems
    for (e, c) in w
        wnew[e] = round(Int32, (c-cmin) / (cmax-cmin) * tmax)
    end
    match = minimum_weight_perfect_matching(g, wnew)
    weight = zero(U)
    for i=1:nv(g)
        j = match.mate[i]
        if j > i
            weight += w[E(i,j)]
        end
    end
    return MatchingResult(weight, match.mate)
end

function minimum_weight_perfect_matching(g::Graph, w::Dict{E,U}) where {U<:Integer, E<:Edge}
    m = BlossomV.Matching(nv(g))
    for (e, c) in w
        BlossomV.add_edge(m, src(e)-1, dst(e)-1, c)
    end
    BlossomV.solve(m)

    mate = fill(-1, nv(g))
    totweight = zero(U)
    for i=1:nv(g)
        j = BlossomV.get_match(m, i-1) + 1
        mate[i] = j <= 0 ? -1 : j
        if i < j
            totweight += w[Edge(i,j)]
        end
    end
    return MatchingResult(totweight, mate)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Use
# const Negative{F} = ComposedFunction{typeof(-),F}
# instead of custom struct Negative
# Move ValueShapes.varshape(f::Negative) to ValueShapes

struct Negative{F<:Base.Callable} <: Function
    orig_f::F
end

@inline (f::Negative)(x) = - f.orig_f(x)

ValueShapes.varshape(f::Negative) = varshape(f.orig_f)
ValueShapes.unshaped(f::Negative) = Negative(unshaped(f.orig_f))

function Base.show(io::IO, f::Negative)
    print(io, Base.typename(typeof(f)).name, "(")
    show(io, f.orig_f)
    print(io, ")")
end

Base.show(io::IO, M::MIME"text/plain", f::Negative) = show(io, f)

function negative end
negative(f::Base.Callable) = Negative(f)
negative(f::Negative) = f.orig_f