text:
using Yields
using Test

@testset "Yields.jl" begin

    @testset "rate types" begin
        rs = Rate.([0.1, 0.02], Yields.Continuous())
        @test rs[1] == Rate(0.1, Yields.Continuous())
        @test rate(rs[1]) == 0.1
    end

    @testset "constructor" begin
        @test Yields.Continuous(0.05) == Rate(0.05, Yields.Continuous())
        @test Yields.Periodic(0.02, 2) == Rate(0.02, Yields.Periodic(2))
        @test Rate(0.02, 2) == Rate(0.02, Yields.Periodic(2))
        @test Rate(0.02, Inf) == Rate(0.02, Yields.Continuous())
    end

    @testset "rate conversions" begin
        m = Rate(0.1, Yields.Periodic(2))
        @test rate(convert(Yields.Continuous(), m)) ≈ rate(Rate(0.09758, Yields.Continuous())) atol = 1e-5
        c = Rate(0.09758, Yields.Continuous())
        @test convert(Yields.Continuous(), c) == c
        @test rate(convert(Yields.Periodic(2), c)) ≈ rate(Rate(0.1, Yields.Periodic(2))) atol = 1e-5
        @test rate(convert(Yields.Periodic(4), m)) ≈ rate(Rate(0.09878030638383972, Yields.Periodic(4))) atol = 1e-5

    end

    @testset "rate equality" begin
        a = Yields.Periodic(0.02, 2)
        b = Yields.Periodic(0.03, 2)
        c = Yields.Continuous(0.02)

        @test a == a
        @test a != b
        @test ~(a ≈ b)
        @test (a ≈ a)
        @test ~(a ≈ c)

    end

    @testset "constant curve and rate -> Constant" begin
        yield = Yields.Constant(0.05)
        rate = Yields.Rate(0.05, Yields.Periodic(1))

        @test Yields.zero(yield, 1) == Rate(0.05, Yields.Periodic(1))
        @test Yields.zero(Yields.Constant(0.05, Yields.Periodic(2)), 10) == Rate(0.05, Yields.Periodic(2))
        @test Yields.zero(yield, 5, Yields.Periodic(2)) == convert(Yields.Periodic(2), Rate(0.05, Yields.Periodic(1)))

        @testset "constant discount time: $time" for time in [0, 0.5, 1, 10]
            @test discount(yield, time) ≈ 1 / (1.05)^time
            @test discount(rate, time) ≈ 1 / (1.05)^time
            @test discount(rate, 0, time) ≈ 1 / (1.05)^time
        end
        @testset "constant discount scalar time: $time" for time in [0, 0.5, 1, 10]
            @test discount(0.05, time) ≈ 1 / (1.05)^time
        end

        @testset "constant accumulation scalar time: $time" for time in [0, 0.5, 1, 10]
            @test accumulation(0.05, time) ≈ 1 * (1.05)^time
        end

        @testset "broadcasting" begin
            @test all(discount.(yield, [1, 2, 3]) .== 1 ./ 1.05 .^ (1:3))
            @test all(accumulation.(yield, [1, 2, 3]) .== 1.05 .^ (1:3))
        end

        @testset "constant accumulation time: $time" for time in [0, 0.5, 1, 10]
            @test accumulation(yield, time) ≈ 1 * 1.05^time
            @test accumulation(rate, time) ≈ 1 * 1.05^time
            @test accumulation(rate, 0, time) ≈ 1 * 1.05^time
        end

        @testset "CompoundingFrequency" begin
            @testset "Continuous" begin
                cnst = Yields.Constant(Yields.Continuous(0.05))
                @test accumulation(cnst, 1) == exp(0.05)
                @test accumulation(cnst, 2) == exp(0.05 * 2)
                @test discount(cnst, 2) == 1 / exp(0.05 * 2)
            end

            @testset "Periodic" begin
                p = Yields.Constant(Rate(0.05, Yields.Periodic(2)))
                @test accumulation(p, 1) == (1 + 0.05 / 2)^(1 * 2)
                @test accumulation(p, 2) == (1 + 0.05 / 2)^(2 * 2)
                @test discount(p, 2) == 1 / (1 + 0.05 / 2)^(2 * 2)

            end
        end

        yield_2x = yield + yield
        yield_add = yield + 0.05
        add_yield = 0.05 + yield
        @testset "constant discount added" for time in [0, 0.5, 1, 10]
            @test discount(yield_2x, time) ≈ 1 / (1.1)^time
            @test discount(yield_add, time) ≈ 1 / (1.1)^time
            @test discount(add_yield, time) ≈ 1 / (1.1)^time
        end

        yield_1bps = yield - Yields.Constant(0.04)
        yield_minus = yield - 0.01
        minus_yield = 0.05 - Yields.Constant(0.01)
        @testset "constant discount subtraction" for time in [0, 0.5, 1, 10]
            @test discount(yield_1bps, time) ≈ 1 / (1.01)^time
            @test discount(yield_minus, time) ≈ 1 / (1.04)^time
            @test discount(minus_yield, time) ≈ 1 / (1.04)^time
        end
    end

    @testset "short curve" begin
        z = Yields.Zero([0.0, 0.05], [1, 2])
        @test rate(zero(z, 1)) ≈ 0.00
        @test discount(z, 1) ≈ 1.00
        @test rate(zero(z, 2)) ≈ 0.05

        # test no times constructor
        z = Yields.Zero([0.0, 0.05])
        @test rate(zero(z, 1)) ≈ 0.00
        @test discount(z, 1) ≈ 1.00
        @test rate(zero(z, 2)) ≈ 0.05
    end

    @testset "Step curve" begin
        y = Yields.Step([0.02, 0.05], [1, 2])

        @test rate(y, 0.5) == 0.02

        @test discount(y, 0.0) ≈ 1
        @test discount(y, 0.5) ≈ 1 / (1.02)^(0.5)
        @test discount(y, 1) ≈ 1 / (1.02)^(1)
        @test rate(y, 1) ≈ 0.02

        @test discount(y, 2) ≈ 1 / (1.02) / 1.05
        @test rate(y, 2) ≈ 0.05
        @test rate(y, 2.5) ≈ 0.05

        @test discount(y, 2) ≈ 1 / (1.02) / 1.05

        @test discount(y, 1.5) ≈ 1 / (1.02) / 1.05^(0.5)

        @testset "broadcasting" begin
            @test all(discount.(y, [1, 2]) .== [1 / 1.02, 1 / 1.02 / 1.05])
            @test all(accumulation.(y, [1, 2]) .== [1.02, 1.02 * 1.05])
        end

        y = Yields.Step([0.02, 0.07])
        @test rate(y, 0.5) ≈ 0.02
        @test rate(y, 1) ≈ 0.02
        @test rate(y, 1.5) ≈ 0.07

    end

    @testset "Salomon Understanding the Yield Curve Pt 1 Figure 9" begin
        maturity = collect(1:10)

        par = [6.0, 8.0, 9.5, 10.5, 11.0, 11.25, 11.38, 11.44, 11.48, 11.5] ./ 100
        spot = [6.0, 8.08, 9.72, 10.86, 11.44, 11.71, 11.83, 11.88, 11.89, 11.89] ./ 100

        # the forwards for 7+ have been adjusted from the table - perhaps rounding issues are exacerbated
        # in the text? forwards for <= 6 matched so reasonably confident that the algo is correct
        # fwd = [6.,10.2,13.07,14.36,13.77,13.1,12.55,12.2,11.97,11.93] ./ 100 # from text
        fwd = [6.0, 10.2, 13.07, 14.36, 13.77, 13.1, 12.61, 12.14, 12.05, 11.84] ./ 100  # modified

        y = Yields.Par(Rate.(par, Yields.Periodic(1)), maturity)

        @testset "UTYC Figure 9 par -> spot : $mat" for mat in maturity
            @test rate(zero(y, mat)) ≈ spot[mat] atol = 0.0001
            @test forward(y, mat - 1) ≈ Yields.Periodic(fwd[mat], 1) atol = 0.0001
        end

    end

    @testset "simple rate and forward" begin
        # Risk Managment and Financial Institutions, 5th ed. Appendix B

        maturity = [0.5, 1.0, 1.5, 2.0]
        zero = [5.0, 5.8, 6.4, 6.8] ./ 100
        curve = Yields.Zero(zero, maturity)

        @test discount(curve, 0) ≈ 1
        @test discount(curve, 1) ≈ 1 / (1 + zero[2])
        @test discount(curve, 2) ≈ 1 / (1 + zero[4])^2

        @test forward(curve, 0.5, 1.0) ≈ Yields.Periodic(6.6 / 100, 1) atol = 0.001
        @test forward(curve, 1.0, 1.5) ≈ Yields.Periodic(7.6 / 100, 1) atol = 0.001
        @test forward(curve, 1.5, 2.0) ≈ Yields.Periodic(8.0 / 100, 1) atol = 0.001

        y = Yields.Zero(zero)

        @test discount(y, 1) ≈ 1 / 1.05
        @test discount(y, 2) ≈ 1 / 1.058^2

        @testset "broadcasting" begin
            @test all(discount.(y, [1, 2]) .≈ [1 / 1.05, 1 / 1.058^2])
            @test all(accumulation.(y, [1, 2]) .≈ [1.05, 1.058^2])
        end

    end

    @testset "Forward Rates" begin
        # Risk Managment and Financial Institutions, 5th ed. Appendix B

        forwards = [0.05, 0.04, 0.03, 0.08]
        curve = Yields.Forward(forwards, [1, 2, 3, 4])

        @testset "discounts: $t" for (t, r) in enumerate(forwards)
            @test discount(curve, t) ≈ reduce((v, r) -> v / (1 + r), forwards[1:t]; init = 1.0)
        end

        # test constructor without times
        curve = Yields.Forward(forwards)

        @testset "discounts: $t" for (t, r) in enumerate(forwards)
            @test discount(curve, t) ≈ reduce((v, r) -> v / (1 + r), forwards[1:t]; init = 1.0)
        end

        @test accumulation(curve, 0, 1) ≈ 1.05
        @test accumulation(curve, 1, 2) ≈ 1.04
        @test accumulation(curve, 0, 2) ≈ 1.04 * 1.05

        @testset "broadcasting" begin
            @test all(accumulation.(curve, [1, 2]) .≈ [1.05, 1.04 * 1.05])
            @test all(discount.(curve, [1, 2]) .≈ 1 ./ [1.05, 1.04 * 1.05])
        end

        # addition / subtraction
        @test discount(curve + 0.1, 1) ≈ 1 / 1.15
        @test discount(curve - 0.03, 1) ≈ 1 / 1.02

        @testset "with specified timepoints" begin
            i = [0.0, 0.05]
            times = [0.5, 1.5]
            y = Yields.Forward(i, times)
            @test discount(y, 0.5) ≈ 1 / 1.0^0.5
            @test discount(y, 1.5) ≈ 1 / 1.0^0.5 / 1.05^1

        end

    end

    @testset "forwardcurve" begin
        maturity = [0.5, 1.0, 1.5, 2.0]
        zero = [5.0, 5.8, 6.4, 6.8] ./ 100
        curve = Yields.Zero(zero, maturity)

        fwd = Yields.ForwardStarting(curve, 1.0)
        @test discount(fwd, 0) ≈ 1
        @test discount(fwd, 0.5) ≈ discount(curve, 1, 1.5)
        @test discount(fwd, 1) ≈ discount(curve, 1, 2)
        @test accumulation(fwd, 1) ≈ accumulation(curve, 1, 2)
    end

    @testset "base + spread" begin
        riskfree_maturities = [0.5, 1.0, 1.5, 2.0]
        riskfree = [5.0, 5.8, 6.4, 6.8] ./ 100 # spot rates

        spread_maturities = [0.5, 1.0, 1.5, 3.0] # different maturities
        spread = [1.0, 1.8, 1.4, 1.8] ./ 100 # spot spreads

        rf_curve = Yields.Zero(riskfree, riskfree_maturities)
        spread_curve = Yields.Zero(spread, spread_maturities)

        yield = rf_curve + spread_curve

        @test discount(yield, 1.0) ≈ 1 / (1 + riskfree[2] + spread[2])^1
        @test discount(yield, 1.5) ≈ 1 / (1 + riskfree[3] + spread[3])^1.5
    end

    @testset "actual cmt treasury" begin
        # Fabozzi 5-5,5-6
        cmt = [5.25, 5.5, 5.75, 6.0, 6.25, 6.5, 6.75, 6.8, 7.0, 7.1, 7.15, 7.2, 7.3, 7.35, 7.4, 7.5, 7.6, 7.6, 7.7, 7.8] ./ 100
        mats = collect(0.5:0.5:10.0)
        curve = Yields.CMT(cmt, mats)
        targets = [5.25, 5.5, 5.76, 6.02, 6.28, 6.55, 6.82, 6.87, 7.09, 7.2, 7.26, 7.31, 7.43, 7.48, 7.54, 7.67, 7.8, 7.79, 7.93, 8.07] ./ 100
        target_periodicity = fill(2, length(mats))
        target_periodicity[2] = 1 # 1 year is a no-coupon, BEY yield, the rest are semiannual BEY
        @testset "Fabozzi bootstrapped rates" for (r, mat, target, tp) in zip(cmt, mats, targets, target_periodicity)
            @test rate(zero(curve, mat, Yields.Periodic(tp))) ≈ target atol = 0.0001
        end

        # Hull, problem 4.34
        adj = ((1 + 0.051813 / 2)^2 - 1) * 100
        cmt = [4.0816, adj, 5.4986, 5.8620] ./ 100
        mats = [0.5, 1.0, 1.5, 2.0]
        curve = Yields.CMT(cmt, mats)
        targets = [4.0405, 5.1293, 5.4429, 5.8085] ./ 100
        @testset "Hull bootstrapped rates" for (r, mat, target) in zip(cmt, mats, targets)
            @test rate(zero(curve, mat, Yields.Continuous())) ≈ target atol = 0.001
        end

        # test that showing the curve doesn't error
        @test isnothing(show(devnull, curve))

        #     # https://www.federalreserve.gov/pubs/feds/2006/200628/200628abs.html
        #     # 2020-04-02 data
        #     cmt = [0.0945,0.2053,0.4431,0.7139,0.9724,1.2002,1.3925,1.5512,1.6805,1.7853,1.8704,1.9399,1.9972,2.045,2.0855,2.1203,2.1509,2.1783,2.2031,2.2261,2.2477,2.2683,2.2881,2.3074,2.3262,2.3447,2.3629,2.3809,2.3987,2.4164] ./ 100
        #     mats = collect(1:30)
        #     curve = Yields.USCMT(cmt,mats)
        #     target = [0.0945,0.2053,0.444,0.7172,0.9802,1.2142,1.4137,1.5797,1.7161,1.8275,1.9183,1.9928,2.0543,2.1056,2.1492,2.1868,2.2198,2.2495,2.2767,2.302,2.3261,2.3494,2.372,2.3944,2.4167,2.439,2.4614,2.4839,2.5067,2.5297] ./ 100

        #     @testset "FRB data" for (t,mat,target) in zip(1:length(mats),mats,target)
        #         @show mat
        #         if mat >= 1
        #             @test rate(zero(curve,mat, Yields.Continuous())) ≈ target[mat] atol=0.001
        #         end
        #     end
    end

    @testset "OIS" begin
        ois = [1.8, 2.0, 2.2, 2.5, 3.0, 4.0] ./ 100
        mats = [1 / 12, 1 / 4, 1 / 2, 1, 2, 5]
        curve = Yields.OIS(ois, mats)
        targets = [0.017987, 0.019950, 0.021880, 0.024693, 0.029994, 0.040401]
        @testset "bootstrapped rates" for (r, mat, target) in zip(ois, mats, targets)
            @test rate(zero(curve, mat, Yields.Continuous())) ≈ target atol = 0.001
        end
    end

    @testset "InstrumentQuotes" begin

        maturities = [1.3, 2.7]
        prices = [1.1, 0.8]
        zcq = Yields.ZeroCouponQuote.(prices, maturities)
        @test isa(first(zcq), Yields.ZeroCouponQuote)

        @test_throws DimensionMismatch Yields.ZeroCouponQuote.([1.3, 2.4, 0.9], maturities)

        rates = [0.4, -0.7]
        swq = Yields.SwapQuote.(rates, maturities, 3)
        @test first(swq).frequency == 3

        @test_throws DimensionMismatch Yields.SwapQuote.([1.3, 2.4, 0.9], maturities, 3)
        @test_throws DomainError Yields.SwapQuote.(rates, maturities, 0)
        @test_throws DomainError Yields.SwapQuote.(rates, maturities, -2)

        rates = [0.4, -0.7]
        bbq = Yields.BulletBondQuote.(rates, prices, maturities, 3)
        @test first(bbq).frequency == 3
        @test first(bbq).yield == first(rates)

        @test_throws DimensionMismatch Yields.BulletBondQuote.([1.3, 2.4, 0.9], prices, maturities, 3)
        @test_throws DimensionMismatch Yields.BulletBondQuote.(rates, prices, [4.3, 5.6, 4.4, 4.4], 3)
        @test first(Yields.BulletBondQuote.(rates, [5.7], maturities, 3)).price == 5.7
        @test_throws DomainError Yields.BulletBondQuote.(rates, prices, maturities, 0)
        @test_throws DomainError Yields.BulletBondQuote.(rates, prices, maturities, -4)

    end

    @testset "SmithWilson" begin

        ufr = 0.03
        α = 0.1
        u = [5.0, 7.0]
        qb = [2.3, -1.2]

        # Basic behaviour
        sw = Yields.SmithWilson(u, qb; ufr = ufr, α = α)
        @test sw.ufr == ufr
        @test sw.α == α
        @test sw.u == u
        @test sw.qb == qb
        @test_throws DomainError Yields.SmithWilson(u, [2.4, -3.4, 8.9], ufr = ufr, α = α)

        # Empty u and Qb should result in a flat yield curve
        # Use this to test methods expected from <:AbstractYieldCurve
        # Only discount and zero are explicitly implemented, so the others should follow automatically
        sw_flat = Yields.SmithWilson(Float64[], Float64[], ufr = ufr, α = α)
        @test discount(sw_flat, 10.0) == exp(-ufr * 10.0)
        @test accumulation(sw_flat, 10.0) ≈ exp(ufr * 10.0)
        @test rate(convert(Yields.Continuous(), zero(sw_flat, 8.0))) ≈ ufr
        @test discount.(sw_flat, [5.0, 10.0]) ≈ exp.(-ufr .* [5.0, 10.0])
        @test rate(convert(Yields.Continuous(), forward(sw_flat, 5.0, 8.0))) ≈ ufr

        # A trivial Qb vector (=0) should result in a flat yield curve
        ufr_curve = Yields.SmithWilson(u, [0.0, 0.0], ufr = ufr, α = α)
        @test discount(ufr_curve, 10.0) == exp(-ufr * 10.0)

        # A single payment at time 4, zero interest
        curve_with_zero_yield = Yields.SmithWilson([4.0], reshape([1.0], 1, 1), [1.0], ufr = ufr, α = α)
        @test discount(curve_with_zero_yield, 4.0) == 1.0

        # In the long end it's still just UFR
        @test rate(convert(Yields.Continuous(), forward(curve_with_zero_yield, 1000.0, 2000.0))) ≈ ufr

        # Three maturities have known discount factors
        times = [1.0, 2.5, 5.6]
        prices = [0.9, 0.7, 0.5]
        cfs = [1 0 0
            0 1 0
            0 0 1]

        curve_three = Yields.SmithWilson(times, cfs, prices, ufr = ufr, α = α)
        @test transpose(cfs) * discount.(curve_three, times) ≈ prices

        # Two cash flows with payments at three times
        prices = [1.0, 0.9]
        cfs = [0.1 0.1
            1.0 0.1
            0.0 1.0]
        curve_nondiag = Yields.SmithWilson(times, cfs, prices, ufr = ufr, α = α)
        @test transpose(cfs) * discount.(curve_nondiag, times) ≈ prices

        # Round-trip zero coupon quotes
        zcq_times = [1.2, 4.5, 5.6]
        zcq_prices = [1.0, 0.9, 1.2]
        zcq = Yields.ZeroCouponQuote.(zcq_prices, zcq_times)
        sw_zcq = Yields.SmithWilson(zcq, ufr = ufr, α = α)
        @testset "ZeroCouponQuotes round-trip" for idx = 1:length(zcq_times)
            @test discount(sw_zcq, zcq_times[idx]) ≈ zcq_prices[idx]
        end

        # uneven frequencies
        swq_maturities = [1.2, 2.5, 3.6]
        swq_interests = [-0.02, 0.3, 0.04]
        frequency = [2, 1, 2]
        swq = Yields.SwapQuote.(swq_interests, swq_maturities, frequency)
        swq_times = 0.5:0.5:3.5   # Maturities are rounded down to multiples of 1/frequency, [1.0, 2.5, 3.5]
        swq_payments = [
            -0.01 0.3 0.02
            0.99 0 0.02
            0 0.3 0.02
            0 0 0.02
            0 1.3 0.02
            0 0 0.02
            0 0 1.02
        ]

        # Round-trip swap quotes
        swq_maturities = [1.2, 2.5, 3.6]
        swq_interests = [-0.02, 0.3, 0.04]
        frequency = 2
        swq = Yields.SwapQuote.(swq_interests, swq_maturities, frequency)
        swq_times = 0.5:0.5:3.5   # Maturities are rounded down to multiples of 1/frequency, [1.0, 2.5, 3.5]
        swq_payments = [-0.01 0.15 0.02
            0.99 0.15 0.02
            0.0 0.15 0.02
            0.0 0.15 0.02
            0.0 1.15 0.02
            0.0 0.0 0.02
            0.0 0.0 1.02]
        sw_swq = Yields.SmithWilson(swq, ufr = ufr, α = α)
        @testset "SwapQuotes round-trip" for swapIdx = 1:length(swq_interests)
            @test sum(discount.(sw_swq, swq_times) .* swq_payments[:, swapIdx]) ≈ 1.0
        end

        # Round-trip bullet bond quotes (reuse data from swap quotes)
        bbq_prices = [1.3, 0.1, 4.5]
        bbq = Yields.BulletBondQuote.(swq_interests, bbq_prices, swq_maturities, frequency)
        sw_bbq = Yields.SmithWilson(bbq, ufr = ufr, α = α)
        @testset "BulletBondQuotes round-trip" for bondIdx = 1:length(swq_interests)
            @test sum(discount.(sw_bbq, swq_times) .* swq_payments[:, bondIdx]) ≈ bbq_prices[bondIdx]
        end

        @testset "SW ForwardStarting" begin
            fwd_time = 1.0
            fwd = Yields.ForwardStarting(sw_swq, fwd_time)

            @test discount(fwd, 3.7) ≈ discount(sw_swq, fwd_time, fwd_time + 3.7)
        end

        # EIOPA risk free rate (no VA), 31 August 2021.
        # https://www.eiopa.europa.eu/sites/default/files/risk_free_interest_rate/eiopa_rfr_20210831.zip
        eiopa_output_qb = [-0.59556534586390800
            -0.07442224713453920
            -0.34193181987682400
            1.54054875814153000
            -2.15552046042343000
            0.73559290752221900
            1.89365225129089000
            -2.75927773116240000
            2.24893737130629000
            -1.51625404117395000
            0.19284859623817400
            1.13410725406271000
            0.00153268224642171
            0.00147942301778158
            -1.85022125156483000
            0.00336230229850928
            0.00324546553910162
            0.00313268874430658
            0.00302383083427276
            1.36047951448615000]
        eiopa_output_u = 1:20
        eiopa_ufr = log(1.036)
        eiopa_α = 0.133394
        sw_eiopa_expected = Yields.SmithWilson(eiopa_output_u, eiopa_output_qb; ufr = eiopa_ufr, α = eiopa_α)

        eiopa_eurswap_maturities = [1:12; 15; 20]
        eiopa_eurswap_rates = [-0.00615, -0.00575, -0.00535, -0.00485, -0.00425, -0.00375, -0.003145,
            -0.00245, -0.00185, -0.00125, -0.000711, -0.00019, 0.00111, 0.00215]   # Reverse engineered from output curve. This is the full precision of market quotes.
        eiopa_eurswap_quotes = Yields.SwapQuote.(eiopa_eurswap_rates, eiopa_eurswap_maturities, 1)
        sw_eiopa_actual = Yields.SmithWilson(eiopa_eurswap_quotes, ufr = eiopa_ufr, α = eiopa_α)

        @testset "Match EIOPA calculation" begin
            @test sw_eiopa_expected.u ≈ sw_eiopa_actual.u
            @test sw_eiopa_expected.qb ≈ sw_eiopa_actual.qb
        end
    end

end


text:
module Operators

export Ax, Ay, Az, volume
export Δzᵃᵃᶜ, Δzᵃᵃᶠ
export Δxᶜᶜᵃ, Δxᶠᶜᵃ, Δxᶠᶠᵃ, Δxᶜᶠᵃ
export Δyᶜᶜᵃ, Δyᶠᶜᵃ, Δyᶠᶠᵃ, Δyᶜᶠᵃ
export Axᵃᵃᶜ, Axᵃᵃᶠ, Axᶜᶜᶜ, Axᶠᶜᶜ, Axᶠᶠᶜ, Axᶠᶜᶠ, Axᶜᶠᶜ, Axᶜᶜᶠ
export Ayᵃᵃᶜ, Ayᵃᵃᶠ, Ayᶜᶜᶜ, Ayᶜᶠᶜ, Ayᶠᶠᶜ, Ayᶜᶠᶠ, Ayᶠᶜᶜ, Ayᶜᶜᶠ
export Azᵃᵃᵃ, Azᶠᶠᵃ, Azᶜᶜᵃ, Azᶠᶜᵃ, Azᶜᶠᵃ
export Vᵃᵃᶜ, Vᵃᵃᶠ, Vᶜᶜᶜ, Vᶠᶜᶜ, Vᶜᶠᶜ, Vᶜᶜᶠ
export δxᶜᵃᵃ, δxᶠᵃᵃ, δyᵃᶜᵃ, δyᵃᶠᵃ, δzᵃᵃᶜ, δzᵃᵃᶠ
export Ax_ψᵃᵃᶠ, Ax_ψᵃᵃᶜ, Ay_ψᵃᵃᶠ, Ay_ψᵃᵃᶜ, Az_ψᵃᵃᵃ
export Ax_uᶠᶜᶜ, Ay_vᶜᶠᶜ, Az_wᶜᶜᵃ
export Ax_ζᶠᶠᶜ, Ay_ζᶠᶠᶜ, Ax_ηᶠᶜᶠ, Az_ηᶠᶜᵃ, Ay_ξᶜᶠᶠ, Az_ξᶜᶠᵃ
export Az_wᶜᶜᵃ, Az_ηᶠᶜᵃ, Az_ξᶜᶠᵃ
export Ax_cᶜᶜᶜ, Ay_cᶜᶜᶜ, Az_cᶜᶜᵃ
export δᴶxᶜᵃᶜ, δᴶxᶜᵃᶠ, δᴶxᶠᵃᶜ, δᴶxᶠᵃᶠ, δᴶyᵃᶜᶜ, δᴶyᵃᶜᶠ, δᴶyᵃᶠᶜ, δᴶyᵃᶠᶠ, δᴶzᵃᵃᶜ, δᴶzᵃᵃᶠ
export ∂xᶜᵃᵃ, ∂xᶠᵃᵃ, ∂yᵃᶜᵃ, ∂yᵃᶠᵃ, ∂zᵃᵃᶜ, ∂zᵃᵃᶠ
export ∂xᶜᶜᵃ, ∂xᶠᶠᵃ, ∂xᶠᶜᵃ, ∂xᶜᶠᵃ, ∂yᶜᶜᵃ, ∂yᶠᶠᵃ, ∂yᶠᶜᵃ, ∂yᶜᶠᵃ
export ∂²xᶜᵃᵃ, ∂²xᶠᵃᵃ, ∂²yᵃᶜᵃ, ∂²yᵃᶠᵃ, ∂²zᵃᵃᶜ, ∂²zᵃᵃᶠ
export ∂³zᵃᵃᶜ, ∂³zᵃᵃᶠ
export ∂⁴xᶜᵃᵃ, ∂⁴xᶠᵃᵃ, ∂⁴yᵃᶜᵃ, ∂⁴yᵃᶠᵃ, ∂⁴zᵃᵃᶜ, ∂⁴zᵃᵃᶠ
export ℑxᶜᵃᵃ, ℑxᶠᵃᵃ, ℑyᵃᶜᵃ, ℑyᵃᶠᵃ, ℑzᵃᵃᶜ, ℑzᵃᵃᶠ
export ℑᴶxᶜᵃᵃ, ℑᴶxᶠᵃᵃ, ℑᴶyᵃᶜᵃ, ℑᴶyᵃᶠᵃ, ℑᴶzᵃᵃᶜ, ℑᴶzᵃᵃᶠ
export ℑxyᶜᶜᵃ, ℑxyᶠᶜᵃ, ℑxyᶠᶠᵃ, ℑxyᶜᶠᵃ, ℑxzᶜᵃᶜ, ℑxzᶠᵃᶜ, ℑxzᶠᵃᶠ, ℑxzᶜᵃᶠ, ℑyzᵃᶜᶜ, ℑyzᵃᶠᶜ, ℑyzᵃᶠᶠ, ℑyzᵃᶜᶠ
export ℑxyzᶜᶜᶠ, ℑxyzᶜᶠᶜ, ℑxyzᶠᶜᶜ, ℑxyzᶜᶠᶠ, ℑxyzᶠᶜᶠ, ℑxyzᶠᶠᶜ
export divᶜᶜᶜ, div_xyᶜᶜᵃ, div_xzᶜᵃᶜ, div_yzᵃᶜᶜ, ζ₃ᶠᶠᵃ
export ∇²ᶜᶜᶜ, ∇²hᶜᶜᶜ, ∇²hᶠᶜᶜ, ∇²hᶜᶠᶜ

using Oceananigans.Grids

#####
##### Convenient aliases
#####

const AG   = AbstractGrid
const ARG  = AbstractRectilinearGrid
const RCG  = RectilinearGrid
const ACG  = AbstractCurvilinearGrid
const AHCG = AbstractHorizontallyCurvilinearGrid

include("difference_operators.jl")
include("interpolation_operators.jl")
include("interpolation_utils.jl")

include("spacings_and_areas_and_volumes.jl")
include("products_between_fields_and_grid_metrics.jl")

include("derivative_operators.jl")
include("divergence_operators.jl")
include("vorticity_operators.jl")
include("laplacian_operators.jl")

end

ext: jl


text:
# Multiples of 3 and 5

LIMIT = 1000

function solve()
    return sum((i % 3 == 0 || i % 5 == 0) ? i : 0 for i in 1:(LIMIT - 1))
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
using JLD2
using DataDeps
using Requires
using Transformers
using Pkg.TOML

export list_pretrains, load_pretrained, @pretrained_str

const pretrained_configs = open(TOML.parse, joinpath(@__DIR__, "pretrains.toml"))

"""
    readckpt(path)

Load weights from a tensorflow checkpoint file into a Dict.
"""
readckpt(path) = error("readckpt require TensorFlow.jl installed. run `Pkg.add(\"TensorFlow\"); using TensorFlow`")

@init @require TensorFlow="1d978283-2c37-5f34-9a8e-e9c0ece82495" begin
    import .TensorFlow
    function readckpt(path)
        weights = Dict{String, Array}()
        TensorFlow.init()
        ckpt = TensorFlow.pywrap_tensorflow.x.NewCheckpointReader(path)
        shapes = ckpt.get_variable_to_shape_map()

        for (name, shape) ∈ shapes
            # Ignore training related data - e.g. learning rates
            # Also ignore scalars and other variables that aren't stored correctly in the checkpoint (shape == Any[])
            (occursin("training", name) || shape == Any[]) && continue
            weight = ckpt.get_tensor(name)
            if length(shape) == 2
                weight = collect(weight')
            end
            weights[name] = weight
        end

        weights
    end
end

"""
    ckpt2_to_jld2(ckptpath::String, ckptname::String, savepath::String)

Loads the pre-trained model weights from a TensorFlow checkpoint and saves to JLD2
"""
function ckpt_to_jld2(ckptpath::String, ckptname::String; savepath::String="./")
    weights = readckpt(joinpath(ckptpath, ckptname))
    jld2name = normpath(joinpath(savepath, ckptname[1:end-5]*".jld2"))
    @info "Saving the model weights to $jld2name"
    JLD2.@save jld2name weights
end

function register_configs(configs)
    for (model_name, config) in pairs(configs)
        model_desc = Transformers.Pretrain.description(config["description"], config["host"], config["link"])
        checksum = config["checksum"]
        url = config["url"]
        dep = DataDep(model_name, model_desc, url, checksum;
                      fetch_method=Transformers.Datasets.download_gdrive)
        DataDeps.register(dep)
    end
end

"""
    list_pretrains()

List all the available pre-trained models.
"""
function list_pretrains()
    println.(keys(pretrained_configs))
    return
end

loading_method(model_name) = throw(error("Loading method for this model is not defined."))

"""
    load_pretrained(path)

Loads a pre-trained Music Transformer model.
"""
function load_pretrained(model_name::String)
    if model_name ∉ keys(pretrained_configs)
        error("""Invalid model.
               Please try list_pretrains() to check the available pre-trained models""")
    end

    model_config = pretrained_configs[model_name]
    loader = loading_method(Val(Symbol(model_config["model_type"])))

    model_path = @datadep_str("$model_name/$model_name.jld2")
    if !endswith(model_path, ".jld2")
        error("""Invalid file. A .jld2 file is required to load the model.
                 If this is a tensorflow checkpoint file, run ckpt_to_jld2 to convert it.""")
    end

    JLD2.@load model_path weights
    loader(weights, model_config)
end

macro pretrained_str(name)
    :(load_pretrained($(esc(name))))
end


text:
export error_if

import Base: show

"""
    error_if(checkFn, errorFn)

Creates an `error_if` operator, which performs a check for every emission on the source Observable with `checkFn`.
If `checkFn` returns `true`, the operator sends an `error` event and unsubscribes from the observable.

# Arguments
- `checkFn`: check function with `(data) -> Bool` signature
- `errorFn`: error object generating function with `(data) -> Any` signature, optional

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([1, 2, 3]) |> error_if((data) -> data > 2, (data) -> "CustomError")

subscription = subscribe!(source, lambda(
    on_next  = (d) -> println("Next: ", d),
    on_error = (e) -> println("Error: ", e),
    on_complete = () -> println("Completed")
))
;

# output
Next: 1
Next: 2
Error: CustomError
```

See also: [`error_if_not`](@ref), [`error_if_empty`](@ref), [`default_if_empty`](@ref), [`lambda`](@ref)
"""
error_if(checkFn::F, errorFn::E = nothing) where { F, E } = ErrorIfOperator{F, E}(checkFn, errorFn)

struct ErrorIfOperator{F, E} <: InferableOperator
    checkFn :: F
    errorFn :: E
end

function on_call!(::Type{L}, ::Type{L}, operator::ErrorIfOperator{F, E}, source) where { L, F, E }
    return proxy(L, source, ErrorIfProxy{F, E}(operator.checkFn, operator.errorFn))
end

operator_right(::ErrorIfOperator, ::Type{L}) where L = L

struct ErrorIfProxy{F, E} <: ActorSourceProxy
    checkFn :: F
    errorFn :: E
end

actor_proxy!(::Type{L}, proxy::ErrorIfProxy{F, E}, actor::A) where { L, A, F, E } = ErrorIfActor{L, A, F, E}(proxy.checkFn, proxy.errorFn, actor, false, voidTeardown)
source_proxy!(::Type{L}, proxy::ErrorIfProxy, source::S) where { L, S } = ErrorIfSource{L, S}(source)

mutable struct ErrorIfActor{L, A, F, E} <: Actor{L}
    checkFn :: F
    errorFn :: E
    actor :: A
    completed::Bool
    subscription
end

error_msg(actor::ErrorIfActor, data) = error_msg(actor, actor.errorFn, data)

error_msg(::ErrorIfActor, ::Nothing, data) = "`error_if` operator check failed for data $(data)"
error_msg(::ErrorIfActor, callback, data)  = callback(data)

function on_next!(actor::ErrorIfActor{L}, data::L) where L
    if !actor.completed
        check = actor.checkFn(data)
        if check
            error!(actor, error_msg(actor, data))
        else
            next!(actor.actor, data)
        end
    end
end

function on_error!(actor::ErrorIfActor, err)
    if !actor.completed
        actor.completed = true
        unsubscribe!(actor.subscription)
        error!(actor.actor, err)
    end
end

function on_complete!(actor::ErrorIfActor)
    if !actor.completed
        actor.completed = true
        complete!(actor.actor)
    end
end

@subscribable struct ErrorIfSource{L, S} <: Subscribable{L}
    source :: S
end

function on_subscribe!(source::ErrorIfSource, actor::ErrorIfActor)
    subscription = subscribe!(source.source, actor)
    if !actor.completed
        actor.subscription = subscription
    end
    return subscription
end

Base.show(io::IO, ::ErrorIfOperator)          = print(io, "ErrorIfOperator()")
Base.show(io::IO, ::ErrorIfProxy)             = print(io, "ErrorIfProxy()")
Base.show(io::IO, ::ErrorIfActor{L}) where L  = print(io, "ErrorIfActor($L)")
Base.show(io::IO, ::ErrorIfSource{L}) where L = print(io, "ErrorIfSource($L)")


text:
# Developed by Jason Ma
# This is a Julia implementation of the strangle model given to us
# This includes the updated version of the finutil_stu file and the
# strangle_stu.py file. My Julia environment had issues importing the
# functions from finutil_stu to strangle_stu, so I merged the files here
# The split between finutil and strangle_stu occurs around line 265.
# Note, I have also included what the files should look like when
# importing works correctly. The only change from this merged file
# to what it should look like is by adding a "fu." prefix to any functions
# from finutil in the strangle_stu.jl file

# This is finutil_stu.jl
# This is a frozen and reduced version of finutil.py. This file is not likely
# to be current with funutil and should not be used in place of finutil (see log).
# This should only be used with models with the _stu appendage.
# These are the utilities that are used throughout the various financial models
# for such things as calculating options prices and volatilities. These are
# designed to be used with Version 3.6 of Python.
# This version 1.2stu, April 11, 2018  NOTE: This version has not been completely
# tested.
# Maintained by Professor Evans
#
#
# daysto calculates the number of days between now and some event, such as days2expiry
# as an integer float.
#
using SpecialFunctions

function daysto(eyear,emonth,eday)
    tnow = Dates.today()
    expiry = Dates.Date(eyear, emonth, eday)
    days2expiry = abs(expiry - tnow)
    return float(days2expiry.days)
end
#
function monthname(mo)
    monthlabel = ["JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG",
                 "SEP", "OCT", "NOV", "DEC" ]
    return monthlabel[mo]
end
#
# csnd integrates a standard normal distribution up to some sigma.
#
function csnd(point)
    return (1.0 + erf(point/sqrt(2.0)))/2.0
end
#
# cnd integrates a Gaussian distribution up to some value.
#
function cnd(center,point,stdev)
    return (1.0 + erf((point - center)/(stdev*sqrt(2.0))))/2.0
end
#
# Calculating the BSM call option price, traditional model. This requires the
# user to provide stock price, strike price, daily volatility, risk-free interest
# rate and days to expiry. (To calculate days use method daysto above).
# This returns the call price, the delta, and duration volatility as a tuple
# array. See popo below for puts.
#
function copo(stock,strike,dayvol,days,rfir)
    d1 = log(stock/strike)+((rfir/365)+(dayvol^2)/2)*days
    durvol = dayvol*sqrt(days)
    delta = csnd(d1/durvol)
    cumd2 = csnd((d1/durvol) - durvol)
    discount = exp(-rfir*days/365)
    callpr = (stock*delta)-(strike*discount*cumd2)
    return [callpr,delta,durvol]
end
#
# An elementary function for calculating the stock price adjusted for drift.
#
function drift(alpha,time)
    return 1.0*exp(alpha*time)
end
#
# An elementary multiplier function for converting daily volatility to duration
# volatility.
#
function durvol(time)
    return 1.0*sqrt(time)
end
#
# An elementary price expected-mean-value adjustment multiplier for log
# distributed prices. The mean of a log-distributed pdf is adjusted by minus
# one-half variance.
#
function lnmeanshift(sigma)
    return 1.0*exp(-1.0*(sigma*sigma/2))
end
#
# A time-discount function to discount the value of a future payment (like an
# option) discounted at the risk-free interest rate. The variable riskfreerate
# is annual and time is in days.
#
function dcount(riskfreerate,time)
    return 1.0*exp(-1.0*(riskfreerate/365)*time)
end
#
# This is an option tranche value calculator function that assumes you have a
# stock price and strike price, adjusted externally (for example, drift is
# adjusted with drift above). This will calculate the strike-price adjusted
# tranche from either -5 sigma to the strike or from the strike to +5 sigma.
# Sigma used here is duration sigma, adjusted outside using durvol above.
# Main program must set call to true if a call, false if a put.
#
function otranche(stock,strike,dursigma,call)
    sspread = (log(strike/stock))/dursigma
    if call
        binborder = linspace(sspread, 5.00, 24)
    else
        binborder = linspace(-5.0, sspread, 24)
    end
    num_bins = 24
    binedgeprob = zeros(num_bins)
#
    for i = 1:num_bins
        binedgeprob[i] = csnd(binborder[i])
    end
    num_bins = num_bins - 1
    binprob = zeros(num_bins)
    binmidprice = zeros(num_bins)
    binvalue = zeros(num_bins)
#
    for i = 1:num_bins
        binprob[i] = binedgeprob[i+1] - binedgeprob[i]
        binmidprice[i] = ((stock*exp(((binborder[i+1]+binborder[i])/2.0)*dursigma))*lnmeanshift(dursigma)) - strike
        binvalue[i] = binmidprice[i]*binprob[i]
    end
    if call
        optionprice = sum(binvalue)
    else
        optionprice = (sum(binvalue))*-1.0
    end
#
    return optionprice
end
#
# This is full tranche value calculator function that assumes you have a stock
# price and reference price, usually a strike price (and still called that
# here). This will calculate the tranche value from either -5 sigma to the
# reference price (left is true) or from the reference to +5 sigma (left
# is false). This is similar to otranche except that it does not subtract the
# strike price and therefore gives the full value of the tranche. It is
# designed to be used primarily by the Aruba Model to calculate the value of
# the remaining stock tranche if the covered call you wrote has been exercised.
#
function ftranche(stock,strike,sigma,left)
    sspread = (log(strike/stock))/sigma
    if left
        binborder = linspace(-5.0, sspread, n=24)
    else
        binborder = linspace(sspread, 5.0, n=24)
    end
    num_bins = 24
    binedgeprob = zeros(24)
#
    for i = 1:num_bins
        binedgeprob[i] = csnd(binborder[i])
    end
    num_bins = num_bins - 1
    binprob = zeros(num_bins)
    binmidprice = np.zeros(num_bins)
    binvalue = np.zeros(num_bins)
#
    for i = 1:num_bins
        binprob[i] = binedgeprob[i+1] - binedgeprob[i]
        binmidprice[i] = ((stock*exp(((binborder[i+1]+binborder[i])/2.0)
        *sigma))*lnmeanshift(sigma))
        binvalue[i] = binmidprice[i]*binprob[i]
    end
#
    trancheprice = sum(binvalue)
    return trancheprice
end
#
# oidv calculates implied daily and duration volatility for a call or a put
# using divide and conquer (the default for most models). Also see oidvnm.
# This uses an iterative process that uses otranche (above) to calculate the
# sigma, here an implied sigma, from the existing option value (ovalue).
# The call variable is True for a call, False for a put. The convergence is
# within the while loop. This function returns a tuple of two values, daily
# IDV and duration IDV.
#
function oidv(stock,strike,ovalue,days,call)
    epsilon = 1e-4
    low = 0.0
    high = 1.0
    daysigma = (high+low)/2
    dursigma = daysigma*durvol(days)
    tempop = otranche(stock,strike,dursigma,call)
    while tempop<=(ovalue-epsilon) || tempop>=(ovalue+epsilon)
        if tempop >= (ovalue+epsilon)
            high = daysigma
        else
            low = daysigma
        end
        daysigma = (high+low)/2
        dursigma = daysigma*durvol(days)
        tempop = otranche(stock,strike,dursigma,call)
    end
    # End of Loop!
    return [daysigma,dursigma]
end
#
# oidvnm calculates implied daily and duration volatility for a call or a put
# using Newton's Method for convergence (default is divide and conquer).
# This uses an iterative process that uses otranche (above) to calculate the
# sigma, here an implied sigma, from the existing option value (ovalue).
# The call variable is True for a call, False for a put. The convergence is
# within the while loop. This function returns a tuple of two values, daily
# IDV and duration IDV.
#
function oidvnm(stock,strike,ovalue,days,call)
    seedsigma = 1e-6
    durseed = seedsigma*durvol(days)
    # NOTE: daysigma is supposed to be set at a reasonable estimate of the
    # actual idv. The Newton method can explode if it is not (especially if you
    # enter option price, strike, and value data that are very unrealistic).
    # Consider setting daysigma at sqrt*time*ln(strike/stock). It was
    # originally set at 0.05
    daysigma = (log(strike/stock))*sqrt(days)
    dursigma = daysigma*durvol(days)
    cutoff = 1e-4
    tempop = 0.00
    #
    # The loop starts here. You start with a test sigma and converge to the
    # actual sigma. The convergence shown here (Newton's method) was designed
    # by Alec Griffith '17
    #
    while abs(tempop - ovalue) > cutoff
        tempop = otranche(stock,strike,dursigma,call)
        price2 = otranche(stock,strike,dursigma+durseed,call)
        deriv = (price2-tempop)/seedsigma
        daysigma -= (tempop-ovalue)/deriv
        dursigma = daysigma*durvol(days)
    end
    # End of Loop!
    return [daysigma,dursigma]
end
#
# Calculating the BSM put option price, traditional model. This requires the
# user to provide stock price, strike price, daily volatility, risk-free interest
# rate and days to expiry. (To calculate days use method daysto above).
# This returns the put price, the delta, and duration volatility as a tuple
# array. See copo below for calls.
#
function popo(stock,strike,dayvol,days,rfir)
    d1 = log(stock/strike)+((rfir/365)+(dayvol^2)/2)*days
    durvol = dayvol*sqrt(days)
    delta = csnd(-d1/durvol)
    cumd2 = csnd(-(d1/durvol - durvol))
    discount = exp(-rfir*days/365)
    putpr = -(stock*delta)+(strike*discount*cumd2)
    return [putpr,delta,durvol]
end
#
# tdfecay adds an extension to the otranche calculator (Taboga model)
# to calculate one-day time decay.
#

function tdecay(stock, strike, daysigma, oprice, days, call)
    days = days - 1.0
    dursigma = daysigma*durvol(days)
    oprice1d = otranche(stock,strike,dursigma,call)
    timedecay = oprice - oprice1d
    return timedecay
end

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

# strangle_stu.py
#
# This is a MOS model: Mudd Open Source only.
# This is a deep strangle model based upon the Taboga concept (student version).
# Much of the documentation can be found in the core Taboga model or finutil.py
# Designed in Python 3.6 on for Economics 136
# Original version developed on May 30, 2017 by Prof Evans
# This is version 2.1stu, dated April 11, 2017
#

# This model draws upon utility model finutil_stu. That model must either be in
# the same directory as this model or the path command below must be used to
# point to its location.
# import sys
# sys.path.append('c:/Users/Prof Gary Evans/Dropbox/PyGo/PyFi')

# This version calculates IDV for the put and call, but also requires that you
# supply historical volatility (252 and 30 day). You can include an earnings
# date if date is relevant.
#
# Load data
aname = "Prof Evans"
stosym = "MSFT"
tnow = Dates.today()
# Set earnings to True if this is an earnings-based strangle.
earnings = true
# If you don't need to calculate the days to expiry and earnings (optional)
# then calcdays is False (default).
calcdays = false
# BEGIN of optional time data. Load below ONLY if calcdays = True
# Expiration date:
exyear = 2018
exmonth = 2
exday = 16
# Earnings date (load only if relevant):
eayear = 2018
eamonth = 1
eaday = 31
# END of optional time data load
#
# This model does not use the alpha.
# stock = 15.15
# callstrike = 16
# callbid = 0.90
# callask = 0.93
# putstrike = 15
# putbid = 1.150
# putask = 1.18
# sigma252 = 0.043534
# sigma30 = 0.025
###########################
stock = 166.86
callstrike = 170
callbid = 3.20
callask = 3.25
putstrike = 165
putbid = 4.50
putask = 4.60
sigma252 = 0.011358
sigma30 = 0.016387
if calcdays
    days = daysto(exyear,exmonth,exday)
    if earnings
        daystoearn = daysto(eayear,eamonth,eaday)
    end
else
    # Here you have to supply the days if you don't calculate them.
    days = 14.0
    if earnings
        daystoearn = 11.0
    end
end
#
moname = monthname(exmonth)
# Set effective options prices (PEG)
peg = 0.6
cspread = callask - callbid
callpr = callbid + (peg*cspread)
pspread = putask - putbid
putpr = putbid + (peg*pspread)
cost = callpr + putpr
# Note: callbe and putbe should be calced based upon the prices they would
# have to be the day after earnings, especially for more distant options.
callbe = stock + cost
putbe = stock - cost
callbegr = log(callbe/stock)*100
putbegr = log(putbe/stock)*100
# Consider adding log percent breakeven to current value
#
# Calculate call and put ideal prices based upon historical volatilities.
hv252 = sigma252*durvol(days)
cop252 = otranche(stock,callstrike,hv252,true)
pop252 = otranche(stock,putstrike,hv252,false)
hv30 = sigma30*durvol(days)
cop30 = otranche(stock,callstrike,hv30,true)
pop30 = otranche(stock,putstrike,hv30,false)
#
# Calculate call and put IDV based upon actual prices
# See the documentation for finutil and tabogacidv to see what these are doing.
# Start with the call:
#
seedsigma = 1e-6
cutoff = 1e-4
tempop = 0.00
price2 = 0.00
derive = 0.00
cdurseed = seedsigma*durvol(days)
cdaysigma = 0.06
cdursigma = cdaysigma*durvol(days)
#
# The convergence shown here (Newton's method) was designed by Alec Griffith '17
#
while abs(tempop - callpr) > cutoff
    tempop = otranche(stock,callstrike,cdursigma,true)
    price2 = otranche(stock,callstrike,cdursigma+cdurseed,true)
    deriv = (price2-tempop)/seedsigma
    cdaysigma -= (tempop-callpr)/deriv
    cdursigma = cdaysigma*durvol(days)
end
# End of Loop!
#
# Repeat for the put:
#
seedsigma = 1e-6
cutoff = 1e-4
tempop = 0.00
price2 = 0.00
derive = 0.00
pdurseed = seedsigma*durvol(days)
pdaysigma = 0.01
pdursigma = pdaysigma*durvol(days)
call = false
#
#
while abs(tempop - putpr) > cutoff
    tempop = otranche(stock,putstrike,pdursigma,false)
    price2 = otranche(stock,putstrike,pdursigma+pdurseed,false)
    deriv = (price2-tempop)/seedsigma
    pdaysigma -= (tempop-putpr)/deriv
    pdursigma = pdaysigma*durvol(days)
end
# End of Loop!
callsigmaratio = cdaysigma/sigma252
putsigmaratio = pdaysigma/sigma252
#
dateString = Dates.format(now(), "e u d y")
println()
println("Strangle Model")
println("Analyst's name: $(aname)")
println( "Date: $(dateString)")
println( "Stock symbol: $(stosym)")
println("Option chain: $(moname), $(exday)")
println("Days to expiry: $(days)")
println( "Stock price: $(stock)")
println( "Call strike price: $(callstrike)")
println( "Call Bid: $(callbid)")
println( "Call Ask: $(callask)")
println( "Call price: $(callpr)")
println( "Put strike price: $(putstrike)")
println( "Put Bid: $(putbid)")
println( "Put Ask: $(putask)")
println( "Put price: $(putpr)")
println( "Position cost: $(cost)")
if earnings
    println("Days to earnings: $(daystoearn)")
println( "Call breakeven price: $(callbe)")
println( "Call breakeven ln percent: $(callbegr)")
println( "Put breakeven price: $(putbe)")
println( "Put breakeven ln percent: $(putbegr)")
println( "Call value based upon 252-day HV ($(sigma252)) : $(cop252)")
println( "Call value based upon 30-day HV ($(sigma30)) : $(cop30)")
println( "Put value based upon 252-day HV ($(sigma252)) : $(pop252)")
println( "Put value based upon 30-day HV ($(sigma30)) : $(pop30)")
println( "Call implied duration volatility: $(cdursigma)")
println( "Put implied duration volatility: $(pdursigma)")
println("CRITICAL")
println( "Call implied daily volatility: $(cdaysigma)")
println( "Put implied daily volatility: $(pdaysigma)")
println( "Call sigma ratio: $(callsigmaratio)")
println( "Put sigma ratio: $(putsigmaratio)")
#
# Calculate one day time decay using our calculated IDV.
#
if days >= 2.0
    days = days - 1.0
    cdursigma = cdaysigma*durvol(days)
    cprice1d = otranche(stock,callstrike,cdursigma,true)
    pdursigma = pdaysigma*durvol(days)
    pprice1d = otranche(stock,putstrike,pdursigma,false)
    timedecay = (callpr - cprice1d) + (putpr - pprice1d)
    println( "One day time decay: $(timedecay)")
else
    println( "No time decay with one day remaining.")
end
end
#


text:
# Dictionary mapping codes to more verbose information #
# See Appendix A of the NI-VISA Programmers Reference
# Dictionary maps integer codes to tuples of (Code, Meaning)
# from: https://github.com/BBN-Q/Instruments.jl
codes = Dict{Int,Tuple{String,String}}(
    VI_SUCCESS => ("VI_SUCCESS", "Operation completed successfully."),
    VI_SUCCESS_EVENT_EN => ("VI_SUCCESS_EVENT_EN", "Specified event is already enabled for at least one of the specified mechanisms."),
    VI_SUCCESS_EVENT_DIS => ("VI_SUCCESS_EVENT_DIS", "Specified event is already disabled for at least one of the specified mechanisms."),
    VI_SUCCESS_QUEUE_EMPTY => ("VI_SUCCESS_QUEUE_EMPTY", "Operation completed successfully, but queue was already empty."),
    VI_SUCCESS_TERM_CHAR => ("VI_SUCCESS_TERM_CHAR", "The specified termination character was read."),
    VI_SUCCESS_MAX_CNT => ("VI_SUCCESS_MAX_CNT", "The number of bytes transferred is equal to the requested input count. More data may be available."),
    VI_SUCCESS_DEV_NPRESENT => ("VI_SUCCESS_DEV_NPRESENT", "Session opened successfully, but the device at the specified address is not responding."),
    VI_SUCCESS_TRIG_MAPPED => ("VI_SUCCESS_TRIG_MAPPED", "The path from trigSrc to trigDest is already mapped."),
    VI_SUCCESS_QUEUE_NEMPTY => ("VI_SUCCESS_QUEUE_NEMPTY", "Wait terminated successfully on receipt of an event notification. There is at least one more event object of the requested type(s) available for this session."),
    VI_SUCCESS_NCHAIN => ("VI_SUCCESS_NCHAIN", "Event handled successfully. Do not invoke any other handlers on this session for this event."),
    VI_SUCCESS_NESTED_SHARED => ("VI_SUCCESS_NESTED_SHARED", "Operation completed successfully, and this session has nested shared locks."),
    VI_SUCCESS_NESTED_EXCLUSIVE=> ("VI_SUCCESS_NESTED_EXCLUSIVE", "Operation completed successfully, and this session has nested exclusive locks."),
    VI_SUCCESS_SYNC => ("VI_SUCCESS_SYNC", "Operation completed successfully, but the operation was actually synchronous rather than asynchronous."),
    VI_WARN_QUEUE_OVERFLOW => ("VI_WARN_QUEUE_OVERFLOW", "VISA received more event information of the specified type than the configured queue size could hold."),
    VI_WARN_CONFIG_NLOADED => ("VI_WARN_CONFIG_NLOADED", "The specified configuration either does not exist or could not be loaded. VISA-specified defaults will be used."),
    VI_WARN_NULL_OBJECT => ("VI_WARN_NULL_OBJECT", "The specified object reference is uninitialized."),
    VI_WARN_NSUP_ATTR_STATE => ("VI_WARN_NSUP_ATTR_STATE", "Although the specified state of the attribute is valid, it is not supported by this implementation."),
    VI_WARN_UNKNOWN_STATUS => ("VI_WARN_UNKNOWN_STATUS", "The status code passed to the operation could not be interpreted."),
    VI_WARN_NSUP_BUF => ("VI_WARN_NSUP_BUF", "The specified I/O buffer type is not supported."),
    VI_WARN_EXT_FUNC_NIMPL => ("VI_WARN_EXT_FUNC_NIMPL", "The operation succeeded, but a lower level driver did not implement the extended functionality."),
    VI_ERROR_SYSTEM_ERROR => ("VI_ERROR_SYSTEM_ERROR", "Unknown system error (miscellaneous error)."),
    VI_ERROR_INV_OBJECT => ("VI_ERROR_INV_OBJECT", "The given session or object reference is invalid."),
    VI_ERROR_RSRC_LOCKED => ("VI_ERROR_RSRC_LOCKED", "Specified type of lock cannot be obtained, or specified operation cannot be performed, because the resource is locked."),
    VI_ERROR_INV_EXPR => ("VI_ERROR_INV_EXPR", "Invalid expression specified for search."),
    VI_ERROR_RSRC_NFOUND => ("VI_ERROR_RSRC_NFOUND", "Insufficient location information or the requested device or resource is not present in the system."),
    VI_ERROR_INV_RSRC_NAME => ("VI_ERROR_INV_RSRC_NAME", "Invalid resource reference specified. Parsing error."),
    VI_ERROR_INV_ACC_MODE => ("VI_ERROR_INV_ACC_MODE", "Invalid access mode."),
    VI_ERROR_TMO => ("VI_ERROR_TMO", "Timeout expired before operation completed."),
    VI_ERROR_CLOSING_FAILED => ("VI_ERROR_CLOSING_FAILED", "The VISA driver failed to properly close the session or object reference. This might be due to an error freeing internal or OS resources, a failed network connection, or a lower-level driver or OS error."),
    VI_ERROR_INV_DEGREE => ("VI_ERROR_INV_DEGREE", "Specified degree is invalid."),
    VI_ERROR_INV_JOB_ID => ("VI_ERROR_INV_JOB_ID", "Specified job identifier is invalid."),
    VI_ERROR_NSUP_ATTR => ("VI_ERROR_NSUP_ATTR", "The specified attribute is not defined or supported by the referenced object."),
    VI_ERROR_NSUP_ATTR_STATE => ("VI_ERROR_NSUP_ATTR_STATE", "The specified state of the attribute is not valid, or is not supported as defined by the object."),
    VI_ERROR_ATTR_READONLY => ("VI_ERROR_ATTR_READONLY", "The specified attribute is read-only."),
    VI_ERROR_INV_LOCK_TYPE => ("VI_ERROR_INV_LOCK_TYPE", "The specified type of lock is not supported by this resource."),
    VI_ERROR_INV_ACCESS_KEY => ("VI_ERROR_INV_ACCESS_KEY", "The access key to the resource associated with the specified session is invalid."),
    VI_ERROR_INV_EVENT => ("VI_ERROR_INV_EVENT", "Specified event type is not supported by the resource."),
    VI_ERROR_INV_MECH => ("VI_ERROR_INV_MECH", "Invalid mechanism specified."),
    VI_ERROR_HNDLR_NINSTALLED => ("VI_ERROR_HNDLR_NINSTALLED", "A handler was not installed."),
    VI_ERROR_INV_HNDLR_REF => ("VI_ERROR_INV_HNDLR_REF", "The given handler reference is either invalid or was not installed."),
    VI_ERROR_INV_CONTEXT => ("VI_ERROR_INV_CONTEXT", "Specified event context is invalid."),
    VI_ERROR_QUEUE_OVERFLOW => ("VI_ERROR_QUEUE_OVERFLOW", "The event queue for the specified type has overflowed (usually due to previous events not having been closed)."),
    VI_ERROR_NENABLED => ("VI_ERROR_NENABLED", "You must be enabled for events of the specified type in order to receive them."),
    VI_ERROR_ABORT => ("VI_ERROR_ABORT", "User abort occurred during transfer."),
    VI_ERROR_RAW_WR_PROT_VIOL => ("VI_ERROR_RAW_WR_PROT_VIOL", "Violation of raw write protocol occurred during transfer."),
    VI_ERROR_RAW_RD_PROT_VIOL => ("VI_ERROR_RAW_RD_PROT_VIOL", "Violation of raw read protocol occurred during transfer."),
    VI_ERROR_OUTP_PROT_VIOL => ("VI_ERROR_OUTP_PROT_VIOL", "Device reported an output protocol error during transfer."),
    VI_ERROR_INP_PROT_VIOL => ("VI_ERROR_INP_PROT_VIOL", "Device reported an input protocol error during transfer."),
    VI_ERROR_BERR => ("VI_ERROR_BERR", "Bus error occurred during transfer."),
    VI_ERROR_IN_PROGRESS => ("VI_ERROR_IN_PROGRESS", "Unable to queue the asynchronous operation because there is already an operation in progress."),
    VI_ERROR_INV_SETUP => ("VI_ERROR_INV_SETUP", "Unable to start operation because setup is invalid (usually due to attributes being set to an inconsistent state)."),
    VI_ERROR_QUEUE_ERROR => ("VI_ERROR_QUEUE_ERROR", "Unable to queue the asynchronous operation (usually due to the I/O completion event not being enabled or insufficient space in the session's queue)."),
    VI_ERROR_ALLOC => ("VI_ERROR_ALLOC", "Insufficient system resources to perform necessary memory allocation."),
    VI_ERROR_INV_MASK => ("VI_ERROR_INV_MASK", "Invalid buffer mask specified."),
    VI_ERROR_IO => ("VI_ERROR_IO", "Could not perform operation because of I/O error."),
    VI_ERROR_INV_FMT => ("VI_ERROR_INV_FMT", "A format specifier in the format string is invalid."),
    VI_ERROR_NSUP_FMT => ("VI_ERROR_NSUP_FMT", "A format specifier in the format string is not supported."),
    VI_ERROR_LINE_IN_USE => ("VI_ERROR_LINE_IN_USE", "The specified trigger line is currently in use."),
    VI_ERROR_NSUP_MODE => ("VI_ERROR_NSUP_MODE", "The specified mode is not supported by this VISA implementation."),
    VI_ERROR_SRQ_NOCCURRED => ("VI_ERROR_SRQ_NOCCURRED", "Service request has not been received for the session."),
    VI_ERROR_INV_SPACE => ("VI_ERROR_INV_SPACE", "Invalid address space specified."),
    VI_ERROR_INV_OFFSET => ("VI_ERROR_INV_OFFSET", "Invalid offset specified."),
    VI_ERROR_INV_WIDTH => ("VI_ERROR_INV_WIDTH", "Invalid access width specified."),
    VI_ERROR_NSUP_OFFSET => ("VI_ERROR_NSUP_OFFSET", "Specified offset is not accessible from this hardware."),
    VI_ERROR_NSUP_VAR_WIDTH => ("VI_ERROR_NSUP_VAR_WIDTH", "Cannot support source and destination widths that are different."),
    VI_ERROR_WINDOW_NMAPPED => ("VI_ERROR_WINDOW_NMAPPED", "The specified session is not currently mapped."),
    VI_ERROR_RESP_PENDING => ("VI_ERROR_RESP_PENDING", "A previous response is still pending, causing a multiple query error."),
    VI_ERROR_NLISTENERS => ("VI_ERROR_NLISTENERS", "No listeners condition is detected (both NRFD and NDAC are deasserted)."),
    VI_ERROR_NCIC => ("VI_ERROR_NCIC", "The interface associated with this session is not currently the controller in charge."),
    VI_ERROR_NSYS_CNTLR => ("VI_ERROR_NSYS_CNTLR", "The interface associated with this session is not the system controller."),
    VI_ERROR_NSUP_OPER => ("VI_ERROR_NSUP_OPER", "The given session or object reference does not support this operation."),
    VI_ERROR_INTR_PENDING => ("VI_ERROR_INTR_PENDING", "An interrupt is still pending from a previous call."),
    VI_ERROR_ASRL_PARITY => ("VI_ERROR_ASRL_PARITY", "A parity error occurred during transfer."),
    VI_ERROR_ASRL_FRAMING => ("VI_ERROR_ASRL_FRAMING", "A framing error occurred during transfer."),
    VI_ERROR_ASRL_OVERRUN => ("VI_ERROR_ASRL_OVERRUN", "An overrun error occurred during transfer. A character was not read from the hardware before the next character arrived."),
    VI_ERROR_TRIG_NMAPPED => ("VI_ERROR_TRIG_NMAPPED", "The path from trigSrc to trigDest is not currently mapped."),
    VI_ERROR_NSUP_ALIGN_OFFSET => ("VI_ERROR_NSUP_ALIGN_OFFSET", "The specified offset is not properly aligned for the access width of the operation."),
    VI_ERROR_USER_BUF => ("VI_ERROR_USER_BUF", "A specified user buffer is not valid or cannot be accessed for the required size."),
    VI_ERROR_RSRC_BUSY => ("VI_ERROR_RSRC_BUSY", "The resource is valid, but VISA cannot currently access it."),
    VI_ERROR_NSUP_WIDTH => ("VI_ERROR_NSUP_WIDTH", "Specified width is not supported by this hardware."),
    VI_ERROR_INV_PARAMETER => ("VI_ERROR_INV_PARAMETER", "The value of some parameter (which parameter is not known) is invalid."),
    VI_ERROR_INV_PROT => ("VI_ERROR_INV_PROT", "The protocol specified is invalid."),
    VI_ERROR_INV_SIZE => ("VI_ERROR_INV_SIZE", "Invalid size of window specified."),
    VI_ERROR_WINDOW_MAPPED => ("VI_ERROR_WINDOW_MAPPED", "The specified session currently contains a mapped window."),
    VI_ERROR_NIMPL_OPER => ("VI_ERROR_NIMPL_OPER", "The given operation is not implemented."),
    VI_ERROR_INV_LENGTH => ("VI_ERROR_INV_LENGTH", "Invalid length specified."),
    VI_ERROR_INV_MODE => ("VI_ERROR_INV_MODE", "Invalid mode specified."),
    VI_ERROR_SESN_NLOCKED => ("VI_ERROR_SESN_NLOCKED", "The current session did not have a lock on the resource."),
    VI_ERROR_MEM_NSHARED => ("VI_ERROR_MEM_NSHARED", "The device does not export any memory."),
    VI_ERROR_LIBRARY_NFOUND => ("VI_ERROR_LIBRARY_NFOUND", "A code library required by VISA could not be located or loaded."),
    VI_ERROR_NSUP_INTR => ("VI_ERROR_NSUP_INTR", "The interface cannot generate an interrupt on the requested level or with the requested statusID value."),
    VI_ERROR_INV_LINE => ("VI_ERROR_INV_LINE", "The value specified by the line parameter is invalid."),
    VI_ERROR_FILE_ACCESS => ("VI_ERROR_FILE_ACCESS", "An error occurred while trying to open the specified file. Possible reasons include an invalid path or lack of access rights."),
    VI_ERROR_FILE_IO => ("VI_ERROR_FILE_IO", "An error occurred while performing I/O on the specified file."),
    VI_ERROR_NSUP_LINE => ("VI_ERROR_NSUP_LINE", "One of the specified lines (trigSrc or trigDest) is not supported by this VISA implementation, or the combination of lines is not a valid mapping."),
    VI_ERROR_NSUP_MECH => ("VI_ERROR_NSUP_MECH", "The specified mechanism is not supported for the given event type."),
    VI_ERROR_INTF_NUM_NCONFIG => ("VI_ERROR_INTF_NUM_NCONFIG", "The interface type is valid but the specified interface number is not configured."),
    VI_ERROR_CONN_LOST => ("VI_ERROR_CONN_LOST", "The connection for the given session has been lost."),
    VI_ERROR_MACHINE_NAVAIL => ("VI_ERROR_MACHINE_NAVAIL", "The remote machine does not exist or is not accepting any connections. If the NI-VISA server is installed and running on the remote machine, it may have an incompatible version or may be listening on a different port."),
    VI_ERROR_NPERMISSION => ("VI_ERROR_NPERMISSION", "Access to the resource or remote machine is denied. This is due to lack of sufficient privileges for the current user or machine")
)


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

""" NURBS segment. """
type NSeg <: AbstractBasis
    order :: Int
    knots :: Vector{Float64}
    weights :: Vector{Float64}
end

function NSeg()
    NSeg(1,
        [-1.0, -1.0, 1.0, 1.0],
        ones(4))
end

function length(basis::NSeg)
    nu = length(basis.knots) - basis.order - 1
    return nu
end

function size(basis::NSeg)
    return (1, length(basis))
end

function eval_basis!(basis::NSeg, N::Matrix, xi)
    pu = basis.order
    tu = basis.knots
    w = basis.weights
    nu = length(tu)-pu-1
    u = xi[1]
    for j=1:nu
        N[j] = w[j]*NURBS(j,pu,u,tu)
    end
    N ./= sum(N)
    return N
end


text:
Flux.trainable(bn::Flux.BatchNorm) = (bn.β, bn.γ, bn.μ, bn.σ²)

struct Shortcut
    shortcut
end
Flux.@functor Shortcut

function (shortcut::Shortcut)(mx::T, x::T) where T
    mx + shortcut.shortcut(x)
end

struct ResNetModel
    entry::Chain
    encoder::Chain
    head::Chain
end
Flux.@functor ResNetModel

function (model::ResNetModel)(x::T) where T
    x |> model.entry |> model.encoder |> model.head
end

function BasicBlock(
    in_filters::Int, out_filters::Int, stride::Int = 1, connection = +,
)::Chain
    layer = Chain(
        Conv((3, 3), in_filters => out_filters, stride=stride, pad=1, bias=Flux.Zeros(Float32)),
        BatchNorm(out_filters, relu),
        Conv((3, 3), out_filters => out_filters, pad=1, bias=Flux.Zeros(Float32)),
        BatchNorm(out_filters),
    )
    Chain(SkipConnection(layer, connection), x -> relu.(x))
end

function Bottleneck(
    in_filters::Int, out_filters::Int, stride::Int = 1, connection = +,
    expansion::Int = 4,
)::Chain
    layer = Chain(
        Conv((1, 1), in_filters => out_filters, bias=Flux.Zeros(Float32)),
        BatchNorm(out_filters, relu),
        Conv((3, 3), out_filters => out_filters, stride=stride, pad=1, bias=Flux.Zeros(Float32)),
        BatchNorm(out_filters, relu),
        Conv((1, 1), out_filters => out_filters * expansion, bias=Flux.Zeros(Float32)),
        BatchNorm(out_filters * expansion),
    )
    Chain(SkipConnection(layer, connection), x -> relu.(x))
end

function make_connection(in_filters::Int, out_filters::Int, stride::Int)
    if stride == 1 && in_filters == out_filters
        return +
    end
    Shortcut(Chain(
        Conv((1, 1), in_filters => out_filters, stride=stride, bias=Flux.Zeros(Float32)),
        BatchNorm(out_filters)
    ))
end

function make_layer(
    in_filters::Int, out_filters::Int, repeat::Int,
    block, expansion::Int, stride::Int = 1,
)::Chain
    connection = make_connection(in_filters, out_filters * expansion, stride)
    layer = [block(in_filters, out_filters, stride, connection)]
    in_filters = out_filters * expansion
    for i in 2:repeat
        push!(layer, block(in_filters, out_filters))
    end
    Chain(layer...)
end

"""
```julia
resnet(size::Int = 18, classes::Int = 1000)
```

Construct ResNet model.

# Parameters
- `size::Int = 18`: Size of the ResNet model.
  Available sizes are `18, 34, 50, 101, 152`.
- `classes::Int = 1000`: Number of classes in the last dense layer.
"""
function resnet(size::Int = 18, classes::Int = 1000)::ResNetModel
    filters = [64, 128, 256, 512]
    strides = [1, 2, 2, 2]
    repeats, block, expansion = Dict(
        18 => ([2, 2, 2, 2], BasicBlock, 1),
        34 => ([3, 4, 6, 3], BasicBlock, 1),
        50 => ([3, 4, 6, 3], Bottleneck, 4),
        101 => ([3, 4, 23, 3], Bottleneck, 4),
        152 => ([3, 8, 36, 3], Bottleneck, 4),
    )[size]

    entry = Chain(
        Conv((7, 7), 3 => 64, pad=(3, 3), stride=(2, 2), bias=Flux.Zeros(Float32)),
        BatchNorm(64, relu),
        MaxPool((3, 3), pad=(1, 1), stride=(2, 2)),
    )
    head = Chain(MeanPool((7, 7)), flatten, Dense(512 * expansion, classes))

    encoder = []
    in_filters = 64
    for (out_filters, repeat, stride) in zip(filters, repeats, strides)
        layer = make_layer(in_filters, out_filters, repeat, block, expansion, stride)
        push!(encoder, layer)
        in_filters = out_filters * expansion
    end
    ResNetModel(entry, Chain(encoder...), head)
end


text:
# startup.jl
#
# Startup script for FFD. This will change as the code progresses

# Includes
push!(LOAD_PATH, "../src/")
push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("SummationByParts"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("MeshMovement"), "src"))

using PdePumiInterface
using SummationByParts
using ODLCommonTools
using ArrayViews
using MPI
using MeshMovement
using FreeFormDeformation

# MPI Declarations
comm = MPI.COMM_WORLD
comm_world = MPI.MPI_COMM_WORLD
comm_self = MPI.COMM_SELF
my_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

opts = PdePumiInterface.get_defaults()
# 2D mesh
opts["numBC"] = 2
opts["BC1"] = [8,11,14,17]
opts["BC1_name"] = "FarField"
opts["BC2"] = [5]
opts["BC2_name"] = "Airfoil"
opts["coloring_distance"] = 2 # 0 For CG Mesh, 2 for DG Mesh

Tmsh = Float64
dmg_name = ".null"
smb_name = "./mesh_files/gvortex1.smb"
order = 1
dofpernode = 1

# SBP & Mesh Parameters
# For SBP Gamma
reorder = false
internal = true
shape_type = 2

# create linear sbp operator
sbp = TriSBP{Float64}(degree=order, reorder=reorder, internal=internal)
ref_verts = [-1. 1 -1; -1 -1 1]
interp_op = SummationByParts.buildinterpolation(sbp, ref_verts)
sbpface = TriFace{Float64}(order, sbp.cub, ref_verts.')
# create linear mesh with 4 dof per node

println("constructing DG mesh")
mesh = PumiMeshDG2{Tmsh}(dmg_name, smb_name, order, sbp, opts, interp_op, sbpface;
       dofpernode=dofpernode, coloring_distance=opts["coloring_distance"],
       shape_type=shape_type)
#
#
# Free Form deformation parameters
ndim = 2
order = [4,4,2]  # Order of B-splines in the 3 directions
nControlPts = [4,4,2]

ffd_map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh)
calcKnot(ffd_map)
# println("ffd_map.edge_knot = \n", ffd_map.edge_knot)

