text:
# Toboggan Trajectory

test_terrain_str = "..##.......
#...#...#..
.#....#..#.
..#.#...#.#
.#...##..#.
..#.##.....
.#.#.#....#
.#........#
#.##...#...
#...##....#
.#..#...#.#"

test_terrain = split(test_terrain_str, '\n')

terrain = readlines("day3_input.txt")

function count_trees(terrain, right, down)
    col = 1
    trees = 0
    len = length(terrain[1])
    for line in terrain[1:down:end]
        if line[col] == '#'
            trees += 1
        end
        col += right
        if col > len
            col = col % len
        end
    end
    return trees
end

println("Found $(count_trees(test_terrain, 3, 1)) trees in the test example")

println("Found $(count_trees(terrain, 3, 1)) trees in the first test")

println()

function count_second_test(terrain)
    mult = 1
    for (right, down) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]
        count = count_trees(terrain, right, down)
        mult *= count
        println("Found $count trees in the test with right = $right and down = $down")
    end
    println()
    println("Multiplying all trees gives us $mult")
end

count_second_test(terrain)

function count_second_test_noprint(terrain)
    mult = 1
    for (right, down) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]
        count = count_trees(terrain, right, down)
        mult *= count
    end
    return mult
end

#=
# I got
```julia
julia> @btime count_second_test_noprint(terrain)
  10.658 μs (7 allocations: 12.38 KiB)
1744787392
```
and no gain in using a single-liner:
```julia
julia> @btime prod([count_trees(terrain, right, down)
                   for (right, down)
                       in ((1,1), (3, 1), (5, 1), (7, 1), (1, 2))])
  10.859 μs (6 allocations: 12.33 KiB)
```
 =#

function count_trees_2(terrain, right, down)
    col = 1
    trees = 0
    len = length(terrain[1])
    for l in 1:down:length(terrain)
        if terrain[l][col] == '#'
            trees += 1
        end
# no gain in speed if we do instead (same time)
#        test = terrain[l][col] == '#' ? trees += 1 : nothing
        col += right
        if col > len
            col = col % len
        end
    end
    return trees
end

function count_second_test_noprint_2(terrain)
    mult = 1
    for (right, down) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]
        count = count_trees_2(terrain, right, down)
        mult *= count
    end
    return mult
end

# Great improvement (half the time) if we index the line
# instead of loop over the lines themselves, surprisingly
#= ```julia
julia> @btime count_second_test_noprint_2(terrain)
  4.534 μs (2 allocations: 176 bytes)
1744787392
```
=#

# Wait! Using iterator is way faster, but it doesn't influence
# much since is it an outside small loop, but in other situations
# it might make a bit difference

function count_second_test_noprint_2b(terrain)
    mult = 1
    for (right, down) in ((1,1), (3, 1), (5, 1), (7, 1), (1, 2))
        count = count_trees_2(terrain, right, down)
        mult *= count
    end
    return mult
end
# We get
#=
```julia
julia> @btime count_second_test_noprint_2(terrain)
  4.528 μs (2 allocations: 176 bytes)
1744787392

julia> @btime count_second_test_noprint_2b(terrain)
  4.465 μs (1 allocation: 16 bytes)
1744787392
```
 =#

###

function count_trees_4(terrain, right, down)
    trees = 0
    len = length(terrain[1])
    for l in 1:down:length(terrain)
        if terrain[l][(div(right*(l-1),down))%length(terrain[1])+1] == '#'
            trees += 1
        end
    end
    return trees
end

function count_second_test_noprint_4(terrain)
    mult = 1
    for (right, down) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]
        mult *= count_trees_4(terrain, right, down)
    end
    return mult
end

# method 4 was way worse:
#= ```julia
julia> @btime count_second_test_noprint_4(terrain)
  51.145 μs (2 allocations: 176 bytes)
1744787392
```
 =#

# I tried with list comprehension but this was not good, probably
# because it had to store too much info and only then compute the result:
# Still, it is interesting to compare the following two codes,
# with the first one needing many more allocations and memory

function count_second_test_listcomprehension_1(terrain)
    mult = 1
    for (right, down) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]
        count = length([c for c in [terrain[l][(div(right*(l-1),down))%length(terrain[1])+1] for l in 1:down:length(terrain)] if c == '#'])
    mult *= count
    end
    return mult
end

function count_second_test_listcomprehension_4(terrain)
    mult = 1
    for (right, down) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]
    mult *= count(c -> c == '#', [terrain[l][(div(right*(l-1),down))%length(terrain[1])+1] for l in 1:down:length(terrain)])
    end
    return mult
end
#=
```julia
julia> @btime count_second_test_listcomprehension_1(terrain)
  65.648 μs (45 allocations: 13.95 KiB)
1744787392

julia> @btime count_second_test_listcomprehension_4(terrain)
  58.912 μs (7 allocations: 6.70 KiB)
1744787392
```
=#

#

function count_trees_5(terrain, right, down)
    trees = 0
    len = length(terrain[1])
    for (l,c) in zip(1:down:length(terrain), 0:right:div(length(terrain)*right, down)-1)
        if terrain[l][c%len+1] == '#'
            trees += 1
        end
    end
    return trees
end

function count_second_test_noprint_5(terrain)
    mult = 1
    for (right, down) in ((1,1), (3, 1), (5, 1), (7, 1), (1, 2))
        mult *= count_trees_5(terrain, right, down)
    end
    return mult
end
# Not bad, but 2 and 2b are still the best
#=
```julia
julia> @btime count_second_test_noprint_5(terrain)
  7.964 μs (1 allocation: 16 bytes)
1744787392
```
 =#

# Now with a one-liner based on the above
count_second_test_noprint_5b(terrain) = prod( sum(1 for (l,c) in zip(1:down:length(terrain), 0:right:div(length(terrain)*right, down)-1) if terrain[l][c%length(terrain[1])+1] == '#') for (right, down) in ((1,1), (3, 1), (5, 1), (7, 1), (1, 2)))
#=
```julia
julia> @btime prod( sum(1 for (l,c) in zip(1:down:length(terrain), 0:right:div(length(terrain)*right, down)-1) if terrain[l][c%length(terrain[1])+1] == '#') for (right, down) in ((1,1), (3, 1), (5, 1), (7, 1), (1, 2)))
  276.830 μs (677 allocations: 11.83 KiB)
1744787392
``` =#

# Using iterator is twice faster but computing cls and mod is not that bad:
#=
```julia
julia> @btime [r * d for (r,d) in [(1,1), (3, 1), (5, 1), (7, 1), (1, 2)]];
  99.231 ns (2 allocations: 288 bytes)

julia> @btime [r * d for (r,d) in ((1,1), (3, 1), (5, 1), (7, 1), (1, 2))];
  52.459 ns (1 allocation: 128 bytes)

julia> @btime [mod(2i-1,8) * cld(i,4) for i in 1:5];
  62.881 ns (1 allocation: 128 bytes)
