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

module ArithmeticTest

using Test
using Dates

@testset "Time arithmetic" begin
    a = DatesPlus.Time(23, 59, 59)
    b = DatesPlus.Time(11, 59, 59)
    @test DatesPlus.CompoundPeriod(a - b) == DatesPlus.Hour(12)
end
@testset "Wrapping arithmetic for Months" begin
    # This ends up being trickier than expected because
    # the user might do 2014-01-01 + Month(-14)
    # monthwrap figures out the resulting month
    # when adding/subtracting months from a date
    @test DatesPlus.monthwrap(1, -14) == 11
    @test DatesPlus.monthwrap(1, -13) == 12
    @test DatesPlus.monthwrap(1, -12) == 1
    @test DatesPlus.monthwrap(1, -11) == 2
    @test DatesPlus.monthwrap(1, -10) == 3
    @test DatesPlus.monthwrap(1, -9) == 4
    @test DatesPlus.monthwrap(1, -8) == 5
    @test DatesPlus.monthwrap(1, -7) == 6
    @test DatesPlus.monthwrap(1, -6) == 7
    @test DatesPlus.monthwrap(1, -5) == 8
    @test DatesPlus.monthwrap(1, -4) == 9
    @test DatesPlus.monthwrap(1, -3) == 10
    @test DatesPlus.monthwrap(1, -2) == 11
    @test DatesPlus.monthwrap(1, -1) == 12
    @test DatesPlus.monthwrap(1, 0) == 1
    @test DatesPlus.monthwrap(1, 1) == 2
    @test DatesPlus.monthwrap(1, 2) == 3
    @test DatesPlus.monthwrap(1, 3) == 4
    @test DatesPlus.monthwrap(1, 4) == 5
    @test DatesPlus.monthwrap(1, 5) == 6
    @test DatesPlus.monthwrap(1, 6) == 7
    @test DatesPlus.monthwrap(1, 7) == 8
    @test DatesPlus.monthwrap(1, 8) == 9
    @test DatesPlus.monthwrap(1, 9) == 10
    @test DatesPlus.monthwrap(1, 10) == 11
    @test DatesPlus.monthwrap(1, 11) == 12
    @test DatesPlus.monthwrap(1, 12) == 1
    @test DatesPlus.monthwrap(1, 13) == 2
    @test DatesPlus.monthwrap(1, 24) == 1
    @test DatesPlus.monthwrap(12, -14) == 10
    @test DatesPlus.monthwrap(12, -13) == 11
    @test DatesPlus.monthwrap(12, -12) == 12
    @test DatesPlus.monthwrap(12, -11) == 1
    @test DatesPlus.monthwrap(12, -2) == 10
    @test DatesPlus.monthwrap(12, -1) == 11
    @test DatesPlus.monthwrap(12, 0) == 12
    @test DatesPlus.monthwrap(12, 1) == 1
    @test DatesPlus.monthwrap(12, 2) == 2
    @test DatesPlus.monthwrap(12, 11) == 11
    @test DatesPlus.monthwrap(12, 12) == 12
    @test DatesPlus.monthwrap(12, 13) == 1
end
@testset "Wrapping arithmetic for years" begin
    # yearwrap figures out the resulting year
    # when adding/subtracting months from a date
    @test DatesPlus.yearwrap(2000, 1, -3600) == 1700
    @test DatesPlus.yearwrap(2000, 1, -37) == 1996
    @test DatesPlus.yearwrap(2000, 1, -36) == 1997
    @test DatesPlus.yearwrap(2000, 1, -35) == 1997
    @test DatesPlus.yearwrap(2000, 1, -25) == 1997
    @test DatesPlus.yearwrap(2000, 1, -24) == 1998
    @test DatesPlus.yearwrap(2000, 1, -23) == 1998
    @test DatesPlus.yearwrap(2000, 1, -14) == 1998
    @test DatesPlus.yearwrap(2000, 1, -13) == 1998
    @test DatesPlus.yearwrap(2000, 1, -12) == 1999
    @test DatesPlus.yearwrap(2000, 1, -11) == 1999
    @test DatesPlus.yearwrap(2000, 1, -2) == 1999
    @test DatesPlus.yearwrap(2000, 1, -1) == 1999
    @test DatesPlus.yearwrap(2000, 1, 0) == 2000
    @test DatesPlus.yearwrap(2000, 1, 1) == 2000
    @test DatesPlus.yearwrap(2000, 1, 11) == 2000
    @test DatesPlus.yearwrap(2000, 1, 12) == 2001
    @test DatesPlus.yearwrap(2000, 1, 13) == 2001
    @test DatesPlus.yearwrap(2000, 1, 23) == 2001
    @test DatesPlus.yearwrap(2000, 1, 24) == 2002
    @test DatesPlus.yearwrap(2000, 1, 25) == 2002
    @test DatesPlus.yearwrap(2000, 1, 36) == 2003
    @test DatesPlus.yearwrap(2000, 1, 3600) == 2300
    @test DatesPlus.yearwrap(2000, 2, -2) == 1999
    @test DatesPlus.yearwrap(2000, 3, 10) == 2001
    @test DatesPlus.yearwrap(2000, 4, -4) == 1999
    @test DatesPlus.yearwrap(2000, 5, 8) == 2001
    @test DatesPlus.yearwrap(2000, 6, -18) == 1998
    @test DatesPlus.yearwrap(2000, 6, -6) == 1999
    @test DatesPlus.yearwrap(2000, 6, 6) == 2000
    @test DatesPlus.yearwrap(2000, 6, 7) == 2001
    @test DatesPlus.yearwrap(2000, 6, 19) == 2002
    @test DatesPlus.yearwrap(2000, 12, -3600) == 1700
    @test DatesPlus.yearwrap(2000, 12, -36) == 1997
    @test DatesPlus.yearwrap(2000, 12, -35) == 1998
    @test DatesPlus.yearwrap(2000, 12, -24) == 1998
    @test DatesPlus.yearwrap(2000, 12, -23) == 1999
    @test DatesPlus.yearwrap(2000, 12, -14) == 1999
    @test DatesPlus.yearwrap(2000, 12, -13) == 1999
    @test DatesPlus.yearwrap(2000, 12, -12) == 1999
    @test DatesPlus.yearwrap(2000, 12, -11) == 2000
    @test DatesPlus.yearwrap(2000, 12, -2) == 2000
    @test DatesPlus.yearwrap(2000, 12, -1) == 2000
    @test DatesPlus.yearwrap(2000, 12, 0) == 2000
    @test DatesPlus.yearwrap(2000, 12, 1) == 2001
    @test DatesPlus.yearwrap(2000, 12, 11) == 2001
    @test DatesPlus.yearwrap(2000, 12, 12) == 2001
    @test DatesPlus.yearwrap(2000, 12, 13) == 2002
    @test DatesPlus.yearwrap(2000, 12, 24) == 2002
    @test DatesPlus.yearwrap(2000, 12, 25) == 2003
    @test DatesPlus.yearwrap(2000, 12, 36) == 2003
    @test DatesPlus.yearwrap(2000, 12, 37) == 2004
    @test DatesPlus.yearwrap(2000, 12, 3600) == 2300
end
@testset "DateTime arithmetic" begin
    # DatesPlus.DateTime arithmetic
    a = DatesPlus.DateTime(2013, 1, 1, 0, 0, 0, 1)
    b = DatesPlus.DateTime(2013, 1, 1, 0, 0, 0, 0)
    @test a - b == DatesPlus.Millisecond(1)
    @test DatesPlus.DateTime(2013, 1, 2) - b == DatesPlus.Millisecond(86400000)

    @testset "DatesPlus.DateTime-Year arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Year(1) == DatesPlus.DateTime(2000, 12, 27)
        @test dt + DatesPlus.Year(100) == DatesPlus.DateTime(2099, 12, 27)
        @test dt + DatesPlus.Year(1000) == DatesPlus.DateTime(2999, 12, 27)
        @test dt - DatesPlus.Year(1) == DatesPlus.DateTime(1998, 12, 27)
        @test dt - DatesPlus.Year(100) == DatesPlus.DateTime(1899, 12, 27)
        @test dt - DatesPlus.Year(1000) == DatesPlus.DateTime(999, 12, 27)
        dt = DatesPlus.DateTime(2000, 2, 29)
        @test dt + DatesPlus.Year(1) == DatesPlus.DateTime(2001, 2, 28)
        @test dt - DatesPlus.Year(1) == DatesPlus.DateTime(1999, 2, 28)
        @test dt + DatesPlus.Year(4) == DatesPlus.DateTime(2004, 2, 29)
        @test dt - DatesPlus.Year(4) == DatesPlus.DateTime(1996, 2, 29)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Year(1) == DatesPlus.DateTime(1973, 6, 30, 23, 59, 59)
        @test dt - DatesPlus.Year(1) == DatesPlus.DateTime(1971, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Year(-1) == DatesPlus.DateTime(1971, 6, 30, 23, 59, 59)
        @test dt - DatesPlus.Year(-1) == DatesPlus.DateTime(1973, 6, 30, 23, 59, 59)

        dt = DatesPlus.DateTime(2000, 1, 1, 12, 30, 45, 500)
        dt2 = dt + DatesPlus.Year(1)
        @test DatesPlus.year(dt2) == 2001
        @test DatesPlus.month(dt2) == 1
        @test DatesPlus.day(dt2) == 1
        @test DatesPlus.hour(dt2) == 12
        @test DatesPlus.minute(dt2) == 30
        @test DatesPlus.second(dt2) == 45
        @test DatesPlus.millisecond(dt2) == 500
    end
    @testset "DateTime-Quarter arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Quarter(1) == DatesPlus.DateTime(2000, 3, 27)
        @test dt + DatesPlus.Quarter(-1) == DatesPlus.DateTime(1999, 9, 27)
    end

    @testset "DateTime-Month arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Month(1) == DatesPlus.DateTime(2000, 1, 27)
        @test dt + DatesPlus.Month(-1) == DatesPlus.DateTime(1999, 11, 27)
        @test dt + DatesPlus.Month(-11) == DatesPlus.DateTime(1999, 1, 27)
        @test dt + DatesPlus.Month(11) == DatesPlus.DateTime(2000, 11, 27)
        @test dt + DatesPlus.Month(-12) == DatesPlus.DateTime(1998, 12, 27)
        @test dt + DatesPlus.Month(12) == DatesPlus.DateTime(2000, 12, 27)
        @test dt + DatesPlus.Month(13) == DatesPlus.DateTime(2001, 1, 27)
        @test dt + DatesPlus.Month(100) == DatesPlus.DateTime(2008, 4, 27)
        @test dt + DatesPlus.Month(1000) == DatesPlus.DateTime(2083, 4, 27)
        @test dt - DatesPlus.Month(1) == DatesPlus.DateTime(1999, 11, 27)
        @test dt - DatesPlus.Month(-1) == DatesPlus.DateTime(2000, 1, 27)
        @test dt - DatesPlus.Month(100) == DatesPlus.DateTime(1991, 8, 27)
        @test dt - DatesPlus.Month(1000) == DatesPlus.DateTime(1916, 8, 27)
        dt = DatesPlus.DateTime(2000, 2, 29)
        @test dt + DatesPlus.Month(1) == DatesPlus.DateTime(2000, 3, 29)
        @test dt - DatesPlus.Month(1) == DatesPlus.DateTime(2000, 1, 29)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Month(1) == DatesPlus.DateTime(1972, 7, 30, 23, 59, 59)
        @test dt - DatesPlus.Month(1) == DatesPlus.DateTime(1972, 5, 30, 23, 59, 59)
        @test dt + DatesPlus.Month(-1) == DatesPlus.DateTime(1972, 5, 30, 23, 59, 59)
    end
    @testset "DateTime-Week arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Week(1) == DatesPlus.DateTime(2000, 1, 3)
        @test dt + DatesPlus.Week(100) == DatesPlus.DateTime(2001, 11, 26)
        @test dt + DatesPlus.Week(1000) == DatesPlus.DateTime(2019, 2, 25)
        @test dt - DatesPlus.Week(1) == DatesPlus.DateTime(1999, 12, 20)
        @test dt - DatesPlus.Week(100) == DatesPlus.DateTime(1998, 1, 26)
        @test dt - DatesPlus.Week(1000) == DatesPlus.DateTime(1980, 10, 27)
        dt = DatesPlus.DateTime(2000, 2, 29)
        @test dt + DatesPlus.Week(1) == DatesPlus.DateTime(2000, 3, 7)
        @test dt - DatesPlus.Week(1) == DatesPlus.DateTime(2000, 2, 22)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Week(1) == DatesPlus.DateTime(1972, 7, 7, 23, 59, 59)
        @test dt - DatesPlus.Week(1) == DatesPlus.DateTime(1972, 6, 23, 23, 59, 59)
        @test dt + DatesPlus.Week(-1) == DatesPlus.DateTime(1972, 6, 23, 23, 59, 59)
    end
    @testset "DateTime-Day arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Day(1) == DatesPlus.DateTime(1999, 12, 28)
        @test dt + DatesPlus.Day(100) == DatesPlus.DateTime(2000, 4, 5)
        @test dt + DatesPlus.Day(1000) == DatesPlus.DateTime(2002, 9, 22)
        @test dt - DatesPlus.Day(1) == DatesPlus.DateTime(1999, 12, 26)
        @test dt - DatesPlus.Day(100) == DatesPlus.DateTime(1999, 9, 18)
        @test dt - DatesPlus.Day(1000) == DatesPlus.DateTime(1997, 4, 1)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Day(1) == DatesPlus.DateTime(1972, 7, 1, 23, 59, 59)
        @test dt - DatesPlus.Day(1) == DatesPlus.DateTime(1972, 6, 29, 23, 59, 59)
        @test dt + DatesPlus.Day(-1) == DatesPlus.DateTime(1972, 6, 29, 23, 59, 59)
    end
    @testset "DateTime-Hour arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Hour(1) == DatesPlus.DateTime(1999, 12, 27, 1)
        @test dt + DatesPlus.Hour(100) == DatesPlus.DateTime(1999, 12, 31, 4)
        @test dt + DatesPlus.Hour(1000) == DatesPlus.DateTime(2000, 2, 6, 16)
        @test dt - DatesPlus.Hour(1) == DatesPlus.DateTime(1999, 12, 26, 23)
        @test dt - DatesPlus.Hour(100) == DatesPlus.DateTime(1999, 12, 22, 20)
        @test dt - DatesPlus.Hour(1000) == DatesPlus.DateTime(1999, 11, 15, 8)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Hour(1) == DatesPlus.DateTime(1972, 7, 1, 0, 59, 59)
        @test dt - DatesPlus.Hour(1) == DatesPlus.DateTime(1972, 6, 30, 22, 59, 59)
        @test dt + DatesPlus.Hour(-1) == DatesPlus.DateTime(1972, 6, 30, 22, 59, 59)
    end
    @testset "DateTime-Minute arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Minute(1) == DatesPlus.DateTime(1999, 12, 27, 0, 1)
        @test dt + DatesPlus.Minute(100) == DatesPlus.DateTime(1999, 12, 27, 1, 40)
        @test dt + DatesPlus.Minute(1000) == DatesPlus.DateTime(1999, 12, 27, 16, 40)
        @test dt - DatesPlus.Minute(1) == DatesPlus.DateTime(1999, 12, 26, 23, 59)
        @test dt - DatesPlus.Minute(100) == DatesPlus.DateTime(1999, 12, 26, 22, 20)
        @test dt - DatesPlus.Minute(1000) == DatesPlus.DateTime(1999, 12, 26, 7, 20)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Minute(1) == DatesPlus.DateTime(1972, 7, 1, 0, 0, 59)
        @test dt - DatesPlus.Minute(1) == DatesPlus.DateTime(1972, 6, 30, 23, 58, 59)
        @test dt + DatesPlus.Minute(-1) == DatesPlus.DateTime(1972, 6, 30, 23, 58, 59)
    end
    @testset "DateTime-Second arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Second(1) == DatesPlus.DateTime(1999, 12, 27, 0, 0, 1)
        @test dt + DatesPlus.Second(100) == DatesPlus.DateTime(1999, 12, 27, 0, 1, 40)
        @test dt + DatesPlus.Second(1000) == DatesPlus.DateTime(1999, 12, 27, 0, 16, 40)
        @test dt - DatesPlus.Second(1) == DatesPlus.DateTime(1999, 12, 26, 23, 59, 59)
        @test dt - DatesPlus.Second(100) == DatesPlus.DateTime(1999, 12, 26, 23, 58, 20)
        @test dt - DatesPlus.Second(1000) == DatesPlus.DateTime(1999, 12, 26, 23, 43, 20)
    end
    @testset "DateTime-Millisecond arithmetic" begin
        dt = DatesPlus.DateTime(1999, 12, 27)
        @test dt + DatesPlus.Millisecond(1) == DatesPlus.DateTime(1999, 12, 27, 0, 0, 0, 1)
        @test dt + DatesPlus.Millisecond(100) == DatesPlus.DateTime(1999, 12, 27, 0, 0, 0, 100)
        @test dt + DatesPlus.Millisecond(1000) == DatesPlus.DateTime(1999, 12, 27, 0, 0, 1)
        @test dt - DatesPlus.Millisecond(1) == DatesPlus.DateTime(1999, 12, 26, 23, 59, 59, 999)
        @test dt - DatesPlus.Millisecond(100) == DatesPlus.DateTime(1999, 12, 26, 23, 59, 59, 900)
        @test dt - DatesPlus.Millisecond(1000) == DatesPlus.DateTime(1999, 12, 26, 23, 59, 59)
        dt = DatesPlus.DateTime(1972, 6, 30, 23, 59, 59)
        @test dt + DatesPlus.Millisecond(1) == DatesPlus.DateTime(1972, 6, 30, 23, 59, 59, 1)
        @test dt - DatesPlus.Millisecond(1) == DatesPlus.DateTime(1972, 6, 30, 23, 59, 58, 999)
        @test dt + DatesPlus.Millisecond(-1) == DatesPlus.DateTime(1972, 6, 30, 23, 59, 58, 999)
    end
end
@testset "Date arithmetic" begin
    @testset "Date-Year arithmetic" begin
        dt = DatesPlus.Date(1999, 12, 27)
        @test dt + DatesPlus.Year(1) == DatesPlus.Date(2000, 12, 27)
        @test dt + DatesPlus.Year(100) == DatesPlus.Date(2099, 12, 27)
        @test dt + DatesPlus.Year(1000) == DatesPlus.Date(2999, 12, 27)
        @test dt - DatesPlus.Year(1) == DatesPlus.Date(1998, 12, 27)
        @test dt - DatesPlus.Year(100) == DatesPlus.Date(1899, 12, 27)
        @test dt - DatesPlus.Year(1000) == DatesPlus.Date(999, 12, 27)
        dt = DatesPlus.Date(2000, 2, 29)
        @test dt + DatesPlus.Year(1) == DatesPlus.Date(2001, 2, 28)
        @test dt - DatesPlus.Year(1) == DatesPlus.Date(1999, 2, 28)
        @test dt + DatesPlus.Year(4) == DatesPlus.Date(2004, 2, 29)
        @test dt - DatesPlus.Year(4) == DatesPlus.Date(1996, 2, 29)
    end
    @testset "Date-Quarter arithmetic" begin
        dt = DatesPlus.Date(1999, 12, 27)
        @test dt + DatesPlus.Quarter(1) == DatesPlus.Date(2000, 3, 27)
        @test dt - DatesPlus.Quarter(1) == DatesPlus.Date(1999, 9, 27)
    end
    @testset "Date-Month arithmetic" begin
        dt = DatesPlus.Date(1999, 12, 27)
        @test dt + DatesPlus.Month(1) == DatesPlus.Date(2000, 1, 27)
        @test dt + DatesPlus.Month(100) == DatesPlus.Date(2008, 4, 27)
        @test dt + DatesPlus.Month(1000) == DatesPlus.Date(2083, 4, 27)
        @test dt - DatesPlus.Month(1) == DatesPlus.Date(1999, 11, 27)
        @test dt - DatesPlus.Month(100) == DatesPlus.Date(1991, 8, 27)
        @test dt - DatesPlus.Month(1000) == DatesPlus.Date(1916, 8, 27)
        dt = DatesPlus.Date(2000, 2, 29)
        @test dt + DatesPlus.Month(1) == DatesPlus.Date(2000, 3, 29)
        @test dt - DatesPlus.Month(1) == DatesPlus.Date(2000, 1, 29)
    end
    @testset "Date-Week arithmetic" begin
        dt = DatesPlus.Date(1999, 12, 27)
        @test dt + DatesPlus.Week(1) == DatesPlus.Date(2000, 1, 3)
        @test dt + DatesPlus.Week(100) == DatesPlus.Date(2001, 11, 26)
        @test dt + DatesPlus.Week(1000) == DatesPlus.Date(2019, 2, 25)
        @test dt - DatesPlus.Week(1) == DatesPlus.Date(1999, 12, 20)
        @test dt - DatesPlus.Week(100) == DatesPlus.Date(1998, 1, 26)
        @test dt - DatesPlus.Week(1000) == DatesPlus.Date(1980, 10, 27)
        dt = DatesPlus.Date(2000, 2, 29)
        @test dt + DatesPlus.Week(1) == DatesPlus.Date(2000, 3, 7)
        @test dt - DatesPlus.Week(1) == DatesPlus.Date(2000, 2, 22)
    end
    @testset "Date-Day arithmetic" begin
        dt = DatesPlus.Date(1999, 12, 27)
        @test dt + DatesPlus.Day(1) == DatesPlus.Date(1999, 12, 28)
        @test dt + DatesPlus.Day(100) == DatesPlus.Date(2000, 4, 5)
        @test dt + DatesPlus.Day(1000) == DatesPlus.Date(2002, 9, 22)
        @test dt - DatesPlus.Day(1) == DatesPlus.Date(1999, 12, 26)
        @test dt - DatesPlus.Day(100) == DatesPlus.Date(1999, 9, 18)
        @test dt - DatesPlus.Day(1000) == DatesPlus.Date(1997, 4, 1)
    end
    @testset "Date-Time arithmetic" begin
        dt = DatesPlus.Date(1999, 12, 27)
        @test dt + DatesPlus.Time(0, 0, 0) == DatesPlus.DateTime(1999, 12, 27, 0, 0, 0)
        @test dt + DatesPlus.Time(1, 0, 0) == DatesPlus.DateTime(1999, 12, 27, 1, 0, 0)
        @test dt + DatesPlus.Time(0, 1, 0) == DatesPlus.DateTime(1999, 12, 27, 0, 1, 0)
        @test dt + DatesPlus.Time(0, 0, 1) == DatesPlus.DateTime(1999, 12, 27, 0, 0, 1)
        @test DatesPlus.Time(0, 0, 1) + dt == DatesPlus.DateTime(1999, 12, 27, 0, 0, 1)

        t = DatesPlus.Time(0, 0, 0) + DatesPlus.Hour(24)
        @test dt + t == DatesPlus.DateTime(1999, 12, 27, 0, 0, 0)
    end
end
@testset "Time-TimePeriod arithmetic" begin
    t = DatesPlus.Time(0)
    @test t + DatesPlus.Hour(1) == DatesPlus.Time(1)
    @test t - DatesPlus.Hour(1) == DatesPlus.Time(23)
    @test t - DatesPlus.Nanosecond(1) == DatesPlus.Time(23, 59, 59, 999, 999, 999)
    @test t + DatesPlus.Nanosecond(-1) == DatesPlus.Time(23, 59, 59, 999, 999, 999)
    @test t + DatesPlus.Hour(24) == t
    @test t + DatesPlus.Nanosecond(86400000000000) == t
    @test t - DatesPlus.Nanosecond(86400000000000) == t
    @test t + DatesPlus.Minute(1) == DatesPlus.Time(0, 1)
    @test t + DatesPlus.Second(1) == DatesPlus.Time(0, 0, 1)
    @test t + DatesPlus.Millisecond(1) == DatesPlus.Time(0, 0, 0, 1)
    @test t + DatesPlus.Microsecond(1) == DatesPlus.Time(0, 0, 0, 0, 1)
    @test_throws MethodError t + DatesPlus.Day(1)
    @testset "Time-TimePeriod arithmetic inequalities" begin
        t = DatesPlus.Time(4, 20)
        @test t - DatesPlus.Nanosecond(1) < t
        @test t + DatesPlus.Nanosecond(1) > t
        @test t + DatesPlus.Hour(24) < typemax(DatesPlus.Time)
        @test t - DatesPlus.Hour(24) > typemin(DatesPlus.Time)
        @test t + DatesPlus.Hour(23) < t
        @test t + DatesPlus.Hour(25) > t
    end
end
@testset "Month arithmetic and non-associativity" begin
    # Month arithmetic minimizes "edit distance", or number of changes
    # needed to get a correct answer
    # This approach results in a few cases of non-associativity
    a = DatesPlus.Date(2012, 1, 29)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 1, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 2, 29)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 3, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 4, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 5, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 6, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 8, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 9, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 10, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
    a = DatesPlus.Date(2012, 11, 30)
    @test (a + DatesPlus.Day(1)) + DatesPlus.Month(1) != (a + DatesPlus.Month(1)) + DatesPlus.Day(1)
end
@testset "Vectorized arithmetic" begin
    # Vectorized Time arithmetic
    a = DatesPlus.Time(1, 1, 1)
    dr = [a, a, a, a, a, a, a, a, a, a]
    b = a + DatesPlus.Hour(1)
    @test dr .+ DatesPlus.Hour(1) == repeat([b], 10)
    b = a + DatesPlus.Second(1)
    @test dr .+ DatesPlus.Second(1) == repeat([b], 10)
    b = a + DatesPlus.Millisecond(1)
    @test dr .+ DatesPlus.Millisecond(1) == repeat([b], 10)
    b = a + DatesPlus.Microsecond(1)
    @test dr .+ DatesPlus.Microsecond(1) == repeat([b], 10)
    b = a + DatesPlus.Nanosecond(1)
    @test dr .+ DatesPlus.Nanosecond(1) == repeat([b], 10)

    b = a - DatesPlus.Hour(1)
    @test dr .- DatesPlus.Hour(1) == repeat([b], 10)
    b = a - DatesPlus.Second(1)
    @test dr .- DatesPlus.Second(1) == repeat([b], 10)
    b = a - DatesPlus.Millisecond(1)
    @test dr .- DatesPlus.Millisecond(1) == repeat([b], 10)
    b = a - DatesPlus.Microsecond(1)
    @test dr .- DatesPlus.Microsecond(1) == repeat([b], 10)
    b = a - DatesPlus.Nanosecond(1)
    @test dr .- DatesPlus.Nanosecond(1) == repeat([b], 10)

    # Vectorized arithmetic
    a = DatesPlus.Date(2014, 1, 1)
    dr = [a, a, a, a, a, a, a, a, a, a]
    b = a + DatesPlus.Year(1)
    @test dr .+ DatesPlus.Year(1) == repeat([b], 10)
    b = a + DatesPlus.Month(1)
    @test dr .+ DatesPlus.Month(1) == repeat([b], 10)
    b = a + DatesPlus.Day(1)
    @test dr .+ DatesPlus.Day(1) == repeat([b], 10)
    b = a - DatesPlus.Year(1)
    @test dr .- DatesPlus.Year(1) == repeat([b], 10)
    b = a - DatesPlus.Month(1)
    @test dr .- DatesPlus.Month(1) == repeat([b], 10)
    b = a - DatesPlus.Day(1)
    @test dr .- DatesPlus.Day(1) == repeat([b], 10)

    # Vectorized arithmetic
    b = a + DatesPlus.Year(1)
    @test dr .+ DatesPlus.Year(1) == repeat([b], 10)
    b = a + DatesPlus.Month(1)
    @test dr .+ DatesPlus.Month(1) == repeat([b], 10)
    b = a + DatesPlus.Day(1)
    @test dr .+ DatesPlus.Day(1) == repeat([b], 10)
    b = a - DatesPlus.Year(1)
    @test dr .- DatesPlus.Year(1) == repeat([b], 10)
    b = a - DatesPlus.Month(1)
    @test dr .- DatesPlus.Month(1) == repeat([b], 10)
    b = a - DatesPlus.Day(1)
    @test dr .- DatesPlus.Day(1) == repeat([b], 10)
    t1 = [DatesPlus.Date(2009, 1, 1) DatesPlus.Date(2009, 1, 2) DatesPlus.Date(2009, 1, 3); DatesPlus.Date(2009, 2, 1) DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3)]
    t2 = [DatesPlus.Date(2009, 1, 2) DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2010, 1, 3); DatesPlus.Date(2010, 2, 1) DatesPlus.Date(2009, 3, 2) DatesPlus.Date(2009, 2, 4)]
    t3 = [DatesPlus.DateTime(2009, 1, 1), DatesPlus.DateTime(2009, 1, 2), DatesPlus.DateTime(2009, 1, 3)]
    t4 = [DatesPlus.DateTime(2009, 1, 1, 0, 0, 1), DatesPlus.DateTime(2009, 1, 2, 0, 1), DatesPlus.DateTime(2009, 1, 3, 1)]
    t5 = [DatesPlus.Time(0, 0, 0) DatesPlus.Time(0, 0, 1) DatesPlus.Time(0, 0, 2); DatesPlus.Time(0, 1, 0) DatesPlus.Time(0, 2, 0) DatesPlus.Time(0, 3, 0)]

    @testset "TimeType, Array{TimeType}" begin
        @test DatesPlus.Date(2010, 1, 1) .- t1 == [DatesPlus.Day(365) DatesPlus.Day(364) DatesPlus.Day(363); DatesPlus.Day(334) DatesPlus.Day(333) DatesPlus.Day(332)]
        @test t1 .- DatesPlus.Date(2010, 1, 1) == [DatesPlus.Day(-365) DatesPlus.Day(-364) DatesPlus.Day(-363); DatesPlus.Day(-334) DatesPlus.Day(-333) DatesPlus.Day(-332)]
        @test DatesPlus.DateTime(2009, 1, 1) .- t3 == [DatesPlus.Millisecond(0), DatesPlus.Millisecond(-86400000), DatesPlus.Millisecond(-172800000)]
        @test t3 .- DatesPlus.DateTime(2009, 1, 1) == [DatesPlus.Millisecond(0), DatesPlus.Millisecond(86400000), DatesPlus.Millisecond(172800000)]
        @test DatesPlus.Time(2) .- t5 == [DatesPlus.Nanosecond(7200000000000) DatesPlus.Nanosecond(7199000000000) DatesPlus.Nanosecond(7198000000000);
                                      DatesPlus.Nanosecond(7140000000000) DatesPlus.Nanosecond(7080000000000) DatesPlus.Nanosecond(7020000000000)]
    end
    @testset "GeneralPeriod, Array{TimeType}" begin
        @test DatesPlus.Day(1) .+ t1 == [DatesPlus.Date(2009, 1, 2) DatesPlus.Date(2009, 1, 3) DatesPlus.Date(2009, 1, 4); DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4)]
        @test DatesPlus.Hour(1) .+ t3 == [DatesPlus.DateTime(2009, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1), DatesPlus.DateTime(2009, 1, 3, 1)]
        @test t1 .+ DatesPlus.Day(1) == [DatesPlus.Date(2009, 1, 2) DatesPlus.Date(2009, 1, 3) DatesPlus.Date(2009, 1, 4); DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4)]
        @test t3 .+ DatesPlus.Hour(1) == [DatesPlus.DateTime(2009, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1), DatesPlus.DateTime(2009, 1, 3, 1)]

        @test (DatesPlus.Month(1) + DatesPlus.Day(1)) .+ t1 == [DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4); DatesPlus.Date(2009, 3, 2) DatesPlus.Date(2009, 3, 3) DatesPlus.Date(2009, 3, 4)]
        @test (DatesPlus.Hour(1) + DatesPlus.Minute(1)) .+ t3 == [DatesPlus.DateTime(2009, 1, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1, 1), DatesPlus.DateTime(2009, 1, 3, 1, 1)]
        @test t1 .+ (DatesPlus.Month(1) + DatesPlus.Day(1)) == [DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4); DatesPlus.Date(2009, 3, 2) DatesPlus.Date(2009, 3, 3) DatesPlus.Date(2009, 3, 4)]
        @test t3 .+ (DatesPlus.Hour(1) + DatesPlus.Minute(1)) == [DatesPlus.DateTime(2009, 1, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1, 1), DatesPlus.DateTime(2009, 1, 3, 1, 1)]

        @test t1 .- DatesPlus.Day(1) == [DatesPlus.Date(2008, 12, 31) DatesPlus.Date(2009, 1, 1) DatesPlus.Date(2009, 1, 2); DatesPlus.Date(2009, 1, 31) DatesPlus.Date(2009, 2, 1) DatesPlus.Date(2009, 2, 2)]
        @test t3 .- DatesPlus.Hour(1) == [DatesPlus.DateTime(2008, 12, 31, 23), DatesPlus.DateTime(2009, 1, 1, 23), DatesPlus.DateTime(2009, 1, 2, 23)]

        @test t1 .- (DatesPlus.Month(1) + DatesPlus.Day(1)) == [DatesPlus.Date(2008, 11, 30) DatesPlus.Date(2008, 12, 1) DatesPlus.Date(2008, 12, 2); DatesPlus.Date(2008, 12, 31) DatesPlus.Date(2009, 1, 1) DatesPlus.Date(2009, 1, 2)]
        @test t3 .- (DatesPlus.Hour(1) + DatesPlus.Minute(1)) == [DatesPlus.DateTime(2008, 12, 31, 22, 59), DatesPlus.DateTime(2009, 1, 1, 22, 59), DatesPlus.DateTime(2009, 1, 2, 22, 59)]
    end
    @testset "#20205" begin
        # ensure commutative subtraction methods are not defined
        @test_throws MethodError DatesPlus.Day(1) .- t1
        @test_throws MethodError DatesPlus.Hour(1) .- t3
        @test_throws MethodError (DatesPlus.Month(1) + DatesPlus.Day(1)) .- t1
        @test_throws MethodError (DatesPlus.Hour(1) + DatesPlus.Minute(1)) .- t3
    end
    @testset "deprecated" begin
        @test DatesPlus.Date(2010, 1, 1) - t1 == [DatesPlus.Day(365) DatesPlus.Day(364) DatesPlus.Day(363); DatesPlus.Day(334) DatesPlus.Day(333) DatesPlus.Day(332)]
        @test t1 - DatesPlus.Date(2010, 1, 1) == [DatesPlus.Day(-365) DatesPlus.Day(-364) DatesPlus.Day(-363); DatesPlus.Day(-334) DatesPlus.Day(-333) DatesPlus.Day(-332)]
        @test DatesPlus.DateTime(2009, 1, 1) - t3 == [DatesPlus.Millisecond(0), DatesPlus.Millisecond(-86400000), DatesPlus.Millisecond(-172800000)]
        @test t3 - DatesPlus.DateTime(2009, 1, 1) == [DatesPlus.Millisecond(0), DatesPlus.Millisecond(86400000), DatesPlus.Millisecond(172800000)]
        @test DatesPlus.Day(1) + t1 == [DatesPlus.Date(2009, 1, 2) DatesPlus.Date(2009, 1, 3) DatesPlus.Date(2009, 1, 4); DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4)]
        @test DatesPlus.Hour(1) + t3 == [DatesPlus.DateTime(2009, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1), DatesPlus.DateTime(2009, 1, 3, 1)]
        @test t1 + DatesPlus.Day(1) == [DatesPlus.Date(2009, 1, 2) DatesPlus.Date(2009, 1, 3) DatesPlus.Date(2009, 1, 4); DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4)]
        @test t3 + DatesPlus.Hour(1) == [DatesPlus.DateTime(2009, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1), DatesPlus.DateTime(2009, 1, 3, 1)]
        @test t5 + DatesPlus.Hour(1) == [DatesPlus.Time(1, 0, 0) DatesPlus.Time(1, 0, 1) DatesPlus.Time(1, 0, 2); DatesPlus.Time(1, 1, 0) DatesPlus.Time(1, 2, 0) DatesPlus.Time(1, 3, 0)]
        @test (DatesPlus.Month(1) + DatesPlus.Day(1)) + t1 == [DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4); DatesPlus.Date(2009, 3, 2) DatesPlus.Date(2009, 3, 3) DatesPlus.Date(2009, 3, 4)]
        @test (DatesPlus.Hour(1) + DatesPlus.Minute(1)) + t3 == [DatesPlus.DateTime(2009, 1, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1, 1), DatesPlus.DateTime(2009, 1, 3, 1, 1)]
        @test t1 + (DatesPlus.Month(1) + DatesPlus.Day(1)) == [DatesPlus.Date(2009, 2, 2) DatesPlus.Date(2009, 2, 3) DatesPlus.Date(2009, 2, 4); DatesPlus.Date(2009, 3, 2) DatesPlus.Date(2009, 3, 3) DatesPlus.Date(2009, 3, 4)]
        @test t3 + (DatesPlus.Hour(1) + DatesPlus.Minute(1)) == [DatesPlus.DateTime(2009, 1, 1, 1, 1), DatesPlus.DateTime(2009, 1, 2, 1, 1), DatesPlus.DateTime(2009, 1, 3, 1, 1)]
        @test t1 - DatesPlus.Day(1) == [DatesPlus.Date(2008, 12, 31) DatesPlus.Date(2009, 1, 1) DatesPlus.Date(2009, 1, 2); DatesPlus.Date(2009, 1, 31) DatesPlus.Date(2009, 2, 1) DatesPlus.Date(2009, 2, 2)]
        @test t3 - DatesPlus.Hour(1) == [DatesPlus.DateTime(2008, 12, 31, 23), DatesPlus.DateTime(2009, 1, 1, 23), DatesPlus.DateTime(2009, 1, 2, 23)]
        @test t1 - (DatesPlus.Month(1) + DatesPlus.Day(1)) == [DatesPlus.Date(2008, 11, 30) DatesPlus.Date(2008, 12, 1) DatesPlus.Date(2008, 12, 2); DatesPlus.Date(2008, 12, 31) DatesPlus.Date(2009, 1, 1) DatesPlus.Date(2009, 1, 2)]
        @test t3 - (DatesPlus.Hour(1) + DatesPlus.Minute(1)) == [DatesPlus.DateTime(2008, 12, 31, 22, 59), DatesPlus.DateTime(2009, 1, 1, 22, 59), DatesPlus.DateTime(2009, 1, 2, 22, 59)]
        @test t2 - t1 == [DatesPlus.Day(1) DatesPlus.Day(31) DatesPlus.Day(365); DatesPlus.Day(365) DatesPlus.Day(28) DatesPlus.Day(1)]
        @test t4 - t3 == [DatesPlus.Millisecond(1000), DatesPlus.Millisecond(60000), DatesPlus.Millisecond(3600000)]
        @test (DatesPlus.Date(2009, 1, 1):DatesPlus.Week(1):DatesPlus.Date(2009, 1, 21)) - (DatesPlus.Date(2009, 1, 1):DatesPlus.Day(1):DatesPlus.Date(2009, 1, 3)) == [DatesPlus.Day(0), DatesPlus.Day(6), DatesPlus.Day(12)]
        @test (DatesPlus.DateTime(2009, 1, 1, 1, 1, 1):DatesPlus.Second(1):DatesPlus.DateTime(2009, 1, 1, 1, 1, 3)) - (DatesPlus.DateTime(2009, 1, 1, 1, 1):DatesPlus.Second(1):DatesPlus.DateTime(2009, 1, 1, 1, 1, 2)) == [DatesPlus.Second(1), DatesPlus.Second(1), DatesPlus.Second(1)]
        @test_throws MethodError DatesPlus.Day(1) - t1
        @test_throws MethodError DatesPlus.Hour(1) - t3
        @test_throws MethodError (DatesPlus.Month(1) + DatesPlus.Day(1)) - t1
        @test_throws MethodError (DatesPlus.Hour(1) + DatesPlus.Minute(1)) - t3
    end
    @testset "TimeZone" begin
        # best we can get in Dates as there is no other tz functionality
        @test ((a, b) -> now(typeof(a))).(UTC(), [1,2,3]) isa Vector{DateTime}
    end
end

@testset "Missing arithmetic" begin
    for t ∈ [Date, Time, Day, Month, Week, Year, Hour, Microsecond, Millisecond, Minute, Nanosecond, Second]
        @test ismissing(t(1) + missing)
        @test ismissing(missing + t(1))
        @test ismissing(t(1) - missing)
        @test ismissing(missing - t(1))
    end
end

@testset "Diff of dates" begin
    for t ∈ [Day, Week, Hour, Minute]
        a = DateTime(2021,1,1):t(1):DateTime(2021,2,1)
        d = diff(a)
        @test d == diff(collect(a))
        @test eltype(d) === typeof(a[1] - a[2])
    end
end

end


text:
function _field_as_vector_space(K::NumField, Q::FlintRationalField)
  BLoverK = absolute_basis(K)
  d = absolute_degree(K)
  m = identity_matrix(Q, d)
  return BLoverK, m
end

function _field_as_vector_space(K::NfRel{nf_elem}, Q::AnticNumberField)
  BLoverK = basis(K)
  d = degree(K)
  m = identity_matrix(Q, d)
  return BLoverK, m
end

function _field_as_vector_space(K::FqNmodFiniteField, Q::GaloisField)
  d = degree(K)
  BLoverK = powers(gen(K), d - 1)
  m = identity_matrix(Q, d)
  return BLoverK, m
end

function _field_as_vector_space(K::FqFiniteField, Q::GaloisFmpzField)
  d = degree(K)
  BLoverK = powers(gen(K), d - 1)
  m = identity_matrix(Q, d)
  return BLoverK, m
end

function _field_as_vector_space(f)
  K = domain(f)
  L = codomain(f)
  d = div(degree(L), degree(K))
  a = absolute_primitive_element(L)
  #BLoverK = elem_type(L)[a^i for i in 0:(d- 1)]
  BLoverK = powers(a, d - 1)
  BK = absolute_basis(K)
  absB = elem_type(L)[]
  for i in 1:d
    for j in 1:length(BK)
      push!(absB, BLoverK[i] * f(BK[j]))
    end
  end

  ad = absolute_degree(L)

  absbasismat = zero_matrix(FlintQQ, ad, ad)

  for i in 1:ad
    for j in 1:ad
      absbasismat[i, j] = absolute_coeff(absB[i], j - 1)
    end
  end

  absbasismatinv = inv(absbasismat)

  m = zero_matrix(K, ad, d)
  for i in 1:d
    for j in 1:degree(K)
      m[(i - 1)*degree(K) + j, i] = BK[j]
    end
  end

  return BLoverK, change_base_ring(K, absbasismatinv) * m
end

mutable struct FldToVecMor{R, S, T, U, V}
  L::R
  K::S
  f::T
  M::U
  B::V
  isone::Bool

  function FldToVecMor(L, f::NfToNfMor)
    return FldToVecMor(f)
  end

  function FldToVecMor(f::NfToNfMor)
    K = domain(f)
    L = codomain(f)
    B, M = _field_as_vector_space(f)
    B = B
    M = M
    z = new{typeof(L), typeof(K), typeof(f), typeof(M), typeof(B)}(L, K, f, M, B)
    z.isone = isone(M)
    return z
  end

  function FldToVecMor(L, Q)
    K = Q
    L = L
    B, M = _field_as_vector_space(L, Q)
    B = B
    M = M
    f = identity
    z = new{typeof(L), typeof(K), typeof(f), typeof(M), typeof(B)}(L, K, f, M, B)
    z.isone = isone(M)
    return z
  end

end

function image(f::FldToVecMor{T, FlintRationalField}, a::NumFieldElem) where {T <: NumField}
  @assert parent(a) == f.L
  L = parent(a)
  d = absolute_degree(L)
  K = f.K
  z = matrix(K, 1, d, elem_type(K)[K(absolute_coeff(a, i)) for i in 0:(d - 1)])
  if f.isone
    v = z
  else
    v = z * f.M
  end
  return elem_type(K)[v[1, i] for i in 1:ncols(v)]
end

function image(f::FldToVecMor{AnticNumberField, AnticNumberField}, a::nf_elem)
  @assert parent(a) == f.L
  L = parent(a)
  d = absolute_degree(L)
  K = f.K
  z = matrix(K, 1, d, elem_type(K)[K(absolute_coeff(a, i)) for i in 0:(d - 1)])
  if f.isone
    v = z
  else
    v = z * f.M
  end
  return elem_type(K)[v[1, i] for i in 1:ncols(v)]
end

function image(f::FldToVecMor{NfRel{nf_elem}, AnticNumberField}, a::NfRelElem{nf_elem})
  @assert parent(a) == f.L
  L = parent(a)
  d = degree(L)
  K = f.K
  z = matrix(K, 1, d, elem_type(K)[coeff(a, i) for i in 0:(d - 1)])
  if f.isone
    v = z
  else
    v = z * f.M
  end
  return elem_type(K)[v[1, i] for i in 1:ncols(v)]
end

function image(f::FldToVecMor{T}, a) where {T <: FinField}
  @assert parent(a) == f.L
  L = parent(a)
  d = degree(L)
  K = f.K
  z = matrix(K, 1, d, elem_type(K)[K(coeff(a, i)) for i in 0:(d - 1)])
  if f.isone
    v = z
  else
    v = z * f.M
  end
  return elem_type(K)[v[1, i] for i in 1:ncols(v)]
end

function preimage(f::FldToVecMor, v::Vector)
  @assert parent(v[1]) == f.K
  return dot(f.B, map(f.f, v))::elem_type(f.L)
end

ext: jl


text:
@testset "Move k backward" begin
    @testset "k=3" begin
        circular = false
        k = 3
        @test squaresub(1, circular, k) == 1
        @test squaresub(2, circular, k) == 1
        @test squaresub(3, circular, k) == 1
        @test squaresub(4, circular, k) == 1

        @test squaresub(11, circular, k) == 1
        @test squaresub(12, circular, k) == 2
        @test squaresub(13, circular, k) == 3
    end
    @testset "k=1" begin
        circular = false
        k = 1
        @test squaresub(1, circular, k) == 1
        @test squaresub(2, circular, k) == 1
        @test squaresub(3, circular, k) == 2
        @test squaresub(4, circular, k) == 3

        @test squaresub(11, circular, k) == 3
        @test squaresub(12, circular, k) == 11
        @test squaresub(13, circular, k) == 12
    end
    @testset "k=10" begin
        circular = false
        k = 10
        @test squaresub(1, circular, k) == 1
        @test squaresub(2, circular, k) == 1
        @test squaresub(3, circular, k) == 1
        @test squaresub(4, circular, k) == 1

        @test squaresub(11, circular, k) == 1
        @test squaresub(12, circular, k) == 1
        @test squaresub(13, circular, k) == 1
    end
end

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
using CartoonLike
points = create_circle(.4)
image, xLoc, yLoc = given_points(points, ninterp=100*100,ε=.001)
heatmap(xLoc, yLoc, image)

points = [-.5  -.9   0.0 .25   .5 .75;
          -.75 -.125 .9  .5  -.25 -.75]'
image, xLoc, yLoc = given_points(points, xSize=samples, ySize=samples)
heatmap(xLoc, yLoc, image')


text:
module day7

function readInput()
    rules = Dict{String,Vector{Tuple{Int,String}}}()
    for line in eachline(joinpath(@__DIR__, "input.txt"))
        tokens = split(line, r"(contain\s)|(,\s)|(\sbags?.?)", keepempty=false)
        key = first(tokens)
        val = Tuple{Int,String}[]
        for str in tokens[2:end]
            m = match(r"(\d+)\s(.*)", str)
            m === nothing && break
            nb, type = m.captures
            push!(val, (parse(Int, nb), type))
        end
        rules[key] = val
    end
    rules
end

function canContain(recipient, target, rules)
    isempty(rules[recipient]) && return false
    any(x -> last(x) == target, rules[recipient]) && return true
    return any(rec -> canContain(rec, target, rules), last.(rules[recipient]))
end

function nbBags(recipient, rules)
    mapreduce(((i, r),) -> (1 + nbBags(r, rules)) * i, +, rules[recipient] ; init=0)
end

part1(rules) = count(canContain.(keys(rules), "shiny gold", Ref(rules)))

part2(rules) = nbBags("shiny gold", rules)

end


text:
int_rules_8_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.1*Error*functions =#
    @apply_utils Antiderivative(erf(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * erf(~a + ~b * ~x)) / ~b + 1 / (~b * sqrt(pi) * ℯ ^ ((~a + ~b * ~x) ^ 2)) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(Erfc(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * Erfc(~a + ~b * ~x)) / ~b - 1 / (~b * sqrt(pi) * ℯ ^ ((~a + ~b * ~x) ^ 2)) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(erfi(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * erfi(~a + ~b * ~x)) / ~b - ℯ ^ ((~a + ~b * ~x) ^ 2) / (~b * sqrt(pi)) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(erf(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * erf(~a + ~b * ~x) ^ 2) / ~b - (4 / sqrt(pi)) * Antiderivative(((~a + ~b * ~x) * erf(~a + ~b * ~x)) / ℯ ^ ((~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(Erfc(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * Erfc(~a + ~b * ~x) ^ 2) / ~b + (4 / sqrt(pi)) * Antiderivative(((~a + ~b * ~x) * Erfc(~a + ~b * ~x)) / ℯ ^ ((~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(erfi(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * erfi(~a + ~b * ~x) ^ 2) / ~b - (4 / sqrt(pi)) * Antiderivative((~a + ~b * ~x) * ℯ ^ ((~a + ~b * ~x) ^ 2) * erfi(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(erf(~(a') + ~(b') * ~x) ^ ~n, ~x) => Unintegrable(erf(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~n], ~x) && (NeQ(~n, 1) && NeQ(~n, 2))
    @apply_utils Antiderivative(Erfc(~(a') + ~(b') * ~x) ^ ~n, ~x) => Unintegrable(Erfc(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~n], ~x) && (NeQ(~n, 1) && NeQ(~n, 2))
    @apply_utils Antiderivative(erfi(~(a') + ~(b') * ~x) ^ ~n, ~x) => Unintegrable(erfi(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~n], ~x) && (NeQ(~n, 1) && NeQ(~n, 2))
    @apply_utils Antiderivative(erf(~(b') * ~x) / ~x, ~x) => ((2 * ~b * ~x) / sqrt(pi)) * HypergeometricPFQ([1 / 2, 1 / 2], [3 / 2, 3 / 2], -((~b) ^ 2) * (~x) ^ 2) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative(Erfc(~(b') * ~x) / ~x, ~x) => log(~x) - Antiderivative(erf(~b * ~x) / ~x, ~x) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative(erfi(~(b') * ~x) / ~x, ~x) => ((2 * ~b * ~x) / sqrt(pi)) * HypergeometricPFQ([1 / 2, 1 / 2], [3 / 2, 3 / 2], (~b) ^ 2 * (~x) ^ 2) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * erf(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * erf(~a + ~b * ~x)) / (~d * (~m + 1)) - ((2 * ~b) / (sqrt(pi) * ~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) / ℯ ^ ((~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * Erfc(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * Erfc(~a + ~b * ~x)) / (~d * (~m + 1)) + ((2 * ~b) / (sqrt(pi) * ~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) / ℯ ^ ((~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * erfi(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * erfi(~a + ~b * ~x)) / (~d * (~m + 1)) - ((2 * ~b) / (sqrt(pi) * ~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) * ℯ ^ ((~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~(m') * erf(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * erf(~b * ~x) ^ 2) / (~m + 1) - ((4 * ~b) / (sqrt(pi) * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ (-((~b) ^ 2) * (~x) ^ 2) * erf(~b * ~x), ~x) <-- FreeQ(~b, ~x) && (IGtQ(~m, 0) || ILtQ((~m + 1) / 2, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * Erfc(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * Erfc(~b * ~x) ^ 2) / (~m + 1) + ((4 * ~b) / (sqrt(pi) * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ (-((~b) ^ 2) * (~x) ^ 2) * Erfc(~b * ~x), ~x) <-- FreeQ(~b, ~x) && (IGtQ(~m, 0) || ILtQ((~m + 1) / 2, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * erfi(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * erfi(~b * ~x) ^ 2) / (~m + 1) - ((4 * ~b) / (sqrt(pi) * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ ((~b) ^ 2 * (~x) ^ 2) * erfi(~b * ~x), ~x) <-- FreeQ(~b, ~x) && (IGtQ(~m, 0) || ILtQ((~m + 1) / 2, 0))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * erf(~a + ~(b') * ~x) ^ 2, ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand(erf(~x) ^ 2, ((~b * ~c - ~a * ~d) + ~d * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * Erfc(~a + ~(b') * ~x) ^ 2, ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand(Erfc(~x) ^ 2, ((~b * ~c - ~a * ~d) + ~d * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * erfi(~a + ~(b') * ~x) ^ 2, ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand(erfi(~x) ^ 2, ((~b * ~c - ~a * ~d) + ~d * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * erf(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * erf(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * Erfc(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * Erfc(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * erfi(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * erfi(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x) ^ ~(n'), ~x) => ((ℯ ^ ~c * sqrt(pi)) / (2 * ~b)) * Subst(Antiderivative((~x) ^ ~n, ~x), ~x, erf(~b * ~x)) <-- FreeQ([~b, ~c, ~d, ~n], ~x) && EqQ(~d, -((~b) ^ 2))
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x) ^ ~(n'), ~x) => ((-(ℯ ^ ~c) * sqrt(pi)) / (2 * ~b)) * Subst(Antiderivative((~x) ^ ~n, ~x), ~x, Erfc(~b * ~x)) <-- FreeQ([~b, ~c, ~d, ~n], ~x) && EqQ(~d, -((~b) ^ 2))
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x) ^ ~(n'), ~x) => ((ℯ ^ ~c * sqrt(pi)) / (2 * ~b)) * Subst(Antiderivative((~x) ^ ~n, ~x), ~x, erfi(~b * ~x)) <-- FreeQ([~b, ~c, ~d, ~n], ~x) && EqQ(~d, (~b) ^ 2)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x), ~x) => ((~b * ℯ ^ ~c * (~x) ^ 2) / sqrt(pi)) * HypergeometricPFQ([1, 1], [3 / 2, 2], (~b) ^ 2 * (~x) ^ 2) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ(~d, (~b) ^ 2)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x), ~x) => Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2), ~x) - Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ(~d, (~b) ^ 2)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x), ~x) => ((~b * ℯ ^ ~c * (~x) ^ 2) / sqrt(pi)) * HypergeometricPFQ([1, 1], [3 / 2, 2], -((~b) ^ 2) * (~x) ^ 2) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ(~d, -((~b) ^ 2))
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(~x * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(a') + ~(b') * ~x), ~x) => (ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x)) / (2 * ~d) - (~b / (~d * sqrt(pi))) * Antiderivative(ℯ ^ (((-((~a) ^ 2) + ~c) - 2 * ~a * ~b * ~x) - ((~b) ^ 2 - ~d) * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative(~x * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(a') + ~(b') * ~x), ~x) => (ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x)) / (2 * ~d) + (~b / (~d * sqrt(pi))) * Antiderivative(ℯ ^ (((-((~a) ^ 2) + ~c) - 2 * ~a * ~b * ~x) - ((~b) ^ 2 - ~d) * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative(~x * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(a') + ~(b') * ~x), ~x) => (ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x)) / (2 * ~d) - (~b / (~d * sqrt(pi))) * Antiderivative(ℯ ^ ((~a) ^ 2 + ~c + 2 * ~a * ~b * ~x + ((~b) ^ 2 + ~d) * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(a') + ~(b') * ~x), ~x) => (((~x) ^ (~m - 1) * ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x)) / (2 * ~d) - (~b / (~d * sqrt(pi))) * Antiderivative((~x) ^ (~m - 1) * ℯ ^ (((-((~a) ^ 2) + ~c) - 2 * ~a * ~b * ~x) - ((~b) ^ 2 - ~d) * (~x) ^ 2), ~x)) - ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 1)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(a') + ~(b') * ~x), ~x) => (((~x) ^ (~m - 1) * ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x)) / (2 * ~d) + (~b / (~d * sqrt(pi))) * Antiderivative((~x) ^ (~m - 1) * ℯ ^ (((-((~a) ^ 2) + ~c) - 2 * ~a * ~b * ~x) - ((~b) ^ 2 - ~d) * (~x) ^ 2), ~x)) - ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 1)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(a') + ~(b') * ~x), ~x) => (((~x) ^ (~m - 1) * ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x)) / (2 * ~d) - (~b / (~d * sqrt(pi))) * Antiderivative((~x) ^ (~m - 1) * ℯ ^ ((~a) ^ 2 + ~c + 2 * ~a * ~b * ~x + ((~b) ^ 2 + ~d) * (~x) ^ 2), ~x)) - ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 1)
    @apply_utils Antiderivative((ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x)) / ~x, ~x) => ((2 * ~b * ℯ ^ ~c * ~x) / sqrt(pi)) * HypergeometricPFQ([1 / 2, 1], [3 / 2, 3 / 2], (~b) ^ 2 * (~x) ^ 2) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ(~d, (~b) ^ 2)
    @apply_utils Antiderivative((ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x)) / ~x, ~x) => Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) / ~x, ~x) - Antiderivative((ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~b * ~x)) / ~x, ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ(~d, (~b) ^ 2)
    @apply_utils Antiderivative((ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x)) / ~x, ~x) => ((2 * ~b * ℯ ^ ~c * ~x) / sqrt(pi)) * HypergeometricPFQ([1 / 2, 1], [3 / 2, 3 / 2], -((~b) ^ 2) * (~x) ^ 2) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ(~d, -((~b) ^ 2))
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(a') + ~(b') * ~x), ~x) => (((~x) ^ (~m + 1) * ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x)) / (~m + 1) - ((2 * ~b) / ((~m + 1) * sqrt(pi))) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ (((-((~a) ^ 2) + ~c) - 2 * ~a * ~b * ~x) - ((~b) ^ 2 - ~d) * (~x) ^ 2), ~x)) - ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(a') + ~(b') * ~x), ~x) => (((~x) ^ (~m + 1) * ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x)) / (~m + 1) + ((2 * ~b) / ((~m + 1) * sqrt(pi))) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ (((-((~a) ^ 2) + ~c) - 2 * ~a * ~b * ~x) - ((~b) ^ 2 - ~d) * (~x) ^ 2), ~x)) - ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(a') + ~(b') * ~x), ~x) => (((~x) ^ (~m + 1) * ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x)) / (~m + 1) - ((2 * ~b) / ((~m + 1) * sqrt(pi))) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ ((~a) ^ 2 + ~c + 2 * ~a * ~b * ~x + ((~b) ^ 2 + ~d) * (~x) ^ 2), ~x)) - ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erf(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * Erfc(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * erfi(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative(erf(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * erf(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ((2 * ~b * ~d * ~n) / sqrt(pi)) * Antiderivative(1 / ℯ ^ ((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(Erfc(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * Erfc(~d * (~a + ~b * log(~c * (~x) ^ ~n))) + ((2 * ~b * ~d * ~n) / sqrt(pi)) * Antiderivative(1 / ℯ ^ ((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(erfi(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * erfi(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ((2 * ~b * ~d * ~n) / sqrt(pi)) * Antiderivative(ℯ ^ ((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~F)(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))) / ~x, ~x) => (1 / ~n) * Subst((~F)(~d * (~a + ~b * ~x)), ~x, log(~c * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && MemberQ([erf, Erfc, erfi], ~F)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * erf(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * erf(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((2 * ~b * ~d * ~n) / (sqrt(pi) * (~m + 1))) * Antiderivative((~e * ~x) ^ ~m / ℯ ^ ((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * Erfc(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * Erfc(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) + ((2 * ~b * ~d * ~n) / (sqrt(pi) * (~m + 1))) * Antiderivative((~e * ~x) ^ ~m / ℯ ^ ((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * erfi(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * erfi(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((2 * ~b * ~d * ~n) / (sqrt(pi) * (~m + 1))) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ ((~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x), ~x) => (im / 2) * Antiderivative(ℯ ^ (-im * ~c - im * ~d * (~x) ^ 2) * erf(~b * ~x), ~x) - (im / 2) * Antiderivative(ℯ ^ (im * ~c + im * ~d * (~x) ^ 2) * erf(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, -((~b) ^ 4))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x), ~x) => (im / 2) * Antiderivative(ℯ ^ (-im * ~c - im * ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) - (im / 2) * Antiderivative(ℯ ^ (im * ~c + im * ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, -((~b) ^ 4))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x), ~x) => (im / 2) * Antiderivative(ℯ ^ (-im * ~c - im * ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) - (im / 2) * Antiderivative(ℯ ^ (im * ~c + im * ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, -((~b) ^ 4))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (-im * ~c - im * ~d * (~x) ^ 2) * erf(~b * ~x), ~x) + (1 / 2) * Antiderivative(ℯ ^ (im * ~c + im * ~d * (~x) ^ 2) * erf(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, -((~b) ^ 4))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (-im * ~c - im * ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) + (1 / 2) * Antiderivative(ℯ ^ (im * ~c + im * ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, -((~b) ^ 4))
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (-im * ~c - im * ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) + (1 / 2) * Antiderivative(ℯ ^ (im * ~c + im * ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, -((~b) ^ 4))
    @apply_utils Antiderivative(sinh(~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~b * ~x), ~x) - (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ 2) * erf(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (~b) ^ 4)
    @apply_utils Antiderivative(sinh(~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) - (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (~b) ^ 4)
    @apply_utils Antiderivative(sinh(~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) - (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (~b) ^ 4)
    @apply_utils Antiderivative(cosh(~(c') + ~(d') * (~x) ^ 2) * erf(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf(~b * ~x), ~x) + (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ 2) * erf(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (~b) ^ 4)
    @apply_utils Antiderivative(cosh(~(c') + ~(d') * (~x) ^ 2) * Erfc(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) + (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ 2) * Erfc(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (~b) ^ 4)
    @apply_utils Antiderivative(cosh(~(c') + ~(d') * (~x) ^ 2) * erfi(~(b') * ~x), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) + (1 / 2) * Antiderivative(ℯ ^ (-(~c) - ~d * (~x) ^ 2) * erfi(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (~b) ^ 4)
    @apply_utils Antiderivative((~F)(~(f') * (~(a') + ~(b') * log(~(c') * (~d + ~(e') * ~x) ^ ~(n')))), ~x) => (1 / ~e) * Subst(Antiderivative((~F)(~f * (~a + ~b * log(~c * (~x) ^ ~n))), ~x), ~x, ~d + ~e * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && MemberQ([erf, Erfc, erfi, FresnelS, FresnelC, ExpIntegralEi, SinIntegral, CosIntegral, SinhIntegral, CoshIntegral], ~F)
    @apply_utils Antiderivative((~g + ~(h') * ~x) ^ ~(m') * (~F)(~(f') * (~(a') + ~(b') * log(~(c') * (~d + ~(e') * ~x) ^ ~(n')))), ~x) => (1 / ~e) * Subst(Antiderivative(((~g * ~x) / ~d) ^ ~m * (~F)(~f * (~a + ~b * log(~c * (~x) ^ ~n))), ~x), ~x, ~d + ~e * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~m, ~n], ~x) && (EqQ(~e * ~f - ~d * ~g, 0) && MemberQ([erf, Erfc, erfi, FresnelS, FresnelC, ExpIntegralEi, SinIntegral, CosIntegral, SinhIntegral, CoshIntegral], ~F))
end

max_stars_repo_path: src/rules/8 Special functions/8.1 Error functions.jl


text:
#=  Type definitions for  site  and  site_collection  and related functions for the infinite sites model.
A site stores statistics for a single site, while a site_collection stores statistics for all of the sites of a trial.
The update_sites1() function is called on every generation to update all sites in the site collection,

GLOBAL VARIABLES:  The global variables declared in this file are used only for correctness checking via the @assert macro.
  There are equivalent fields in site_collection, and results depend only on these fields, not on global variables.
  Thus, all global variables and global variable computations can be removed without affecting functionality.
=#
export site_type, site,
    site_collection, sc_push!, update_sites!, update_neutral, update_selected,
    print_summary, average_time_to_extinction, average_time_to_fixation,
    fixed_fraction, average_fitness_fixed, average_fitness_extinct, average_fitness_all,
    sites_per_gen, heterozygosity_per_gen

#=
Stores the properties of an "site" (possibly deleterious or advantageous).
=#

mutable struct site_type
  identifier::Int64   # integer for this site, not sure if this is needed
  start_gen::Int64    # the generation (time step) when the site was generated
  final_gen::Int64    # the generation when the site went extinct.  Should be zero while site is evolving
  fitness_coefficient::Float64    # The selection coefficient of the site (allele)
  sum_gens::Int64     # sum of the number of generations where counts and heterozygosity are accumulated
  sum_counts::Int64   # the current accumulated sum of the counts of the mutant allele per generation for infinite sites
  sum_sq_counts::Int64   # the current accumulated sum of squares of the counts per generation for infinite sites
  sum_heteroz::Float64  # the current accumulated sum of heterozygosities for infinite sites
  sum_sq_heteroz::Float64  # the current accumulated sum of squares ofheterozygosities for infinite sites
  previous_allele_freq::Int64  # Allele frequence of the previous generation
end

@doc """ function site()
  Initializes the properties of the site.
  Initial values for sum_counts, sum_sq_counts, sum_heteroz, sum_sq_heteroz are added.
  This means that the number of generations is one larger than the number of increments in update_innovations!(),
    and 1 is added to sum_gens in sites_per_gen() and heterozygosity_per_gen().
"""
function site( id::Int64, N::Int64, start_gen::Int64, fitness_coef::Float64=1.0 )
  # See the comments on global variables at the beginning of this file.
  global sum_counts += 1
  global sum_sq_counts += 1
  initial_heteroz = 1.0 - watterson_homozygosity([N-1,1])
  global sum_heteroz += initial_heteroz
  global sum_sq_heteroz += initial_heteroz^2
  global sum_gens += 1
  return site_type( id, start_gen, 0, fitness_coef, 1, 1, 1, initial_heteroz, initial_heteroz^2, 1  )
end

#=
Stores a collection of sites.
Sites are partitioned into 3 subsets: active, fixed, and extinct.
=#
mutable struct site_collection
  N::Int64                # popsize
  list::Dict{Int64,site_type} #    Stores a collection of sites
  active::BitSet   # indices of actively evolving sites
  fixed::BitSet   # indices of sites that have fixed
  extinct::BitSet   # indices of sites that have gone extinct
  fix_minimum::Float64  # minimum fraction of popsize for fixation for infinite alleles
  sum_gens::Int64     # sum of counts of generations
  sum_counts::Int64   # sum of counts of mutant alleles for infinite sites
  sum_sq_counts::Float64   # sum of counts of mutant alleles for infinite sites
  sum_heteroz::Float64  # sum of heterozygosities for infinite sites
  sum_sq_heteroz::Float64  # sum of heterozygosities squared for infinite sites
  sum_generations::Int64  # Total number of generations corresponding to sum_counts and sum_heteroz
  count_fixed_del::Int64  # Number of fixed mutations that are deleterious (fit_coef < 1.0)
  count_fixed_adv::Int64  # Number of fixed mutations that are advantageous (fit_coef > 1.0)
  in_use::Bool      # If false, not used
end

# Constructor for a new empty site collection with a value for fix_minimum
function site_collection( N::Int64, fix_min::Float64, in_use::Bool=true )
  global sum_counts = 0
  global sum_sq_counts = 0
  global sum_heteroz = 0.0
  global sum_sq_heteroz = 0.0
  global sum_gens = 0
  site_collection( N, Dict{Int64,site_type}(), BitSet(), BitSet(), BitSet(), fix_min, 0, 0, 0.0, 0.0, 0.0, 0,0,0, in_use )
end

@doc """ sc_push!()
 Adds an site to the collection.
"""
function sc_push!( innov_collection::site_collection, innov::site_type )
  if !innov_collection.in_use
    return
  end
  innov_collection.list[innov.identifier] = innov
  if innov.final_gen == 0
    Base.push!( innov_collection.active, innov.identifier )
  elseif innov.history[end] > 0
    Base.push!( innov_collection.fixed, innov.identifier )
  else
    Base.push!( innov_collection.extinct, innov.identifier )
  end
  innov_collection
end

@doc """ fix_test()
The most basic test for fixation, namely that the allele frequency is the population size.
This is used in the infitite sites model.
"""
function fix_test( N::Int64, new_allele_freq::Int64 )
  return N == new_allele_freq
end

@doc """ function update_sites!( )
 Update all active sites, and make some extinct and fixed
 This version is used only by infsites.jl.
"""
function update_sites!( sc::site_collection, g::Int64, N::Int64 )
  if !sc.in_use
    return
  end
  for index in sc.active  # updates sites to the next generation
    # update_selected does the equivalent of proportional selection
    new_allele_freq = update_selected( index, N, sc.list[index].previous_allele_freq, sc.list[index].fitness_coefficient )
    innov = sc.list[index]
    if new_allele_freq > 0 && new_allele_freq < N
      global sum_counts += new_allele_freq
      global sum_sq_counts += new_allele_freq^2
      innov.sum_gens += 1
      innov.sum_counts += new_allele_freq
      innov.sum_sq_counts += new_allele_freq^2
      heteroz =  1.0 - watterson_homozygosity([N-new_allele_freq, new_allele_freq])
      global sum_heteroz += heteroz
      global sum_sq_heteroz += heteroz^2
      global sum_gens +=1
      innov.sum_heteroz += heteroz
      innov.sum_sq_heteroz += heteroz^2
    end
    innov.previous_allele_freq = new_allele_freq  # save for the next call to update_selected.
    if new_allele_freq == 0 ||  new_allele_freq >= N
      sc.sum_gens += innov.sum_gens
      sc.sum_counts += innov.sum_counts
      sc.sum_sq_counts += innov.sum_sq_counts
      sc.sum_heteroz += innov.sum_heteroz
      sc.sum_sq_heteroz += innov.sum_sq_heteroz
      sc.sum_generations += innov.sum_gens
      Base.pop!( sc.active,index)
      innov.final_gen = g
      if new_allele_freq == 0  # extinction
        Base.push!( sc.extinct,index)
      elseif new_allele_freq >= N # fixation
        Base.push!( sc.fixed,index)
      end
    end
  end
end

@doc """ function update_neutral()
  Updates site according to the Wright-Fisher model of drift.
  Since mutation has already happened by the time a new site is established, there is no mutation.
"""
function update_neutral( site::Int64, N::Int64, old_allele_freq::Int64 )
  p = Float64(old_allele_freq)/N
  new_allele_freq = rand(Binomial(N,p))
  return new_allele_freq
end

@doc """ function update_selected()
  Updates site according to the Wright-Fisher model of drift with selection.
  Since mutation has already happened by the time a new site is established, there is no mutation.
  This function is called by update_sites!() in site_collection.jl.

"""
function update_selected( site::Int64, N::Int64, old_allele_freq::Int64, select_coef::Float64 )
  p = min(1.0,(select_coef*Float64(old_allele_freq)/N))
  new_allele_freq = rand(Binomial(N,p))
  return new_allele_freq
end

function N_inf_sites( sc::site_collection )
  if !sc.in_use
    return
  end
  sum_N = 0
  for index in sc.active
    sum_N += sc.list[index].history[end]
  end
  sum_N
end

@doc """ check_globals_against_site_collection()
  Check that the values of global variables approximately agree with the corresponding variables in sc.
"""
function check_globals_against_site_collection( sc::site_collection )
  global sum_counts
  global sum_sq_counts
  global sum_heteroz
  global sum_sq_heteroz
  global sum_gens
  @assert isapprox(sc.sum_gens,sum_gens)
  @assert isapprox(sc.sum_counts,sum_counts)
  @assert isapprox(sc.sum_sq_counts,sum_sq_counts)
  @assert isapprox(sc.sum_heteroz,sum_heteroz)
  @assert isapprox(sc.sum_sq_heteroz,sum_sq_heteroz)
end

function print_summary( sc::site_collection; print_lists::Bool=false )
  println("print summary of site collection: ")
  global sum_counts
  global sum_sq_counts
  global sum_heteroz
  global sum_sq_heteroz
  global sum_gens

  #global sum_counts
  #println("psum sum_counts: ",sum_counts)
  if !sc.in_use
    return
  end
  if print_lists
    println("active list: ",sc.active)
    println("fixed list: ",sc.fixed)
    println("extinct list: ",sc.extinct)
  end
  println("number active: ",length(sc.active))
  println("number fixed: ",length(sc.fixed))
  println("number extinct: ",length(sc.extinct))
  println("fixed fraction: ",fixed_fraction(sc))
  println("sum_counts: ",sc.sum_counts,"  glb: ",sum_counts)
  println("sum_sq_counts: ",sc.sum_sq_counts,"  glb: ",sum_sq_counts)
  println("sum_heteroz: ",sc.sum_heteroz,"  glb: ",sum_heteroz)
  println("sum_sq_heteroz: ",sc.sum_sq_heteroz,"  glb: ",sum_sq_heteroz)
  println("sum_gens: ",sc.sum_generations,"  glb: ",sum_gens)
  (avg,stderr) = sites_per_gen(sc)
  println("avg sites per generation: ",avg)
  println("stderr sites per generation: ",stderr)
  (avg,stderr) = heterozygosity_per_gen(sc)
  println("avg per generation heterozygosity: ",avg)
  println("stderr per generation heterozygosity: ",stderr)
  println("avg time to fixation: ",average_time_to_fixation(sc))
  println("avg time to extinction: ",average_time_to_extinction(sc))
  println("avg fitness fixed: ",average_fitness_fixed(sc))
  println("avg fitness extinct: ",average_fitness_extinct(sc))
  println("avg fitness all: ",average_fitness_all(sc))
end

function average_time_to_extinction( innov_collection::site_collection )
  if !innov_collection.in_use
    return 0.0
  end
  if length(innov_collection.extinct) == 0
    println("no extinct sites")
    return -1.0
  end
  #print("average_time_to_extinction:  ")
  sum = 0
  for i in innov_collection.extinct
    sum += innov_collection.list[i].final_gen - innov_collection.list[i].start_gen
    #println("id: ",innov_collection.list[i].identifier,"  time: ",(innov_collection.list[i].final_gen-innov_collection.list[i].start_gen),"  sum: ",sum)
  end
  return Float64(sum)/length(innov_collection.extinct)
end

function average_time_to_fixation( innov_collection::site_collection )
  if !innov_collection.in_use
    return 0.0
  end
  if length(innov_collection.fixed) == 0
    println("no fixed sites")
    return -1.0
  end
  sum = 0
  for i in innov_collection.fixed
    sum += innov_collection.list[i].final_gen - innov_collection.list[i].start_gen
  end
  return Float64(sum)/length(innov_collection.fixed)
end

# Fraction of fixed sites out of fixed plus extinct sites
function fixed_fraction( innov_collection::site_collection )
  if !innov_collection.in_use
    return 0.0
  end
  if  length(innov_collection.extinct) == 0 && length(innov_collection.fixed) == 0
    println("no extinct and no fixed sites")
    return -1.0
  end
  return length(innov_collection.fixed)/Float64(length(innov_collection.fixed)+length(innov_collection.extinct))
end

function average_fitness_fixed( sc::site_collection )
  if !sc.in_use
    return 0.0
  end
  if  length(sc.fixed) == 0
    return -1.0
  end
  sum = 0.0
  for i in sc.fixed
    sum += sc.list[i].fitness_coefficient
  end
  return sum/length(sc.fixed)
end

function average_fitness_extinct( sc::site_collection )
  if !sc.in_use
    return 0.0
  end
  if  length(sc.extinct) == 0
    return -1.0
  end
  sum = 0.0
  for i in sc.extinct
    sum += sc.list[i].fitness_coefficient
  end
  return sum/length(sc.extinct)
end

function average_fitness_all( sc::site_collection )
  if !sc.in_use
    return 0.0
  end
  if  length(sc.extinct) == 0
    return -1.0
  end
  sum = 0.0
  for i in sc.extinct
    sum += sc.list[i].fitness_coefficient
  end
  for i in sc.fixed
    sum += sc.list[i].fitness_coefficient
  end
  return sum/(length(sc.extinct) + length(sc.fixed))
end

@doc """ function count_adv_del_fixed( )
  Counts the number of fixed sites that are advantageous, and the number that are deleterious at the end of the run.
  Stores these in site collection sc, and also returns them as a pair.
"""
function count_adv_del_fixed( sc::site_collection )
  #println("count_adv_del_fixed")
  if !sc.in_use
    return 0.0
  end
  if  length(sc.fixed) == 0
    return 0, 0
  end
  count_adv = 0
  count_del = 0
  for i in sc.fixed
    fit_coef = sc.list[i].fitness_coefficient
    #println("i: ",i,"  fit_coef: ",fit_coef)
    if fit_coef > 1.0
      count_adv += 1
    end
    if fit_coef < 1.0
      count_del += 1
    end
  end
  sc.count_fixed_adv = count_adv
  sc.count_fixed_del = count_del
  return count_adv, count_del
end

function sites_per_gen( sc::site_collection )
  global sum_gens
  global sum_counts
  global sum_sq_counts
  global counts_var = sqrt(sum_sq_counts/(sum_gens-1) - sum_counts^2/sum_gens/(sum_gens-1))/sum_gens
  sumgens = sc.sum_generations
  average = sc.sum_counts/sumgens
  coef_var = sqrt(sc.sum_sq_counts/(sumgens-1) - sc.sum_counts^2/sumgens/(sumgens-1))/sumgens
  @assert isapprox( coef_var, counts_var )
  if coef_var >= 0
    stderr = coef_var
  else
    stderr = -1.0
  end
  return (average, stderr)
end

function heterozygosity_per_gen( sc::site_collection )
  global sum_gens
  global sum_heteroz
  global sum_sq_heteroz
  heteroz_var = sqrt(sum_sq_heteroz/(sum_gens-1) - sum_heteroz^2/sum_gens/(sum_gens-1))/sum_gens
  sumgens = sc.sum_generations
  average = sc.sum_heteroz/sumgens
  coef_var = sqrt(sc.sum_sq_heteroz/(sumgens-1) - sc.sum_heteroz^2/sumgens/(sumgens-1))/sumgens
  @assert isapprox( coef_var, heteroz_var )
  if coef_var >= 0
    stderr = coef_var
  else
    stderr = -1.0
  end
  return (average, stderr)
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
#Kacper Kulczak 279079
#Program do zadania 2.12 na pracownie z przedmiotu Analiza Numeryczna M

#Odwracanie listy

for i = 5:-1:0
  println(i)
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
function createQuadraticModel128(qpdata; name = "qp_pb")
  return QuadraticModel(
    convert(Array{Float128}, qpdata.c),
    qpdata.qrows,
    qpdata.qcols,
    convert(Array{Float128}, qpdata.qvals),
    Arows = qpdata.arows,
    Acols = qpdata.acols,
    Avals = convert(Array{Float128}, qpdata.avals),
    lcon = convert(Array{Float128}, qpdata.lcon),
    ucon = convert(Array{Float128}, qpdata.ucon),
    lvar = convert(Array{Float128}, qpdata.lvar),
    uvar = convert(Array{Float128}, qpdata.uvar),
    c0 = Float128(qpdata.c0),
    x0 = zeros(Float128, length(qpdata.c)),
    name = name,
  )
end

@testset "multi_mode" begin
  stats1 = ripqp(QuadraticModel(qps1), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  stats2 = ripqp(QuadraticModel(qps2), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable

  stats3 = ripqp(QuadraticModel(qps3), iconf = InputConfig(mode = :multi), display = false)
  @test isapprox(stats3.objective, 5.32664756, atol = 1e-2)
  @test stats3.status == :acceptable
end

@testset "Float128" begin
  qm128_1 = createQuadraticModel128(qps1)
  stats1 = ripqp(
    qm128_1,
    itol = InputTol(ϵ_rb32 = 0.1, ϵ_rb64 = 0.01),
    iconf = InputConfig(mode = :multi, normalize_rtol = false),
    display = false,
  )
  @test isapprox(stats1.objective, -1.59078179, atol = 1e-2)
  @test stats1.status == :acceptable

  qm128_2 = createQuadraticModel128(qps2)
  stats2 = ripqp(qm128_2, display = false)
  @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-2)
  @test stats2.status == :acceptable
end


text:
# Command: Status
"""
	status(wd::WebDriver)::Bool
Returns `true` if the server is running and ready to accept new sessions.
"""
function status(wd::RemoteWebDriver)::Bool
    @unpack addr = wd
    response = HTTP.get("$addr/status", [("Content-Type" => "application/json")])
    @assert response.status == 200
    json = JSON3.read(response.body).value
    json.ready && json.message == "Server is running"
end


text:
@testset "get_coefficients.jl" begin

end


text:
### Exports ###

export Layer, GenericLayer

### Types ###

struct Layer{S,N} <: AbstractNode
  synapses::S
  neurons::N
end

const GenericLayer = Layer{GenericSynapses,GenericNeurons}
function GenericLayer(sz::Int)
  Layer(GenericSynapses(sz), GenericNeurons(sz))
end

### Methods ###

addedge!(layer::Layer, dstcont, dst, op, conf) =
  addedge!(layer.synapses, dstcont, dst, op, conf)
deledge!(layer::Layer, dst, op) =
  deledge!(layer.synapses, dst, op)
deledges!(layer::Layer) =
  deledges!(layer.synapses)

function reinit!(layer::Layer)
  reinit!(layer.synapses)
  reinit!(layer.neurons)
end

Base.size(layer::Layer) = size(layer.neurons)

Base.getindex(layer::Layer, idx...) =
  getindex(layer.neurons, idx...)

# ==(layer1::Layer, layer2::Layer) =
#   (layer1.synapses == layer2.synapses) && (layer1.neurons == layer2.neurons)

const ConnWrapper = Tuple{D1, CPUContainer{D2}, Dict{Symbol,Any}} where {D1,D2}

function eforward!(lcont::CPUContainer{Layer{S,N}}, args) where {S,N}
  layer = root(lcont)
  synapses, neurons = layer.synapses, layer.neurons

  # TODO: Sort conns and args by UUID for speed?
  conns = ConnWrapper[]
  for conn in connections(synapses)
    idx = findfirst(arg->arg[2]==conn.uuid, args)
    @assert idx !== nothing
    cont = args[idx][3]
    inputs = cont[:]
    state = Dict{Symbol,Any}()
    state[:inputs] = inputs
    push!(conns, (conn.data, cont, state))
  end

  # Reset run-local variables
  reset_run!(synapses)

  # Initialize connection-local states
  for conn in conns
    initialize_state!(synapses, conn)
  end

  # Get inputs and shift frontends
  for conn in conns
    shift_frontend!(synapses, conn)
  end

  # Early modulation
  for conn in conns
    early_modulate!(synapses, conn, conns)
  end

  # Calculate outputs
  for conn in conns
    calculate_outputs!(synapses, conn)
  end

  # Late modulation
  for conn in conns
    late_modulate!(synapses, conn, conns)
  end

  # Cycle neurons
  for idx in 1:size(neurons)
    # FIXME: Use dispatch properly
    neurons.state.I[idx] = synapses[idx]
  end
  cycle_neurons!(neurons)

  # TODO: Learning modulation?

  # Learn weights
  for conn in conns
    learn_weights!(synapses, conn, neurons)
  end
end

initialize_state!(synapses::GenericSynapses, conn::ConnWrapper) =
  initialize_state!(synapses, conn[1], conn[3])

shift_frontend!(synapses::GenericSynapses, conn::ConnWrapper) =
  shift_frontend!(synapses, conn[1], conn[3][:inputs])

early_modulate!(synapses::GenericSynapses, conn::ConnWrapper, conns) =
  early_modulate!(synapses, conn[1].modulator, conn, conns)

calculate_outputs!(synapses::GenericSynapses, conn::ConnWrapper) =
  calculate_outputs!(synapses, conn[1], conn[3][:inputs])

late_modulate!(synapses::GenericSynapses, conn::ConnWrapper, conns) =
  late_modulate!(synapses, conn[1].modulator, conn, conns)


text:
# To disable debug assertions, change DASSERT_LEVEL to `0` in your environment
# and re-compile.  To enable more aggressive debug assertions, change
# DASSERT_LEVEL to `2` (or higher).

function _assert_level()
    if haskey(ENV, "DASSERT_LEVEL")
        return parse(Int, ENV["DASSERT_LEVEL"])
    else
        return 1
    end
end

const DASSERT_LEVEL = _assert_level()

macro dassert1(ex, msgs...)
    if DASSERT_LEVEL ≥ 1
        esc(:(@assert($ex, $(msgs...))))
    else
        :(nothing)
    end
end

export @dassert1


text:
using Stheno: GPC, EQ, Exp

@timedtestset "compose" begin
    @timedtestset "general" begin
        rng, N, N′, gpc = MersenneTwister(123456), 5, 3, GPC()
        x, x′ = randn(rng, N), randn(rng, N′)
        f, g, h = GP(sin, EQ(), gpc), cos, GP(exp, Exp(), gpc)
        fg = f ∘ g

        # Check marginals statistics inductively.
        @test mean(fg(x)) == mean(f(map(g, x)))
        @test cov(fg(x)) == cov(f(map(g, x)))

        # Check cross covariance between transformed process and original inductively.
        @test cov(fg(x), fg(x′)) == cov(f(map(g, x)), f(map(g, x′)))
        @test cov(fg(x), f(x′)) == cov(f(map(g, x)), f(x′))
        @test cov(f(x), fg(x′)) == cov(f(x), f(map(g, x′)))

        # Check cross covariance between transformed process and independent absolutely.
        @test cov(fg(x), h(x′)) == zeros(length(x), length(x′))
        @test cov(h(x), fg(x′)) == zeros(length(x), length(x′))

        @timedtestset "Consistency Tests" begin
            P, Q = 4, 3
            x0, x1, x2, x3 = randn(rng, P), randn(rng, Q), randn(rng, Q), randn(rng, P)
            abstractgp_interface_tests(fg, f, x0, x1, x2, x3)
            abstractgp_interface_tests(stretch(f, 0.1), f, x0, x1, x2, x3)

            # f = GP(EQ(), GPC())
            # abstractgp_interface_tests(periodic(f, 0.1), f, x0, x1, x2, x3)
        end
        @timedtestset "Diff Tests" begin
            standard_1D_tests(
                MersenneTwister(123456),
                Dict(:σ=>0.5),
                θ->begin
                    f = θ[:σ] * GP(sin, EQ(), GPC())
                    return stretch(f, 0.5), f
                end,
                collect(range(-2.0, 2.0; length=N)),
                collect(range(-1.5, 2.2; length=N′)),
            )
        end
    end
    @timedtestset "Stretch" begin
        @timedtestset "scalar stretch" begin
            rng, N, λ = MersenneTwister(123456), 3, 0.51
            x = randn(rng)
            f = GP(1.3, EQ(), GPC())
            g = stretch(f, λ)

            @timedtestset "scalar input" begin
                @test first(cov(f, g, [0.0], [0.0])) == 1.0
                @test first(cov(f, g, [λ * x], [x])) == 1.0
                standard_1D_tests(
                    MersenneTwister(123456),
                    Dict(:σ=>0.5, :l=>0.32),
                    θ->begin
                        f_ = θ[:σ] * GP(sin, EQ(), GPC())
                        return stretch(f_, θ[:l]), f_
                    end,
                    collect(range(-2.0, 2.0; length=N)),
                    collect(range(-1.5, 2.2; length=5)),
                )
            end
            @timedtestset "vector input" begin
                D = 11
                X = randn(rng, D, 1)
                @test first(cov(f, g, ColVecs(zeros(D, 1)), ColVecs(zeros(D, 1)))) == 1.0
                @test first(cov(f, g, ColVecs(λ .* X), ColVecs(X))) == 1.0
            end
        end
        @timedtestset "Vector stretch" begin
            rng, N, D = MersenneTwister(123456), 3, 7
            λ = randn(rng, D)
            f = GP(1.0, EQ(), GPC())
            g = stretch(f, λ)

            X = randn(rng, D, 1)
            @test first(cov(f, g, ColVecs(zeros(D, 1)), ColVecs(zeros(D, 1)))) == 1.0
            @test first(cov(f, g, ColVecs(Diagonal(λ) * X), ColVecs(X))) == 1.0
        end
        @timedtestset "Matrix stretch" begin
            rng, N, D = MersenneTwister(123456), 3, 7
            A = randn(rng, D, D)
            f = GP(1.0, EQ(), GPC())
            g = stretch(f, A)

            X = randn(rng, D, 1)
            @test first(cov(f, g, ColVecs(zeros(D, 1)), ColVecs(zeros(D, 1)))) == 1.0
            @test first(cov(f, g, ColVecs(A * X), ColVecs(X))) == 1.0
        end
    end
    @timedtestset "Select" begin
        rng, N, D = MersenneTwister(123456), 3, 6
        idx = [1, 3]
        f = GP(1.3, EQ(), GPC())
        g = select(f, idx)

        X = randn(rng, D, N)
        X_f = ColVecs(X[idx, :])
        X_g = ColVecs(X)
        @test cov(f, g, X_f, X_g) ≈ cov(f, X_f, X_f)
        @test cov(f, g, X_f, X_g) ≈ cov(g, X_g, X_g)
    end
    @timedtestset "Shift" begin
        @timedtestset "Shift{Float64}" begin
            rng, N, D = MersenneTwister(123456), 3, 6
            a = randn(rng)
            f = GP(1.3, EQ(), GPC())
            g = shift(f, a)

            x = randn(rng, N)
            x_f = x .- a
            x_g = x
            @test cov(f, g, x_f, x_g) ≈ cov(f, x_f, x_f)
            @test cov(f, g, x_f, x_g) ≈ cov(g, x_g, x_g)

            X = randn(rng, D, N)
            X_f = ColVecs(X .- a)
            X_g = ColVecs(X)
            @test cov(f, g, X_f, X_g) ≈ cov(f, X_f, X_f)
            @test cov(f, g, X_f, X_g) ≈ cov(g, X_g, X_g)
        end
        @timedtestset "Shift{Vector{Float64}}" begin
            rng, N, D = MersenneTwister(123456), 3, 6
            a = randn(rng, D)
            f = GP(1.3, EQ(), GPC())
            g = shift(f, a)

            X = randn(rng, D, N)
            X_f = ColVecs(X .- a)
            X_g = ColVecs(X)
            @test cov(f, g, X_f, X_g) ≈ cov(f, X_f, X_f)
            @test cov(f, g, X_f, X_g) ≈ cov(g, X_g, X_g)
        end
    end
end