# Create Bounding box
offset = [0., 0., 0.5]
ffd_box = PumiBoundingBox{Tmsh}(ffd_map, mesh, sbp, offset)

# println("ffd_box.origin = ", ffd_box.origin)

controlPoint(ffd_map, ffd_box)
#=
for k = 1:ffd_map.nctl[3]
  for j = 1:ffd_map.nctl[2]
    for i = 1:ffd_map.nctl[1]
      println("cp_xyz[3,$i,$j,$k] = ", ffd_map.cp_xyz[3,i,j,k])
    end
    println('\n')
  end
end
=#

# Create Linear Mapping
calcParametricMappingNonlinear(ffd_map, ffd_box, mesh)
#println("ffd_map.xi = \n", ffd_map.xi)

# Translate control points along x & y by  5 units
ffd_map.cp_xyz[1,:,:,:] += 2
ffd_map.cp_xyz[2,:,:,:] += 3
#=
# Check angular rotation
# Rotate control points by 90 degrees about Z axis
gamma = 0.5*pi # angle in radians by which ctls are rotated
rotx = [1. 0. 0.;0. cos(gamma) -sin(gamma);0. sin(gamma) cos(gamma)]
roty = [cos(gamma) 0. sin(gamma);0. 1. 0.;-sin(gamma) 0. cos(gamma)]
rotz = [cos(gamma) -sin(gamma) 0;sin(gamma) cos(gamma) 0;0 0 1]
# Perform the rotation
for k = 1:ffd_map.nctl[3]
  for j = 1:ffd_map.nctl[2]
    for i = 1:ffd_map.nctl[1]
      ffd_map.cp_xyz[:,i,j,k] = rotz*ffd_map.cp_xyz[:,i,j,k]
    end
  end
end
=#
evalVolume(ffd_map, mesh)

# println("mesh.coords = \n", mesh.coords)

for i = 1:mesh.numEl
  update_coords(mesh, i, mesh.vert_coords[:,:,i])
end
commit_coords(mesh)
PumiInterface.writeVtkFiles("Translation_noUpdate", mesh.m_ptr)

# geom_bounds = zeros(2,3)
# FreeFormDeformation.calcGeomBounds(mesh.coords, geom_bounds)

#=
map = Mapping(ndim, order, nControlPts, nnodes)
calcKnot(map)  # Create knot vectors

for i = 1:3
  println("map.edge_knot[$i] = ", map.edge_knot[i])
end

# Create Bounding Box
offset = [0.5,0.5,0.5]  # offset for the bounding box
geom_bounds = [1. 1. 1.;3. 3. 3.]
box = BoundingBox(ndim, geom_bounds, offset)
println("origin = $(box.origin)")

# Unit vectors
# box.unitVector = eye(Float64, 3)

println("box.lengths = $(box.lengths)")
println("box.box_bound = $(box.box_bound)\n\n")

#------------------------------------
# Representative obeject
# nnodes = [5,3,4]
nodes_xyz = zeros(nnodes[1], nnodes[2], nnodes[3], 3)
origin = [1,1,1]
incx = 0.0
incy = 0.0
incz = 0.0
for k = 1:nnodes[3]
  incy = 0.0
  for j = 1:nnodes[2]
    incx = 0.0
    for i = 1:nnodes[1]
      nodes_xyz[i,j,k,1] = origin[1] + incx
      nodes_xyz[i,j,k,2] = origin[2] + incy
      nodes_xyz[i,j,k,3] = origin[3] + incz
      incx += 1
    end
    incy += 1
  end
  incz += 1
end

# Output the values to check for errors
for k = 1:nnodes[3]
  for j = 1:nnodes[2]
    for i = 1:nnodes[1]
      println("nodes_xyz[$i,$j,$k,:] = ", nodes_xyz[i,j,k,:])
    end
    println('\n')
  end
end
#----------------------------------------

# Define the control points
controlPoint(map,box)

for k = 1:map.nctl[3]
  for j = 1:map.nctl[2]
    for i = 1:map.nctl[1]
      println("cp_xyz[$i,$j,$k,:] = ", round(map.cp_xyz[i,j,k,:],2) )
    end
    println('\n')
  end
end

# calcParametricMappingLinear(map, box, nodes_xyz)
calcParametricMappingNonlinear(map, box, nodes_xyz)

for k = 1:nnodes[3]
  for j = 1:nnodes[2]
    for i = 1:nnodes[1]
      println("nodes_stu[$i,$j,$k,:] = ", map.xi[i,j,k,:])
    end
  end
end

Vol = zeros(nodes_xyz)
evalVolume(map, Vol)

for k = 1:nnodes[3]
  for j = 1:nnodes[2]
    for i = 1:nnodes[1]
      println("nodes_xyz_err[$i,$j,$k,:] = ", round(Vol[i,j,k,:]-nodes_xyz[i,j,k,:],2))
    end
    println('\n')
  end
  println('\n')
end

# Check linear scaling (Translation)
fill!(Vol, 0.0)
# move the x coordinates of all the control by 4 units
map.cp_xyz[:,:,:,3] += 2
evalVolume(map, Vol)
# Output the values to check for errors
for k = 1:nnodes[3]
  for j = 1:nnodes[2]
    for i = 1:nnodes[1]
      println( "Vol_err[$i,$j,$k,:] = ", round(Vol[i,j,k,:]-nodes_xyz[i,j,k,:],2) )
    end
    println('\n')
  end
  println('\n')
end

# Check angular rotation
# Rotate control points by 90 degrees about Z axis
gamma = 0.5*pi # angle in radians by which ctls are rotated
rotx = [1. 0. 0.;0. cos(gamma) -sin(gamma);0. sin(gamma) cos(gamma)]
roty = [cos(gamma) 0. sin(gamma);0. 1. 0.;-sin(gamma) 0. cos(gamma)]
rotz = [cos(gamma) -sin(gamma) 0;sin(gamma) cos(gamma) 0;0 0 1]
controlPoint(map,box) # Recreate unperturbed control points
# Perform the rotation
for k = 1:map.nctl[3]
  for j = 1:map.nctl[2]
    for i = 1:map.nctl[1]
      map.cp_xyz[i,j,k,:] = rotx*map.cp_xyz[i,j,k,:]
    end
  end
end

println("Rotated coordinates\n")
for k = 1:map.nctl[3]
  for j = 1:map.nctl[2]
    for i = 1:map.nctl[1]
      println("cp_xyz[$i,$j,$k,:] = ", round(map.cp_xyz[i,j,k,:],2) )
    end
    println('\n')
  end
end

evalVolume(map, Vol)
# Check Output
for k = 1:nnodes[3]
  for j = 1:nnodes[2]
    for i = 1:nnodes[1]
      println( "Vol[$i,$j,$k,:] = ", round(Vol[i,j,k,:],2) )
    end
    println('\n')
  end
  println('\n')
end
=#


text:
module ProBeTune

import Base.@kwdef

using OrdinaryDiffEq
using DiffEqFlux
using Statistics
using Random

include("samplers.jl")

export PBTProblemParameters
"""
    PBTProblemParameters(; kw...)

## Keywords:
- `N_samples::Int`: Number of samples to use.
- `size_p_spec::Int`: Size of the parameter space of the spec.
- `size_p_sys::Int`: Size of the parameter space of the system.
- `output_metric`: Metric that calculates how far the outputs of system and spec are apart.
    Called as `output_metric(sol_system, sol_spec)`, where sol is a DifferentialEquations solution object.
"""
Base.@kwdef mutable struct PBTProblemParameters
    N_samples::Int
    size_p_spec::Int
    size_p_sys::Int
    output_metric
end

export AbstractPBTProblem
"""
    AbstractPBTProblem

Supertype type to define PBTProblems.

Each concrete subtype has to have the fields
- p::PBTProblemParameters
- input_sampler::AbstractInputSampler

Each concrete subtype has to implement

    solve_sys_spec(pbt::ConcretePBTProblem, input, p_sys, p_spec) -> sol_sys, sol_spec

which returns a tuple of two DiffEq solution objects: one for the sys, one for the spec.
"""
abstract type AbstractPBTProblem end

###################
# Definition of the default 'rhs' PBT problem
###################
export PBTProblem
"""
    PBTProblem <: AbstractPBTProblem

Characterize problems by their rhs.
"""
mutable struct PBTProblem <: AbstractPBTProblem
    p::PBTProblemParameters
    input_sampler::ContinousInputSampler
    f_spec
    f_sys
    y0_spec
    y0_sys
    tsteps
    t_span
    solver
end

"""
    PBTProblem(kw...)

Keyword constructor for default PBT Problem.

## Keywords:
- `N_samples::Int`: Number of samples to use.
- `size_p_spec::Int`: Size of the parameter space of the spec.
- `size_p_sys::Int`: Size of the parameter space of the system.
- `output_metric`: Metric that calculates how far the outputs of system and spec are apart.
    Called as `output_metric(sol_system, sol_spec)`, where sol is a DifferentialEquations solution object.
- `input_sampler`: `ContinousInputSampler`
- `f_spec`: Right hand side function of the system, called as f(dy, y, i, p, t).
- `f_sys`: Right hand side function of the specification, called as f(dy, y, i, p, t).
- `y0_spec`: Initial conditions for integrating the spec.
- `y0_sys`: Initial conditions for integrating the system.
- `tsteps`: Times are which to save the solution
- `t_span`: Start and end time of the integration
- `solver = Tsit5()`: Solver for integrating spec and system.
"""
function PBTProblem(; N_samples, size_p_spec, size_p_sys, output_metric,
                    input_sampler, f_spec, f_sys, y0_spec, y0_sys, tsteps, t_span, solver = Tsit5())
    p = PBTProblemParameters(; N_samples, size_p_spec, size_p_sys, output_metric)
    PBTProblem(p, input_sampler, f_spec, f_sys, y0_spec, y0_sys, tsteps, t_span, solver)
end

"""
    resample!(pbt::AbstractPBTProblem; n = pbt.p.N_samples)

Resample the inputs.
"""
function resample!(pbt::AbstractPBTProblem; n = pbt.p.N_samples)
    pbt.p.N_samples = n
    resample!(pbt.input_sampler)
end

###################
# Generating solutions
# The following functions generate solutions
###################

export solve_sys_spec
"""
    solve_sys_spec(pbt::PBTProblem, i, p_sys, p_spec; solver_options...)

Solve system and specification for given parameters `p_sys` and `p_spec` and input `i`.
Wraps the `f_sys` and `f_spec` rhs functions as an ODEProblem function. `solver_options` kw arguments get
passed to both solve calls.
"""
function solve_sys_spec(pbt::PBTProblem, i, p_sys, p_spec; solver_options...)
    sol_sys = solve(ODEProblem((dy,y,p,t) -> pbt.f_sys(dy, y, i(t), p, t), pbt.y0_sys, pbt.t_span, p_sys), pbt.solver; saveat=pbt.tsteps, solver_options...)
    sol_spec = solve(ODEProblem((dy,y,p,t) -> pbt.f_spec(dy, y, i(t), p, t), pbt.y0_spec, pbt.t_span, p_spec), pbt.solver; saveat=pbt.tsteps, solver_options...)
    sol_sys, sol_spec
end

export solve_sys_spec_n
"""
    solve_sys_spec_n(pbt::AbstractPBTProblem, n::Int, p; solver_options...)

Solve system and spec for the n-th input sample given a stacked parameter array p.
Parameters:
- `pbt`: existing AbstractPBTProblem instance
- `n`: number of input in the array
- `p`: combined array of sys and spec parameters
"""
function solve_sys_spec_n(pbt::AbstractPBTProblem, n::Int, p; solver_options...)
    p_sys = view(p,1:pbt.p.size_p_sys)
    p_spec = view(p,(pbt.p.size_p_sys + 1 + (n - 1) * pbt.p.size_p_spec):(pbt.p.size_p_sys + n * pbt.p.size_p_spec))

    i = pbt.input_sampler[n]

    solve_sys_spec(pbt, i, p_sys, p_spec; solver_options...)
end

###################
# Loss functions
# The following functions are loss functions actually used to
# tune the system and evaluate behavioural distances.
###################

export pbt_loss

"""
    pbt_loss(pbt::AbstractPBTProblem, p; solver_options...)

pbt_loss calculates the loss function for the Probabilistic Tuning Problem.
It will add up the loss provided by the `output_metric` for all samples.
"""
function pbt_loss(pbt::AbstractPBTProblem, p; solver_options...)
    loss = 0.

        for n in 1:pbt.p.N_samples
            sol_sys, sol_spec = solve_sys_spec_n(pbt, n, p; solver_options...)
            loss += pbt.p.output_metric(sol_sys, sol_spec)
        end

    loss / pbt.p.N_samples
end

"""
    pbt_individual_loss(pbt::AbstractPBTProblem, n, p_sys, p_spec; solver_options...)

pbt_individual_loss provides the contribution to the loss function for the Probabilistic Tuning Problem for a single sample.
"""
function pbt_individual_loss(pbt::AbstractPBTProblem, n, p_sys, p_spec; solver_options...)
  # loss function for sample n only
  i = pbt.input_sampler[n]
  pbt.p.output_metric(solve_sys_spec(pbt, i, p_sys, p_spec; solver_options...)...)
end

###################
# Behavioural Distance
# For a fixed system we can provide a given p_sys.
###################

export behavioural_distance

"""
    behavioural_distance(pbt::AbstractPBTProblem, p; verbose=true, optimizer=DiffEqFlux.ADAM(0.01),
                         optimizer_options=(:maxiters => 1000,), solver_options...)

Tune the specification parameters only to calculate the behvaioural distance of the system
to the specification.

# Arguments:
- `pbt`: PBT problem
- `p`: stacked array of system parameters and initial parameters for the specifications
- `optimizer`: choose optimization algorithm (default `DiffEqFlux.ADAM(0.01)`)
"""
function behavioural_distance(pbt::AbstractPBTProblem, p; verbose=true, optimizer=DiffEqFlux.ADAM(0.01),
                              optimizer_options=(:maxiters => 1000,), solver_options...)
  # Optimize the specs only

  p_tuned = copy(p)

  @views begin
      p_sys = p_tuned[1:pbt.p.size_p_sys]
      p_specs = [p_tuned[pbt.p.size_p_sys + 1 + (n - 1) * pbt.p.size_p_spec:pbt.p.size_p_sys + n * pbt.p.size_p_spec] for n in 1:pbt.p.N_samples]
  end

  verbose && print("0 out of $(pbt.p.N_samples) samples evaluated")

  distances = zeros(pbt.p.N_samples)

  for n in 1:pbt.p.N_samples
    res = DiffEqFlux.sciml_train(
        x -> pbt_individual_loss(pbt, n, p_sys, x; solver_options...),
        Array(p_specs[n]),
        optimizer;
        optimizer_options...
        )
    verbose && print("\r$n out of $(pbt.p.N_samples) samples evaluated")
    p_specs[n] .= res.minimizer
    distances[n] = res.minimum
  end
  print("\n")
  mean(distances), p_tuned, distances
end

###################
# Tuning the system
###################

export pbt_tuning

"""
    pbt_tuning(pbt::AbstractPBTProblem, p; optimizer=DiffEqFlux.ADAM(0.01), optimizer_options=(:maxiters => 100,), solver_options...)

Tune the system to the specification.

# Arguments:
- `pbt`: PBT problem
- `p`: stacked array of initial system and specification parameters
- `optimizer`: choose optimization algorithm (default `DiffEqFlux.ADAM(0.01)`)
- `optimizer_options`: choose optimization options (default `(:maxiters => 100,)`)
- `solver_options...` all further options are passed through to the differential equations solver
"""
function pbt_tuning(pbt::AbstractPBTProblem, p; optimizer=DiffEqFlux.ADAM(0.01), optimizer_options=(:maxiters => 100,), solver_options...)
  DiffEqFlux.sciml_train(
    x -> pbt_loss(pbt, x; solver_options...),
    p,
    optimizer;
    optimizer_options...
    )
end

#############################
#############################
#############################
#############################
# Library of further things

export PBTLibrary
include("PBTLibrary.jl")

export PBTExample
include("PBTExample.jl")

end


text:
struct WeekAheadIndices <: AbstractModelIndices
    hours::Vector{Int}
    plants::Vector{Plant}
    segments::Vector{Int}
end
plants(indices::WeekAheadIndices) = indices.plants

function HydroModels.modelindices(data::WeekAheadData, horizon::Horizon; areas::Vector{Area} = [0], rivers::Vector{River} = [:All])
    hours = collect(1:nhours(horizon))
    plants = plants_in_areas_and_rivers(hydrodata(data),areas,rivers)
    if isempty(plants)
        error("No plants in given set of price areas and rivers")
    end
    segments = collect(1:2)
    return WeekAheadIndices(hours, plants, segments)
end


text:
function handle_callbacks!(integrator)
    discrete_callbacks = integrator.opts.callback.discrete_callbacks
    continuous_callbacks = integrator.opts.callback.continuous_callbacks
    atleast_one_callback = false

    continuous_modified = false
    discrete_modified = false
    saved_in_cb = false
    if !(typeof(continuous_callbacks) <: Tuple{})
        time, upcrossing, event_occured, event_idx, idx, counter =
            DiffEqBase.find_first_continuous_callback(integrator, continuous_callbacks...)
        if event_occured
            integrator.event_last_time = idx
            integrator.vector_event_last_time = event_idx
            continuous_modified, saved_in_cb = DiffEqBase.apply_callback!(
                integrator,
                continuous_callbacks[idx],
                time,
                upcrossing,
                event_idx,
            )
        else
            integrator.event_last_time = 0
            integrator.vector_event_last_time = 1
        end
    end
    if !(typeof(discrete_callbacks) <: Tuple{})
        discrete_modified, saved_in_cb =
            DiffEqBase.apply_discrete_callback!(integrator, discrete_callbacks...)
    end

    integrator.u_modified = continuous_modified || discrete_modified
    if integrator.u_modified
        handle_callback_modifiers!(integrator)
    end

    if !saved_in_cb
        savevalues!(integrator)
    end

    integrator.u_modified = false
end

function DiffEqBase.savevalues!(
    integrator::AbstractSundialsIntegrator,
    force_save = false,
)::Tuple{Bool, Bool}
    saved, savedexactly = false, false
    !integrator.opts.save_on && return saved, savedexactly
    uType = typeof(integrator.sol.prob.u0)
    # The call to first is an overload of Base.first implemented in DataStructures
    while !isempty(integrator.opts.saveat) &&
        integrator.tdir * first(integrator.opts.saveat) < integrator.tdir * integrator.t
        saved = true
        curt = pop!(integrator.opts.saveat)

        tmp = integrator(curt)
        save_value!(integrator.sol.u, tmp, uType, integrator.sizeu, integrator.opts.save_idxs, Val{false})
        push!(integrator.sol.t, curt)
        if integrator.opts.dense
            tmp = integrator(curt, Val{1})
            save_value!(integrator.sol.interp.du, tmp, uType, integrator.sizeu, integrator.opts.save_idxs, Val{false})
        end
    end

    if force_save || integrator.opts.save_everystep || (integrator.opts.save_everystep && (isempty(integrator.sol.t) ||
                                         (integrator.t !== integrator.sol.t[end])))
        saved = true
        save_value!(integrator.sol.u, integrator.u, uType, integrator.sizeu, integrator.opts.save_idxs)
        push!(integrator.sol.t, integrator.t)
        if integrator.opts.dense
            tmp = integrator(integrator.t, Val{1})
            save_value!(integrator.sol.interp.du, tmp, uType, integrator.sizeu, integrator.opts.save_idxs)
        end
    end
    savedexactly = !isempty(integrator.sol.t) && last(integrator.sol.t) == integrator.t
    return saved, savedexactly