```
=#


text:
cases = ["B-splines","Cheb","Lin"]

# construct basis
holder = (BasisMatrices.Basis(BasisMatrices.SplineParams(15,-1,1,1),BasisMatrices.SplineParams(20,-5,2,3)),
    BasisMatrices.Basis(BasisMatrices.ChebParams(15,-1,1),BasisMatrices.ChebParams(20,-5,2)),
    BasisMatrices.Basis(BasisMatrices.LinParams(15,-1,1),BasisMatrices.LinParams(20,-5,2)))

@testset for (i, case) in enumerate(cases)
    basis = holder[i]

    # get nodes
    X, x12 = BasisMatrices.nodes(basis)

    # function to interpolate
    f(x1, x2) = cos(x1) ./ exp(x2)
    f(X::Matrix) = f(X[:, 1], X[:, 2])

    # function at nodes
    y = f(X)

    # benchmark coefficients
    c_direct, bs_direct = BasisMatrices.funfitxy(basis, X, y)

    @testset "test funfitxy for tensor and direct agree on coefs" begin
        c_tensor, bs_tensor = BasisMatrices.funfitxy(basis, x12, y)
        @test maxabs(c_tensor -  c_direct) <=  1e-12
    end

    @testset "test funfitf" begin
        c = BasisMatrices.funfitf(basis,f)
        @test maxabs(c -  c_direct) <=  1e-12
    end

    @testset "test funeval methods" begin
        # single point
        sp = BasisMatrices.funeval(c_direct,basis,X[5:5,:])[1]
        @test maxabs(sp -  y[5]) <= 1e-12

        # multiple points using tensor directly
        mp = BasisMatrices.funeval(c_direct,basis,x12)
        @test maxabs(mp -  y) <=  1e-12

        # multiple points using direct
        mp = BasisMatrices.funeval(c_direct,basis,X)
        @test maxabs(mp -  y) <=  1e-12

        # multiple points giving basis in direct form
        mpd = BasisMatrices.funeval(c_direct,bs_direct)
        @test maxabs(mpd -  y) <=  1e-12

        # multiple points giving basis in expanded form
        Phiexp = Base.convert(BasisMatrices.Expanded,bs_direct)
        mpe = BasisMatrices.funeval(c_direct,Phiexp)
        @test maxabs(mpe -  y) <=  1e-12

    end

    @testset "test interpoland methods" begin
        # (Basis,BasisMatrix,..)
        intp1 = BasisMatrices.Interpoland(basis,bs_direct,y)
        @test maxabs(BasisMatrices.evaluate(intp1,X) - y) <= 1e-12
        # (Basis,Array,..)
        intp2 = BasisMatrices.Interpoland(basis,X,y)
        @test maxabs(BasisMatrices.evaluate(intp2,X) - y) <= 1e-12
        # (BasisParams,Function)
        intp3 = BasisMatrices.Interpoland(basis,f)
        @test maxabs(BasisMatrices.evaluate(intp3,X) - y) <= 1e-12
    end

    @testset "Printing" begin
        iob = IOBuffer()
        show(iob, BasisMatrices.Interpoland(basis, bs_direct, y))
    end

    # TODO: call show on an interpoland instance to get coverage for writemime

end # testset

ext: jl


text:
using Yields
using Test

@testset "Yields.jl" begin

    @testset "rate types" begin
        rs = Rate.([0.1,.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(.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 "constant curve and rate -> Constant" begin
        yield = Yields.Constant(0.05)
        rate = Yields.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.,8.,9.5,10.5,11.0,11.25,11.38,11.44,11.48,11.5] ./ 100
        spot = [6.,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.,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) ≈ fwd[mat] 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) ≈ 6.6 / 100 atol = 0.001
        @test forward(curve, 1.0, 1.5) ≈ 7.6 / 100 atol = 0.001
        @test forward(curve, 1.5, 2.0) ≈ 8.0 / 100 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 "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.)
        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 + .051813 / 2)^2 - 1) * 100
        cmt  = [4.0816,adj,5.4986,5.8620] ./ 100
        mats =  [.5,1.,1.5,2.]
        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

    #     # 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.ZeroCouponQuotes(prices, maturities)
        @test zcq.prices == prices
        @test zcq.maturities == maturities

        @test_throws DomainError Yields.ZeroCouponQuotes([1.3, 2.4, 0.9], maturities)

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

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

        rates = [0.4, -0.7]
        bbq = Yields.BulletBondQuotes(rates, maturities, prices, 3)
        @test bbq.interests == rates
        @test bbq.maturities == maturities
        @test bbq.prices == prices
        @test bbq.frequency == 3

        @test_throws DomainError Yields.BulletBondQuotes([1.3, 2.4, 0.9], maturities, prices, 3)
        @test_throws DomainError Yields.BulletBondQuotes(rates, [4.3, 5.6, 4.4, 4.4], prices, 3)
        @test_throws DomainError Yields.BulletBondQuotes(rates, maturities, [5.7], 3)
        @test_throws DomainError Yields.BulletBondQuotes(rates, maturities, prices, 0)
        @test_throws DomainError Yields.BulletBondQuotes(rates, maturities, prices, -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(), Rate(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(), Rate(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 in 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 in 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 in 1:length(swq_interests)
            @test sum(discount.(sw_bbq, swq_times) .* swq_payments[:, bondIdx]) ≈ bbq_prices[bondIdx]
        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:
#include(Pkg.dir("SDE","src", "Schauder.jl"))
#using Schauder
using SDE.Schauder
using Base.Test

th1 = sin(linspace(0,10,2^5-1))
th2 = copy(th1)
pickup!(th2)
drop!(th2)
@test norm(th1 - th2) < 2.0*eps()

@test [level(ones(i)) for i in 1:15] ==  [1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4]

@test all([vectoroflevels(5,0),  5, vectoroflevels(5,0)] .== vectoroflevels(6,0))

L = 4
@test [number(L,0), 2^L-1+(1:5)] == number(L, 5)

@test unique(diff(vectoroflevels(L,0)[sortperm(number(L,0))])) == [0,1]

@test_approx_eq norm(fe_transf(sin, 0., 2pi, 3).^2 + fe_transf(sin, pi/2, 2.5pi, 3).^2 -1.0) 0

@test norm(fe_transfB1(sin, 0, 2pi, L) - [fe_transf(sin, 0, 2pi, L), 0]) < 4eps()

@test norm(fe_transfB2(sin, 0, 2pi, L) - [fe_transf(sin, 0, 2pi, L),0,0]) < 4eps()

@test dropB1([0,2,0,1]) == [1,3,1,1]

@test dropB2([0.,2.,1.,1.,1.]) == [1., 1., 3., 2., 1.]


text:
function dict end
macro dict(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    Base.getindex(p::$ProxyType, args...)	= Base.getindex(ProxyInterfaces.dict(p), args...)
    Base.length(p::$ProxyType) = Base.length(ProxyInterfaces.dict(p))
    Base.in(a, p::$ProxyType) = Base.in(a, ProxyInterfaces.dict(p))
    Base.haskey(p::$ProxyType, key) = Base.haskey(ProxyInterfaces.dict(p), key)
    Base.keys(p::$ProxyType) = Base.keys(ProxyInterfaces.dict(p))
    Base.values(p::$ProxyType) = Base.values(ProxyInterfaces.dict(p))
    Base.merge(p1::$ProxyType, p2::$ProxyType) = $ProxyType(Base.merge(
      ProxyInterfaces.dict(p1), ProxyInterfaces.dict(p2)))
  end
end

macro dict_mutable(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    @dict($ProxyType)
    Base.setindex!(p::$ProxyType, args...) = Base.setindex!(ProxyInterfaces.dict(p), args...)
    Base.empty!(p::$ProxyType) = Base.empty!(ProxyInterfaces.dict(p))
  end
end


text:
using Documenter, ExercismGenerators

makedocs(;
    modules=[ExercismGenerators],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/SaschaMann/ExercismGenerators.jl/blob/{commit}{path}#L{line}",
    sitename="ExercismGenerators.jl",
    authors="Sascha Mann",
    assets=String[],
)


text:
module VersusMatch
using BenchmarkTools
using Statistics
using Gadfly
using MLStyle
using DataFrames
import Match
using ..ArbitrarySampler
using ..Utils

sym = @spec ::Symbol
scall = @spec :($$sym($$sym, $$sym{0,5}...)) || :($$sym())
atom = @spec $sym || $scall
xcall = @spec :($$atom($$atom, $$atom{0,5}...)) || :($$atom()) || $atom
specs = [
    :s1 => @spec(:(function $$xcall($$xcall, $$xcall)
        $$xcall + $ $xcall
    end)),
    :s2 => @spec(:(function $$xcall($$xcall, $($xcall{0,5}...))
        $$xcall
    end)),
    :s3 => @spec(Expr(:call, Expr(:., $xcall, QuoteNode($sym)), $xcall{0,3}...)),
    :s4 => @spec(:(struct $$sym <: $$xcall
        $$sym::$ $sym
        $$sym::$ $xcall
    end)),
    :s5 => @spec(:(const $$sym = $$scall)),
    :s6 => @spec(:($$xcall = $$scall + $ $sym)),
]

implementations = [
    Symbol(:MLStyle, " Expr-pattern") => let
        function extract_name(e)
            MLStyle.@match e begin
                ::Symbol => e
                Expr(:<:, a, _) => extract_name(a)
                Expr(:struct, _, name, _) => extract_name(name)
                Expr(:call, f, _...) => extract_name(f)
                Expr(:., subject, attr, _...) => extract_name(subject)
                Expr(:function, sig, _...) => extract_name(sig)
                Expr(:const, assn, _...) => extract_name(assn)
                Expr(:(=), fn, body, _...) => extract_name(fn)
                Expr(expr_type, _...) => error(
                    "Can't extract name from ",
                    expr_type,
                    " expression:\n",
                    "    $e\n",
                )
            end
        end
        extract_name
    end,
    Symbol(:MLStyle, " AST-pattern") => let
        function extract_name_homoiconic(e)
            MLStyle.@match e begin
                ::Symbol => e
                :($a <: $_) => extract_name_homoiconic(a)
                :(struct $name <: $_
                    $(_...)
                end) => name
                :($f($(_...))) => extract_name_homoiconic(f)
                :($subject.$_) => extract_name_homoiconic(subject)
                :(function $name($(_...))
                    $(_...)
                end) => extract_name_homoiconic(name)
                :(const $assn = $_) => extract_name_homoiconic(assn)
                :($fn = $_) => extract_name_homoiconic(fn)
                Expr(expr_type, _...) => error(
                    "Can't extract name from ",
                    expr_type,
                    " expression:\n",
                    "    $e\n",
                )
            end
        end
        extract_name_homoiconic
    end,
    Symbol("Match.jl") => let
        extract_name(e::Symbol) = e
        function extract_name(e::Expr)
            Match.@match e begin
                Expr(:<:, [a, b]) => extract_name(a)
                Expr(:struct, [_, name, _]) => extract_name(name)
                Expr(:call, [f, _...]) => extract_name(f)
                Expr(:., [subject, attr, _...]) => extract_name(subject)
                Expr(:function, [sig, _...]) => extract_name(sig)
                Expr(:const, [assn, _...]) => extract_name(assn)
                Expr(:(=), [fn, body, _...]) => extract_name(fn)
                Expr(expr_type, _...) => error(
                    "Can't extract name from ",
                    expr_type,
                    " expression:\n",
                    "    $e\n",
                )
            end
        end
        extract_name
    end,
]

records = NamedTuple{(:time_mean, :implementation, :case)}[]
for (spec_id, spec) in specs
    for (impl_id, impl_fn) in implementations
        bench′ =
            @benchmark $impl_fn(sample) setup = (sample = $generate($spec)) samples = 2000
        time′ = mean(bench′.times)
        @info :bench (spec_id, impl_id, time′)
        push!(records, (time_mean = time′, implementation = impl_id, case = spec_id))
    end
end

df = DataFrame(records)
@info df

report_meantime, df_time =
    report(df, Guide.title("Example from Match.jl Documentation"); benchfield = :time_mean)

open("stats/bench-versus-match.txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/bench-versus-match.svg", 14inch, 6inch), report_meantime)

end


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

include("spatialvolume.jl")
include("parambounds.jl")
include("named_tuple_dist.jl")
include("parameter_mapping.jl")
include("density_sample.jl")


text:
module MathTools

using StaticArrays
using Base.Iterators
using Base: product
using Interpolations

struct Point2D{T <: Real} <: FieldVector{2, T}
    x :: T
    y :: T
end

x(p :: Point2D) = p.x
y(p :: Point2D) = p.y

struct Point3D{T <: Real} <: FieldVector{3, T}
    x :: T
    y :: T
    z :: T
end

x(p :: Point3D) = p.x
y(p :: Point3D) = p.y
z(p :: Point3D) = p.z

# Copying NumPy's linspace function
linspace(min, max, step) = min:(max - min)/step:max
columns(M) = tuple((view(M, :, i) for i in 1:size(M, 2))...)

## Haskell Master Race
#===========================================================================#

span(pred, iter) = (takewhile(pred, iter), dropwhile(pred, iter))
splitat(n, xs) = @views (xs[1:n,:], xs[n+1:end,:])

lisa(pred, iter) = span(!pred, iter)

Base.Iterators.partition(pred, xs, f, g) = f.(filter(pred, xs)), g.(filter(!pred, xs))
Base.Iterators.partition(pred, xs, f = identity) = partition(pred, xs, f, f)

## Renaming math operations
#===========================================================================#

"""
"Lenses" to access subfields on lists of objects.
"""
|>(obj, fields :: Array{Symbol}) = foldl(getproperty, fields, init = obj)
|>(list_objs :: Array{T}, fields :: Array{Symbol}) where T <: Any = list_objs .|> [fields]

field << obj = getfield(obj, field)

# Convert homogeneous struct entries to lists
structtolist(x) = [ name << x for name ∈ (fieldnames ∘ typeof)(x) ]

## Renaming math operations
#===========================================================================#

⊗(A, B)    = kron(A, B)

×(xs, ys)   = product(xs, ys)
dot(V₁, V₂) = sum(V₁ .* V₂)
# ×(xs, ys) = (collect ∘ zip)(xs' ⊗ (ones ∘ length)(ys), (ones ∘ length)(xs)' ⊗ ys)

## Basic transformations/geometry
#===========================================================================#

# Transforms (x, y) to the coordinate system with (x_s, y_s) as origin oriented at α_s.
affine_2D(x, y, x_s, y_s, α_s) = rotation(x - x_s, y - y_s, α_s)
inverse_rotation(x, y, angle) = SVector(x * cos(angle) - y * sin(angle), x * sin(angle) + y * cos(angle))
rotation(x, y, angle) = SVector(x * cos(angle) + y * sin(angle), -x * sin(angle) + y * cos(angle))

# Matrix versions
rotation(θ)         = [ cos(θ) sin(θ) ;
                       -sin(θ) cos(θ) ]
inverse_rotation(θ) = rotation(-θ)

# Cartesian-polar coordinates
magnitude(U, W) = sqrt(U^2 + W^2)
angle(U, W) = atan(W, U)

slope(x1, y1, x2, y2) = (y2 - y1) / (x2 - x1)

## Array conversions
#===========================================================================#

tuparray(xs)  = tuple.(eachcol(xs)...)
vectarray(xs) = SVector.(eachcol(xs)...)

extend_yz(coords) = @views [ coords[:,1] zero(coords[:,1]) coords[:,2] ]

reflect_mapper(f, xs) = @views [ f(xs[:,end:-1:1]) xs ]

## Difference operations
#===========================================================================#

forward_difference_matrix(n) = [ I zeros(n) ] - [ zeros(n) I ]

forward_sum(xs) = @views @. xs[2:end] + xs[1:end-1]
forward_difference(xs) = @views @. xs[2:end] - xs[1:end-1]
forward_division(xs) = @views @. xs[2:end] / xs[1:end-1]
ord2diff(xs) = @views @. xs[3:end] - 2 * xs[2:end-1] + xs[1:end-2]

adj3(xs) = @views zip(xs[1:end-2,:], xs[2:end-1,:], xs[3:end,:])

# Central differencing schema for pairs except at endpoints
midpair_map(f :: H, xs) where {H} =
    @views  [        f.(xs[1,:], xs[2,:])'        ;
               f.(xs[1:end-2,:], xs[3:end,:]) / 2 ;
                 f.(xs[end-1,:], xs[end,:])'      ]

# stencil(xs, n) = [ xs[n+1:end] xs[1:length(xs) - n] ]
# parts(xs) = let adj = stencil(xs, 1); adj[1,:], adj[end,:] end

# Lazy? versions
# stencil(xs, n) = zip(xs[n+1:end], xs[1:length(xs) - n])
# parts(xs) = (first ∘ stencil)(xs, 1), (last ∘ stencil)(xs, 1)

# function midgrad(xs)
#     first_two_pairs, last_two_pairs = permutedims.(parts(xs))
#     central_diff_pairs = stencil(xs, 2)

#     [first_two_pairs; central_diff_pairs; last_two_pairs]
# end

## Spacing formulas
#===========================================================================#

uniform_spacing(x1, x2, n) = range(x1, x2, length = n)
linear_spacing(x_center, len, n :: Integer) = @. x_center + len * 0:1/(n-1):1
cosine_spacing(x_center, diameter, n :: Integer = 40) = @. x_center + (diameter / 2) * cos(-π:π/(n-1):0)

function sine_spacing(x1, x2, n :: Integer = 40)
    d = x2 - x1
    if n < 0
       @. x2 - d * sin(π/2 * (1. - (1:1/(n+1):0)))[end:-1:1]
    else
       @. x1 + d * sin(π/2 * (0:1/(n-1):1))
    end
end

function cosine_interp(xs, n :: Integer = 40)
    d = maximum(xs) - minimum(xs)
    x_center = (maximum(xs) + minimum(xs)) / 2
    x_circ = cosine_spacing(x_center, d, n)
end

## Iterator methods
#===========================================================================#

# Need to improve this via recursion; the for loop seems really unnecessary
function accumap(f, n, xs)
    data = [ xs ]
    for i = 1:n
        ys = map(f, xs)
        data = [ reduce(vcat, data); ys ]
        xs = ys
    end
    return reduce(hcat, data)
end

## Helper functions
#===========================================================================#

"""
    weighted(x1, x2, μ)

Compute the weighted value between two values ``x_1`` and ``x_2`` with weight ``\\mu \\in [0,1]``.
"""
weighted(x1, x2, μ) = (1 - μ) * x1 + μ * x2

"""
Compute the weighted average (μ) of two vectors.
"""
weighted_vector(x1, x2, μ) = weighted.(x1, x2, μ)

## Macros
#===========================================================================#

# macro getter(obj)
#     return :((function name(x)
#                 getfield(name, x)
#             end))
# end

reshape_array(arr, inds, sizes) = @views [ reshape(arr[i1+1:i2], size) for (i1, i2, size) in zip(inds[1:end-1], inds[2:end], sizes) ]

end


text:
############################################################################
# Find the ground state of the Ising Hamitlonian and get the magnetization #
############################################################################
using MatrixProductStates

let
    # Number of Spins
    N = 20
    # Bond Dimension
    D = 20
    # Physical dimension
    d = 2
    # Desired Accuracy
    acc = 1E-10
    # The number of time steps we take
    nsteps = 200
    # The size of each time step
    dt = 5E-2

    # Coupling between nearest neigbours
    J = 1.0
    # Coupling to external field
    lambda = 1.0
    # Generate Ising Hamiltonian H = -J * sum_{i=1}^{N-1} X^i X^i+1 - lambda * sum_{i=1}^N Z^i.
    H = getIsingMPO(N, J, lambda)
    # Get the MPO for the total spin
    Sz = getTotalSpinMPO(N)

    #=
    In the following we prepare the time evolution operators in MPO form using an
    odd-even decomposition for the evolution operator. Since the Ising Hamiltonian is
    of the form H = \sum_i h_{i,i+1} we can approximate the matrix exponential using
    a Suzuki-Trotter decomposition as follows
    exp(-dt * H) ≈ \prod_{i even} exp(-dt * h_{i,i+1}) \prod_{i odd} exp(-dt * h_{i,i+1}).
    Since all terms starting at an odd (even) site commute, we can find an MPO
    representation by decomping one of the local terms exp(-dt * h_{i,i+1}) into an MPO,
    and using the same tensors at every site. Special care has to be taken on the
    boundaries.
    =#

    # The required Pauli matrices for building the evolution MPOs
    Id, X, _, Z = getPauliMatrices()
    # Prepare the matrix representation for the different exponentials appearing in the odd-even decomposition
    Uloc_left_boundary = exp(-dt * (-J * kron(X, X) - lambda * kron(Z, Id) - 0.5 * lambda * kron(Id, Z)))
    Uloc_center = exp(-dt * (-J * kron(X, X) - 0.5 * lambda * (kron(Z, Id) + kron(Id, Z))))
    Uloc_right_boundary = exp(-dt * (-J * kron(X, X) - 0.5 * lambda * kron(Z, Id) - lambda * kron(Id, Z)))
    # The identity operator as MPO tensor
    id_tensor = reshape(Id, (1, 1, 2, 2))
    # Find the mpo representation for the local terms
    mpo_Uloc_left_boundary = decompose_into_mpo(Uloc_left_boundary, d)
    mpo_Uloc_center = decompose_into_mpo(Uloc_center, d)
    mpo_Uloc_right_boundary = decompose_into_mpo(Uloc_right_boundary, d)
    # Build the MPO for evolving the odd sites
    Uodd = MPO{Float64}(undef, N)
    for i = 1:2:N-1
        if i == 1
            Uodd[i] = mpo_Uloc_left_boundary[1]
            Uodd[i+1] = mpo_Uloc_left_boundary[2]
        elseif i == N - 1
            Uodd[i] = mpo_Uloc_right_boundary[1]
            Uodd[i+1] = mpo_Uloc_right_boundary[2]
        else
            Uodd[i] = mpo_Uloc_center[1]
            Uodd[i+1] = mpo_Uloc_center[2]
        end
    end
    if isodd(N)
        Uodd[N] = id_tensor
    end
    # Build the MPO for evolving the even sites
    Ueven = MPO{Float64}(undef, N)
    Ueven[1] = id_tensor
    for i = 2:2:N-1
        if i == N - 1
            Ueven[i] = mpo_Uloc_right_boundary[1]
            Ueven[i+1] = mpo_Uloc_right_boundary[2]
        else
            Ueven[i] = mpo_Uloc_center[1]
            Ueven[i+1] = mpo_Uloc_center[2]
        end
    end
    if iseven(N)
        Ueven[N] = id_tensor
    end

    # Evolve in imaginary time starting from a random, normalized state
    psi = random_mps_obc(N, D, d, Float64)
    gaugeMPS!(psi, :left, true)
    for i = 1:nsteps
        print("\rStep: ", i)
        # Apply the evolution operators
        psi = apply_operator(Uodd, psi)
        psi = apply_operator(Ueven, psi)
        # Compress and renormalize
        psi = svd_compress_mps(psi, D)
        gaugeMPS!(psi, :left, true)
    end
    E0 = expectation_value(psi, H)

    # Print the results
    println(" ")
    println("Groundstate energy:            ", E0)
    println("Groundstate energy density:    ", E0 / N)
    println("Total spin:                    ", expectation_value(psi, Sz))
    println("Entropy along center bond      ", compute_entropy(psi))

    nothing
end


text:
# using Unitful: K, °C, m, mm, s, kg
# using Unitful: @u_str, ustrip, uconvert

@testset "Interface" begin

# test that load return a ClimGrid type
file1 = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")
file2 = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")
files = [file1, file2]
C = load(files, "tas")
filenc = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")
C = load(filenc, "tas")
@test load(filenc, "tas", data_units = "Celsius")[2] == "°C"
@test load(filenc, "pr", data_units = "mm")[2] == "mm"

# With a shapefile
filename = joinpath(dirname(@__FILE__), "data", "SudQC_GCM.shp")
polyshp = read(filename,Shapefile.Handle)
P = shapefile_coords_poly(polyshp.shapes[1])
C = load(filenc, "tas", poly=P)
@test size(C[1]) == (23, 12, 1)
@test C[1][1,1,1] == 294.6609f0

# Dataset functions
ds = Dataset(filenc)
@test ClimateTools.longridname(ds) == "lon"
@test ClimateTools.latgridname(ds) == "lat"
@test getdim_lon(ds) == ("lon", false)
@test getdim_lat(ds) == ("lat", false)

# Meridian polygon
P = [NaN 10 -10 -10 10 10;NaN 20 20 30 30 20]
C = load(filenc, "tas", poly=P)
@test size(C[1]) == (256, 7, 1)
@test findmax(C, skipnan=true) == (303.45667f0, CartesianIndex(124, 4, 1))

fileorog = joinpath(dirname(@__FILE__), "data", "orog_fx_GFDL-ESM2G_historicalMisc_r0i0p0.nc")
orog = load2D(fileorog, "orog")
@test size(orog[1]) == (144, 90)
@test orog.typeofvar == "orog"

filename = joinpath(dirname(@__FILE__), "data", "SudQC_GCM.shp")
polyshp = read(filename,Shapefile.Handle)
P = shapefile_coords_poly(polyshp.shapes[1])
orog = load2D(fileorog, "orog", poly = P)
@test size(orog[1]) == (13, 8)
@test orog.typeofvar == "orog"

P = [NaN 10 -10 -10 10 10;NaN 20 20 30 30 20]
orog = load2D(fileorog, "orog", poly=P)
@test size(orog[1]) == (144, 5)
@test findmax(orog, skipnan=true) == (1116.9403076171875, CartesianIndex(76, 3))

# INTERFACE
C = load(filenc, "tas")
P = shapefile_coords_poly(polyshp.shapes[1])
# B = vcat(C, C)
# @test size(B.data) == (256, 128, 2) # vcat does not look at dimensions
B = merge(C, C)
@test size(B.data) == (256, 128, 1) # C being similar, they should not add up, as opposed to vcat
# Operators +, -, *, /
B = C + C; @test B[1].data[1, 1, 1] == 438.4457f0
B = C * C; @test B[1].data[1, 1, 1] == 48058.66f0
B = C / C; @test B[1].data[1, 1, 1] == 1.0f0
B = C - C; @test B[1].data[1, 1, 1] == 0.0f0
B = C - 1.0; @test B[1].data[1, 1, 1] == 218.2228546142578
B = C - 1; @test B[1].data[1, 1, 1] == 218.22285f0
B = C / 2; @test B[1].data[1, 1, 1] == 109.61143f0
B = C / 2; @test B[1].data[1, 1, 1] == 109.61143f0
B = C / 2.2; @test B[1].data[1, 1, 1] == 99.6467520973899
B = C * 2; @test B[1].data[1, 1, 1] == 438.4457f0
B = C * 2; @test B[1].data[1, 1, 1] == 438.4457f0
B = C * 2.2; @test B[1].data[1, 1, 1] == 482.2902801513672

@test ClimateTools.mean(C) == 278.6421f0
@test ClimateTools.maximum(C) == 309.09613f0
@test ClimateTools.minimum(C) == 205.24321f0
@test ClimateTools.std(C) == 21.92836f0
@test round(ClimateTools.var(C), digits=3) == 480.853f0

# @test typeof(show(C)) == Dict{Any, Any}
# @test typeof(C[1].data) == Array{Unitful.Quantity{Float32,Unitful.Dimensions{(Unitful.Dimension{:Temperature}(1//1),)},Unitful.FreeUnits{(Unitful.Unit{:Kelvin,Unitful.Dimensions{(Unitful.Dimension{:Temperature}(1//1),)}}(0, 1//1),),Unitful.Dimensions{(Unitful.Dimension{:Temperature}(1//1),)}}},3}

@test C[2] == "K"
@test C[3] == "N/A"
@test C[4] == "720 ppm stabilization experiment (SRESA1B)"
@test C[5] == "N/A"
@test C[6] == "degrees_east"
@test C[7] == "degrees_north"
@test C[8] == joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")
@test C[9] == "tas"
@test annualmax(C)[9] == "annualmax"
@test C[10] == "tas"
@test C[11] == "noleap"
@test typeof(C[12]) == Dict{String, Any}
@test C[12]["project_id"] == "IPCC Fourth Assessment"
@test_throws ErrorException C[13]
@test annualmax(C)[10] == "tas"
@test size(C) == (22, )
@test size(C, 1) == 22
@test length(C) == 22
# @test endof(C) == 21
@test_throws MethodError C[end]
@test ndims(C) == 1

# Spatial subset
C = load(filenc, "tas")
Csub = spatialsubset(C, P)
@test size(Csub[1]) == (23, 12, 1)
@test Csub[1][1, 1, 1] == 294.6609f0
Csub = spatialsubset(C, P')
@test size(Csub[1]) == (23, 12, 1)
@test Csub[1][1, 1, 1] == 294.6609f0
C = load(filenc, "ua")
Csub = spatialsubset(C, P)
@test size(Csub[1]) == (23, 12, 17, 1)
@test Csub[1][12, 1, 1, 1] == 6.658482f0
@test isnan(Csub[1][1, 1, 1, 1])

poly= [[NaN 10 -10 -10 10 10];[NaN -10 -20 10 10 -10]] # meridian test
C = load(filenc, "tas", poly=poly)

# Spatial subset
C = load(filenc, "tas")
Csub = temporalsubset(C, (2000, 05, 16, 12, 0, 0), (2000, 05, 16, 12, 0, 0))
@test Csub[1][1, 1, 1] == 219.22285f0
@test Csub[1][Axis{:time}][1] == DateTimeNoLeap(2000, 05, 16, 12)
B = load(filenc, "tas", start_date=(2000, 05, 16, 12), end_date=(2000, 05, 16, 12))
@test B[1] == C[1]

# Time resolution
timevec = [1, 2, 3]
@test ClimateTools.timeresolution(timevec) == "24h"
@test ClimateTools.daymean_factor(ClimateTools.timeresolution(timevec)) == 1

timevec = [1.0, 1.5, 2.0]
@test ClimateTools.timeresolution(timevec) == "12h"
@test ClimateTools.daymean_factor(ClimateTools.timeresolution(timevec)) == 2

timevec = [1.25, 1.5, 1.75]
@test ClimateTools.timeresolution(timevec) == "6h"
@test ClimateTools.daymean_factor(ClimateTools.timeresolution(timevec)) == 4

timevec = [1.125, 1.25, 1.375]
@test ClimateTools.timeresolution(timevec) == "3h"
@test ClimateTools.daymean_factor(ClimateTools.timeresolution(timevec)) == 8

timevec = [1.0, 1.04167, 1.08334]
@test ClimateTools.timeresolution(timevec) == "1h"
@test ClimateTools.daymean_factor(ClimateTools.timeresolution(timevec)) == 24

timevec = [0.0, 365.0, 365.0*2]
@test ClimateTools.timeresolution(timevec) == "Yearly"

timevec = NetCDF.ncread(filenc, "time")
@test ClimateTools.timeresolution(timevec) == "N/A"
@test ClimateTools.daymean_factor(ClimateTools.timeresolution(timevec)) == 1

# MESHGRID
YV = [1, 2, 3]
XV = [4, 5, 6]
@test meshgrid(XV, YV) == ([4 5 6; 4 5 6; 4 5 6], [1 1 1; 2 2 2; 3 3 3])
@test meshgrid(XV) == ([4 5 6; 4 5 6; 4 5 6], [4 4 4; 5 5 5; 6 6 6])
Q = Array{Int64, 3}(undef, 3, 3, 3)
R = Array{Int64, 3}(undef, 3, 3, 3)
S = Array{Int64, 3}(undef, 3, 3, 3)
Q[:, :, 1] .= [4 5 6; 4 5 6; 4 5 6]
Q[:, :, 2] .= [4 5 6; 4 5 6; 4 5 6]
Q[:, :, 3] .= [4 5 6; 4 5 6; 4 5 6]
R[:, :, 1] .= [1 1 1; 2 2 2; 3 3 3]
R[:, :, 2] .= [1 1 1; 2 2 2; 3 3 3]
R[:, :, 3] .= [1 1 1; 2 2 2; 3 3 3]
S[:, :, 1] .= [4 4 4; 4 4 4; 4 4 4]
S[:, :, 2] .= [5 5 5; 5 5 5; 5 5 5]
S[:, :, 3] .= [6 6 6; 6 6 6; 6 6 6]

@test meshgrid(XV, YV, XV) == (Q, R, S)

@test ndgrid(XV, YV) == ([4 4 4; 5 5 5; 6 6 6], [1 2 3; 1 2 3; 1 2 3])
@test ndgrid([XV, YV, XV]) ==   [[4, 5, 6], [1, 2, 3], [4, 5, 6]]
@test ndgrid(XV) == [4, 5, 6]

# isdefined
C = 1
@test ClimateTools.@isdefined C
@test (ClimateTools.@isdefined T) == false

## INPOLY
@test ClimateTools.leftorright(0.5,0.5, 1,0,1,1) == -1
@test ClimateTools.leftorright(1.5,.5, 1,0,1,1) == 1
@test ClimateTools.leftorright(1,0.5, 1,0,1,1) == 0

poly = Float64[0 0
               0 1
               1 1
               1 0
               0 0]'
p1 = [0.5, 0.5]
p2 = [0.5, 0.99]
p22 = [0.5, 1] # on edge
p23 = [0.5, 0] # on edge
p24 = [0, 0]   # on corner
p25 = [0, .4]   # on edge
p3 = [0.5, 1.1]

@test inpoly(p1, poly)
@test inpoly(p2, poly)
@test inpoly(p22, poly)
@test inpoly(p23, poly)
@test inpoly(p24, poly)
@test inpoly(p25, poly)
@test !inpoly(p3, poly)

# clockwise poly
poly = Float64[0 0
               1 0
               1 1
               0 1
               0 0]'

@test inpoly(p1, poly)
@test inpoly(p2, poly)
@test inpoly(p22, poly)
@test inpoly(p23, poly)
@test inpoly(p24, poly)
@test inpoly(p25, poly)
@test !inpoly(p3, poly)

# cross-over poly
poly = Float64[0 0
               1 0
               0 1
               1 1
               0 0]'
eval(:(@test_broken inpoly(p1, poly) )) # should be true
@test inpoly(p2, poly)
@test inpoly(p22, poly)
@test inpoly(p23, poly)
@test inpoly(p24, poly)
@test !inpoly(p25, poly) # different
@test !inpoly(p3, poly)

# with interior region
poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside interior poly: i.e. labeled as outside
@test !inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               # in-interior
               0.4 0.4
               0.4 0.2
               0.2 0.2
               0.2 0.4
               0.4 0.4
               # interior
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside in-interior poly
@test inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               # in-interior
               0.4 0.4
               0.2 0.4
               0.2 0.2
               0.4 0.2
               0.4 0.4
               # interior
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside in-interior poly
@test inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior #1
               0.1 0.1
               0.1 0.6
               0.4 0.6
               0.4 0.6
               0.4 0.1
               0.1 0.1
               0 0
               # interior #2
               0.6 0.4
               0.6 0.6
               0.8 0.6
               0.8 0.4
               0.6 0.4
               0 0
               # exterior
               0 1
               1 1
               1 0
               0 0]'
@test !inpoly([0.2,0.4], poly)
@test !inpoly([0.3,0.15], poly)
@test inpoly([0.5,0.4], poly)
@test inpoly([0.5,0.2], poly)
@test !inpoly([0.7,0.5], poly)

# Test Interpolation/regridding

C = load(filenc, "tas")
# Get lat lon vector
lat = Float32.(C[1][Axis{:lat}][:])
lon = Float32.(C[1][Axis{:lon}][:])
latgrid = Float32.(C.latgrid)
longrid = Float32.(C.longrid)
# Shift longitude by 1
lon .+= Float32(1.0)
longrid .+= Float32(1.0)
axisdata = AxisArray(C[1].data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(C[1][Axis{:time}][:]))
C2 = ClimGrid(axisdata, variable = "tas", longrid=longrid, latgrid=latgrid, msk=C.msk)
@test regrid(C, C2)[1].data[1, 1, 1] == 219.24705505371094 #219.28975099839346#219.2400638156467
@test regrid(C, C2, min=0.0, max=0.0)[1].data[1, 1, 1] == 0.0
@test regrid(C, lon, lat)[1].data[1, 1, 1] == 219.24705505371094#219.28975099839346 #219.2400638156467
@test regrid(C, lon, lat, min=0.0, max=0.0)[1].data[1, 1, 1] == 0.0

# Test applymask
# 1-D data
data = randn(3)
mask = [NaN; 1.;1.]
@test isnan(applymask(data, mask)[1])
@test applymask(data, mask)[2] == data[2]
@test applymask(data, mask)[3] == data[3]

# 2-D data
data = randn(3, 2)
mask = [[NaN; 1; 1] [1.; NaN;1.]]
@test isnan(applymask(data, mask)[1, 1]) && isnan(applymask(data, mask)[2, 2])
@test applymask(data, mask)[2, 1] == data[2, 1]
@test applymask(data, mask)[1, 2] == data[1, 2]
@test applymask(data, mask)[3, 1] == data[3, 1]
@test applymask(data, mask)[3, 2] == data[3, 2]

# 3-D data
data = randn(3, 2, 3)
mask = [[NaN; 1; 1] [1.; NaN;1.]]
@test isnan(applymask(data, mask)[1, 1, 1]) && isnan(applymask(data, mask)[2,2,1]) && isnan(applymask(data, mask)[1, 1, 2]) && isnan(applymask(data, mask)[2, 2, 2]) && isnan(applymask(data, mask)[1,1,3]) && isnan(applymask(data, mask)[2, 2, 3])

for i = 1:size(data, 1)
    @test applymask(data, mask)[2,1,i] == data[2,1,i]
    @test applymask(data, mask)[1,2,i] == data[1,2,i]
    @test applymask(data, mask)[3,1,i] == data[3,1,i]
    @test applymask(data, mask)[3,2,i] == data[3,2,i]
end

# 4-D data
data = randn(3,2,1,3)
mask = [[NaN; 1; 1] [1.; NaN;1.]]
@test isnan(applymask(data, mask)[1, 1, 1, 1]) && isnan(applymask(data, mask)[2,2,1,1]) && isnan(applymask(data, mask)[1, 1, 1,2]) && isnan(applymask(data, mask)[2, 2, 1,2]) && isnan(applymask(data, mask)[1, 1, 1,3]) && isnan(applymask(data, mask)[2, 2, 1,3])

for i = 1:size(data, 1)
    @test applymask(data, mask)[2,1,1,i] == data[2,1,1,i]
    @test applymask(data, mask)[1,2,1,i] == data[1,2,1,i]
    @test applymask(data, mask)[3,1,1,i] == data[3,1,1,i]
    @test applymask(data, mask)[3,2,1,i] == data[3,2,1,i]
end

# Test timeresolution and pr_timefactor
filename = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")
timevec = NetCDF.ncread(filename, "time")
@test ClimateTools.pr_timefactor(ClimateTools.timeresolution(timevec)) == 1.0
@test ClimateTools.pr_timefactor("24h") == 86400.0
@test ClimateTools.pr_timefactor("12h") == 43200.0
@test ClimateTools.pr_timefactor("6h") == 21600.0
@test ClimateTools.pr_timefactor("3h") == 10800.0
@test ClimateTools.pr_timefactor("1h") == 3600.0

# timeindex tests
d = collect(DateTime(2000, 01, 01):Day(1):DateTime(2005, 12, 31))
T = typeof(d[1])

start_date = (Inf,)
end_date = (Inf,)
idxtimebeg, idxtimeend = ClimateTools.timeindex(d, start_date, end_date, T)
@test idxtimebeg == 1
@test idxtimeend == length(d)

start_date = (2001, 01, 01)
end_date = (Inf,)
idxtimebeg, idxtimeend = ClimateTools.timeindex(d, start_date, end_date, T)
@test idxtimebeg == 367
@test idxtimeend == length(d)

start_date = (2001, 01, 01)
end_date = (2003, 12, 31)
idxtimebeg, idxtimeend = ClimateTools.timeindex(d, start_date, end_date, T)
@test idxtimebeg == 367
@test idxtimeend == 1461

@test_throws ArgumentError ClimateTools.timeindex(d, end_date, start_date, T)

start_date = (2001,)
end_date = (2004,)
idxtimebeg, idxtimeend = ClimateTools.timeindex(d, start_date, end_date, T)
@test idxtimebeg == 367
@test idxtimeend == 1462

start_date = (2001, 01, 01, 0, 0)
end_date = (2003, 12, 31, 0, 0)
idxtimebeg, idxtimeend = ClimateTools.timeindex(d, start_date, end_date, T)
@test idxtimebeg == 367
@test idxtimeend == 1461

end


text:
# Some utility functions for testing Ferrite.jl

#####################################
# Volume for the reference elements #
#####################################
reference_volume(::Interpolation{dim, RefCube}) where {dim} = 2^dim
reference_volume(::Interpolation{dim, RefTetrahedron}) where {dim} = 1 / factorial(dim)
# For faces
reference_volume(fs::Interpolation, ::Int) = reference_volume(Ferrite.getlowerdim(fs))
reference_volume(fs::Interpolation{2, RefTetrahedron}, face::Int) = face == 1 ? sqrt(2) : 1.0
reference_volume(fs::Interpolation{3, RefTetrahedron}, face::Int) = face == 3 ? sqrt(2 * 1.5) / 2.0 : 0.5

######################################################
# Coordinates and normals for the reference elements #
######################################################

# Lagrange{1, RefCube}
function reference_normals(::Lagrange{1, RefCube})
    return [Vec{1, Float64}((-1.0,)),
            Vec{1, Float64}(( 1.0,))]
end

# Lagrange{2, RefCube}
function reference_normals(::Lagrange{2, RefCube})
    return [Vec{2, Float64}(( 0.0, -1.0)),
            Vec{2, Float64}(( 1.0,  0.0)),
            Vec{2, Float64}(( 0.0,  1.0,)),
            Vec{2, Float64}((-1.0,  0.0,))]
end

# Lagrange{2, RefTetrahedron}
function reference_normals(::Lagrange{2, RefTetrahedron})
    return [Vec{2, Float64}((1/√2, 1/√2)),
            Vec{2, Float64}((-1.0, 0.0)),
            Vec{2, Float64}((0.0, -1.0))]
end

# Lagrange{3, RefTetrahedron}
function reference_normals(::Lagrange{3, RefTetrahedron})
    return [Vec{3, Float64}((0.0, 0.0, -1.0)),
            Vec{3, Float64}((0.0 ,-1.0, 0.0)),
            Vec{3, Float64}((1/√3, 1/√3, 1/√3)),
            Vec{3, Float64}((-1.0, 0.0, 0.0))]
end

# Lagrange{3, Cube}
function reference_normals(::Lagrange{3, RefCube, 1})
    return [Vec{3, Float64}(( 0.0,  0.0, -1.0)),
            Vec{3, Float64}(( 0.0, -1.0,  0.0)),
            Vec{3, Float64}(( 1.0,  0.0,  0.0)),
            Vec{3, Float64}(( 0.0,  1.0,  0.0)),
            Vec{3, Float64}((-1.0,  0.0,  0.0)),
            Vec{3, Float64}(( 0.0,  0.0,  1.0))]
end

# Serendipity{2, RefCube}
reference_normals(::Serendipity{2, RefCube, 2}) = reference_normals(Lagrange{2, RefCube,1}())

##################################
# Valid coordinates by expanding #
# and rotating reference shape   #
##################################
function rotmat(dim, θ=π/6)
    if dim == 1
        R = Tensor{2,1}((cos(θ),))
        return R
    elseif dim == 2
        R = Tensor{2,2}((cos(θ), sin(θ), -sin(θ), cos(θ)))
        return R
    else
        u = Vec{3}((1.0, 2.0, 3.0)); u /= norm(u)
        ux = Tensor{2,3}((0.0, u[3], -u[2], -u[3], 0.0, u[1], u[2], -u[1], 0.0))
        R = cos(θ)*one(Tensor{2,3}) + sin(θ)*ux + (1-cos(θ))*u⊗u
        return R
    end
end

function valid_coordinates_and_normals(fs::Interpolation{dim, shape, order}) where {dim, shape, order}
    x = Ferrite.reference_coordinates(fs)
    n = reference_normals(fs)
    R = rotmat(dim)
    return [2.0 * (R ⋅ x[i]) for i in 1:length(x)] , [(R ⋅ n[i]) / norm((R ⋅ n[i])) for i in 1:length(n)]
end

#######################################
# Volume of cells (with planar edges) #
#######################################
function calculate_volume(::Lagrange{1, RefCube, 1}, x::Vector{Vec{dim, T}}) where {T, dim}
    vol = norm(x[2] - x[1])
    return vol
end

function calculate_volume(::Lagrange{1, RefCube, 2}, x::Vector{Vec{dim, T}}) where {T, dim}
    vol = norm(x[3] - x[1]) + norm(x[2]-x[3])
    return vol
end

function calculate_volume(::Lagrange{2, RefCube, 1}, x::Vector{Vec{dim, T}}) where {T, dim}
    vol = norm((x[4] - x[1]) × (x[2] - x[1])) * 0.5 +
          norm((x[4] - x[3]) × (x[2] - x[3])) * 0.5
    return vol
end

function calculate_volume(::Lagrange{2, RefCube, 2}, x::Vector{Vec{dim, T}}) where {T, dim}
    vol = norm((x[8] - x[1]) × (x[5] - x[1])) * 0.5 +
          norm((x[8] - x[9]) × (x[5] - x[9])) * 0.5 +
          norm((x[5] - x[2]) × (x[6] - x[2])) * 0.5 +
          norm((x[5] - x[9]) × (x[6] - x[9])) * 0.5 +
          norm((x[6] - x[3]) × (x[7] - x[3])) * 0.5 +
          norm((x[6] - x[9]) × (x[7] - x[9])) * 0.5 +
          norm((x[7] - x[4]) × (x[8] - x[4])) * 0.5 +
          norm((x[7] - x[9]) × (x[8] - x[9])) * 0.5
    return vol
end

function calculate_volume(::Lagrange{2, RefTetrahedron, 1}, x::Vector{Vec{dim, T}}) where {T, dim}
    vol = norm((x[1] - x[3]) × (x[2] - x[3])) * 0.5
    return vol
end

function calculate_volume(::Lagrange{2, RefTetrahedron, 2}, x::Vector{Vec{dim, T}}) where {T, dim}
    vol = norm((x[6] - x[3]) × (x[5] - x[3])) * 0.5 +
          norm((x[6] - x[4]) × (x[5] - x[4])) * 0.5 +
          norm((x[1] - x[6]) × (x[4] - x[6])) * 0.5 +
          norm((x[4] - x[5]) × (x[2] - x[5])) * 0.5
    return vol
end

function calculate_volume(::Lagrange{3, RefTetrahedron, order}, x::Vector{Vec{3, T}}) where {T, order}
    vol = norm((x[2] - x[1]) ⋅ ((x[3] - x[1]) × (x[4] - x[1]))) / 6.0
    return vol
end

function calculate_volume(::Lagrange{3, RefCube, 1}, x::Vector{Vec{3, T}}) where T
    vol = norm((x[1] - x[5]) ⋅ ((x[2] - x[5]) × (x[4] - x[5]))) / 6.0 +
          norm((x[2] - x[7]) ⋅ ((x[3] - x[7]) × (x[4] - x[7]))) / 6.0 +
          norm((x[2] - x[7]) ⋅ ((x[4] - x[7]) × (x[5] - x[7]))) / 6.0 +
          norm((x[2] - x[7]) ⋅ ((x[5] - x[7]) × (x[6] - x[7]))) / 6.0 +
          norm((x[4] - x[8]) ⋅ ((x[5] - x[8]) × (x[7] - x[8]))) / 6.0
    return vol
end

function calculate_volume(::Serendipity{2, RefCube, 2}, x::Vector{Vec{2, T}}) where T
    vol = norm((x[5] - x[1]) × (x[8] - x[1])) * 0.5 +
          norm((x[6] - x[2]) × (x[5] - x[2])) * 0.5 +
          norm((x[7] - x[3]) × (x[6] - x[3])) * 0.5 +
          norm((x[8] - x[4]) × (x[7] - x[4])) * 0.5 +
          norm((x[6] - x[5]) × (x[8] - x[5])) * 0.5 +
          norm((x[6] - x[7]) × (x[8] - x[7])) * 0.5
    return vol
end

# For faces
function calculate_volume(::Lagrange{0, RefCube, order}, ::Vector{Vec{1, T}}) where {order, T}
    return one(T)
end

getnfaces(::Interpolation{dim, RefCube}) where {dim} = 2*dim
getnfaces(::Interpolation{2, RefTetrahedron}) = 3
getnfaces(::Interpolation{3, RefTetrahedron}) = 4

coords_on_faces(x, ::Lagrange{1, RefCube, 1}) = ([x[1]], [x[2]])
coords_on_faces(x, ::Lagrange{1, RefCube, 2}) = ([x[1]], [x[2]])
coords_on_faces(x, ::Lagrange{2, RefCube, 1}) =
    ([x[1],x[2]], [x[2],x[3]], [x[3],x[4]], [x[1],x[4]])
coords_on_faces(x, ::Lagrange{2, RefCube, 2}) =
    ([x[1],x[2],x[5]], [x[2],x[3],x[6]], [x[3],x[4],x[7]], [x[1],x[4],x[8]])
coords_on_faces(x, ::Lagrange{2, RefTetrahedron, 1}) =
    ([x[1],x[2]], [x[2],x[3]], [x[1],x[3]])
coords_on_faces(x, ::Lagrange{2, RefTetrahedron, 2}) =
    ([x[1],x[2],x[4]], [x[2],x[3],x[5]], [x[1],x[3],x[6]])
coords_on_faces(x, ::Lagrange{3, RefTetrahedron, 1}) =
    ([x[1],x[2],x[3]], [x[1],x[2],x[4]], [x[2],x[3],x[4]], [x[1],x[3],x[4]])
coords_on_faces(x, ::Lagrange{3, RefTetrahedron, 2}) =
    ([x[1],x[2],x[3],x[5],x[6],x[7]], [x[1],x[2],x[4],x[5],x[8],x[9]], [x[2],x[3],x[4],x[6],x[9],x[10]], [x[1],x[3],x[4],x[7],x[8],x[10]])
coords_on_faces(x, ::Lagrange{3, RefCube, 1}) =
    ([x[1],x[2],x[3],x[4]], [x[1],x[2],x[5],x[6]], [x[2],x[3],x[6],x[7]],[x[3],x[4],x[7],x[8]],[x[1],x[4],x[5],x[8]],[x[5],x[6],x[7],x[8]])
coords_on_faces(x, ::Serendipity{2, RefCube, 2}) =
    ([x[1],x[2],x[5]], [x[2],x[3],x[6]], [x[3],x[4],x[7]], [x[1],x[4],x[8]])


text:
using CUDA
using BenchmarkTools
using GemmKernels
using Statistics
using Printf

function summarise_bench(b)
    mu = mean(b.times)
    sigma = std(b.times)

    if mu < 1e3
        mu, sigma, unit = mu, sigma, "ns"
    elseif mu < 1e6
        mu, sigma, unit = mu / 1e3, sigma / 1e3, "μs"
    elseif mu < 1e9
        mu, sigma, unit = mu / 1e6, sigma / 1e6, "ms"
    else
        mu, sigma, unit = mu / 1e9, sigma / 1e9, "s"
    end

    println(string(@sprintf("(%.2f ± %.2f)", mu, sigma), " ", unit))
end

M = 4096
N = 4096
K = 4096

transpose_a = false
transpose_b = false

f(x) = max(x, 0)
g(x) = x + 1

a_h = rand(Float16, (M, K)) / sqrt(Float16(K))
b_h = rand(Float16, (K, N)) / sqrt(Float16(K))
c_h = rand(Float32, (M, N))

a   = CuArray(a_h)
b   = CuArray(b_h)
c   = CuArray(c_h)

bias = CuArray(rand(Float32, (1, N)))

function bench_cublas(a, b, c, M, N, K, transpose_a, transpose_b)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : 100
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
        end
    end
end

function bench_cublas_relu(a, b, c, M, N, K, transpose_a, transpose_b)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : 100
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
            c = f.(c)
        end
    end
end

function bench_cublas_bias(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : 100
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
            c = c .+ bias
        end
    end
end

function bench_cublas_biasrelu(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : 100
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
            c = f.(c .+ bias)
        end
    end
end

function bench_cublas_biasrelutwice(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : 100
            c = f.(c)
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
            c = f.(c .+ bias)
        end
    end
end

function bench_cublas_biasrelutwice_ab_elop(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUBLAS.cublasSetMathMode(CUBLAS.handle(), CUBLAS.CUBLAS_TENSOR_OP_MATH)

    CUDA.@sync begin
        for i = 1 : 100
            a = g.(a)
            b = g.(b)
            c = f.(c)
            CUDA.CUBLAS.gemmEx!(!transpose_a ? 'N' : 'T', !transpose_b ? 'N' : 'T', Float32(1), a, b, Float32(1), c)
            c = f.(c .+ bias)
        end
    end
end

function bench_gemmkernels(a, b, c, M, N, K, transpose_a, transpose_b)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                                      gemm_shape = (M = M, N = N, K = K),
        operator = Operator.WMMAOp{16, 16, 16, Float32},
        global_a_layout = transpose_a ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},
        global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

        global_c_layout = Layout.AlignedColMajor{Float32},
        global_d_layout = Layout.AlignedColMajor{Float32},

        is_a_col_major = !transpose_a,
        is_b_col_major = !transpose_b,
       )

        for i = 1 : 100
            GemmKernels.matmul(a, b, c, c, conf;
                               kernel = Kernel.matmul_pipelined)
        end
    end
end

function bench_gemmkernels_relu(a, b, c, M, N, K, transpose_a, transpose_b)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                                      gemm_shape = (M = M, N = N, K = K),
        operator = Operator.WMMAOp{16, 16, 16, Float32},
        global_a_layout = transpose_a ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},
        global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

        global_c_layout = Layout.AlignedColMajor{Float32},
        global_d_layout = Layout.AlignedColMajor{Float32},

        is_a_col_major = !transpose_a,
        is_b_col_major = !transpose_b,
       )

        for i = 1 : 100
            GemmKernels.matmul(a, b, c, c, conf;
                               transform_regs_to_shared_d = Transform.Elementwise(f),
                               kernel = Kernel.matmul_pipelined)
        end
    end
end

function bench_gemmkernels_bias(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                                      gemm_shape = (M = M, N = N, K = K),
        operator = Operator.WMMAOp{16, 16, 16, Float32},
        global_a_layout = transpose_a ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},
        global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

        global_c_layout = Layout.AlignedColMajor{Float32},
        global_d_layout = Layout.AlignedColMajor{Float32},

        is_a_col_major = !transpose_a,
        is_b_col_major = !transpose_b,
       )

        for i = 1 : 100
            GemmKernels.matmul(a, b, c, c, conf;
                               epilogue = Epilogue.Bias(pointer(bias)),
                               kernel = Kernel.matmul_pipelined)
        end
    end
end

function bench_gemmkernels_biasrelu(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                                      gemm_shape = (M = M, N = N, K = K),
        operator = Operator.WMMAOp{16, 16, 16, Float32},
        global_a_layout = transpose_a ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},
        global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

        global_c_layout = Layout.AlignedColMajor{Float32},
        global_d_layout = Layout.AlignedColMajor{Float32},

        is_a_col_major = !transpose_a,
        is_b_col_major = !transpose_b,
       )

        for i = 1 : 100
            GemmKernels.matmul(a, b, c, c, conf;
                               transform_regs_to_shared_d = Transform.Elementwise(f),
                               epilogue = Epilogue.Bias(pointer(bias)),
                               kernel = Kernel.matmul_pipelined)
        end
    end
end

function bench_gemmkernels_biasrelutwice(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                                      gemm_shape = (M = M, N = N, K = K),
        operator = Operator.WMMAOp{16, 16, 16, Float32},
        global_a_layout = transpose_a ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},
        global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

        global_c_layout = Layout.AlignedColMajor{Float32},
        global_d_layout = Layout.AlignedColMajor{Float32},

        is_a_col_major = !transpose_a,
        is_b_col_major = !transpose_b,
       )

        for i = 1 : 100
            GemmKernels.matmul(a, b, c, c, conf;
                               transform_shared_to_regs_c = Transform.Elementwise(f),
                               transform_regs_to_shared_d = Transform.Elementwise(f),
                               epilogue = Epilogue.Bias(pointer(bias)),
                               kernel = Kernel.matmul_pipelined)
        end
    end
end

function bench_gemmkernels_biasrelutwice_ab_elop(a, b, c, bias, M, N, K, transpose_a, transpose_b)
    CUDA.@sync begin
        conf = GemmKernels.get_config(
                                      gemm_shape = (M = M, N = N, K = K),
        operator = Operator.WMMAOp{16, 16, 16, Float32},
        global_a_layout = transpose_a ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},
        global_b_layout = transpose_b ? Layout.AlignedRowMajor{Float16} : Layout.AlignedColMajor{Float16},

        global_c_layout = Layout.AlignedColMajor{Float32},
        global_d_layout = Layout.AlignedColMajor{Float32},

        is_a_col_major = !transpose_a,
        is_b_col_major = !transpose_b,
       )

        for i = 1 : 100
            GemmKernels.matmul(a, b, c, c, conf;
                               transform_shared_to_regs_a = Transform.Elementwise(g),
                               transform_shared_to_regs_b = Transform.Elementwise(g),
                               transform_shared_to_regs_c = Transform.Elementwise(f),
                               transform_regs_to_shared_d = Transform.Elementwise(f),
                               epilogue = Epilogue.Bias(pointer(bias)),
                               kernel = Kernel.matmul_pipelined)
        end
    end
end


text:
@testset "Basic Vec3 Functionality" begin
    # Check the constructors
    a = Vec3(1.0)

    @test a.x[] == 1.0 && a.y[] == 1.0 && a.z[] == 1.0

    a = Vec3([1.0, 1.0])

    @test a.x == [1.0, 1.0] && a.y == [1.0, 1.0] && a.z == [1.0, 1.0]

    # Operators +, -, * | Both Vec3
    for op in (:+, :-, :*)
        len = rand(1:25)

        a1 = rand(len)
        b1 = rand(len)
        c1 = rand(len)

        a2 = rand(len)
        b2 = rand(len)
        c2 = rand(len)

        vec3_1 = Vec3(a1, b1, c1)
        vec3_2 = Vec3(a2, b2, c2)

        res = @eval $(op)($vec3_1, $vec3_2)
        res_x = @eval broadcast($(op), $a1, $a2)
        res_y = @eval broadcast($(op), $b1, $b2)
        res_z = @eval broadcast($(op), $c1, $c2)

        @test res.x == res_x && res.y == res_y && res.z == res_z
    end

    # Operators +, -, *, /, % | Only 1 Vec3
    for op in (:+, :-, :*, :/, :%)
        len = rand(1:25)

        a = rand(len)
        b = rand(len)
        c = rand(len)

        num = rand()

        vec3_1 = Vec3(a, b, c)

        res1 = @eval $(op)($vec3_1, $num)
        res1_x = @eval broadcast($(op), $a, $num)
        res1_y = @eval broadcast($(op), $b, $num)
        res1_z = @eval broadcast($(op), $c, $num)

        @test res1.x == res1_x && res1.y == res1_y && res1.z == res1_z

        res2 = @eval $(op)($num, $vec3_1)
        res2_x = @eval broadcast($(op), $a, $num)
        res2_y = @eval broadcast($(op), $b, $num)
        res2_z = @eval broadcast($(op), $c, $num)

        @test res2.x == res2_x && res2.y == res2_y && res2.z == res2_z
    end

    # Other Operators
    a_vec = Vec3(1.0, 2.0, 3.0)
    b_vec = Vec3(4.0, 5.0, 6.0)
    neg_a_vec = - a_vec
    cross_prod = RayTracer.cross(a_vec, b_vec)
    b_clamped = clamp(b_vec, 4.5, 5.5)

    @test neg_a_vec.x[] == -1.0 && neg_a_vec.y[] == -2.0 && neg_a_vec.z[] == -3.0
    @test RayTracer.dot(a_vec, b_vec)[] == 32.0
    @test RayTracer.l2norm(a_vec)[] == 14.0
    @test RayTracer.l2norm(RayTracer.normalize(a_vec))[] == 1.0
    @test maximum(a_vec) == 3.0
    @test minimum(a_vec) == 1.0
    @test size(a_vec) == (1,)
    @test getindex(a_vec, 1) == (x = 1.0, y = 2.0, z = 3.0)
    @test cross_prod.x[] == -3.0 && cross_prod.y[] == 6.0 && cross_prod.z[] == -3.0
    @test b_clamped.x[] == 4.5 && b_clamped.y[] == 5.0 && b_clamped.z[] == 5.5

end


text:
using LinLogQuantization
using Test

@testset "minpos" begin
    @test minpos([0,0,0,0]) == 0
    @test minpos([0,0,1,0]) == 1

    A = randn(100)
    @test minpos(A) > 0
    @test minpos(A) == minimum(A[A.>0])

    A = rand(100) .+ 1
    @test minpos(A) == minimum(A)
end

@testset "Linear quantization" begin
    for T in [Float64,Float32,Float16]
        for s in [(100,),
                    (10,20),
                    (13,14,15),
                    (23,17,12,5)]

            A = rand(T,s...)

            for LinQ in [LinQuant8Array,
                        LinQuant16Array,
                        LinQuant24Array,
                        LinQuant32Array]

                # initial conversion is not reversible
                # due to rounding errors
                A2 = Array{T}(LinQ(A))

                # then test whether back&forth conversion is reversible
                @test A2 == Array{T}(LinQ(A2))
            end
        end
    end
end

@testset "Log quantization" begin
    for T in [Float64,Float32,Float16]
        for s in [(100,),
                    (10,20),
                    (13,14,15),
                    (23,17,12,5)]

            A = rand(T,s...)

            for LogQ in [LogQuant8Array,
                        LogQuant16Array,
                        LogQuant24Array,
                        LogQuant32Array]

                for rn in [:linspace,:logspace]

                    # initial conversion is not reversible
                    # due to rounding errors
                    A2 = Array{T}(LogQ(A,rn))

                    # then test whether back&forth conversion is reversible
                    @test A2 == Array{T}(LogQ(A2,rn))
                end
            end
        end
    end
end

@testset "LinQuant along dimension" begin
    A = rand(Float32,10,20,30,40)
    Q = LinQuant32Array(A,4)
    @test A ≈ Array{Float32}(Q)

    Q = LinQuant24Array(A,4)
    @test A ≈ Array{Float32}(Q)

    Q = LinQuant16Array(A,4)
    @test A ≈ Array{Float32}(Q)

    Q = LinQuant8Array(A,4)
    @test all(isapprox.(A,Array{Float32}(Q),atol=1e-1))
end

@testset "LogQuant along dimension" begin
    A = rand(Float32,10,20,30,40)
    Q = LogQuant32Array(A,4)
    @test A ≈ Array{Float32}(Q)

    Q = LogQuant24Array(A,4)
    @test A ≈ Array{Float32}(Q)

    Q = LogQuant16Array(A,4)
    @test A ≈ Array{Float32}(Q)

    Q = LogQuant8Array(A,4)
    @test all(isapprox.(A,Array{Float32}(Q),atol=1e-1))
end

@testset "Linear quantization will all 0" begin
    for T in [Float64,Float32,Float16]
        for s in [(100,),
                    (10,20),
                    (3,4,5),
                    (6,7,8,9)]

            A = zeros(T,s...)

            for LinQ in [LinQuant8Array,
                        LinQuant16Array,
                        LinQuant24Array,
                        LinQuant32Array]

                # initial conversion is not reversible
                # due to rounding errors
                A2 = Array{T}(LinQ(A))

                # then test whether back&forth conversion is reversible
                @test A2 == Array{T}(LinQ(A2))
            end
        end
    end
end

@testset "Log quantization will all 0" begin
    for T in [Float64,Float32,Float16]
        for s in [(100,),
                    (10,20),
                    (3,4,5),
                    (6,7,8,9)]

            A = zeros(T,s...)

            for LogQ in [LogQuant8Array,
                LogQuant16Array,
                LogQuant24Array,
                LogQuant32Array]

                # initial conversion is not reversible
                # due to rounding errors
                A2 = Array{T}(LogQ(A))

                # then test whether back&forth conversion is reversible
                @test A2 == Array{T}(LogQ(A2))
            end
        end
    end
end

@testset "Linear quantization will all c" begin
    for T in [Float64,Float32,Float16]
        for s in [(100,),
                    (10,20),
                    (3,4,5),
                    (6,7,8,9)]

            A = zeros(T,s...) .+ T(randn())

            for LinQ in [LinQuant8Array,
                        LinQuant16Array,
                        LinQuant24Array,
                        LinQuant32Array]

                # initial conversion is not reversible
                # due to rounding errors
                A2 = Array{T}(LinQ(A))

                # then test whether back&forth conversion is reversible
                @test A2 == Array{T}(LinQ(A2))
            end
        end
    end
end

@testset "Log quantization will all c" begin
    for T in [Float64,Float32,Float16]
        for s in [(100,),
                    (10,20),
                    (3,4,5),
                    (6,7,8,9)]

            A = zeros(T,s...) .+ T(rand())

            for LogQ in [LogQuant8Array,
                LogQuant16Array,
                LogQuant24Array,
                LogQuant32Array]

                # initial conversion is not reversible
                # due to rounding errors
                A2 = Array{T}(LogQ(A))

                # then test whether back&forth conversion is reversible
                @test A2 == Array{T}(LogQ(A2))
            end
        end
    end
end


text:
using ForwardDiff
import Libxc: evaluate!, Functional

function energy_per_particle(::Val{identifier}, args...; kwargs...) where {identifier}
    error("Fallback functional for $identifier not implemented.")
end

include("lda_x.jl")
include("lda_c_vwn.jl")

# Function that always dispatches to the DFTK fallback implementations
function xc_fallback!(func::Functional, ::Val{:lda}, ρ::AbstractArray;
                      zk=nothing, vrho=nothing, v2rho2=nothing)
    func.n_spin == 1  || error("Fallback functionals only for $(func.n_spin) == 1")
    zk = reshape(zk, size(ρ))

    fE = ρ -> energy_per_particle(Val(func.identifier), ρ)
    if zk !== nothing
        zk .= fE.(ρ)
    end

    fV(ρ) = ForwardDiff.derivative(ρ -> ρ*fE(ρ), ρ)
    if vrho !== nothing
        vrho .= fV.(ρ)
    end

    fV2(ρ) = ForwardDiff.derivative(fV, ρ)
    if v2rho2 !== nothing
        v2rho2 .= fV2.(ρ)
    end
end

# For cases where the Array type is not a plain Julia array and the Floating point
# type is not Float64, use xc_fallback! to evaluate the functional.
# Note: This is type piracy on the evaluate! function from Libxc.jl
Libxc.evaluate!(args...; kwargs...) = xc_fallback!(args...; kwargs...)


text:
include("operators.jl")

### Terms
# - A Term is something that, given a state, returns a named tuple (E, hams) with an energy
#   and a list of RFO (for each kpoint).
# - Each term must overload
#     `ene_ops(t, ψ, occ; kwargs...)` -> (E::Real, ops::Vector{RealFourierOperator}).
# - Note that terms are allowed to hold on to references to ψ (eg Fock term),
#   so ψ should not mutated after ene_ops

# The Hamiltonian is defined as half of the gradient of the energy
# with respect to the density matrix sum_n fn |ψn><ψn|.
# In particular, dE/dψn = 2 fn |Hψn> (plus weighting for k-point sampling)
abstract type Term end

### Builders are objects X that store the term parameters, and produce a
# XTerm <: Term when instantiated with a `basis`

"""
A term with a constant zero energy.
"""
struct NoopTerm <: Term
    basis::PlaneWaveBasis
end
function ene_ops(term::NoopTerm, ψ, occ; kwargs...)
    (E=zero(eltype(term.basis)), ops=[NoopOperator(term.basis, kpoint)
                                      for kpoint in term.basis.kpoints])
end

include("Hamiltonian.jl")

# breaks_symmetries on a term builder answers true if this term breaks
# the symmetries of the lattice/atoms (in which case k-point reduction
# is invalid)
breaks_symmetries(term_builder::Any) = false

include("kinetic.jl")

include("local.jl")
breaks_symmetries(term_builder::ExternalFromReal) = true
breaks_symmetries(term_builder::ExternalFromFourier) = true

include("nonlocal.jl")
include("hartree.jl")
include("power_nonlinearity.jl")
include("xc.jl")
include("ewald.jl")
include("psp_correction.jl")
include("entropy.jl")

include("magnetic.jl")
breaks_symmetries(term_builder::Magnetic) = true

include("anyonic.jl")
breaks_symmetries(term_builder::Anyonic) = true

# forces computes either nothing or an array forces[el][at][α]
compute_forces(term::Term, ψ, occ; kwargs...) = nothing  # by default, no force

@doc raw"""
    compute_kernel(basis::PlaneWaveBasis; kwargs...)

Computes a matrix representation of the full response kernel
(derivative of potential with respect to density) in real space.
For non-spin-polarized calculations the matrix dimension is
`prod(basis.fft_size)` × `prod(basis.fft_size)` and
for collinear spin-polarized cases it is
`2prod(basis.fft_size)` × `2prod(basis.fft_size)`.
In this case the matrix has effectively 4 blocks
```math
\left(\begin{array}{cc}
    K_{αα} & K_{αβ}\\
    K_{βα} & K_{ββ}
\end{array}\right)
```
"""
@timing function compute_kernel(basis::PlaneWaveBasis{T}; kwargs...) where {T}
    n_spin = basis.model.n_spin_components
    kernel = zeros(T, n_spin * prod(basis.fft_size), n_spin * prod(basis.fft_size))
    for term in basis.terms
        isnothing(term) && continue
        kernel .+= compute_kernel(term; kwargs...)
    end
    kernel
end
compute_kernel(::Term; kwargs...) = nothing  # By default no kernel

"""
    apply_kernel(basis::PlaneWaveBasis, δρ; kwargs...)

Computes the potential response to a perturbation δρ in real space,
as a 4D (i,j,k,σ) array.
"""
@timing function apply_kernel(basis::PlaneWaveBasis, δρ;
                              RPA=false, kwargs...)
    n_spin = basis.model.n_spin_components
    @assert 1 ≤ n_spin ≤ 2

    δV = zero(δρ)
    for term in basis.terms
        # Skip XC term if RPA is selected
        RPA && term isa TermXc && continue

        δV_term = apply_kernel(term, δρ; kwargs...)
        if !isnothing(δV_term)
            δV .+= δV_term
        end
    end
    δV
end
apply_kernel(::Term, δρ; kwargs...) = nothing  # by default, no kernel


text:
# ellipse to multiple forward simulated shapes
# only final configurations observed

using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD2
using NPZ

workdir = @__DIR__
cd(workdir)

cardiac = npzread("cardiac.npy")  # heart data (left ventricles, the one we used in https://arxiv.org/abs/1705.10943
cardiacx = cardiac[:,:,1]  # x-coordinates of landmarks
cardiacy = cardiac[:,:,2]  # y-coordinates of landmarks

landmarksset = 1:6:66
nshapes = 14

n = length(landmarksset)
xobsT = fill(zeros(PointF,66),nshapes)
for i in 1:nshapes # for each image
    xobsT[i] = [PointF(cardiacx[i,j], cardiacy[i,j]) for j in landmarksset ]
end

xobs0 = 0

JLD2.@save "data_cardiac.jld2" xobs0 xobsT n nshapes


text:
#!/usr/bin/julia
# FizzBuzz in one line
# Author: @joaogui1

map(println, collect(i % 15 == 0 ? "FizzBuzz" : i % 5 == 0 ? "Buzz" : i % 3 == 0 ? "Fizz" : i for i in 1:100))


text:
"""
Return the version of the loaded module.

Source: https://stackoverflow.com/questions/60587041/julia-getting-the-version-number-of-my-module
"""
loadedversion(m::Module) =  VersionNumber(Pkg.TOML.parsefile(abspath(string(first(methods(m.eval)).file), "..", "..", "Project.toml"))["version"])

"Parse MIToS scripts command line arguments."
function parse_commandline(args...; description::AbstractString="Made with MIToS",
                                    output::AbstractString=".mitos.",
                                    stdout::Bool=true, mitos_version="")
    settings = ArgParseSettings(description = description,
                                version = "MIToS $mitos_version",
                                add_version = true,
                                epilog =    """
                                            \n
                                            MIToS $mitos_version\n
                                            \n
                                            Bioinformatics Unit\n
                                            Leloir Institute Foundation\n
                                            Av. Patricias Argentinas 435, CP C1405BWE, Buenos Aires, Argentina
                                            """
                                )

    add_arg_table!(settings,
                  "FILE",
                  Dict(
                      :help => "File name. If it is not used, the script reads from STDIN.",
                      :required => false
                      ),
                  ["--list", "-l"],
                  Dict(
                      :help => "The input is a list of file names. If -p is used, files will be processed in parallel.",
                      :action => :store_true
                      ),
                  ["--output", "-o"],
                  Dict(
                      :help => string( """Name of the output file. Output will be gzip if the extension is ".gz".
                      If it starts with a dot, the name is used as a suffix or extension of the input filename.
                      If it ends with a dot, is used as a prefix.
                      If the output name starts and ends with dots, it's used as an interfix before the extension.""",
                      stdout ? """If a single file is used and there is not a file name (STDIN), the output will be print into
                      STDOUT, unless a output filename is used. You can use "STDOUT" to force print into STDOUT.
                      STDOUT can not be use with --list.""" : ""),
                      :arg_type => AbstractString,
                      :default => output
                      ),
                  ["--parallel", "-p"],
                  Dict(
                      :help => "Number of worker processes.",
                      :arg_type => Int,
                      :default => 1
                      ),

                  args...)

    return parse_args(settings)
end

function _generate_output_name(file, output)
    begins = startswith(output, '.')
    ends = endswith(output,'.')

    if file !== nothing && file !== stdin

        if begins && ends
            parts = split(file, '.')
            return( string(join(parts[1:end-1],'.'), output, parts[end]) )
        elseif begins
            return( string(file, output) )
        elseif ends
            return( string(output, file) )
        else
            return(output)
        end

    elseif !begins && !ends

        return(output)

    else

        return("STDOUT")

    end
end

"Opens the output file or returns STDOUT."
function open_output(file, output)
    output_name = _generate_output_name(file, output)
    if output_name != "STDOUT"
        fh = open(output_name, "w")
        if endswith(output_name, ".gz")
            return(GzipCompressorStream(fh))
        end
        return(fh)
    else
        return(stdout)
    end
end

"Close output (check if output is STDOUT)."
function close_output(fh_out)
    fh_out === stdout ? nothing : close(fh_out)
    nothing
end

"Adds the needed number of workers."
function set_parallel(parallel)
    N = nprocs()
    if N < parallel
        addprocs(parallel - N + 1)
    end
    nothing
end

script(x, y, z) = throw("Define your script function!")

"""
Opens and closes the output stream, tries to run the defined script.
"""
function run_single_script(input, args)
    fh_out = open_output(input, args["output"])
    try
        script(input, args, fh_out)
    catch err
        println("ERROR: ", input)
        println(err)
    finally
        close(fh_out)
    end
end

"""
If FILE is not used, calls run_single_script with STDIN. Otherwise calls run_single_script with the FILE.
If `list` is `true`, for loop or `pmap` is used over eachline of the input.
"""
function runscript(args)
    file = args["FILE"]
    if args["list"]
        fh_in = file === nothing ? stdin : open(file, "r")
        if nprocs() != 1
            pmap(eachline(fh_in)) do line
                run_single_script(line, args)
            end
        else
            for line in eachline(fh_in)
                run_single_script(line, args)
            end
        end
    else
        run_single_script(file === nothing ? stdin : file, args)
    end
end

"Decides if read or parse, uses parse with STDIN"
readorparse(input::AbstractString, args...; kargs...) = read(input, args...; kargs...)
readorparse(input::Base.LibuvStream, args...; kargs...) = parse(input, args...; kargs...)


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>5.0, :alg=>INT(δ=0.01), :fixpoint_check=>false),
    "X0"   => [(1, Interval(68., 69.))],
    "info" => Dict(
        "name" => "thermostat",
        "dim" => 1,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary