text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module RangesTest

using Test
using Dates

using InteractiveUtils: subtypes

let
    for T in (DatesPlus.Date, DatesPlus.DateTime)
        f1 = T(2014); l1 = T(2013, 12, 31)
        f2 = T(2014); l2 = T(2014)
        f3 = T(-2000); l3 = T(2000)
        f4 = typemin(T); l4 = typemax(T)

        for P in subtypes(DatesPlus.DatePeriod)
            for pos_step in (P(1), P(2), P(50), P(2048), P(10000))
                # empty range
                dr = f1:pos_step:l1
                len = length(dr)
                @test len == 0
                @test isa(len, Int64)
                @test isempty(dr)
                @test first(dr) == f1
                @test last(dr) < f1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr), dr) == Int64[]
                @test [dr;] == T[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) < f1
                @test last(reverse(dr)) >= f1
                @test issorted(dr)
                @test sortperm(dr) === StepRange{Int64,Int}(1:1:0)
                @test !(f1 in dr)
                @test !(l1 in dr)
                @test !(f1 - pos_step in dr)
                @test !(l1 + pos_step in dr)
                @test dr == []
                @test hash(dr) == hash([])

                for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                    dr = f:pos_step:l
                    len = length(dr)
                    @test len > 0
                    @test isa(len, Int64)
                    @test !isempty(dr)
                    @test first(dr) == f
                    @test last(dr) <= l
                    @test minimum(dr) == first(dr)
                    @test maximum(dr) == last(dr)
                    @test dr[1] == f
                    @test dr[end] <= l
                    @test iterate(dr) == (first(dr), (length(dr), 1))

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr), dr) == [1:len;]
                        @test length([dr;]) == len
                        @test dr == dr1
                        @test hash(dr) == hash(dr1)
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test last(reverse(dr)) == f
                    @test issorted(dr)
                    @test f in dr

                end
            end
            for neg_step in (P(-1), P(-2), P(-50), P(-2048), P(-10000))
                # empty range
                dr = l1:neg_step:f1
                len = length(dr)
                @test len == 0
                @test isa(len, Int64)
                @test isempty(dr)
                @test first(dr) == l1
                @test last(dr) > l1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr), dr) == Int64[]
                @test [dr;] == T[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) > l1
                @test last(reverse(dr)) <= l1
                @test issorted(dr)
                @test sortperm(dr) === StepRange{Int64,Int}(1:1:0)
                @test !(l1 in dr)
                @test !(l1 in dr)
                @test !(l1 - neg_step in dr)
                @test !(l1 + neg_step in dr)
                @test dr == []
                @test hash(dr) == hash([])

                for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                    dr = l:neg_step:f
                    len = length(dr)
                    @test len > 0
                    @test isa(len, Int64)
                    @test !isempty(dr)
                    @test first(dr) == l
                    @test last(dr) >= f
                    @test minimum(dr) == last(dr)
                    @test maximum(dr) == first(dr)
                    @test dr[1] == l
                    @test dr[end] >= f
                    @test iterate(dr) == (first(dr), (length(dr), 1))

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr), dr) == [1:len;]
                        @test length([dr;]) == len
                        @test dr == dr1
                        @test hash(dr) == hash(dr1)
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test issorted(dr) == (len <= 1)
                    @test l in dr
                end
            end
        end
        if T == DatesPlus.DateTime
            for P in subtypes(DatesPlus.TimePeriod)
                P in (DatesPlus.Microsecond, DatesPlus.Nanosecond) && continue
                for pos_step in (P(1), P(2), P(50), P(2048), P(10000))
                    # empty range
                    dr = f1:pos_step:l1
                    len = length(dr)
                    @test len == 0
                    @test isa(len, Int64)
                    @test isempty(dr)
                    @test first(dr) == f1
                    @test last(dr) < f1
                    @test length([i for i in dr]) == 0
                    @test_throws ArgumentError minimum(dr)
                    @test_throws ArgumentError maximum(dr)
                    @test_throws BoundsError dr[1]
                    @test findall(in(dr), dr) == Int64[]
                    @test [dr;] == T[]
                    @test isempty(reverse(dr))
                    @test length(reverse(dr)) == 0
                    @test first(reverse(dr)) < f1
                    @test last(reverse(dr)) >= f1
                    @test issorted(dr)
                    @test sortperm(dr) === StepRange{Int64,Int}(1:1:0)
                    @test !(f1 in dr)
                    @test !(l1 in dr)
                    @test !(f1 - pos_step in dr)
                    @test !(l1 + pos_step in dr)
                    @test dr == []
                    @test hash(dr) == hash([])

                    for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                        dr = f:pos_step:l
                        len = length(dr)
                        @test len > 0
                        @test isa(len, Int64)
                        @test !isempty(dr)
                        @test first(dr) == f
                        @test last(dr) <= l
                        @test minimum(dr) == first(dr)
                        @test maximum(dr) == last(dr)
                        @test dr[1] == f
                        @test dr[end] <= l
                        @test iterate(dr) == (first(dr), (length(dr), 1))

                        if len < 10000
                            dr1 = [i for i in dr]
                            @test length(dr1) == len
                            @test findall(in(dr), dr) == [1:len;]
                            @test length([dr;]) == len
                            @test dr == dr1
                            @test hash(dr) == hash(dr1)
                        end
                        @test !isempty(reverse(dr))
                        @test length(reverse(dr)) == len
                        @test last(reverse(dr)) == f
                        @test issorted(dr)
                        @test f in dr

                    end
                end
                for neg_step in (P(-1), P(-2), P(-50), P(-2048), P(-10000))
                    # empty range
                    dr = l1:neg_step:f1
                    len = length(dr)
                    @test len == 0
                    @test isa(len, Int64)
                    @test isempty(dr)
                    @test first(dr) == l1
                    @test last(dr) > l1
                    @test length([i for i in dr]) == 0
                    @test_throws ArgumentError minimum(dr)
                    @test_throws ArgumentError maximum(dr)
                    @test_throws BoundsError dr[1]
                    @test findall(in(dr), dr) == Int64[]
                    @test [dr;] == T[]
                    @test isempty(reverse(dr))
                    @test length(reverse(dr)) == 0
                    @test first(reverse(dr)) > l1
                    @test last(reverse(dr)) <= l1
                    @test issorted(dr)
                    @test sortperm(dr) === StepRange{Int64,Int}(1:1:0)
                    @test !(l1 in dr)
                    @test !(l1 in dr)
                    @test !(l1 - neg_step in dr)
                    @test !(l1 + neg_step in dr)
                    @test dr == []
                    @test hash(dr) == hash([])

                    for (f, l) in ((f2, l2), (f3, l3), (f4, l4))
                        dr = l:neg_step:f
                        len = length(dr)
                        @test len > 0
                        @test isa(len, Int64)
                        @test !isempty(dr)
                        @test first(dr) == l
                        @test last(dr) >= f
                        @test minimum(dr) == last(dr)
                        @test maximum(dr) == first(dr)
                        @test dr[1] == l
                        @test dr[end] >= f
                        @test iterate(dr) == (first(dr), (length(dr), 1))

                        if len < 10000
                            dr1 = [i for i in dr]
                            @test length(dr1) == len
                            @test findall(in(dr), dr) == [1:len;]
                            @test length([dr;]) == len
                            @test dr == dr1
                            @test hash(dr) == hash(dr1)
                        end
                        @test !isempty(reverse(dr))
                        @test length(reverse(dr)) == len
                        @test last(reverse(dr)) >= l
                        @test issorted(dr) == (len <= 1)
                        @test l in dr

                    end
                end
            end
        end
    end
end

# Dates are physical units, and ranges should require an explicit step.
# See #19896 and https://discourse.julialang.org/t/type-restriction-on-unitrange/6557/12
if VERSION >= v"1.0.0"
    @test_throws MethodError DatesPlus.DateTime(2013, 1, 1):DatesPlus.DateTime(2013, 2, 1)
end

# All the range representations we want to test
# Date ranges
dr  = DatesPlus.DateTime(2013, 1, 1):DatesPlus.Day(1):DatesPlus.DateTime(2013, 2, 1)
dr1 = DatesPlus.DateTime(2013, 1, 1):DatesPlus.Day(1):DatesPlus.DateTime(2013, 1, 1)
dr2 = DatesPlus.DateTime(2013, 1, 1):DatesPlus.Day(1):DatesPlus.DateTime(2012, 2, 1) # empty range
dr3 = DatesPlus.DateTime(2013, 1, 1):DatesPlus.Day(-1):DatesPlus.DateTime(2012) # negative step
# Big ranges
dr4 = DatesPlus.DateTime(0):DatesPlus.Day(1):DatesPlus.DateTime(20000, 1, 1)
dr5 = DatesPlus.DateTime(0):DatesPlus.Day(1):DatesPlus.DateTime(200000, 1, 1)
dr6 = DatesPlus.DateTime(0):DatesPlus.Day(1):DatesPlus.DateTime(2000000, 1, 1)
dr7 = DatesPlus.DateTime(0):DatesPlus.Day(1):DatesPlus.DateTime(20000000, 1, 1)
dr8 = DatesPlus.DateTime(0):DatesPlus.Day(1):DatesPlus.DateTime(200000000, 1, 1)
dr9 = typemin(DatesPlus.DateTime):DatesPlus.Day(1):typemax(DatesPlus.DateTime)
# Other steps
dr10 = typemax(DatesPlus.DateTime):DatesPlus.Day(-1):typemin(DatesPlus.DateTime)
dr11 = typemin(DatesPlus.DateTime):DatesPlus.Week(1):typemax(DatesPlus.DateTime)

dr12 = typemin(DatesPlus.DateTime):DatesPlus.Month(1):typemax(DatesPlus.DateTime)
dr13 = typemin(DatesPlus.DateTime):DatesPlus.Year(1):typemax(DatesPlus.DateTime)

dr14 = typemin(DatesPlus.DateTime):DatesPlus.Week(10):typemax(DatesPlus.DateTime)
dr15 = typemin(DatesPlus.DateTime):DatesPlus.Month(100):typemax(DatesPlus.DateTime)
dr16 = typemin(DatesPlus.DateTime):DatesPlus.Year(1000):typemax(DatesPlus.DateTime)
dr17 = typemax(DatesPlus.DateTime):DatesPlus.Week(-10000):typemin(DatesPlus.DateTime)
dr18 = typemax(DatesPlus.DateTime):DatesPlus.Month(-100000):typemin(DatesPlus.DateTime)
dr19 = typemax(DatesPlus.DateTime):DatesPlus.Year(-1000000):typemin(DatesPlus.DateTime)
dr20 = typemin(DatesPlus.DateTime):DatesPlus.Day(2):typemax(DatesPlus.DateTime)

drs = Any[dr, dr1, dr2, dr3, dr4, dr5, dr6, dr7, dr8, dr9, dr10,
          dr11, dr12, dr13, dr14, dr15, dr16, dr17, dr18, dr19, dr20]
drs2 = map(x->DatesPlus.Date(first(x)):step(x):DatesPlus.Date(last(x)), drs)

@test map(length, drs) == map(x->size(x)[1], drs)
@test map(length, drs) == map(x->length(DatesPlus.Date(first(x)):step(x):DatesPlus.Date(last(x))), drs)
@test map(length, drs) == map(x->length(reverse(x)), drs)
@test all(x->findall(in(x), x)==[1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == range(last(x), step=-step(x), length=length(x)), drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test all(drs[1:3]) do dd
    for (i, d) in enumerate(dd)
        @test d == (first(dd) + DatesPlus.Day(i - 1))
    end
    true
end
@test_throws MethodError dr .+ 1
a = DatesPlus.DateTime(2013, 1, 1)
b = DatesPlus.DateTime(2013, 2, 1)
@test map!(x->x + DatesPlus.Day(1), Vector{DatesPlus.DateTime}(undef, 32), dr) == [(a + DatesPlus.Day(1)):DatesPlus.Day(1):(b + DatesPlus.Day(1));]
@test map(x->x + DatesPlus.Day(1), dr) == [(a + DatesPlus.Day(1)):DatesPlus.Day(1):(b + DatesPlus.Day(1));]

@test map(x->a in x, drs[1:4]) == [true, true, false, true]
@test a in dr
@test b in dr
@test DatesPlus.DateTime(2013, 1, 3) in dr
@test DatesPlus.DateTime(2013, 1, 15) in dr
@test DatesPlus.DateTime(2013, 1, 26) in dr
@test !(DatesPlus.DateTime(2012, 1, 1) in dr)

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test length(b:DatesPlus.Day(-1):a) == 32
@test length(b:DatesPlus.Day(1):a) == 0
@test length(b:DatesPlus.Day(1):a) == 0
@test length(a:DatesPlus.Day(2):b) == 16
@test last(a:DatesPlus.Day(2):b) == DatesPlus.DateTime(2013, 1, 31)
@test length(a:DatesPlus.Day(7):b) == 5
@test last(a:DatesPlus.Day(7):b) == DatesPlus.DateTime(2013, 1, 29)
@test length(a:DatesPlus.Day(32):b) == 1
@test last(a:DatesPlus.Day(32):b) == DatesPlus.DateTime(2013, 1, 1)
@test (a:DatesPlus.Day(1):b)[1] == DatesPlus.DateTime(2013, 1, 1)
@test (a:DatesPlus.Day(1):b)[2] == DatesPlus.DateTime(2013, 1, 2)
@test (a:DatesPlus.Day(1):b)[7] == DatesPlus.DateTime(2013, 1, 7)
@test (a:DatesPlus.Day(1):b)[end] == b
@test first(a:DatesPlus.Day(1):DatesPlus.DateTime(20000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.DateTime(200000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.DateTime(2000000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.DateTime(20000000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.DateTime(200000000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):typemax(DatesPlus.DateTime)) == a
@test first(typemin(DatesPlus.DateTime):DatesPlus.Day(1):typemax(DatesPlus.DateTime)) == typemin(DatesPlus.DateTime)

# Date ranges
dr  = DatesPlus.Date(2013, 1, 1):DatesPlus.Day(1):DatesPlus.Date(2013, 2, 1)
dr1 = DatesPlus.Date(2013, 1, 1):DatesPlus.Day(1):DatesPlus.Date(2013, 1, 1)
dr2 = DatesPlus.Date(2013, 1, 1):DatesPlus.Day(1):DatesPlus.Date(2012, 2, 1) # empty range
dr3 = DatesPlus.Date(2013, 1, 1):DatesPlus.Day(-1):DatesPlus.Date(2012, 1, 1) # negative step
# Big ranges
dr4 = DatesPlus.Date(0):DatesPlus.Day(1):DatesPlus.Date(20000, 1, 1)
dr5 = DatesPlus.Date(0):DatesPlus.Day(1):DatesPlus.Date(200000, 1, 1)
dr6 = DatesPlus.Date(0):DatesPlus.Day(1):DatesPlus.Date(2000000, 1, 1)
dr7 = DatesPlus.Date(0):DatesPlus.Day(1):DatesPlus.Date(20000000, 1, 1)
dr8 = DatesPlus.Date(0):DatesPlus.Day(1):DatesPlus.Date(200000000, 1, 1)
dr9 = typemin(DatesPlus.Date):DatesPlus.Day(1):typemax(DatesPlus.Date)
# Other steps
dr10 = typemax(DatesPlus.Date):DatesPlus.Day(-1):typemin(DatesPlus.Date)
dr11 = typemin(DatesPlus.Date):DatesPlus.Week(1):typemax(DatesPlus.Date)
dr12 = typemin(DatesPlus.Date):DatesPlus.Month(1):typemax(DatesPlus.Date)
dr13 = typemin(DatesPlus.Date):DatesPlus.Year(1):typemax(DatesPlus.Date)
dr14 = typemin(DatesPlus.Date):DatesPlus.Week(10):typemax(DatesPlus.Date)
dr15 = typemin(DatesPlus.Date):DatesPlus.Month(100):typemax(DatesPlus.Date)
dr16 = typemin(DatesPlus.Date):DatesPlus.Year(1000):typemax(DatesPlus.Date)
dr17 = typemax(DatesPlus.Date):DatesPlus.Week(-10000):typemin(DatesPlus.Date)
dr18 = typemax(DatesPlus.Date):DatesPlus.Month(-100000):typemin(DatesPlus.Date)
dr19 = typemax(DatesPlus.Date):DatesPlus.Year(-1000000):typemin(DatesPlus.Date)
dr20 = typemin(DatesPlus.Date):DatesPlus.Day(2):typemax(DatesPlus.Date)

drs = Any[dr, dr1, dr2, dr3, dr4, dr5, dr6, dr7, dr8, dr9, dr10,
          dr11, dr12, dr13, dr14, dr15, dr16, dr17, dr18, dr19, dr20]

@test map(length, drs) == map(x->size(x)[1], drs)
@test all(x->findall(in(x), x) == [1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == last(x): - step(x):first(x), drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test all(drs[1:3]) do dd
    for (i, d) in enumerate(dd)
        @test d == (first(dd) + DatesPlus.Day(i - 1))
    end
    true
end
@test_throws MethodError dr .+ 1
a = DatesPlus.Date(2013, 1, 1)
b = DatesPlus.Date(2013, 2, 1)
@test map!(x->x + DatesPlus.Day(1), Vector{DatesPlus.Date}(undef, 32), dr) == [(a + DatesPlus.Day(1)):DatesPlus.Day(1):(b + DatesPlus.Day(1));]
@test map(x->x + DatesPlus.Day(1), dr) == [(a + DatesPlus.Day(1)):DatesPlus.Day(1):(b + DatesPlus.Day(1));]

@test map(x->a in x, drs[1:4]) == [true, true, false, true]
@test a in dr
@test b in dr
@test DatesPlus.Date(2013, 1, 3) in dr
@test DatesPlus.Date(2013, 1, 15) in dr
@test DatesPlus.Date(2013, 1, 26) in dr
@test !(DatesPlus.Date(2012, 1, 1) in dr)

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test length(b:DatesPlus.Day(-1):a) == 32
@test length(b:DatesPlus.Day(1):a) == 0
@test length(a:DatesPlus.Day(2):b) == 16
@test last(a:DatesPlus.Day(2):b) == DatesPlus.Date(2013, 1, 31)
@test length(a:DatesPlus.Day(7):b) == 5
@test last(a:DatesPlus.Day(7):b) == DatesPlus.Date(2013, 1, 29)
@test length(a:DatesPlus.Day(32):b) == 1
@test last(a:DatesPlus.Day(32):b) == DatesPlus.Date(2013, 1, 1)
@test (a:DatesPlus.Day(1):b)[1] == DatesPlus.Date(2013, 1, 1)
@test (a:DatesPlus.Day(1):b)[2] == DatesPlus.Date(2013, 1, 2)
@test (a:DatesPlus.Day(1):b)[7] == DatesPlus.Date(2013, 1, 7)
@test (a:DatesPlus.Day(1):b)[end] == b
@test first(a:DatesPlus.Day(1):DatesPlus.Date(20000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.Date(200000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.Date(2000000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.Date(20000000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):DatesPlus.Date(200000000, 1, 1)) == a
@test first(a:DatesPlus.Day(1):typemax(DatesPlus.Date)) == a
@test first(typemin(DatesPlus.Date):DatesPlus.Day(1):typemax(DatesPlus.Date)) == typemin(DatesPlus.Date)

@test length(typemin(DatesPlus.Date):DatesPlus.Week(1):typemax(DatesPlus.Date)) == 26351950414948059
# Big Month/Year ranges
@test length(typemin(DatesPlus.Date):DatesPlus.Month(1):typemax(DatesPlus.Date)) == 6060531933867600
@test length(typemin(DatesPlus.Date):DatesPlus.Year(1):typemax(DatesPlus.Date)) == 505044327822300
@test length(typemin(DatesPlus.DateTime):DatesPlus.Month(1):typemax(DatesPlus.DateTime)) == 3507324288
@test length(typemin(DatesPlus.DateTime):DatesPlus.Year(1):typemax(DatesPlus.DateTime)) == 292277024

@test length(typemin(DatesPlus.DateTime):DatesPlus.Week(1):typemax(DatesPlus.DateTime)) == 15250284420
@test length(typemin(DatesPlus.DateTime):DatesPlus.Day(1):typemax(DatesPlus.DateTime)) == 106751990938
@test length(typemin(DatesPlus.DateTime):DatesPlus.Hour(1):typemax(DatesPlus.DateTime)) == 2562047782512
@test length(typemin(DatesPlus.DateTime):DatesPlus.Minute(1):typemax(DatesPlus.DateTime)) == 153722866950720
@test length(typemin(DatesPlus.DateTime):DatesPlus.Second(1):typemax(DatesPlus.DateTime)) == 9223372017043200
@test length(typemin(DateTime):DatesPlus.Millisecond(1):typemax(DateTime)) == 9223372017043199001

c = DatesPlus.Date(2013, 6, 1)
@test length(a:DatesPlus.Month(1):c) == 6
@test [a:DatesPlus.Month(1):c;] == [a + DatesPlus.Month(1)*i for i in 0:5]
@test [a:DatesPlus.Month(2):DatesPlus.Date(2013, 1, 2);] == [a]
@test [c:DatesPlus.Month(-1):a;] == reverse([a:DatesPlus.Month(1):c;])

@test length(range(Date(2000), step=DatesPlus.Day(1), length=366)) == 366
let n = 100000
    local a = DatesPlus.Date(2000)
    for i = 1:n
        @test length(range(a, step=DatesPlus.Day(1), length=i)) == i
    end
    return a + DatesPlus.Day(n)
end

@test typeof(step(DatesPlus.DateTime(2000):DatesPlus.Day(1):DatesPlus.DateTime(2001))) == DatesPlus.Day

a = DatesPlus.Date(2013, 1, 1)
b = DatesPlus.Date(2013, 2, 1)
d = DatesPlus.Date(2020, 1, 1)
@test length(a:DatesPlus.Year(1):d) == 8
@test first(a:DatesPlus.Year(1):d) == a
@test last(a:DatesPlus.Year(1):d) == d
@test length(a:DatesPlus.Month(12):d) == 8
@test first(a:DatesPlus.Month(12):d) == a
@test last(a:DatesPlus.Month(12):d) == d
@test length(a:DatesPlus.Week(52):d) == 8
@test first(a:DatesPlus.Week(52):d) == a
@test last(a:DatesPlus.Week(52):d) == DatesPlus.Date(2019, 12, 24)
@test length(a:DatesPlus.Day(365):d) == 8
@test first(a:DatesPlus.Day(365):d) == a
@test last(a:DatesPlus.Day(365):d) == DatesPlus.Date(2019, 12, 31)

a = DatesPlus.Date(2013, 1, 1)
b = DatesPlus.Date(2013, 2, 1)
@test length(a:DatesPlus.Year(1):DatesPlus.Date(2020, 2, 1)) == 8
@test length(a:DatesPlus.Year(1):DatesPlus.Date(2020, 6, 1)) == 8
@test length(a:DatesPlus.Year(1):DatesPlus.Date(2020, 11, 1)) == 8
@test length(a:DatesPlus.Year(1):DatesPlus.Date(2020, 12, 31)) == 8
@test length(a:DatesPlus.Year(1):DatesPlus.Date(2021, 1, 1)) == 9
@test length(DatesPlus.Date(2000):DatesPlus.Year(-10):DatesPlus.Date(1900)) == 11
@test length(DatesPlus.Date(2000, 6, 23):DatesPlus.Year(-10):DatesPlus.Date(1900, 2, 28)) == 11
@test length(DatesPlus.Date(2000, 1, 1):DatesPlus.Year(1):DatesPlus.Date(2000, 2, 1)) == 1

let n = 100000
    local a, b
    a= b = DatesPlus.Date(0)
    for i = 1:n
        @test length(a:DatesPlus.Year(1):b) == i
        b += DatesPlus.Year(1)
    end
end

let n = 10000,
    a = DatesPlus.Date(1985, 12, 5),
    b = DatesPlus.Date(1986, 12, 27),
    c = DatesPlus.DateTime(1985, 12, 5),
    d = DatesPlus.DateTime(1986, 12, 27)
    for i = 1:n
        @test length(a:DatesPlus.Month(1):b) == 13
        @test length(a:DatesPlus.Year(1):b) == 2
        @test length(c:DatesPlus.Month(1):d) == 13
        @test length(c:DatesPlus.Year(1):d) == 2
        a += DatesPlus.Day(1)
        b += DatesPlus.Day(1)
    end
end

let n = 100000
    local a, b
    a = b = DatesPlus.Date(2000)
    for i = 1:n
        @test length(a:DatesPlus.Month(1):b) == i
        b += DatesPlus.Month(1)
    end
end

@test length(DatesPlus.Year(1):DatesPlus.Year(1):DatesPlus.Year(10)) == 10
@test length(DatesPlus.Year(10):DatesPlus.Year(-1):DatesPlus.Year(1)) == 10
@test length(DatesPlus.Year(10):DatesPlus.Year(-2):DatesPlus.Year(1)) == 5
@test length(typemin(DatesPlus.Year):DatesPlus.Year(1):typemax(DatesPlus.Year)) == 0 # overflow
@test_throws MethodError DatesPlus.Date(0):DatesPlus.DateTime(2000)
@test_throws MethodError DatesPlus.Date(0):DatesPlus.Year(10)
@test length(range(DatesPlus.Date(2000), step=DatesPlus.Day(1), length=366)) == 366
@test last(range(DatesPlus.Date(2000), step=DatesPlus.Day(1), length=366)) == DatesPlus.Date(2000, 12, 31)
@test last(range(DatesPlus.Date(2001), step=DatesPlus.Day(1), length=365)) == DatesPlus.Date(2001, 12, 31)
@test last(range(DatesPlus.Date(2000), step=DatesPlus.Day(1), length=367)) == last(range(DatesPlus.Date(2000), step=DatesPlus.Month(12), length=2)) == last(range(DatesPlus.Date(2000), step=DatesPlus.Year(1), length=2))
@test last(range(DatesPlus.DateTime(2000), step=DatesPlus.Day(366), length=2)) == last(range(DatesPlus.DateTime(2000), step=DatesPlus.Hour(8784), length=2))

# Issue 5
lastdaysofmonth = [DatesPlus.Date(2014, i, DatesPlus.daysinmonth(2014, i)) for i=1:12]
@test [Date(2014, 1, 31):DatesPlus.Month(1):Date(2015);] == lastdaysofmonth

# Range addition/subtraction:
let d = DatesPlus.Day(1)
    @test (DatesPlus.Date(2000):d:DatesPlus.Date(2001)) + d == (DatesPlus.Date(2000) + d:d:DatesPlus.Date(2001) + d)
    @test (DatesPlus.Date(2000):d:DatesPlus.Date(2001)) - d == (DatesPlus.Date(2000) - d:d:DatesPlus.Date(2001) - d)
end

# Time ranges
dr  = DatesPlus.Time(23, 1, 1):DatesPlus.Second(1):DatesPlus.Time(23, 2, 1)
dr1 = DatesPlus.Time(23, 1, 1):DatesPlus.Second(1):DatesPlus.Time(23, 1, 1)
dr2 = DatesPlus.Time(23, 1, 1):DatesPlus.Second(1):DatesPlus.Time(22, 2, 1) # empty range
dr3 = DatesPlus.Time(23, 1, 1):DatesPlus.Minute(-1):DatesPlus.Time(22, 1, 1) # negative step
# Big ranges
dr8 = typemin(DatesPlus.Time):DatesPlus.Second(1):typemax(DatesPlus.Time)
dr9 = typemin(DatesPlus.Time):DatesPlus.Nanosecond(1):typemax(DatesPlus.Time)
# Other steps
dr10 = typemax(DatesPlus.Time):DatesPlus.Microsecond(-1):typemin(DatesPlus.Time)
dr11 = typemin(DatesPlus.Time):DatesPlus.Millisecond(1):typemax(DatesPlus.Time)
dr12 = typemin(DatesPlus.Time):DatesPlus.Minute(1):typemax(DatesPlus.Time)
dr13 = typemin(DatesPlus.Time):DatesPlus.Hour(1):typemax(DatesPlus.Time)
dr14 = typemin(DatesPlus.Time):DatesPlus.Millisecond(10):typemax(DatesPlus.Time)
dr15 = typemin(DatesPlus.Time):DatesPlus.Minute(100):typemax(DatesPlus.Time)
dr16 = typemin(DatesPlus.Time):DatesPlus.Hour(1000):typemax(DatesPlus.Time)
dr17 = typemax(DatesPlus.Time):DatesPlus.Millisecond(-10000):typemin(DatesPlus.Time)
dr18 = typemax(DatesPlus.Time):DatesPlus.Minute(-100):typemin(DatesPlus.Time)
dr19 = typemax(DatesPlus.Time):DatesPlus.Hour(-10):typemin(DatesPlus.Time)
dr20 = typemin(DatesPlus.Time):DatesPlus.Microsecond(2):typemax(DatesPlus.Time)

drs = Any[dr, dr1, dr2, dr3, dr8, dr9, dr10,
          dr11, dr12, dr13, dr14, dr15, dr16, dr17, dr18, dr19, dr20]

@test map(length, drs) == map(x->size(x)[1], drs)
@test all(x->findall(in(x), x) == [1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == last(x): - step(x):first(x), drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)), drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)), drs[4:end])
@test_throws MethodError dr .+ 1

a = DatesPlus.Time(23, 1, 1)
@test map(x->a in x, drs[1:4]) == [true, true, false, true]
@test a in dr

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x), drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x), drs)

@test !(1 in Date(2017, 01, 01):DatesPlus.Day(1):Date(2017, 01, 05))
@test !(Complex(1, 0) in Date(2017, 01, 01):DatesPlus.Day(1):Date(2017, 01, 05))
@test !(π in Date(2017, 01, 01):DatesPlus.Day(1):Date(2017, 01, 05))
@test !("a" in Date(2017, 01, 01):DatesPlus.Day(1):Date(2017, 01, 05))

@test hash(Any[Date("2018-1-03"), Date("2018-1-04"), Date("2018-1-05")]) ==
      hash([Date("2018-1-03"), Date("2018-1-04"), Date("2018-1-05")]) ==
      hash(Date("2018-1-03"):Day(1):Date("2018-1-05"))

@testset "range overflow" begin
    # DateTime ranges interactions with overflow. If not handled correctly `DatesPlus.len` could
    # infinite loop
    @test length(DateTime(0):typemax(Millisecond):DateTime(0)) == 1
    @test length(typemax(DateTime):typemax(Millisecond):typemax(DateTime)) == 1

    # Overflow interaction is easier to comprehend with using UTM extremes
    utm_typemin = DateTime(DatesPlus.UTM(typemin(Int64)))
    utm_typemax = DateTime(DatesPlus.UTM(typemax(Int64)))

    @test length(utm_typemax:Millisecond(1):utm_typemax) == 1
    @test length(utm_typemin:-Millisecond(1):utm_typemin) == 1
end

end  # RangesTest module


text:
# testing performance of AD
# conclusion is that it scales linearly with dimension;
# however it is a factor 100 slower than the objective.

using BenchmarkTools
using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile
import ForwardDiff

ϕ(r) = exp( - 6 * (abs(r)-1)) - 2 * exp( - 3 * (abs(r)-1))
dϕ(r) = ForwardDiff.derivative(ϕ, r)
ρ(r) = exp( - 4 * (abs(r)-1) )
dρ(r) = ForwardDiff.derivative(ρ, r)
ψ(t) = sqrt(t)
dψ(t) = ForwardDiff.derivative(ψ, t)

function get_rho(x)
   N = length(x)
   rho = zeros(eltype(x), N)
   for n = 2:N
      t = ρ(x[n]-x[n-1])
      rho[n-1] += t
      rho[n] += 1
   end
   for n = 3:N
      t = ρ(x[n] - x[n-2])
      rho[n] += t
      rho[n-2] += t
   end
   return rho
end

function F(x)
   N = length(x)
   E = 0.0
   for j = 1:2, n = (1+j):N
      E += ϕ(x[n]-x[n-j])
   end
   return E + sum(ψ, get_rho(x))
end

function dF(x)
   N = length(x)
   dE = zeros(N)
   F = ψ.(get_rho(x))
   for j = 1:2, n = (1+j):N
      dϕ_ = dϕ(x[n]-x[n-j])
      dρ_ = dρ(x[n]-x[n-j])
      a = dϕ_ + (F[n] - F[n-j]) * dρ_
      dE[n] += a
      dE[n-j] -= a
   end
   return dE
end

x = gresult = compiled_f_tape = nothing
for N in (100, 400, 1600)
   # pre-record a GradientTape for `F` using inputs with Float64 elements
   f_tape = GradientTape(F, rand(N))
   # compile `f_tape` into a more optimized representation
   compiled_f_tape = compile(f_tape)
   # some inputs and work buffers to play around with
   gresult = zeros(N)
   x = collect(1:N) + rand(N)
   println("N = $N")
   print("                 Time for F: ")
   @btime F($x)
   print("  Time for DF (ForwardDiff): ")
   @btime ForwardDiff.gradient($F, $x)
   print("  Time for DF (ReverseDiff): ")
   @btime gradient!($gresult, $compiled_f_tape, $x)
   print("        Time for DF(manual): ")
   @btime dF($x)
end

# RESULTS 03/04/2018, MacBook Pro, Julia v0.6.2,
# - ForwardDiff                   0.7.4
# - ReverseDiff                   0.2.0
#
# N = 8000
#                Time for F:
#   564.137 μs (2 allocations: 62.58 KiB)
#   Time for DF (ForwardDiff):
# N = 100
#                  Time for F:   6.191 μs (1 allocation: 896 bytes)
#   Time for DF (ForwardDiff):   409.391 μs (6043 allocations: 659.25 KiB)
#   Time for DF (ReverseDiff):   451.058 μs (0 allocations: 0 bytes)
#         Time for DF(manual):   17.426 μs (3 allocations: 2.63 KiB)
# N = 400
#                  Time for F:   26.206 μs (1 allocation: 3.25 KiB)
#   Time for DF (ForwardDiff):   6.538 ms (96204 allocations: 10.17 MiB)
#   Time for DF (ReverseDiff):   1.999 ms (0 allocations: 0 bytes)
#         Time for DF(manual):   63.400 μs (3 allocations: 9.75 KiB)
# N = 1600
#                  Time for F:   109.590 μs (1 allocation: 12.63 KiB)
#   Time for DF (ForwardDiff):   155.275 ms (1536804 allocations: 162.25 MiB)
#   Time for DF (ReverseDiff):   15.896 ms (0 allocations: 0 bytes)
#         Time for DF(manual):   260.063 μs (3 allocations: 37.88 KiB)

ext: jl


text:
export OffPolicy, VectorWSARTTrajectory

Base.@kwdef struct OffPolicy{P,B} <: AbstractPolicy
    π_target::P
    π_behavior::B
end

(π::OffPolicy)(env) = π.π_behavior(env)

const VectorWSARTTrajectory = Trajectory{<:NamedTuple{(:weight, SART...)}}

function VectorWSARTTrajectory(;weight=Float64, state=Int, action=Int, reward=Float32, terminal=Bool)
    VectorTrajectory(;weight=Float64, state=state, action=action, reward=reward, terminal=terminal)
end

function RLBase.update!(
    p::OffPolicy,
    t::VectorWSARTTrajectory,
    e::AbstractEnv,
    s::AbstractStage
)
    update!(p.π_target, t, e, s)
end

function RLBase.update!(
    t::VectorWSARTTrajectory,
    p::OffPolicy,
    env::AbstractEnv,
    ::PreActStage,
    a
)
    s = state(env)
    push!(t[:state], s)
    push!(t[:action], a)

    w = prob(p.π_target, env, a) / prob(p.π_behavior, env, a)
    push!(t[:weight], w)
end

function RLBase.update!(
    t::VectorWSARTTrajectory,
    p::Union{OffPolicy{<:QBasedPolicy{<:TDLearner}}, OffPolicy{<:VBasedPolicy{<:MonteCarloLearner}}},
    env::AbstractEnv,
    s::PreEpisodeStage,
)
    empty!(t)
end

function RLBase.update!(
    t::VectorWSARTTrajectory,
    p::Union{OffPolicy{<:QBasedPolicy{<:TDLearner}}, OffPolicy{<:VBasedPolicy{<:MonteCarloLearner}}},
    env::AbstractEnv,
    s::PostEpisodeStage,
)
    action = rand(action_space(env))

    push!(t[:state], state(env))
    push!(t[:action], action)
    push!(t[:weight], 1.0)
end


text:
module Curry

export @curry, curry
export ldump, sldump
export generic_expr_walker, remove_lines!

"""
   @curry <function definition>

Define function methods with multiple argument lists.
The definition
    `@curry f(x)(y) = anything`
is equivalent with
    `f(x) = (y) -> anything`.
Each argument list supports all features of argument lists as keyword arguments
and ellipse.
"""
macro curry(ex)
    curry!(__source__, ex)
end

"""
    curry(expr::Expr[, line=LineNumberNode())

User interface to implementation of Expr transformation for [`@curry`](@ref).
The input expression is not modified.
"""
function curry(ex::Expr; line::LineNumberNode=LineNumberNode(0))
    curry!(line, deepcopy(ex))
end

# implementation.
function curry!(__source__::LineNumberNode, ex::Expr)
    ass = findex(isin([:(=),:function]), ex)
    ass != nothing && length(ass.args) == 2 || error("no function definition")
    cal = ass.args[1] # the lhs of = or f(..) part after function

    # if not :call, unmodified return
    cal isa Expr && cal.head == :call || return finish(__source__, ex)
    fna = cal.args[1] # potentially second :call for multiple arg lists

    # if not a second :call, unmodified return
    fna isa Expr && fna.head == :call || return finish(__source__, ex)

    # extract function arguments
    far = length(cal.args) <= 2 ? cal.args[2] : :(($(cal.args[2:end]...),))
    blk = ass.args[2] # execution block of rhs
    ass.args[1] = fna # put second :call at top
    ass.args[2] = Expr(:block, Expr(:(->), far, blk)) # generate -> expr
    curry!(__source__, ex) # recursively build blocks for further :call s
end

isin(a) = x -> x in a

# add line info into first block and surround with esc.
function finish(__source__::LineNumberNode, ex::Expr)
    blk = findex(isequal(:block), ex)
    blk != nothing || return esc(ex)
    pushfirst!(blk.args, __source__)
    esc(ex)
end

# find expression with symbol as head in expression tree
function findex(pred::Function, ex::Expr)
    pred(ex.head) && return ex
    for ex in ex.args
        r = findex(pred, ex)
        r == nothing || return r
    end
    nothing
end
findex(::Function, ::Any) = nothing

"""
    ldump(:IO, ::Expr[, indent=string, delta=string])

Dump expression in a more readable list form.
Each expression is started with a new line `( head`
Each argument starts with a new line afetr indentation.
The keyword argument `indent` is the initial indentation.
The keyword argument `delta` is the additional indentation after the current.
"""
function ldump(io::IO, ex::Expr; indent="", delta="  ")
    indent = string(indent)
    println(io, indent, '(', ex.head)
    for a in ex.args
        delta = string(delta)
        ldump(io, a, indent=indent*delta, delta=delta)
    end
    println(io, indent, ')')
end

ldump(io::IO, ex::LineNumberNode; indent="", delta="") = nothing
ldump(io::IO, ex::Any; indent="", delta="") = println(io, indent, typeof(ex), ' ', ex)
ldump(ex::Any; indent="", delta="  ") = ldump(stdout, ex, indent=indent, delta=delta)

"""
    sldump(::Expr[, indent=string, delta=string])

    List dump of an expression. See [`ldump`](@ref).
Store output in string.
"""
function sldump(ex::Any; indent="", delta="  ")
    io = IOBuffer()
    ldump(io, ex, indent=indent, delta=delta)
    String(take!(io))
end

"""
    Expression tree walker
"""
function generic_expr_walker(a::Function, ex::Expr)
    i = 1
    j = 0
    while i <= length(ex.args)
        arg = ex.args[i]
        argn = a(arg)
        if argn != nothing
            j += 1
            argn = something(argn)
            if j != i || argn != arg
                ex.args[j] = argn
            end
        end
        i += 1
    end
    resize!(ex.args, j)
    ex
end

identity_walker(ex::Expr) = generic_expr_walker(example_walker, ex)
example_walker(x) = Some(x)

remove_lines!(ex::Expr) = generic_expr_walker(remove_lines!, ex)
remove_lines!(::LineNumberNode) = nothing
remove_lines!(x) = Some(x)

end # module


text:
# This table represents the optimal "power tree"
# based on Knuth's "TAOCP vol 2: Seminumerical Algorithms",
# third edition, section 4.6.3, Fig. 15.  It was
# transcribed into this array form in the gcc source
# code (tree-ssa-math-opts.c), but since this is just
# a table of mathematical facts it should not be copyrightable.
#
# To compute x^q in a minimal number of multiplications
# for 1 ≤ q ≤ 255, you compute x^r * x^s for
# r = power_tree[i] and s = i-q, recursively
# (memoizing each power computation), and we implicitly
# define power_tree[0] = 0.
const power_tree = [
      1,   1,   2,   2,   3,   3,   4,       #   1 -   7
      4,   6,   5,   6,   6,  10,   7,   9,  #   8 -  15
      8,  16,   9,  16,  10,  12,  11,  13,  #  16 -  23
     12,  17,  13,  18,  14,  24,  15,  26,  #  24 -  31
     16,  17,  17,  19,  18,  33,  19,  26,  #  32 -  39
     20,  25,  21,  40,  22,  27,  23,  44,  #  40 -  47
     24,  32,  25,  34,  26,  29,  27,  44,  #  48 -  55
     28,  31,  29,  34,  30,  60,  31,  36,  #  56 -  63
     32,  64,  33,  34,  34,  46,  35,  37,  #  64 -  71
     36,  65,  37,  50,  38,  48,  39,  69,  #  72 -  79
     40,  49,  41,  43,  42,  51,  43,  58,  #  80 -  87
     44,  64,  45,  47,  46,  59,  47,  76,  #  88 -  95
     48,  65,  49,  66,  50,  67,  51,  66,  #  96 - 103
     52,  70,  53,  74,  54, 104,  55,  74,  # 104 - 111
     56,  64,  57,  69,  58,  78,  59,  68,  # 112 - 119
     60,  61,  61,  80,  62,  75,  63,  68,  # 120 - 127
     64,  65,  65, 128,  66, 129,  67,  90,  # 128 - 135
     68,  73,  69, 131,  70,  94,  71,  88,  # 136 - 143
     72, 128,  73,  98,  74, 132,  75, 121,  # 144 - 151
     76, 102,  77, 124,  78, 132,  79, 106,  # 152 - 159
     80,  97,  81, 160,  82,  99,  83, 134,  # 160 - 167
     84,  86,  85,  95,  86, 160,  87, 100,  # 168 - 175
     88, 113,  89,  98,  90, 107,  91, 122,  # 176 - 183
     92, 111,  93, 102,  94, 126,  95, 150,  # 184 - 191
     96, 128,  97, 130,  98, 133,  99, 195,  # 192 - 199
    100, 128, 101, 123, 102, 164, 103, 138,  # 200 - 207
    104, 145, 105, 146, 106, 109, 107, 149,  # 208 - 215
    108, 200, 109, 146, 110, 170, 111, 157,  # 216 - 223
    112, 128, 113, 130, 114, 182, 115, 132,  # 224 - 231
    116, 200, 117, 132, 118, 158, 119, 206,  # 232 - 239
    120, 240, 121, 162, 122, 147, 123, 152,  # 240 - 247
    124, 166, 125, 214, 126, 138, 127, 153,  # 248 - 255
]

# return the inlined/unrolled expression to compute x^p
# by repeated multiplications (using an optimal addition
# chain for |p|<256 and power-by-squaring thereafter.
function addition_chain_pow(x::Symbol, p::Int)
    p == 0 && return :(one($x))
    ex = Union{Expr,Symbol}[] # expressions to compute intermediate powers
    if p < 0
        x′ = gensym()
        push!(ex, :($x′ = inv($x)))
        x = x′
        p = -p
    end
    pows_computed = Dict(1 => x) # powers q => variable storing x^q
    pows_to_compute = [p]
    while !isempty(pows_to_compute)
        q = pop!(pows_to_compute)
        if !haskey(pows_computed, q)
            if q ≤ length(power_tree) # use optimal addition chain
                q1 = power_tree[q]
                q2 = q-q1
                @assert q1 > 0 && q2 > 0
                has1 = haskey(pows_computed, q1)
                has2 = haskey(pows_computed, q2)
                if has1 && has2
                    xq = gensym()
                    push!(ex, :($xq = $(pows_computed[q1]) * $(pows_computed[q2])))
                    pows_computed[q] = xq
                else
                    push!(pows_to_compute, q) # try again after computing q1, q2
                    has1 || push!(pows_to_compute, q1)
                    has2 || push!(pows_to_compute, q2)
                end
            else # q too big, use power-by-squaring algorithm
                q1 = q >> 1
                if haskey(pows_computed, q1)
                    xq = gensym()
                    xq1 = pows_computed[q1]
                    push!(ex, iseven(q) ? :($xq = $xq1 * $xq1) : :($xq = $xq1 * $xq1 * $x))
                    pows_computed[q] = xq
                else
                    push!(pows_to_compute, q) # try again after computing q1
                    push!(pows_to_compute, q1)
                end
            end
        end
    end
    push!(ex, pows_computed[p]) # the variable for the final result x^p
    return Expr(:block, ex...)
end
addition_chain_pow(x::Symbol, p::Integer) = addition_chain_pow(x, Int(p))

function addition_chain_pow(expr, p::Integer)
    x = gensym()
    return Expr(:let, Expr(:(=), x, expr), Expr(:block, addition_chain_pow(x, p)))
end


text:
## SARS-CoV-2 /COVID19
## Affan Shoukat, 2020

using DifferentialEquations, Plots, Parameters, DataFrames, CSV, LinearAlgebra
using StatsPlots, Query, Distributions, Statistics, Random, DelimitedFiles, ProgressMeter

∑(x) = sum(x) # synctactic sugar
heaviside(x) = x <= 0 ? 0 : 1

@with_kw mutable struct ModelParameters
    # default parameter values, (fixed: fixed parameters, sampled: sampled at time of run, input: given as input,i.e. scenarios)
    ## parameters for transmission dynamics.
    β::Float64 = 0.037 # (input) 0.037: R0 ~ 2.5, next generation method
    κ::Float64 = 0.5 # (fixed) reduction in transmission due to mild symptoms.
    σ::Float64 = 1/5.2 # (sampled)incubation period 5.2 days mean, LogNormal
    q::NTuple{4, Float64} = (0.05, 0.05, 0.05, 0.05) # (fixed) proportion of self-quarantine
    h::NTuple{4, Float64} = (0.02075, 0.02140, 0.025, 0.03885) ## (sampled), default values mean of distribution proporition of TOTAL infections going to hospital.
    fₐ::Float64 = 0.05 ##
    fᵢ::Float64 = 0.80 ##
    c::NTuple{4, Float64} = (0.0129, 0.03875, 0.0705, 0.15) ## (sampled), default values mean of distribution, mean is average of Chiense and US CDC.
    τₐ::Float64 = 0.5   ##
    τᵢ::Float64 = 1   ##
    γ::Float64 = 1/4.6 # (fixed) symptom onset to recovery, assumed fixed, based on serial interval... sampling creates a problem negative numbers
    δ::Float64 = 1/3.5 # (sampled), default value mean of distribution, symptom onset to hospitalization
    θ::NTuple{4, Float64} = (0.8, 0.8, 0.4, 0.2) ## (fixed)

    ## recovery and mortality
    mH::Float64 = 0.2296  ## prob of death in hospital
    μH::Float64 = 1/9.7   ## length of hospital stay before death (average of CDC reports)
    ψH::Float64 = 1/10    ## length of hospital stay before recovery (JAMA paper, send PR to CDC GitHub jama paper is median)
    mC::Float64 = 0.1396 ## prob of death in ICU
    μC::Float64 = 1/7     ## length of ICU stay before death Lancet Res "Clinical course and outcomes of critically ill Yang"
    ψC::Float64 = 1/13.25 ## length of ICU before recovery (Chad's ventilation calculation)

    ## internal parameters
    pop::NTuple{4, Float64} = (81982665,129596376,63157200,52431193)

end

function contact_matrix()
    ## contact matrix for general population and in household.
    M = ones(Float64, 4, 4)
    M[1, :] = [9.76190295449833 3.76760818611057 1.50823242972404 0.603972723406184]
    M[2, :] = [3.76760818611057 9.42964041341912 3.05467747979113 0.702942024329397]
    M[3, :] = [1.50823242972404 3.05467747979113 2.95716998311431 0.760587605942274]
    M[4, :] = [0.603972723406184 0.702942024329397 0.760587605942274 1.24948301635913]
    M̃ = ones(Float64, 4, 4)
    M̃[1, :] = [2.03651666746364 1.56405653862347 0.504076940943109 0.380678379808072]
    M̃[2, :] = [1.56405653862347 1.50876370897303 0.445015376096172 0.239252424523331]
    M̃[3, :] = [0.504076940943109 0.445015376096172 1.03746276893954 0.191098354935966]
    M̃[4, :] = [0.380678379808072 0.239252424523331 0.191098354935966 0.641792954732621]
    return M, M̃
end

function Model!(du, u, p, t)
    # model v4, with age structure, feb 20
    # 4 age groups: 0 - 18, 19 - 49, 50 - 65, 65+
    S₁, S₂, S₃, S₄,
    E₁, E₂, E₃, E₄,
    F₁, F₂, F₃, F₄,
    Iₙ₁, Iₙ₂, Iₙ₃, Iₙ₄,
    Qₙ₁, Qₙ₂, Qₙ₃, Qₙ₄,
    Iₕ₁, Iₕ₂, Iₕ₃, Iₕ₄,
    Qₕ₁, Qₕ₂, Qₕ₃, Qₕ₄,
    H₁, H₂, H₃, H₄,
    C₁, C₂, C₃, C₄,
    N₁, N₂, N₃, N₄,

    # internal incidence equations
    Z₁, Z₂, Z₃, Z₄,
    CX₁, CX₂, CX₃, CX₄,
    CY₁, CY₂, CY₃, CY₄,
    DX₁, DX₂, DX₃, DX₄,
    DY₁, DY₂, DY₃, DY₄,

    ## asymptomatic classes
    Aₙ₁, Aₙ₂, Aₙ₃, Aₙ₄,
    Aₛ₁, Aₛ₂, Aₛ₃, Aₛ₄,

    # set up the vectors for syntactic sugar
    S = (S₁, S₂, S₃, S₄)
    E = (E₁, E₂, E₃, E₄)
    F = (F₁, F₂, F₃, F₄)
    Iₙ = (Iₙ₁, Iₙ₂, Iₙ₃, Iₙ₄)
    Qₙ = (Qₙ₁, Qₙ₂, Qₙ₃, Qₙ₄)
    Iₕ = (Iₕ₁, Iₕ₂, Iₕ₃, Iₕ₄)
    Qₕ = (Qₕ₁, Qₕ₂, Qₕ₃, Qₕ₄)
    Aₙ = (Aₙ₁, Aₙ₂, Aₙ₃, Aₙ₄)
    Aₛ = (Aₛ₁, Aₛ₂, Aₛ₃, Aₛ₄)
    H = (H₁, H₂, H₃, H₄)
    C = (C₁, C₂, C₃, C₄)
    N = (N₁, N₂, N₃, N₄)

    # internal incidence euqations
    Z = (Z₁, Z₂, Z₃, Z₄)
    CX = (CX₁, CX₂, CX₃, CX₄)
    CY = (CY₁, CY₂, CY₃, CY₄)
    DX = (DX₁, DX₂, DX₃, DX₄)
    DY = (DY₁, DY₂, DY₃, DY₄)

    #get the contact matrix
    M, M̃ = contact_matrix()

    # constants
    pop = p.pop  ## fixed population
    Nᵥ = p.strat

    ## working in new model

    @unpack β, κ, ξ, ν, σ, q, h, fₐ, fᵢ, τₐ, τᵢ, γ, δ, θ, c, mH, μH, ψH, mC, μC, ψC = p
    for a = 1:4
        # sus S
        du[a+0] = -β*S[a]*(dot(M[a, :], Iₙ./pop) + dot(M[a, :], Iₕ./pop) + κ*dot(M[a, :], Aₙ./pop)  -
                   β*S[a]*(dot(M̃[a, :], Qₙ./pop) + dot(M̃[a, :], Qₕ./pop) + κ*dot(M̃[a, :], Aₛ./pop)
        # exposed E
        du[a+4]  = β*S[a]*(dot(M[a, :], Iₙ./pop) + dot(M[a, :], Iₕ./pop) + κ*dot(M[a, :], Aₙ./pop) +
                   β*S[a]*(dot(M̃[a, :], Qₙ./pop) + dot(M̃[a, :], Qₕ./pop) + κ*dot(M̃[a, :], Aₛ./pop) -
                   σ*E[a]
        # In class
        du[a+8] = (1 - θ[a])*(1 - q[a])*(1 - h[a])*σ*E[a]  - (1 - fᵢ)*γ*Iₙ[a] - fᵢ*τᵢ*Iₙ[a]
        # Qn class
        du[a+12] = (1 - θ[a])*q[a]*(1 - h[a])*σ*E[a] + fᵢ*τᵢ*Iₙ[a] - γ*Qₙ[a]
        # Ih class
        du[a+16] = (1 - θ[a])*(1 - q[a])*h[a]*σ*E[a] - (1 - fᵢ)*δ*Iₕ[a] - fᵢ*τᵢ*Iₕ[a]
        # Qh class
        du[a+20] = (1 - θ[a])*q[a]*h[a]*σ*E[a] + fᵢ*τᵢ*Iₕ[a] - δ*Qₕ[a]
        # Ha class
        du[a+24] = (1 - c[a])*(1 - fᵢ)*δ*Iₕ[a] + (1 - c[a])*δ*Qₕ[a] - (mH*μH + (1 - mH)*ψH)*H[a]
        # Ca class
        du[a+28] = c[a]*(1 - fᵢ)*δ*Iₕ[a] + c[a]*δ*Qₕ[a] - (mC*μC + (1 - mC)*ψC)*C[a]
        # Na class
        du[a+32] = -mC*μC*C[a] - mH*μH*H[a]

        # Z class to collect cumulative incidence (only from suspectibles)
        du[a+36] =  β*S[a]*(dot(M[a, :], Iₙ./pop) + dot(M[a, :], Iₕ./pop) + κ*dot(M[a, :], Aₙ./pop)  +
                    β*S[a]*(dot(M̃[a, :], Qₙ./pop) + dot(M̃[a, :], Qₕ./pop) + κ*dot(M̃[a, :], Aₛ./pop)
        ## collect cumulative incidence of hospitalization and ICU class 69 to 84
        du[a+40] = (1 - c[a])*(1 - fᵢ)*δ*Iₕ[a] + (1 - c[a])*δ*Qₕ[a]
        du[a+44] = c[a]*(1 - fᵢ)*δ*Iₕ[a] + c[a]*δ*Qₕ[a]
        du[a+48] = (mH*μH)*H[a]  ## DX
        du[a+52] = (mC*μC)*C[a]  ## DY

        ## Asymptomatic classes
        # Aₙ₁, Aₙ₂, Aₙ₃, Aₙ₄,
        du[a+56] = θ[a]*σ*(F[a] + E[a]) - (1 - fₐ)*γ*Aₙ[a] - fₐ*τₐ*Aₙ[a]
        # Aₛ₁, Aₛ₂, Aₛ₃, Aₛ₄,
        du[a+60] = fₐ*τₐ*Aₙ[a] - γ*Aₛ[a]
    end
end

function run_model(p::ModelParameters, nsims=1)
    ## set up ODE time and initial conditions
    tspan = (0.0, 1000.0)
    u0 = zeros(Float64, 64) ## total compartments
    sols = []
    for mc = 1:nsims
        ## reset the initial conditions
        u0[1] = u0[33] = p.pop[1]
        u0[2] = u0[34] = p.pop[2]
        u0[3] = u0[35] = p.pop[3]
        u0[4] = u0[36] = p.pop[4]
        # initial infected person
        u0[13] = 1

        # sample the parameters needed per simulation
        p.δ = 1/(rand(Uniform(2, 5)))
        p.σ = 1/(rand(LogNormal(log(5.2), 0.1)))
        p.h = (rand(Uniform(0.0085, 0.033)), rand(Uniform(0.0088, 0.034)), rand(Uniform(0.01, 0.042)), rand(Uniform(0.017, 0.066))) ./ (1 .- p.θ)
        p.c = (rand(Uniform(0.013, 0.015)), rand(Uniform(0.04, 0.044)), rand(Uniform(0.05, 0.1)), rand(Uniform(0.10, 0.20))) ## c dosnt change

        ## solve the ODE model
        #print("...sim: $mc, params: $(p.τ), $(p.f), $(p.nva), $(p.ν) \r")
        prob = ODEProblem(Model!, u0, tspan, p)
        sol = solve(prob, Rodas4(autodiff=false), dt=1.0, adaptive=false, callback=cb)  ## WORKS
        #sol = solve(prob, Rosenbrock23(autodiff=false),  dt=0.1, adaptive=false, callback=cbset)  ## WORKS
        push!(sols, sol)
    end
    #println("\n simulation scenario finished")
    return sols
end

function run_single()
    # A function that runs a single simulation for testing/calibrating purposes.
    ## setup parameters (overwrite some default ones if needed)
    p = ModelParameters()
    p.β = getβ("2.0")
    p.fₐ = 0.20
    p.τₐ = 1
    p.fᵢ = 0.80
    p.τᵢ = 1
    dump(p)
    sol = run_model(p, 1)[1] ## the [1] is required since runsims returns an array of solutions
    return sol
end

## hospital capacity: https://hcup-us.ahrq.gov/reports/statbriefs/sb185-Hospital-Intensive-Care-Units-2011.jsp
## https://www.aha.org/system/files/media/file/2020/01/2020-aha-hospital-fast-facts-new-Jan-2020.pdf
## https://www.sccm.org/Communications/Critical-Care-Statistics


text:
using BinDeps

@BinDeps.setup

libtweetnacl = library_dependency("tweetnacl", aliases=["libtweetnacl"])

const libtweetnacl_dir = joinpath(BinDeps.depsdir(libtweetnacl), "src")

provides(SimpleBuild,
         (@build_steps begin
              @build_steps begin
                  ChangeDirectory(libtweetnacl_dir)
                  `make`
                  `mv libtweetnacl.so ../usr/lib`
              end

          end), libtweetnacl, os=:Unix)

@BinDeps.install Dict([(:tweetnacl, :tweetnacl)])

max_stars_repo_name: LachlanGunn/TweetNaCl

max_issues_repo_name: LachlanGunn/TweetNaCl

max_forks_repo_name: LachlanGunn/TweetNaCl


text:
# Test factor analysis implementation using a test similar to the one in
# scikit-learn, so we can directly compare implementations if needed.

using MultivariateStats

n_samples = 20
n_features = 5
n_components = 3

X = readdlm("X_fa.csv", ',')'
fa, L = fit(FactorAnalysis, X, n_components)

X_t = transform(fa, X)

X_rec = fa.Vbar .+ fa.F*X_t .+ fa.Ψ.*randn(n_features,n_samples)


text:
module TubeScreamer

using ACME: ACME, DiscreteModel, @circuit, capacitor, diode, opamp, resistor, run!,
    voltagesource, voltageprobe
using Interpolations: BSpline, Cubic, Line, OnGrid, interpolate, scale
using NumericalIntegration: cumul_integrate
using StaticArrays: SVector

export runsim, runsimaaout, runsimaa

# circuit description
const circ = @circuit begin
    # power supply
    j3 = voltagesource(9), [+] ⟷ vcc, [-] ⟷ gnd # 9V power supply

    # simplified input stage
    j1 = voltagesource(), [-] ⟷ gnd # input
    c2 = capacitor(18e-9), [1] ⟷ j1[+]
    r5 = resistor(100e3), [1] ⟷ c2[2], [2] ⟷ gnd

    # distortion stage
    ic1a = opamp(), ["in+"] ⟷ c2[2], ["out-"] ⟷ gnd
    d1 = diode(is=4e-9, η=2), [-] ⟷ ic1a["out+"], [+] ⟷ ic1a["in-"] # 1N914
    d2 = diode(is=3e-9, η=2), [-] ⟷ ic1a["in-"] # 1N914
    d3 = diode(is=5e-9, η=2), [+] ⟷ ic1a["out+"], [-] ⟷ d2[+] # 1N914
    r6 = resistor(500e3), [1] ⟷ ic1a["in-"], [2] ⟷ ic1a["out+"]
    c99 = capacitor(51e-12), [1] ⟷ ic1a["in-"], [2] ⟷ ic1a["out+"]
    c4 = capacitor(47e-9), [1] ⟷ ic1a["in-"]
    r7 = resistor(4.7e3), [1] ⟷ c4[2], [2] ⟷ gnd

    # output stage
    c10 = capacitor(1e-6), [1] ⟷ ic1a["out+"]
    r16 = resistor(100e3), [1] ⟷ c10[2], [2] ⟷ gnd
    j2 = voltageprobe(), [+] ⟷ c10[2], [-] ⟷ gnd # output
end

# precompute antiderivate table
function makeantiderivtable(model)
    pmax = 1e-3 # guesstimate of relevant range
    ps = range(-pmax, pmax, length=1024)
    zs = [SVector{ACME.nn(model,1)}(ACME.solve(model.solvers[1], [p])) for p in ps]
    zsint = cumul_integrate(ps, zs)
    return scale(interpolate(zsint, BSpline(Cubic(Line(OnGrid())))), ps)
end

# simulation without any antialiasing
function runsim(fs, u)
    model = DiscreteModel(circ, 1//fs)
    permutedims(run!(model, permutedims(u)))
end

# simulation with antialiasing of the output only
function runsimaaout(fs, u)
    model = DiscreteModel(circ, 1//fs)
    y = similar(u)
    adtab = makeantiderivtable(model)
    last_x = copy(model.x)
    last_p = 0.0
    last_u = 0.0
    last_ad = adtab(last_p)
    for n in eachindex(u)
        p = (model.dqs[1]*model.x + model.eqs[1]*[u[n]])[1]
        ad = adtab(p)
        if abs(p-last_p) > 1e-15
            z = (ad - last_ad) / (p - last_p)
        else
            z = ACME.solve(model.solvers[1], [0.5*(p+last_p)])
        end
        y[n] = (0.5*model.dy*(model.x+last_x) + 0.5*model.ey*[u[n]+last_u] + model.fy*z + model.y0)[1]
        last_p = p
        last_ad = ad
        last_u = u[n]
        copyto!(last_x, model.x)
        copyto!(model.x, model.a * model.x + model.b * [u[n]] + model.c * ACME.solve(model.solvers[1], [p]) + model.x0)
    end
    return y
end

# simulation with full antialiasing
function runsimaa(fs, u)
    fs = fs/1.5
    model = DiscreteModel(circ, 1/fs)
    y = similar(u)
    adtab = makeantiderivtable(model)
    last_x = copy(model.x)
    old_x = copy(model.x)
    last_p = 0.0
    last_u = 0.0
    last_ad = adtab(last_p)
    for n in eachindex(u)
        p = (model.dqs[1]*model.x + model.eqs[1]*[u[n]])[1]
        ad = adtab(p)
        if abs(p-last_p) > 1e-15
            z = (ad - last_ad) / (p - last_p)
        else
            z = ACME.solve(model.solvers[1], [0.5*(p+last_p)])
        end
        y[n] = (0.5*model.dy*(model.x+last_x) + 0.5*model.ey*[u[n]+last_u] + model.fy*z + model.y0)[1]
        copyto!(old_x, model.x)
        copyto!(model.x, 0.5*model.a*(model.x+last_x) + 0.5*model.b*[u[n]+last_u] + model.c*z + model.x0)
        copyto!(last_x, old_x)
        last_p = p
        last_ad = ad
        last_u = u[n]
    end
    return y
end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the nutation and EO computation of IAU-2006.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Wallace, P. T., Capitaine, N (2006). Precession-nutation procedures
#       consistent with IAU 2006 resolutions. Astronomy & Astrophysics.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/iau2006/nutation_eo.jl
# ==================================================

# Function nutation_eo_iau2006
# ----------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Numerical example at Appendix of [1]
#
# According to this example, using:
#
#   UTC    = 2006 January 15, 21h 24.375
#   JD_TT  = 2400000.5+53750.892855138888889
#   JD_UT1 = 2400000.5+53750.892104561342593
#
# one gets:
#
#   mϵ_2000 = +84378.576696215 arcsec (called ϵ_A in [1])
#   Δϵ_2000 =   +8.656841020   arcsec (called Δϵ in [1])
#   ΔΨ_2000 =   –1.071332969   arcsec (called ΔΨ in [1])
#   EO      = -277.646996035   arcsec
#
################################################################################

@testset "Function nutation_eo_iau2006" begin
    JD_TT  = 2400000.5+53750.892855138888889
    JD_UT1 = 2400000.5+53750.892104561342593

    r2a = 180/π*3600

    mϵ_2000, Δϵ_2000, ΔΨ_2000, EO = nutation_eo_iau2006(JD_TT)

    @test mϵ_2000*r2a ≈ 84378.576696215 atol = 1e-9
    @test Δϵ_2000*r2a ≈     8.656841020 atol = 1e-7
    @test ΔΨ_2000*r2a ≈    -1.071332969 atol = 1e-7
    @test EO*r2a      ≈  -277.646996035 atol = 1e-6
end


text:
module day16

include("inputs.jl")

ex = inputs.ex_d16()
ex2 = inputs.ex2_d16()
input = inputs.input_d16()

function parse_input(input)
    lines = split(input, '\n', keepempty=false)

    fields = Dict{String, Array{Int}}()
    rgx = r"^([a-zA-Z]+\s*[a-zA-Z]+): (\d+)-(\d+) or (\d+)-(\d+)$"
    for line in lines
        if contains(line, ": ")
            m = match(rgx, line).captures
            r = parse.(Int, m[2:5])
            rng = [r[1]:r[2]..., r[3]:r[4]...]
            fields[m[1]] = rng
        end
    end

    your_ticket = findfirst(line->line=="your ticket:", lines) + 1 |>
                    x->getindex(lines, x) |>
                    x->split(x, ",") |>
                    x->parse.(Int, x)

    ntline = findfirst(line->line=="nearby tickets:", lines) + 1
    nearby_tickets = []
    for line in lines[ntline:end]
        push!(nearby_tickets, split(line, ",") |> x->parse.(Int, x))
    end

    return fields, your_ticket, nearby_tickets
end

# Part 1
function find_invalid_items(ticket, fields)
    # index of ticket entries that *are* found somewhere in fields' values:
    idx = [any([item ∈ value for value in values(fields)]) for item in ticket]
    return ticket[.!idx]
end

function part1(input)
    fields, _, nearby_tickets = parse_input(input)
    sum(sum.([find_invalid_items(ticket, fields) for ticket in nearby_tickets]))
end

@assert part1(input) == 21081
@info "Day 16, Part 1 answer: $(part1(input))"

# Part 2
isvalid(ticket, fields) = all([any([item ∈ value for value in values(fields)]) for item in ticket])

function part2(input)
    """Get the valid tickets"""
    f, yt, nts = parse_input(input)
    valid_nearbys = nts[map(x->isvalid(x, f), nts)]
    tickets = permutedims(hcat(valid_nearbys...))

    """Find the possible associations between fields and ticket positions."""
    assoc = Dict{String, Array{Int}}()
    for (field, range) in f # iterate over range of fields
        for (i, position) in enumerate(eachcol(tickets))
            """Check if each "pos" (a column in the tickets) is entirely contained
            in the range of a field. If so, that pos could be associated with that
            field."""

            if Set(intersect(position, range)) == Set(position)
                # println("Field ", field, " could be assoc. with pos ", i)
                haskey(assoc, field) ? push!(assoc[field], i) : assoc[field] = [i]
            end
        end
    end

    """Identify which pairs have only one valid assocation, remove those values from the other pairs and repeat."""
    while !all(length.(values(assoc)) .== 1)  # while all the values have >1 one el
        for (field, positions) in assoc
            if length(assoc[field]) == 1
                # println("Associate ", field, " with ", assoc[field])
                for (f, p) in assoc
                    if f != field
                        assoc[f] = filter(x->x!=assoc[field][1], assoc[f])
                    end
                end
            end
        end
    end

    """Puzzle answer asks to multiply your ticket where the identified positions starts with departure."""
    prod = 1
    for field in keys(assoc)
        if startswith(field, "departure")
            prod *= yt[assoc[field][1]]
        end
    end

    return prod
end

@assert part2(input) == 314360510573
@info "Day 16, Part 2 answer: $(part2(input))"

end


text:
export Problem

"""
    Problem( model, initial, param ; verbose=true)
or
    Problem( model, initial, param, solver ; verbose=true)

Builds an initial-value problem which can then be solved (integrated in time) through `solve!( problem )`

# Arguments
- `model   :: AbstractModel`,  the system of equation solved.
May be built, e.g., by `WaterWaves(param)`;
- `initial :: InitialData`, the initial data.
May be buit, e.g., by `Init(η,v)` where
`η` is the surface deformation and `v` the derivative of the trace of the velocity potential at the surface;
- `param   :: NamedTuple`, must contain values for
    - `N`, the number of collocation points of the spatial grid
    - `L`, the half-length of the spatial grid
    - `T`, the final time of integration
    - `dt`, the timestep
    - additionally, it may contain `Ns` the number of computed data or `ns` for storing data every `ns` computation steps (by default, every computed data is stored).

- `solver  :: TimeSolver`, the solver for time integration (optional, default is explicit Runge-Kutta fourth order solver).
May be built, e.g., by `RK4(model)` or `RK4_naive()`.

Information are not printed if keyword `verbose = false` (default is `true`).

"""
mutable struct Problem

    model   :: AbstractModel
    initial :: InitialData
    param   :: NamedTuple
    solver  :: TimeSolver
    times   :: Times
    mesh    :: Mesh
    data    :: Data

    function Problem(model   :: AbstractModel,
                     initial :: InitialData,
                     param   :: NamedTuple;
                     solver = RK4(model)  :: TimeSolver,
                     verbose = true :: Bool)

        if verbose == true
            @info string("\nBuild the initial-value problem for model ",model.label,"\n",
                         "with parameters\n",param)
        end
        if in(:Ns,keys(param))
            times = Times(param.dt, param.T; Ns = param.Ns)
        elseif in(:ns,keys(param))
            times = Times(param.dt, param.T; ns = param.ns)
        else
            times = Times(param.dt, param.T)
        end

        mesh  = Mesh(param)
        data  = Data(model.mapto(initial))

        # A basic check
        try
            step!(solver, model, copy(last(data.U)), 1)
        catch
            @warn "The model and the solver are incompatible. solve! will not work."
        end

        new(model, initial, param, solver, times, mesh, data)

    end

end

export solve!

"""
    solve!( problem; verbose=true )

Solves (i.e. integrates in time) an initial-value problem

The argument `problem` should be of type `Problem`.
It may be buit, e.g., by `Problem(model, initial, param)`

Information are not printed if keyword `verbose = false` (default is `true`).

"""
function solve!(problem :: Problem; verbose=true::Bool)

    ci = get(ENV, "CI", nothing) == "true"

    if verbose == true
        @info string("\nNow solving the initial-value problem for model ",problem.model.label,"\n",
            "with parameters\n",problem.param)
    end

    U = copy(last(problem.data.U))

    dt     = problem.times.dt
    solver = problem.solver
    model  = problem.model
    data   = problem.data.U

    if problem.times.tc == problem.times.ts
        pbar = Progress(problem.times.Ns-1; enabled = !ci)
        for j in 1:problem.times.Ns-1
            step!(solver, model, U, dt)
            push!(data,copy(U))
            next!(pbar)
        end

    elseif length(problem.times.ts) > 25
        pbar = Progress(problem.times.Ns-1; enabled = !ci)
        for j in 1:problem.times.Ns-1
            for l in 1:problem.times.ns[j]
                step!(solver, model, U, dt)
            end
            push!(data,copy(U))
            next!(pbar)
        end
    else
        for j in 1:problem.times.Ns-1
            pbar = Progress(problem.times.ns[j], desc=string("Step ",j,"/",problem.times.Ns-1,"...") ; enabled = !ci)
            for l in 1:problem.times.ns[j]
                step!(solver, model, U, dt)
                next!(pbar)
            end
            push!(data,copy(U))
            println()

        end
    end

    println()

end

using Base.Threads

"""
    solve!( problems; verbose=true )

Solves (i.e. integrates in time) a collection of initial-value problems.

The argument `problems` should be a collection (list, array...) of elements of type `Problem`.

"""
function solve!(problems; verbose=true::Bool)
    U=[];nsteps=0
    for problem in problems
        push!(U,copy(last(problem.data.U)))
        nsteps+=problem.times.Ns-1
    end
    pg=Progress(nsteps;dt=1)
    @threads for i in 1:length(problems)

        if problems[i].times.ns == 1

            for j in 1:problems[i].times.Ns-1
                step!(problems[i].solver, problems[i].model, U[i], problems[i].times.dt)
                push!(problems[i].data.U,copy(U[i]))
                next!(pg)
            end

        else
            for j in 1:problems[i].times.Ns-1
                for l in 1:problems[i].times.ns[j]
                    step!(problems[i].solver, problems[i].model, U[i], problems[i].times.dt)
                end
                push!(problems[i].data.U,copy(U[i]))
                next!(pg)
            end
        end
        if verbose == true
            @info string("\nDone solving the model ",problems[i].model.label,"\n",
                "with parameters\n",problems[i].param)
        end

    end

    println()

end


text:
# ## Read Data

# you can read data using vega dataset

using VegaDatasets
data = dataset("iris")
print(data)

## Plot data

# using graph, it become much easier to find pattern from data.

# You can plot using VegaLite.
# In Julia, `:` means symbolifying such as `:sepalLength` which is not variable, it is a name of variable.(I am sure but it seems to work like a pointer.)

using VegaLite
p = data |> @vlplot(:point, x=:sepalLength, y=:sepalWidth, color=:species)

save(joinpath(@OUTPUT, "plot.svg"), p)  # hide

# \fig{plot}

# Plotting pairs is much easier because you can see all data at once.

p = data |>
    @vlplot(
        repeat={
            row=[:sepalLength,:sepalWidth,:petalLength, :petalWidth],
            column=[:sepalLength,:sepalWidth,:petalLength, :petalWidth]
        }
    ) +
    @vlplot(
        :point,
        x={field={repeat=:column}, type=:quantitative},
        y={field={repeat=:row}, type=:quantitative},
        color={field=:species},
    )

save(joinpath(@OUTPUT, "plot2.svg"), p)  # hide

# \fig{plot2}

# From the graph above, you can see 'Setosa' is pretty easy to identify but 'versicolor' and 'virginica' are overlapped.

# You can do all the process much easier with "DataVoyager".


text:
# @generated function Base.zero{U<:AbstractUnum}(::Type{U})
#   c = unumConstants(U)
#   :($(c.zero))
# end

# # convert any floating point number to a unum
# # lets worry about the base-2 case for now, then generalize later:
# @generated function call{E,UINT, FLOAT<:AbstractFloat}(::Type{Unum{E,UINT}}, x::FLOAT)
#   B = 2
#   c = unumConstants(B, E, UINT)
#   f = FloatInfo(FLOAT)

#   # this is the actual conversion function:
#   quote
#     println($(f.uintType))

#     ival = reinterpret($(f.uintType), x)
#     exponent = (ival & $(f.emask)) >> $(f.fpos)
#     fraction = (ival & $(f.fmask))
#     sign = ival >> $(f.nbits-1)

#     if exponent == 0
#       if fraction == 0
#         # zero (exact)
#         return _i2u($B,$E,$(c.zero))
#       end
#     end

#     # for x in (ival, exponent, fraction, sign)
#     #   println(bits(x))
#     # end

#   end
# end

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

type FloatInfo{F<:AbstractFloat, I<:Integer}
  nbits::Int
  esize::Int
  fsize::Int
  emask::I
  fmask::I
  smask::I
  UINT::DataType
  INT::DataType
  FloatInfo() = new()
end

ispositive(x::AbstractFloat) = x >= 0.
isnegative(x::AbstractFloat) = x < 0.

const _floatinfo = Dict(
    Float16 => (UInt16, Int16, 5,  10),
    Float32 => (UInt32, Int32, 8,  23),
    Float64 => (UInt64, Int64, 11, 52),
  )

# keep a cache for given parameter sets so we don't keep rebuilding the constants
const floatConstCache = Dict{DataType, FloatInfo}()

# expect this to be called from a generated function, so we're being passed type params
function floatConstants(F::DataType)
  get!(floatConstCache, F) do
    UINT, INT, ES, FS = _floatinfo[F]
    info = FloatInfo{F,INT}()
    info.nbits = numbits(F)
    info.esize = ES
    info.fsize = FS
    info.fmask = createmask(INT, FS,    FS)
    info.emask = createmask(INT, FS+ES, ES)
    info.smask = createmask(INT, info.nbits)
    info.UINT = UINT
    info.INT = INT
    info
  end
end

function Base.show{F,INT}(io::IO, info::FloatInfo{F,INT})
  println("FloatInfo{$F,$INT}:")
  for fn in fieldnames(info)[1:3]
    println(@sprintf("  %8s %6d", fn, getfield(info, fn)))
  end
  for fn in fieldnames(info)[4:end-2]
    println(@sprintf("  %8s %s", fn, bits(getfield(info, fn))))
  end
end

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

# NOTE: these are not ideal conversions... just want to get a ballpark answer
function Base.float(base, esize, fsize, e, f)
  if e == 0
    return base ^ (e - 2^(esize-1)) * f / (2^fsize)
  else
    return base ^ (e - 2^(esize-1) - 1) * (1 + (base-1) * f / (2^fsize))
  end
end
Base.float{B,ESS,FSS}(u::AbstractUnum{B,ESS,FSS}) = (isnegative(u) ? -1 : 1) * float(Float64(B), esize(u), fsize(u), exponent(u), significand(u))

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

@generated function Base.convert{U<:AbstractUnum, F<:AbstractFloat}(::Type{F}, u::U)
end

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

# this is missing?
# Base.issubnormal(x::Float16) = (x != zero(Float16) && reinterpret(UInt16,x) & 0x7c00) == zero(UInt16))
# Base.issubnormal(x::Float16) = (abs(x) < $(box(Float32,unbox(UInt32,0x00800000)))) & (x!=0)
@eval Base.issubnormal(x::Float16) = (abs(x) < $(box(Float16, unbox(UInt16, 0x0400)))) & (x != zero(Float16))

# convert a float to a unum
# credit to @ityonemo... I used his function as reference
@generated function Base.convert{U<:AbstractUnum, F<:AbstractFloat}(::Type{U}, x::F)
  c = unumConstants(U)
  fc = floatConstants(F)

  quote
    # #some checks for special values
    # (isnan(x)) && return nan(Unum{ESS,FSS})
    # (isinf(x)) && return ((x < 0) ? neg_inf(Unum{ESS,FSS}) : pos_inf(Unum{ESS,FSS}))

    # nan, inf or zero
    isnan(x) && return $(c.nan)
    isinf(x) && return x<0.0 ? $(c.neginf) : $(c.posinf)

    # #convert the floating point x to its integer equivalent
    # I = fp.intequiv                 #the integer type of the same width
    # _esize = fp.esize               #how many bits in the exponent
    # _fsize = fp.fsize               #how many bits in the fraction
    # _bits = _esize + _fsize + 1     #how many total bits
    # _ebias = 1 << (_esize - 1) - 1   #exponent bias (= _emax)
    # _emin = -(_ebias) + 1           #minimum exponent

    # exponent bias and minimum
    # ebias = 1 << ($(fc.fsize) - 1) - 1
    ebias = $(createmask(c.INT, fc.fsize-1, fc.fsize-1))
    # emin = -ebias + 1
    # TODO: emin can be hardcoded

    # ibits = uint64(reinterpret(I, x)[1])

    # get the raw bits
    ibits = reinterpret($(c.INT), x)

    # fraction = ibits & mask(_fsize) << (64 - _fsize)
    # #make some changes to the data for subnormal numbers.
    # (x == 0) && return zero(Unum{ESS,FSS})
    x == 0.0 && return $(c.zero)

    # #grab the sign
    # flags = (ibits & (one(I) << (_esize + _fsize))) != 0 ? UNUM_SIGN_MASK : z16
    u_signmask = x > 0 ? $(c.signbitmask) : $(c.zero)

    # #grab the exponent part
    # biased_exp::Int16 = ibits & mask(_fsize:(_fsize + _esize - 1)) >> _fsize

    # get the float fraction
    x_fraction = ibits & $(fc.fmask)

    # #generate the unbiased exponent and remember to take frac_move into account.
    # unbiased_exp::Int16 = biased_exp - _ebias + ((biased_exp == 0) ? 1 : 0)

    # get the biased and unbiased exponent
    x_exponent_biased = (ibits & $(fc.emask)) >> fc.fsize
    x_exponent_unbiased = x_exponent_biased - ebias + (biased_exp == 0 ? 1 : 0)

    # if issubnormal(x)
    #   #keeping in mind that the fraction bits are now left-aligned, calculate
    #   #how much further we have to push the fraction bits.
    #   frac_move::Int16 = clz(fraction) + 1
    #   fraction = fraction << frac_move
    #   unbiased_exp -= frac_move
    # end

    if issubnormal(x)
      # he shifts the fraction to far left, then counts leading zeros
      # TODO: i think there's a lot of busy work here that's unnecessary... maybe start from scratch...
    end

    # #grab the fraction part

    # #check to see if the exponent is too low.
    # if (unbiased_exp < min_exponent(ESS))
    #   #right shift the fraction by the requisite amount.
    #   shift = min_exponent(ESS) - unbiased_exp
    #   #make sure we don't have any bits in the shifted segment.
    #   #first, are there more bits in the shift than the width?
    #   if (shift > 64)
    #     ((fraction != 0) || ((shift > 65) && (unbiased_exp != 0))) && (flags |= UNUM_UBIT_MASK)
    #   else
    #     ((fraction & mask(shift)) == 0) || (flags |= UNUM_UBIT_MASK)
    #   end
    #   #shift fraction by the amount.
    #   fraction = fraction >> shift
    #   #punch in the one
    #   fraction |= ((biased_exp == 0) ? 0 : t64 >> (shift - 1))
    #   #set to subnormal settings.
    #   esize = uint16((1 << ESS) - 1)
    #   exponent = z64
    # elseif (unbiased_exp > max_exponent(ESS))
    #   return mmr(Unum{ESS,FSS}, flags & UNUM_SIGN_MASK)
    # else
    #   (esize, exponent) = encode_exp(unbiased_exp)
    # end

    # #for really large FSS fractions pad some zeroes in front.
    # (__frac_cells(FSS) > 1) && (fraction = [zeros(Uint64, __frac_cells(FSS) - 1),fraction])

    # r = unum(Unum{ESS,FSS}, min(_fsize, max_fsize(FSS)), esize, flags, fraction, exponent)
    # #check for the "infinity hack" where we "accidentally" create inf.
    # is_inf(r) ? mmr(Unum{ESS,FSS}, flags & UNUM_SIGN_MASK) : r
  end
end

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

# to extract a value:

# # this stores some key sizes and masks for doing float conversions
# type FloatInfo{F<:AbstractFloat, U<:AbstractUnum}
#   nbits::Int
#   epos::Int
#   fpos::Int
#   esize::Int
#   fsize::Int
#   signbitmask::U
#   emask::U
#   fmask::U
#   fType::DataType
#   uType::DataType
# end

# function Base.show{F,U}(io::IO, info::FloatInfo{F,U})
#   println("FloatInfo{$F,$U}:")
#   for fn in fieldnames(info)[1:5]
#     println(@sprintf("  %15s %6d", fn, getfield(info, fn)))
#   end
#   for fn in fieldnames(info)[6:8]
#     println(@sprintf("  %15s %s", fn, bits(getfield(info, fn))))
#   end
# end

# FloatInfo(::Type{Float16}) = FloatInfo(Float16, UInt16, 16, 5)
# FloatInfo(::Type{Float32}) = FloatInfo(Float32, UInt32, 32, 8)
# FloatInfo(::Type{Float64}) = FloatInfo(Float64, UInt64, 64, 11)
# # FloatInfo(::Type{Float128}) = FloatInfo(Float128, UInt128, 128, 15, 113)

# function FloatInfo{F,U}(::Type{F}, ::Type{U}, nbits::Int, esize::Int)
#   fsize = nbits - esize - 1
#   FloatInfo{F,U}(nbits,
#                     nbits-1,
#                     fsize,
#                     esize,
#                     fsize,
#                     createmask(U, nbits, 1),
#                     createmask(U, nbits-1, esize),
#                     createmask(U, fsize, fsize),
#                     F,
#                     U)
# end


text:
"""
$(TYPEDSIGNATURES)

Create Grid from Triangle input data.

See [`default_options`](@ref) for available `kwargs`.
"""
function ExtendableGrids.simplexgrid(::Type{TriangulateType},Triangulate, input; kwargs...)

    opts=blendoptions!(default_options();kwargs...)

    flags=makeflags(opts,:triangle)

    if opts[:verbose]
        @show flags
    end
    if !isnothing(opts[:unsuitable])
        Triangulate.triunsuitable(opts[:unsuitable])
    end

    triout,vorout=Triangulate.triangulate(flags,input)

    pointlist=triout.pointlist

    trianglelist=triout.trianglelist

    if size(triout.triangleattributelist,2)==0
        # Add default for cellregions if that was not created
        cellregions=ones(Int32,size(trianglelist,2))
    else
        cellregions=Vector{Int32}(vec(triout.triangleattributelist))
    end

    segmentlist=triout.segmentlist

    segmentmarkerlist=triout.segmentmarkerlist

    ExtendableGrids.simplexgrid(pointlist,trianglelist,cellregions,segmentlist,segmentmarkerlist)
end

"""
$(TYPEDSIGNATURES)
Create a TriangulateIO structure
from a number of input arrays.
The 2D input arrays are transposed if necessary and converted to
the proper data types for Triangulate.

This conversion is not performed if the data types are those
indicated in the defaults and the leading dimension of 2D arrays
corresponds to the space dimension.

"""
function triangulateio(Triangulate;
                       points=Array{Cdouble,2}(undef,0,0),
                       bfaces=Array{Cint,2}(undef,0,0),
                       bfaceregions=Array{Cint,1}(undef,0),
                       regionpoints=Array{Cdouble,2}(undef,0,0),
                       regionnumbers=Array{Cint,1}(undef,0),
                       regionvolumes=Array{Cdouble,1}(undef,0)
                       )

    ndims(points)==2 || throw(DimensionMismatch("Wrong space dimension"))
    if size(points,2)==2
        points=transpose(points)
    end
    if typeof(points)!=Array{Cdouble,2}
        points=Array{Cdouble,2}(points)
    end
    size(points,2)>2  || throw(ErrorException("Expected more than 2 input points"))

    @assert ndims(bfaces)==2
    if size(bfaces,2)==2
        bfaces=transpose(bfaces)
    end
    if typeof(bfaces)!=Array{Cint,2}
        bfaces=Array{Cint,2}(bfaces)
    end
    # @assert(size(bfaces,2)>0)

    ndims(bfaceregions)==1 || throw(DimensionMismatch("bfaceregions must be vector"))
    size(bfaceregions,1) == size(bfaces,2) || throw(DimensionMismatch("size(bfaceregions,1) != size(bfaces,2)"))
    if typeof(bfaceregions)!=Array{Cint,1}
        bfaceregions=Array{Cint,1}(bfaceregions)
    end

    ndims(regionpoints)==2 || throw(DimensionMismatch("Region point must be 2D"))
    if size(regionpoints,1)!=2
        regionpoints=transpose(regionpoints)
    end
    if typeof(regionpoints)!=Array{Cdouble,2}
        regionpoints=Array{Cdouble,2}(regionpoints)
    end
    # @assert(size(regionpoints,2)>0)

    @assert ndims(regionnumbers) ==1 || throw(DimensionMismatch("regionnumbers  must be vector"))
    @assert ndims(regionvolumes) ==1 || throw(DimensionMismatch("regionvolumes  must be vector"))
    @assert size(regionnumbers,1) == size(regionpoints,2) || throw(DimensionMismatch("size(regionnumbers,1) != size(regionpoints,2)"))
    @assert size(regionvolumes,1) == size(regionpoints,2) || throw(DimensionMismatch("size(regionvolumes,1) !== size(regionpoints,2)"))

    nholes=0
    nregions=0
    for i=1:length(regionnumbers)
        if regionnumbers[i]==0
            nholes+=1
        else
            nregions+=1
        end
    end

    regionlist=Array{Cdouble,2}(undef,4,nregions)
    holelist=Array{Cdouble,2}(undef,2,nholes)

    ihole=1
    iregion=1
    for i=1:length(regionnumbers)
        if regionnumbers[i]==0
            holelist[1,ihole]=regionpoints[1,i]
            holelist[2,ihole]=regionpoints[2,i]
            ihole+=1
        else
            regionlist[1,iregion]=regionpoints[1,i]
            regionlist[2,iregion]=regionpoints[2,i]
            regionlist[3,iregion]=regionnumbers[i]
            regionlist[4,iregion]=regionvolumes[i]
            iregion+=1
        end
    end

    tio=Triangulate.TriangulateIO()
    tio.pointlist=points

    if size(bfaces,2)>0
        tio.segmentlist=bfaces
    end

    if size(bfaceregions,1)>0
        tio.segmentmarkerlist=bfaceregions
    end

    if size(regionlist,2)>0
        tio.regionlist=regionlist
    end

    if size(holelist,2)>0
        tio.holelist=holelist
    end

    tio
end

"""
$(TYPEDSIGNATURES)

Create triangle input from the current state of the builder.
"""
function triangulateio(this::SimplexGridBuilder)
    dim_space(this)==2 || throw(error("dimension !=2 not implemented"))
    facets=Array{Cint,2}(undef,2,length(this.facets))
    for i=1:length(this.facets)
        facets[1,i]=this.facets[i][1]
        facets[2,i]=this.facets[i][2]
    end

    triangulateio(this.Generator,
                  points=this.pointlist.points,
                  bfaces=facets,
                  bfaceregions=this.facetregions,
                  regionpoints=this.regionpoints,
                  regionnumbers=this.regionnumbers,
                  regionvolumes=this.regionvolumes)

end


text:
using RDatasets
using Random
using Statistics
using LinearAlgebra
using Plots

iris = dataset("datasets", "iris")

X = Matrix(iris[:, 1:4])
y = iris.Species

# Exercise

function split(X, y::AbstractVector; dims=1, ratio_train=0.8, kwargs...)
    n = length(y)
    size(X, dims) == n || throw(DimensionMismatch("..."))

    n_train = round(Int, ratio_train*n)
    i_rand = randperm(n)
    i_train = i_rand[1:n_train]
    i_test = i_rand[n_train+1:end]

    return selectdim(X, dims, i_train), y[i_train], selectdim(X, dims, i_test), y[i_test]
end

X_train, y_train, X_test, y_test = split(X, y)

# Exercise

function normalize(X_train, X_test; dims=1, kwargs...)
    col_mean = mean(X_train; dims)
    col_std = std(X_train; dims)

    return (X_train .- col_mean) ./ col_std, (X_test .- col_mean) ./ col_std
end

X_train, X_test = normalize(X_train, X_test)

# Exercise

function onehot(y, classes)
    y_onehot = falses(length(classes), length(y))
    for (i, class) in enumerate(classes)
        y_onehot[i, y .== class] .= 1
    end
    return y_onehot
end

onecold(y, classes) = [classes[argmax(y_col)] for y_col in eachcol(y)]

classes = unique(y)

isequal(onecold(onehot(y, classes), classes), y)

# Preparing data

function prepare_data(X, y; do_normal=true, do_onehot=true, kwargs...)
    X_train, y_train, X_test, y_test = split(X, y; kwargs...)

    if do_normal
        X_train, X_test = normalize(X_train, X_test; kwargs...)
    end

    classes = unique(y)

    if do_onehot
        y_train = onehot(y_train, classes)
        y_test = onehot(y_test, classes)
    end

    return X_train, y_train, X_test, y_test, classes
end

Random.seed!(666)

iris = dataset("datasets", "iris")

X = Matrix(iris[:, 1:4])
y = iris.Species

X_train, y_train, X_test, y_test, classes = prepare_data(X', y; dims=2)

# Testing dims in prepare_data

Random.seed!(666)
aux1 = prepare_data(X, y; dims=1)
Random.seed!(666)
aux2 = prepare_data(X', y; dims=2)

norm(aux1[1] - aux2[1]')

# SimpleNet

struct SimpleNet{T<:Real}
    W1::Matrix{T}
    b1::Vector{T}
    W2::Matrix{T}
    b2::Vector{T}
end

# Exercise

SimpleNet(n1, n2, n3) = SimpleNet(randn(n2, n1), randn(n2), randn(n3, n2), randn(n3))

# Initialize SimpleNet

Random.seed!(666)

m = SimpleNet(size(X_train,1), 5, size(y_train,1))

# Exercise

function (m::SimpleNet)(x)
    z1 = m.W1*x .+ m.b1
    a1 = max.(z1, 0)
    z2 = m.W2*a1 .+ m.b2
    return exp.(z2) ./ sum(exp.(z2), dims=1)
end

# Evaluate SimpleNet

m(X_train[:,1:2])

# Gradient computation

function grad(m::SimpleNet, x::AbstractVector, y; ϵ=1e-10)
    z1 = m.W1*x .+ m.b1
    a1 = max.(z1, 0)
    z2 = m.W2*a1 .+ m.b2
    a2 = exp.(z2) ./ sum(exp.(z2), dims=1)
    l = -sum(y .* log.(a2 .+ ϵ))

    e_z2 = exp.(z2)
    l_part = (- e_z2 * e_z2' + Diagonal(e_z2 .* sum(e_z2))) / sum(e_z2)^2

    l_a2 = - y ./ (a2 .+ ϵ)
    l_z2 = l_part * l_a2
    l_a1 = m.W2' * l_z2
    l_z1 = l_a1 .* (a1 .> 0)
    l_x = m.W1' * l_z1

    l_W2 = l_z2 * a1'
    l_b2 = l_z2
    l_W1 = l_z1 * x'
    l_b1 = l_z1

    return l, l_W1, l_b1, l_W2, l_b2
end

# Few comments about the gradient

g_all = [grad(m, X_train[:,k], y_train[:,k]) for k in 1:size(X_train,2)]

typeof(g_all)

mean_tuple(d::AbstractArray{<:Tuple}) = Tuple([mean([d[k][i] for k in 1:length(d)]) for i in 1:length(d[1])])

g_mean = mean_tuple(g_all)

typeof(g_mean)

# Exercise

α = 1e-1
max_iter = 200
L = zeros(max_iter)
for iter in 1:max_iter
    grad_all = [grad(m, X_train[:,k], y_train[:,k]) for k in 1:size(X_train,2)]
    grad_mean = mean_tuple(grad_all)

    L[iter] = grad_mean[1]

    m.W1 .-= α*grad_mean[2]
    m.b1 .-= α*grad_mean[3]
    m.W2 .-= α*grad_mean[4]
    m.b2 .-= α*grad_mean[5]
end

plot(L, xlabel="Iteration", ylabel="Loss function", label="", title="Loss function on the training set")

# Exercise

predict(X) = m(X)
accuracy(X, y) = mean(onecold(predict(X), classes) .== onecold(y, classes))

println("Train accuracy = ", accuracy(X_train, y_train))
println("Test accuracy = ", accuracy(X_test, y_test))

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

plot(x=rand(100), y=rand(100), size=rand([1mm, 2mm, 3mm], 100), Geom.point)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SCIP"
version = v"0.2.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://scipopt.org/download/release/scipoptsuite-8.0.0.tgz", "b6ada57f14a728198cc9bc5d22bc6e61dc18cb750a17846df9d853ff6380fbd1"),
]

# Bash recipe for building across all platforms
script = raw"""
cd scipoptsuite*
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix\
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}\
  -DCMAKE_BUILD_TYPE=Release\
  -DPAPILO=0\
  -DZIMPL=OFF\
  -DGCG=0\
  -DUG=0\
  -DAMPL=0\
  -DBOOST=off\
  -DSYM=bliss\
  -DIPOPT_DIR=${prefix} -DIPOPT_LIBRARIES=${libdir} ..
make -j${nproc}
make install

mkdir -p ${prefix}/share/licenses/SCIP
for dir in papilo scip soplex; do
    cp $WORKSPACE/srcdir/scipoptsuite*/${dir}/COPYING ${prefix}/share/licenses/SCIP/LICENSE_${dir}
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libscip", :libscip),
    ExecutableProduct("scip", :scipexe),
    LibraryProduct("libsoplexshared", :libsoplex),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="bliss_jll", uuid="508c9074-7a14-5c94-9582-3d4bc1871065")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.2.0"),
    Dependency(PackageSpec(name="Ipopt_jll", uuid="9cc047cb-c261-5740-88fc-0cf96f7bdcc7")),
    Dependency(PackageSpec(name="Readline_jll", uuid="05236dd9-4125-5232-aa7c-9ec0c9b2c25a")),
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"6", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Round 1A 2014
# Problem B: ???

function main()
	n = int(readline())
	for i = 1:n
	end
end

main()


text:
###This file intends to demonstrate the rendering and basically follows basic.jl
import MuJoCo
import GLFW

using MuJoCo
using MuJoCo.time

mj_activate()
xmlfile = joinpath(@__DIR__, "cartpoleBalance.xml")
model = MJModel(xmlfile)
data = MJData(model)

window = GLFW.CreateWindow(640, 480, "demo")
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)

renderer = MJRenderer(data)

width = Int32(0)
height = Int32(0)
mj_step(model, data)
while !GLFW.WindowShouldClose(window)

    #simulate
    simstart = time(data)
    while (time(data) - simstart < 1.0/60.0 )
        mj_step(model, data)
    end

    #Render
    width, height = GLFW.GetFramebufferSize(window);
    mjv_updateScene(renderer)
    mjr_render(renderer, width, height)

    # Swap front and back buffers
	GLFW.SwapBuffers(window)

	# Poll for and process events
	GLFW.PollEvents()
end

GLFW.Terminate()


text:
# Load ImplicitEquations

## meta these
preds = [(:Lt, :≪, :<), # \ll
         (:Le, :≦, :<=), # \leqq
         (:Eq, :⩵, :(==)), # \Equal
         (:Ge, :≧, :>=), # \gegg
         (:Gt, :≫, :>) # \gg
         ]

#import SymPy: Lt, Le, Eq, Ge, Gt, ≪, ≦, ⩵, ≧, ≫
for (fn, uop, op) in preds
    fnname =  string(fn)

    @eval begin
        ($fn)(f::Function, x::Real) = Pred(f, $op, x)
        ($uop)(f::Function, x::Real) = ($fn)(f, x)
        ($fn)(f::Function, g::Function) = $(fn)((x,y) -> f(x,y) - g(x,y), 0)
        ($uop)(f::Function, g::Function) = ($fn)(f, g)
    end
    eval(Expr(:export, fn))
    eval(Expr(:export, uop))
end
import ImplicitEquations: Neq
export Neq


text:
# single bars

# Jl Plots
using StatsBase, Random
Random.seed!(2018)
s = ["Poor", "Fair", "Good", "Excellent"]
data = sample(s, Weights([1,1,2,2]), 100)
datamap = countmap(data)
bar((x -> datamap[x]).(s), xticks=(1:4, s))

# grouped bars

# Jl Plots
using StatsPlots

measles = [38556, 24472, 14556, 18060, 19549, 8122, 28541, 7880, 3283, 4135, 7953, 1884]
mumps = [20178, 23536, 34561, 37395, 36072, 32237, 18597, 9408, 6005, 6268, 8963, 13882]
chickenPox = [37140, 32169, 37533, 39103, 33244, 23269, 16737, 5411, 3435, 6052, 12825, 23332]

ticklabel = string.(collect('A':'L'))
groupedbar([measles mumps chickenPox], bar_position = :dodge, bar_width=0.7, xticks=(1:12, ticklabel))

# Makie
using CarioMakie

colors = Makie.wong_colors()

# Figure and Axis
fig = Figure(backgroundcolor = :honeydew)
ax = Axis(fig[1,1], xticks = (1:3, ["left", "middle", "right"]),
        title = "Dodged bars with legend")

# Plot
barplot!(ax, tbl.x, tbl.height,
        dodge = tbl.grp,
        color = colors[tbl.grp])

# Legend
labels = ["group 1", "group 2", "group 3"]
elements = [PolyElement(polycolor = colors[i]) for i in 1:length(labels)]
title = "Groups"

Legend(fig[1,2], elements, labels, title)

# Horizontal grouped
barplot(
    tbl.x, tbl.height,
    dodge = tbl.grp,
    color = tbl.grp,
    bar_labels = :y,
    axis = (xticks = (1:3, ["left", "middle", "right"]),
            title = "Dodged bars horizontal with labels"),
    colormap = [:lightsalmon, :mistyrose, :indianred3],
    color_over_background=:red,
    color_over_bar=:white,
    flip_labels_at=0.85,
    direction=:x,
)

# stacked bars

# Jl Plots
using StatsPlots

measles = [38556, 24472, 14556, 18060, 19549, 8122, 28541, 7880, 3283, 4135, 7953, 1884]
mumps = [20178, 23536, 34561, 37395, 36072, 32237, 18597, 9408, 6005, 6268, 8963, 13882]
chickenPox = [37140, 32169, 37533, 39103, 33244, 23269, 16737, 5411, 3435, 6052, 12825, 23332]

ticklabel = string.(collect('A':'L'))
groupedbar([measles mumps chickenPox],
        bar_position = :stack,
        bar_width=0.7,
        xticks=(1:12, ticklabel),
        label=["measles" "mumps" "chickenPox"])

# Makie
# makie bars

# stacked
using CairoMakie

tbl = (x = [1, 1, 1, 2, 2, 2, 3, 3, 3],
       height = 0.1:0.1:0.9,
       grp = [1, 2, 3, 1, 2, 3, 1, 2, 3],
       grp1 = [1, 2, 2, 1, 1, 2, 1, 1, 2],
       grp2 = [1, 1, 2, 1, 2, 1, 1, 2, 1]
       )

barplot(tbl.x, tbl.height,
        stack = tbl.grp,
        color = tbl.grp,
        axis = (xticks = (1:3, ["left", "middle", "right"]),
                title = "Stacked bars"),
        )

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using Test
using ADVUMPS
using ADVUMPS: energy, num_grad, diaglocal
using OMEinsum, Zygote, Random
using Optim, LineSearches
using LinearAlgebra: svd, norm

@testset "non-interacting" begin
    Random.seed!(100)
    model = diaglocal([1,-1.0])
    h = hamiltonian(model)
    as = (rand(3,3,3,3,2) for _ in 1:10)
    @test all(a -> -1 < energy(h, model, SquareIPEPS(a); χ=5, tol=1e-10, maxiter=20)/2 < 1, as)

    a = zeros(2,2,2,2,2) .+ 1e-12 * rand(2,2,2,2,2)
    a[1,1,1,1,2] = randn()
    @test energy(h, model, SquareIPEPS(a); χ=4, tol=1e-10, maxiter=20)/2 ≈ -1

    a = zeros(2,2,2,2,2) .+ 1e-12 * rand(2,2,2,2,2)
    a[1,1,1,1,1] = randn()
    @test energy(h, model, SquareIPEPS(a); χ=4, tol=1e-10, maxiter=20)/2 ≈ 1

    a = zeros(2,2,2,2,2) .+ 1e-12 * rand(2,2,2,2,2)
    a[1,1,1,1,2] = a[1,1,1,1,1] = randn()
    @test abs(energy(h,model,SquareIPEPS(a); χ=4, tol=1e-10, maxiter=20)) < 1e-9

    grad = let energy = x -> real(energy(h, model, SquareIPEPS(x); χ=4, tol=1e-10, maxiter=20))
        res = optimize(energy,
            Δ -> Zygote.gradient(energy,Δ)[1], a, LBFGS(m=20), inplace = false)
    end
    @test grad != Nothing

    Random.seed!(100)
    hdiag = [0.3,-0.43]
    model = diaglocal(hdiag)
    ipeps, key = init_ipeps(model; D=2, χ=4, tol=1e-10, maxiter=20)
    res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    e = minimum(res)/2
    @test isapprox(e, minimum(hdiag), atol=1e-3)
end

@testset "gradient" begin
    Random.seed!(0)
    model = TFIsing(1.0)
    h = hamiltonian(model)
    ipeps, key = init_ipeps(model; D=2, χ=4, tol=1e-10, maxiter=20)
    gradzygote = first(Zygote.gradient(ipeps) do x
        energy(h,model,x; χ=4, tol=1e-10, maxiter=20)
    end).bulk
    gradnum = num_grad(ipeps.bulk, δ=1e-3) do x
        energy(h,model,SquareIPEPS(x); χ=4, tol=1e-10, maxiter=20)
    end
    @test isapprox(gradzygote, gradnum, atol=1e-3)

    # Random.seed!(3)
    # model = Heisenberg()
    # h = hamiltonian(model)
    # ipeps, key = init_ipeps( model; D=2, χ=4, tol=1e-10, maxiter=20)
    # gradzygote = first(Zygote.gradient(ipeps) do x
    #     energy(h,model,x; χ=4, tol=1e-10, maxiter=20)
    # end).bulk
    # gradnum = num_grad(ipeps.bulk, δ=1e-3) do x
    #     energy(h,model, SquareIPEPS(x); χ=4, tol=1e-10, maxiter=20)
    # end
    # @test isapprox(gradzygote , gradnum, atol=1e-3)
end

@testset "TFIsing" begin
    # comparison with results from https://github.com/wangleiphy/tensorgrad
    Random.seed!(3)
    model = TFIsing(1.0)
    ipeps, key = init_ipeps(model; D=2, χ=20, tol=1e-10, maxiter=20)
    res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    e = minimum(res)
    @test isapprox(e, -2.12566, atol = 1e-2)

    # Random.seed!(3)
    # model = TFIsing(0.5)
    # ipeps, key = init_ipeps(model; D=2, χ=20, tol=1e-10, maxiter=20)
    # res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -2.0312, atol = 1e-2)

    # Random.seed!(3)
    # model = TFIsing(2.0)
    # ipeps, key = init_ipeps(model; D=2, χ=20, tol=1e-10, maxiter=20)
    # res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -2.5113, atol = 1e-2)
end

@testset "heisenberg" begin
    # comparison with results from https://github.com/wangleiphy/tensorgrad
    Random.seed!(100)
    model = Heisenberg(1.0,1.0,1.0)
    ipeps, key = init_ipeps(model; D=2, χ=20, tol=1e-10, maxiter=10)
    res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    e = minimum(res)
    @test isapprox(e, -0.66023, atol = 1e-4)

    # Random.seed!(100)
    # model = Heisenberg(1.0,2.0,2.0)
    # ipeps, key = init_ipeps(model; D=2, χ=20, tol=1e-10, maxiter=20)
    # res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -1.190, atol = 1e-3)

    # Random.seed!(100)
    # model = Heisenberg(2.0,0.5,0.5)
    # ipeps, key = init_ipeps(model; D=2, χ=20, tol=1e-10, maxiter=20)
    # res = optimiseipeps(ipeps, key; f_tol = 1e-6)
    # e = minimum(res)
    # @test isapprox(e, -1.0208, atol = 1e-3)
end


text:
"""
    x=cvp(y,H,infinite=Val{true},Umax=1)

Solve the problem `argmin_x ||y-Hx||` for integer x using the technique from
the paper below.  The input vector `y` is of length `n`, with `H` of
dimension `n` by `n`, and the returned vector `x` of length `n`. If
`infinite==Val{true}` then we search the (infinite) lattice, otherwise we
search integers in `[-Umax,Umax]`.  At present cvp does not handle complex
numbers.

Uses alg from "Faster Recursions in Sphere Decoding" Arash Ghasemmehdi, Erik
Agrell, IEEE Transactions on Information Theory, vol 57, issue 6 , June 2011.

# Examples
```jldoctest
julia> H=[1 2; 3 4]; Q,R=qr(H); uhat = cvp(Q'*[0,2],R)
2-element Array{Float64,1}:
  2.0
 -1.0

julia> n=100;H=randn(n,n);Q,R=qr(H);

julia> u=Int.(rand(0:1e10,n));y=H*u+rand(n)/100;

julia> uhat=cvp(Q'*y,R); sum(abs.(u-uhat))
0.0

julia> n=500;H=randn(n,n);Q,R=qr(H);

julia> u=Int.(rand(-1:1,n));y=H*u+rand(n)/10;

julia> uhat=cvp(Q'*y,R,Val{false},1); sum(abs.(u-uhat))
0.0
```
"""
function cvp(r::AbstractArray{Td,1},G::AbstractArray{Td,2},
               infinite=Val{true},Umax=1) where {Td<:Number}

    roundGA(x) = roundFinite(x,Td(Umax))
# nx=1
    n = length(r)
# mxNx=Int(ceil(log2(n)*100000))
    C = Inf
    i = n+1
    d = ones(Int,n)*n
    λ = zeros(n+1)
    F = zeros(Td,n,n)
    F[:,n] = r
    p = zeros(n)
    u = NaN*zeros(Int,n)
    û = NaN*zeros(Int,n)
    Δ = zeros(Int,n)

    begin
        @label LOOP

#        nx+=1
        while λ[i]<C
            if i != 1
                i-=1
                for j=d[i]:-1:i+1
                    F[i,j-1] = F[i,j] - G[i,j]*u[j]
                end
                p[i] = F[i,i]/G[i,i]
                if infinite==Val{true}
                    u[i] = round(p[i])
                else
                    u[i] = roundGA(p[i])
                end
                y = (p[i]-u[i])*G[i,i]
                Δ[i] = signGA(y)
                λ[i] = λ[i+1]+y*y
            else
                û .= u
                C = λ[1]
            end
        end
        m = i

        while λ[i]≥C
#            if i==n || nx>mxNx
            if i==n
                # if nx>mxNx
                #     @warn("terminated search after $(nx) iterations.")
                # end
                return û
            else
                i+=1
                if infinite ≠ Val{true}
                    y = typemax(Td)
                end
                u[i]+=Δ[i]
                Δ[i]=-Δ[i]-signGA(Δ[i])
                if infinite ≠ Val{true}
                    if -Umax ≤ u[i] ≤ Umax
                        y = (p[i]-u[i])*G[i,i]
                    else
                        u[i]+=Δ[i]
                        Δ[i]=-Δ[i]-signGA(Δ[i])
                        if -Umax ≤ u[i] ≤ Umax
                            y = (p[i]-u[i])*G[i,i]
                        end
                    end
                else
                    y = (p[i]-u[i])*G[i,i]
                end
                λ[i] = λ[i+1]+y*y
            end
        end
        d[m:i-1] .= i
        for j=m-1:-1:1
            if d[j]<i
                d[j]=i
            else
                @goto LOOP
            end
        end
        @goto LOOP
    end
end

signGA(x) = x<=0 ? -one(x) : one(x)

function roundFinite(x,Umax)
    y=round(x)
    return  abs(y)>Umax ?  Umax*signGA(y) : y
end

"""
    b = svp(B)

Find the shortest basis vector `b` for the lattice formed by the matrix
`B`. This solves the 'shortest vector problem' (SVP).

We use the [`cvp`](@ref) function in the library, which is not necessarily
the fastest SVP solver.  Roughly follows the CVP-to-SVP reduction in
http://web.eecs.umich.edu/~cpeikert/lic15/lec06.pdf

# Examples
```jldoctest
julia> H=[1 2; 3 4]; svp(H)
2-element Array{Int64,1}:
 -1
  1

julia> H= BigFloat.([2.5 2; 3 4]); svp(H)
2-element Array{BigFloat,1}:
  0.50
 -1.0

```
"""
function svp(B::AbstractArray{Td,2}) where Td
    m,n= size(B)
    V = zeros(Td,m,n)
    c = zeros(Td,n)

    for i = 1:n
        Bi = copy(B)
        Bi[:,i] *=2
        Q,R = qr(Bi)
        vc = cvp(Q'*B[:,i],R)
        V[:,i] = Bi*vc-B[:,i]
        c[i] = V[:,i]'*V[:,i]
    end
    idx = sortperm(c)
    return V[:,idx[1]]
end