end

function save_value!(
    save_array,
    val,
    ::Type{T},
    sizeu, save_idxs,
    make_copy::Type{Val{bool}} = Val{true},
) where {T <: Number, bool}
    push!(save_array, first(val))
end
function save_value!(
    save_array,
    val,
    ::Type{T},
    sizeu, save_idxs,
    make_copy::Type{Val{bool}} = Val{true},
) where {T <: Vector, bool}
    save = if save_idxs !== nothing
        val[save_idxs]
    else
        bool ? copy(val) : val
    end
    push!(save_array, save)
end
function save_value!(
    save_array,
    val,
    ::Type{T},
    sizeu, save_idxs,
    make_copy::Type{Val{bool}} = Val{true},
) where {T <: AbstractArray, bool}
    save = if save_idxs !== nothing
        reshape(val, sizeu)[save_idxs]
    else
        x = bool ? copy(val) : val
        reshape(x, sizeu)
    end
    push!(save_array, save)
end

function handle_callback_modifiers!(integrator::CVODEIntegrator)
    CVodeReInit(integrator.mem, integrator.t, integrator.u)
end

function handle_callback_modifiers!(integrator::ARKODEIntegrator)
    ARKStepReInit(integrator.mem, integrator.t, integrator.u)
end

function handle_callback_modifiers!(integrator::IDAIntegrator)
    IDAReInit(integrator.mem, integrator.t, integrator.u, integrator.du)
    DiffEqBase.initialize_dae!(integrator)
end

function DiffEqBase.add_tstop!(integrator::AbstractSundialsIntegrator, t)
    integrator.tdir * (t - integrator.t) < 0 &&
        error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
    push!(integrator.opts.tstops, t)
end

function DiffEqBase.add_saveat!(integrator::AbstractSundialsIntegrator, t)
    integrator.tdir * (t - integrator.t) < 0 &&
        error("Tried to add a saveat that is behind the current time. This is strictly forbidden")
    push!(integrator.opts.saveat, t)
end

DiffEqBase.get_tmp_cache(integrator::AbstractSundialsIntegrator) = (integrator.tmp,)

@inline function DiffEqBase.u_modified!(integrator::AbstractSundialsIntegrator, bool::Bool)
    integrator.u_modified = bool
end

function DiffEqBase.terminate!(
    integrator::AbstractSundialsIntegrator,
    retcode = :Terminated,
)
    integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, retcode)
    integrator.opts.tstops.valtree = typeof(integrator.opts.tstops.valtree)()
end

@inline function DiffEqBase.get_du(integrator::CVODEIntegrator)
    integrator(integrator.t, Val{1})
end

@inline function DiffEqBase.get_du!(out, integrator::CVODEIntegrator)
    integrator(out, integrator.t, Val{1})
end

@inline function DiffEqBase.get_du(integrator::IDAIntegrator)
    reshape(integrator.du, integrator.sizedu)
end

@inline function DiffEqBase.get_du!(out, integrator::IDAIntegrator)
    out .= reshape(integrator.du, integrator.sizedu)
end

function DiffEqBase.change_t_via_interpolation!(integrator::AbstractSundialsIntegrator, t)
    integrator.t = t
    integrator(integrator.u, integrator.t)
    return nothing
end

@inline function Base.getproperty(integrator::AbstractSundialsIntegrator, sym::Symbol)
    if sym == :dt
        return integrator.t - integrator.tprev
    else
        return getfield(integrator, sym)
    end
end

DiffEqBase.reeval_internals_due_to_modification!(integrator::AbstractSundialsIntegrator) =
    nothing
DiffEqBase.reeval_internals_due_to_modification!(integrator::IDAIntegrator) =
    handle_callback_modifiers!(integrator::IDAIntegrator)

# Required for callbacks
DiffEqBase.set_proposed_dt!(i::AbstractSundialsIntegrator, dt) = nothing

DiffEqBase.initialize_dae!(integrator::AbstractSundialsIntegrator) = nothing
function DiffEqBase.initialize_dae!(integrator::IDAIntegrator)
    integrator.f(integrator.tmp, integrator.du, integrator.u, integrator.p, integrator.t)
    if any(abs.(integrator.tmp) .>= integrator.opts.reltol)
        if integrator.sol.prob.differential_vars === nothing && !integrator.alg.init_all
            error("Must supply differential_vars argument to DAEProblem constructor to use IDA initial value solver.")
        end
        if integrator.alg.init_all
            init_type = IDA_Y_INIT
        else
            init_type = IDA_YA_YDP_INIT
            integrator.flag =
                IDASetId(integrator.mem, integrator.sol.prob.differential_vars)
        end
        integrator.flag = IDACalcIC(integrator.mem, init_type, integrator.dt)
    end
end

DiffEqBase.has_reinit(integrator::AbstractSundialsIntegrator) = true
function DiffEqBase.reinit!(integrator::AbstractSundialsIntegrator,u0 = integrator.sol.prob.u0;
  t0 = integrator.sol.prob.tspan[1], tf = integrator.sol.prob.tspan[2],
  erase_sol = true,
  tstops = integrator.opts.tstops_cache,
  saveat = integrator.opts.saveat_cache,
  reinit_callbacks = true, initialize_save = true,
  reinit_cache = true)

  if isinplace(integrator.sol.prob)
    copyto!(integrator.u,u0)
    copyto!(integrator.uprev,integrator.u)
  else
    integrator.u = u0
    integrator.uprev = integrator.u
  end

  integrator.t = t0
  integrator.tprev = t0

  tType = typeof(integrator.t)
  tspan = (tType(t0), tType(tf))
  tdir = sign(tspan[2] - tspan[1])

  tstops_internal, saveat_internal =
      tstop_saveat_disc_handling(tstops, saveat, tdir, tspan, eltype(integrator.t))

  integrator.opts.tstops = tstops_internal
  integrator.opts.saveat = saveat_internal

  if erase_sol
    if integrator.opts.save_start
      resize_start = 1
    else
      resize_start = 0
    end
    resize!(integrator.sol.u,resize_start)
    resize!(integrator.sol.t,resize_start)
    if integrator.opts.save_start
      integrator.sol.t[1] = t0
      if integrator.opts.save_idxs === nothing
        integrator.sol.u[1] = u0
      else
        u_initial = u0[integrator.opts.save_idxs]
        integrator.sol.u[1] = u_initial
      end
    end
    if integrator.sol.u_analytic !== nothing
      resize!(integrator.sol.u_analytic,0)
    end
  end

  integrator.u_modified = false
  handle_callback_modifiers!(integrator)

  if reinit_callbacks
    initialize_callbacks!(integrator, initialize_save)
  end

  nothing
end


text:
using MixedModels, Test

@testset "GHnorm" begin
    gh2 = GHnorm(2)
    @test gh2.z  == [-1.0, 1.0]
    @test gh2.w == [0.5, 0.5]
    @test GHnorm(2) === gh2
    gh9 = GHnorm(9.0)
    @test sum(gh9.w) ≈ 1
end


text:
# ---
# title: 212. Word Search II
# id: problem212
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Backtracking, Trie
# link: <https://leetcode.com/problems/word-search-ii/description/>
# hidden: true
# ---
#
# Given an `m x n` `board` of characters and a list of strings `words`, return
# _all words on the board_.
#
# Each word must be constructed from letters of sequentially adjacent cells,
# where **adjacent cells** are horizontally or vertically neighboring. The same
# letter cell may not be used more than once in a word.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/07/search1.jpg)
#
#
#
#     Input: board = [["o","a","a","n"],["e","t","a","e"],["i","h","k","r"],["i","f","l","v"]], words = ["oath","pea","eat","rain"]
#     Output: ["eat","oath"]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/07/search2.jpg)
#
#
#
#     Input: board = [["a","b"],["c","d"]], words = ["abcb"]
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `m == board.length`
#   * `n == board[i].length`
#   * `1 <= m, n <= 12`
#   * `board[i][j]` is a lowercase English letter.
#   * `1 <= words.length <= 3 * 104`
#   * `1 <= words[i].length <= 10`
#   * `words[i]` consists of lowercase English letters.
#   * All the strings of `words` are unique.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using HDF5

function write_hdf5!(group, dict::AbstractDict)
    for (name, value) in dict
        name = String(name)
        write_value_hdf5!(group, name, value)
    end
end

function write_hdf5!(group, df::AbstractDataFrame)
    for (name, value) in zip(names(df), eachcol(df))
        write_value_hdf5!(group, String(name), value)
    end
end

function write_value_hdf5!(group, name::String, value)
    # for unknown datatypes write a string representation of the value as an attribute
    attributes(group)[name] = repr(value)
end

function write_value_hdf5!(group, name::String, value::Union{AbstractDict, AbstractDataFrame})
    # if we can't write the value directly as a dataset we fall back
    # to creating a new group
    newgroup = create_group(group, name)
    write_hdf5!(newgroup, value)
end

# non array values are written as attributes
function write_value_hdf5!(group, name::String, value::Union{String,Number})
    attributes(group)[name] = value
end

function write_value_hdf5!(group, name::String, value::AbstractArray{<:Number})
    group[name] = value
end

function write_value_hdf5!(group, name::String, value::AbstractVector{<:Array{T,N}}) where {T,N}
    outer_len = length(value)
    if outer_len < 1
        group[name] = zeros(T, 0)
        return
    end
    inner_size = size(value[1])
    dset = create_dataset(group, name, datatype(T), dataspace(inner_size..., outer_len), chunk=(inner_size..., 1))
    for (i, subarray) in enumerate(value)
        dset[axes(subarray)..., i] = subarray
    end
end


text:
@testset "unconstrained" begin

    using STMO.Unconstrained

    P = [10 1; 1 5]
    q = [-5, 7]
    r = 2.0

    f(x) = 0.5x' * P * x + q' * x + r
    ∇f(x) = P * x .+ q
    ∇²f(x) = P

    xstar = - P \ q
    x0 = zero(xstar)

    @test isapprox(Unconstrained.gradient_descent(f, copy(x0), ∇f, ν=1e-6), xstar, atol=1e-5)
    @test isapprox(coordinate_descent(f, copy(x0), ∇f, ν=1e-6), xstar, atol=1e-5)
    @test isapprox(newtons_method(f, copy(x0), ∇f, ∇²f, ϵ=1e-6), xstar, atol=1e-5)
end


text:
__precompile__()
module Dagger

using Compat

import Base.collect

include("lib/util.jl")
include("lib/logging.jl")

# Distributed data
include("processor.jl")
include("thunk.jl")
include("domain.jl")
include("chunks.jl")

# Task scheduling
include("compute.jl")

# Array computations
include("array/darray.jl")
include("array/alloc.jl")
include("array/map-reduce.jl")

# File IO
include("file-io.jl")

include("array/operators.jl")
include("array/getindex.jl")
include("array/setindex.jl")
include("array/matrix.jl")
include("array/sparse_partition.jl")
include("array/sort.jl")

include("array/show.jl")

include("ui/graph.jl")

end # module


text:
using Test
using Laplacians
using LinearAlgebra
using Random
using SparseArrays
using LinearAlgebra
using Statistics
using RandomV06

L = Laplacians
include("testIncludes.jl")
include("test_graph_generators.jl")

include("testByFile.jl")

include("testByExport.jl")

include("testDrawing.jl")

include("testSolvers.jl")

include("testPCG.jl")


text:
using GeometryBasics
using Makie
using FractalTools
using Random
Random.seed!(0)

# Define the triangle
trig = Triangle(
    Point(BigFloat(-1.), BigFloat(-1.)),
    Point(BigFloat(1.), BigFloat(-1.)),
    Point(BigFloat(0.), BigFloat(2))
)

# Get valid points in the domain
pnts2d = disperse(trig, 10)

# Define function
f(x, y) = x^2 + y^2

# Evaluate function
pnts3d = [Point(pnt[1], pnt[2], f(pnt[1], pnt[2])) for pnt in pnts2d]  # Interpolation data

# Check transformations
_, msh3 = triangulate(pnts3d)
outtrig = findouttriangle(pnts3d)
outtrigpnts2d = project(outtrig.points)

# Add observables for interactivity
i = Node(39)

ttrig = @lift(msh3[$i])

ttrigpnts2d = lift(ttrig) do ttrig
    project(ttrig.points)
end

mtrig = lift(ttrig) do ttrig
    tfm = gettransform(outtrig, ttrig, 0.01)
    L(x) = tfm.A[1:2, 1:2] * x + tfm.b[1:2]
    Triangle(map(pnt -> Point(L(pnt)...), outtrigpnts2d)...)
end

# Plot figures
fig = Figure()
ls = LScene(fig[1, 1])
buttongrid = fig[1, 2] =  GridLayout(tellheight=false)
button1 = buttongrid[1, 1] = Button(fig, label = "Inc")
button2 = buttongrid[2, 1] = Button(fig, label = "Dec")
label = buttongrid[3, 1] = Label(fig, "Tri idx:$(string(i[]))")
mesh!(ls, outtrig)
wireframe!(outtrigpnts2d, color=:black, linewidth=3)
mesh!(ttrig, color=:blue)
wireframe!(ttrigpnts2d, color=:blue, linewidth=3)
mesh!(mtrig, color=:green)

# Button intrections
on(button1.clicks) do n
    i[] += 1
    label.text[] = "Tri idx: " * string(i[])
end

on(button2.clicks) do n
    i[] -= 1
    label.text[] = "Tri idx: " * string(i[])
end

display(fig)