text:
# using Base.Test
# using MIToS.MSA
# using MIToS.Information

let s = res"ARNDCQEGHILKMFPSTWYV-",
  r = reverse(s),
  g = res"GGGGGGGGGGGGGGGGGGGG",
  Pg = probabilities(BigFloat, g),
  Pgg = probabilities(BigFloat, g,g),
  Pggg = probabilities(BigFloat, g,g,g),
  Ps = probabilities(BigFloat, s),
  Pss = probabilities(BigFloat, s,s),
  Psr = probabilities(BigFloat, s,r),
  Psss = probabilities(BigFloat, s,s,s),
  prob_random = probabilities(BigFloat, AdditiveSmoothing(1.0), Residue[], Residue[]),
  count_random = count(AdditiveSmoothing(1.0), Residue[], Residue[])

  print("""

Tests for Information Measures
==============================
""")

  print("""

Entropy
-------
""")

  print("""
Entropy: H(X)
""")

  @test estimate(Entropy{BigFloat}(), Pg) == zero(BigFloat)
  @test_approx_eq estimate(Entropy{BigFloat}(), Ps) log(big"20")
  @test_approx_eq estimate(Entropy{BigFloat}(2), Ps) log(2, big"20")

  print("""
Joint Entropy: H(X,Y)
""")

  @test estimate(Entropy{BigFloat}(), Pgg) == zero(BigFloat)
  @test_approx_eq estimate(Entropy{BigFloat}(), Pss) log(big"20")
  @test_approx_eq estimate(Entropy{BigFloat}(), Psr) log(big"19")
  @test_approx_eq estimate(Entropy{BigFloat}(2), Pss) log(2, big"20")
  @test_approx_eq estimate(Entropy{BigFloat}(), count_random) log(big"400")
  @test_approx_eq estimate(Entropy{BigFloat}(),  prob_random) log(big"400")

  print("""
Joint Entropy: H(X,Y,Z)
""")

  @test estimate(Entropy{BigFloat}(), Pggg) == zero(BigFloat)
  @test_approx_eq estimate(Entropy{BigFloat}(), Psss) log(big"20")

  print("""
Entropy using ResidueCount
""")
  @test estimate(Entropy{BigFloat}(), count(g,g)) == estimate(Entropy{BigFloat}(), Pgg)
  @test estimate(Entropy{BigFloat}(), count(s,s)) == estimate(Entropy{BigFloat}(), Pss)
  @test_approx_eq estimate(Entropy{BigFloat}(), count(s,r)) estimate(Entropy{BigFloat}(), Psr)

  print("""

Kullback-Leibler
----------------
""")

  @test estimate(Information.KullbackLeibler{Float64}(Float64[1/20 for i in 1:20]), Ps) == 0.0
  @test estimate(Information.KullbackLeibler{Float64}(BLOSUM62_Pi), Ps) != 0.0
  @test_approx_eq estimate(Information.KullbackLeibler{Float64}(BLOSUM62_Pi), Ps) sum( map(i -> 0.05 * log(0.05 / BLOSUM62_Pi[i]), 1:20) )

  @test estimate(Information.KullbackLeibler{BigFloat}(Ps), Ps) == zero(BigFloat)
  @test estimate(Information.KullbackLeibler{BigFloat}(Psr), Psr) == zero(BigFloat)

  print("""

Mutual Information
------------------
""")

  @test estimate(MutualInformation{BigFloat}(), Pgg) == zero(BigFloat)
  @test estimate(MutualInformation{BigFloat}(), probabilities(BigFloat, g,s)) == zero(BigFloat)
  @test_approx_eq estimate(MutualInformation{BigFloat}(),  prob_random) zero(BigFloat) # diff for Float64 is 1e-15

  print("""
MI(X,Y) = H(X) + H(Y) - H(X,Y)
""")
  @test estimate(MutualInformation{BigFloat}(), Psr) == estimate_on_marginal(Entropy{BigFloat}(), Psr, 1) + estimate_on_marginal(Entropy{BigFloat}(), Psr, 2) - estimate(Entropy{BigFloat}(), Psr)
  @test estimate(MutualInformation{BigFloat}(), Pss) == estimate_on_marginal(Entropy{BigFloat}(), Pss, 1) + estimate_on_marginal(Entropy{BigFloat}(), Pss, 2) - estimate(Entropy{BigFloat}(), Pss)
  @test estimate(MutualInformation{BigFloat}(2), Psr) == estimate_on_marginal(Entropy{BigFloat}(2), Psr, 1) + estimate_on_marginal(Entropy{BigFloat}(2), Psr, 2) - estimate(Entropy{BigFloat}(2), Psr)
  @test estimate(MutualInformation{BigFloat}(20), Pss) == estimate_on_marginal(Entropy{BigFloat}(20), Pss, 1) + estimate_on_marginal(Entropy{BigFloat}(20), Pss, 2) - estimate(Entropy{BigFloat}(20), Pss)
  @test_approx_eq_eps estimate(MutualInformation{BigFloat}(), prob_random) ( estimate_on_marginal(Entropy{BigFloat}(), prob_random, 1) +
                                                                              estimate_on_marginal(Entropy{BigFloat}(),  prob_random, 2) -
                                                                              estimate(Entropy{BigFloat}(), prob_random) ) 1e-74 # diff for Float64 is 1e-15

  print("""
MI using ResidueCount
""")
  @test_approx_eq estimate(MutualInformation{BigFloat}(), count_random) zero(BigFloat)
  @test estimate(MutualInformation{BigFloat}(), count(g,g)) == estimate(MutualInformation{BigFloat}(), Pgg)
  @test estimate(MutualInformation{BigFloat}(), count(s,s)) == estimate(MutualInformation{BigFloat}(), Pss)
  @test_approx_eq estimate(MutualInformation{BigFloat}(), count(s,r)) estimate(MutualInformation{BigFloat}(), Psr)

  print("""
MI(X,Y,Z)
""")
  @test estimate(MutualInformation{BigFloat}(), Pggg) == zero(BigFloat)
  @test estimate(MutualInformation{BigFloat}(), probabilities(BigFloat, g,s,r)) == zero(BigFloat)
  @test_approx_eq estimate(MutualInformation{BigFloat}(), probabilities(BigFloat, s,s,r)) estimate(MutualInformation{BigFloat}(), count(s,s,r))

  print("""
MI(X,Y,Z) = H(X) + H(Y) + H(Z) - H(X,Y) - H(X,Z) - H(Y,Z) + H(X,Y,Z)
""")
  @test estimate(MutualInformation{BigFloat}(), Psss) == ( estimate_on_marginal(Entropy{BigFloat}(), Psss, 1) + estimate_on_marginal(Entropy{BigFloat}(), Psss, 2) + estimate_on_marginal(Entropy{BigFloat}(), Psss, 3) -
                                                            estimate(Entropy{BigFloat}(), Pss) - estimate(Entropy{BigFloat}(), Pss) - estimate(Entropy{BigFloat}(), Pss) + estimate(Entropy{BigFloat}(), Psss) )

  print("""
MI(X,Y,Z) <= min{ H(X,Y), H(X,Z), H(Y,Z) }
""")
  @test_approx_eq estimate(MutualInformation{BigFloat}(), Psss) estimate(MutualInformation{BigFloat}(),Pss)

  print("""

Pairwise Gap Percentage
-----------------------
""")

  @test estimate(GapUnionPercentage{Float64}(),        count(res"AA--", res"--AA", usegap=true)) == 100.0
  @test estimate(GapIntersectionPercentage{Float64}(), count(res"AA--", res"--AA", usegap=true)) == 0.0

  @test estimate(GapUnionPercentage{Float64}(),        count(res"AA--", res"--AA", usegap=true, weight=Float64[.25, .25, .25, .25])) == 100.0
  @test estimate(GapIntersectionPercentage{Float64}(), count(res"AA--", res"--AA", usegap=true, weight=Float64[.25, .25, .25, .25])) == 0.0

  @test estimate(GapUnionPercentage{Float64}(),        count(res"AAA-", res"AA--", usegap=true)) == 50.0
  @test estimate(GapIntersectionPercentage{Float64}(), count(res"AAA-", res"AA--", usegap=true)) == 25.0

  @test estimate(GapUnionPercentage{Float64}(),        count(res"AAA-", res"AA--", usegap=true, weight=Float64[.2, .2, .2, .4])) == 60.0
  @test estimate(GapIntersectionPercentage{Float64}(), count(res"AAA-", res"AA--", usegap=true, weight=Float64[.2, .2, .2, .4])) == 40.0
end


text:
"""
    Basis

Abstract supertype of various basis sets that vector spherical harmonics may be decomposed in.
"""
abstract type Basis end
@doc raw"""
    SphericalCovariant <: Basis

The spherical covariant basis ``\chi_\mu`` for ``\mu\in\{-1,0,1\}``
"""
struct SphericalCovariant <: Basis end
@doc raw"""
    Polar <: Basis

The spherical polar basis ``\hat{r}``, ``\hat{\theta}``, ``\hat{\phi}``.
"""
struct Polar <: Basis end
@doc raw"""
    HelicityCovariant  <: Basis

The helicity basis ``\mathbf{e}_\mu`` for ``\mu\in\{-1,0,1\}``
"""
struct HelicityCovariant <: Basis end
@doc raw"""
    Cartesian  <: Basis

The Cartesian basis ``\hat{x}``, ``\hat{y}``, ``\hat{z}``
"""
struct Cartesian <: Basis end

Base.broadcastable(B::Basis) = Ref(B)

_basisinds(::Union{HelicityCovariant, SphericalCovariant}) = -1:1
_basisinds(::Union{Cartesian, Polar}) = 1:3

#= Define matrices that convert between components C1 and C2 in bases B1 and B2 as
    M * C1 = C2

The corresponding functions are named B1_B2_conversion.
These matrices are conjugates of those that convert between the bases, so conj.(M) * B1 = B2
=#
const HelicitySphericalPolarConversionMatrix = SMatrix{3,3}([
                                        0       1   0
                                        1/√2    0   -1/√2
                                        -im/√2   0    -im/√2
                                        ])

helicity_polar_conversion(θ, ϕ) = HelicitySphericalPolarConversionMatrix
polar_helicity_conversion(θ, ϕ) = helicity_polar_conversion(θ, ϕ)'

function helicity_spherical_conversion(θ, ϕ)
    cisϕ = cis(ϕ)
    sinθ, cosθ = sincos(θ)
    invsqrt2 = 1/√2
    sin²θby2 = (1 - cosθ)/2
    cos²θby2 = (1 + cosθ)/2
    SMatrix{3,3}((
        cos²θby2 * cisϕ,
        -sinθ * invsqrt2,
        sin²θby2 * conj(cisϕ),

        sinθ * invsqrt2 * cisϕ,
        cosθ,
        -sinθ * invsqrt2 * conj(cisϕ),

        sin²θby2 * cisϕ,
        sinθ * invsqrt2,
        cos²θby2 * conj(cisϕ),
        ))
end
spherical_helicity_conversion(θ, ϕ) = helicity_spherical_conversion(θ, ϕ)'

function helicity_cartesian_conversion(θ, ϕ)
    sinθ, cosθ = sincos(θ)
    sinϕ, cosϕ = sincos(ϕ)
    invsqrt2 = 1/√2

    SMatrix{3,3}((
        (cosθ * cosϕ + im * sinϕ) * invsqrt2,
        (cosθ * sinϕ - im * cosϕ) * invsqrt2,
        -sinθ * invsqrt2,

        sinθ * cosϕ,
        sinθ * sinϕ,
        cosθ,

        -(cosθ * cosϕ - im * sinϕ) * invsqrt2,
        -(cosθ * sinϕ + im * cosϕ) * invsqrt2,
        sinθ * invsqrt2,
        ))
end
cartesian_helicity_conversion(θ, ϕ) = helicity_cartesian_conversion(θ, ϕ)'

function spherical_polar_conversion(θ, ϕ)
    invsqrt2 = 1/√2
    invsqrt2cisϕ = invsqrt2 * cis(ϕ)
    sinθ, cosθ = sincos(θ)
    normsinθcisϕ = sinθ * invsqrt2cisϕ
    normcosθcisϕ = cosθ * invsqrt2cisϕ
    SMatrix{3,3}((
        conj(normsinθcisϕ),
        conj(normcosθcisϕ),
        -im * conj(invsqrt2cisϕ),

        cosθ,
        -sinθ,
        0,

        -normsinθcisϕ,
        -normcosθcisϕ,
        -im * invsqrt2cisϕ,
        ))
end
polar_spherical_conversion(θ, ϕ) = spherical_polar_conversion(θ, ϕ)'

const SphericalCartesianConversionMatrix = SMatrix{3,3}((
        1/√2,
        -im/√2,
        0,

        0,
        0,
        1,

        -1/√2,
        -im/√2,
        0,
        ))

spherical_cartesian_conversion(θ, ϕ) = SphericalCartesianConversionMatrix
cartesian_spherical_conversion(θ, ϕ) = spherical_cartesian_conversion(θ, ϕ)'

function cartesian_polar_conversion(θ, ϕ)
    sinθ, cosθ = sincos(θ)
    sinϕ, cosϕ = sincos(ϕ)
    SMatrix{3,3}((
        sinθ*cosϕ,
        cosθ*cosϕ,
        -sinϕ,

        sinθ*sinϕ,
        cosθ*sinϕ,
        cosϕ,

        cosθ,
        -sinθ,
        0,
        ))
end
polar_cartesian_conversion(θ, ϕ) = cartesian_polar_conversion(θ, ϕ)'

"""
    basisconversionmatrix(B1, B2, θ, ϕ)

Return the matrix that converts the components of a vector from the basis `B1` to the basis `B2`
at the point `(θ,ϕ)`.

!!! note
    For the complex bases `SphericalCovariant` and `HelicityCovariant`, the matrix transforms the
    contravariant components between bases.

# Examples
```jldoctest
julia> v = [1,0,0] # x
3-element Vector{Int64}:
 1
 0
 0

julia> θ, ϕ = pi/2, 0
(1.5707963267948966, 0)

julia> M = VectorSphericalHarmonics.basisconversionmatrix(Cartesian(), Polar(), θ, ϕ)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  1.0          0.0   6.12323e-17
  6.12323e-17  0.0  -1.0
 -0.0          1.0   0.0

julia> M * v # polar coordinates along r, θ, ϕ
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 1.0
 6.123233995736766e-17
 0.0
```
"""
basisconversionmatrix(::T, ::T, θ, ϕ) where {T<:Basis} = I
basisconversionmatrix(::HelicityCovariant, ::SphericalCovariant, θ, ϕ) = helicity_spherical_conversion(θ, ϕ)
basisconversionmatrix(::SphericalCovariant, ::HelicityCovariant, θ, ϕ) = spherical_helicity_conversion(θ, ϕ)
basisconversionmatrix(::HelicityCovariant, ::Polar, θ, ϕ) = helicity_polar_conversion(θ, ϕ)
basisconversionmatrix(::Polar, ::HelicityCovariant, θ, ϕ) = polar_helicity_conversion(θ, ϕ)
basisconversionmatrix(::HelicityCovariant, ::Cartesian, θ, ϕ) = helicity_cartesian_conversion(θ, ϕ)
basisconversionmatrix(::Cartesian, ::HelicityCovariant, θ, ϕ) = cartesian_helicity_conversion(θ, ϕ)

basisconversionmatrix(::SphericalCovariant, ::Polar, θ, ϕ) = spherical_polar_conversion(θ, ϕ)
basisconversionmatrix(::Polar, ::SphericalCovariant, θ, ϕ) = polar_spherical_conversion(θ, ϕ)
basisconversionmatrix(::SphericalCovariant, ::Cartesian, θ, ϕ) = spherical_cartesian_conversion(θ, ϕ)
basisconversionmatrix(::Cartesian, ::SphericalCovariant, θ, ϕ) = cartesian_spherical_conversion(θ, ϕ)

basisconversionmatrix(::Cartesian, ::Polar, θ, ϕ) = cartesian_polar_conversion(θ, ϕ)
basisconversionmatrix(::Polar, ::Cartesian, θ, ϕ) = polar_cartesian_conversion(θ, ϕ)

ext: jl


text:
module Decapodes

using Requires

include("Diagrams.jl")
include("OpenDiagrams.jl")
include("Schedules.jl")
include("Simulations.jl")
include("Examples.jl")

function __init__()
  @require CairoMakie="13f3f980-e62b-5c42-98c6-ff1f3baf88f0" include("Debug.jl")
  @require AlgebraicPetri="4f99eebe-17bf-4e98-b6a1-2c4f205a959b" include("PetriNets.jl")
end
end


text:
# BitIndex
# --------
#
# Utils for indexing bits in a vector of 64-bit integers (internal use only).
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

struct BitIndex
    val::Int64
end

BitIndex(index, nbits) = BitIndex((index - 1) << trailing_zeros(nbits))

#         index(i)-1        index(i)        index(i)+1
# ....|................|..X.............|................|....
#                          |<-offset(i)-|
#                      |<--- 64 bits -->|

index(i::BitIndex) = (i.val >> 6) + 1
offset(i::BitIndex) = i.val & 0b111111

Base.:+(i::BitIndex, n::Int) = BitIndex(i.val + n)
Base.:-(i::BitIndex, n::Int) = BitIndex(i.val - n)
Base.:-(i1::BitIndex, i2::BitIndex) = i1.val - i2.val
Base.:(==)(i1::BitIndex, i2::BitIndex) = i1.val == i2.val
Base.isless(i1::BitIndex, i2::BitIndex) = isless(i1.val, i2.val)
Base.cmp(i1::BitIndex, i2::BitIndex) = cmp(i1.val, i2.val)

function Base.iterate(i::BitIndex, s=1)
    if s == 1
        return (index(i), 2)
    elseif s == 2
        return (offset(i), 3)
    else
        return nothing
    end
end

Base.show(io::IO, i::BitIndex) = print(io, '(', index(i), ", ", offset(i), ')')

# Create a bit mask that fills least significant `n` bits (`n` must be a
# non-negative integer).
bitmask(n::Integer) = bitmask(UInt64, n)
bitmask(::Type{T}, n::Integer) where {T} = (one(T) << n) - one(T)


text:
#".\\src\\testing\\withWords.TextGrid"
#input = raw".\\src\\testing\\withWords.TextGrid"
input = raw"C:\Users\hemad\Desktop\Master\GoogleSpeech\Example1\file2\A032_Sarah_edited.TextGrid"
f = open(input, "r")
#CHECK FILE TYPE
occursin("ooTextFile", readline(f)) ? print("Successuflly opened ooTextFile\n") : print("File type NOT: \"ooTextFile\n")#line#1
#CHECK OBJECT CLASS
if !(readline(f) == "Object class = \"TextGrid\"")#line#2
    throw(UndefVarError(:NotTextGridFile))
end

readline(f)#SKIP BLANK LINE (line#3)
#READ START TIME
xmin = parse(Float64,readline(f)[8:end])#line#4
#READ END TIME
xmax = parse(Float64,readline(f)[8:end])#line#5

#READ TIERS INFO
tiers = occursin("tiers? <exists>", readline(f))#line#6
tierSize = parse(Int64,readline(f)[8:end])#line#7

#CREATE VECTOR VARIABLES
class = Array{String}(undef, 1, tierSize)
name = Array{String}(undef, 1, tierSize)
xmin = Array{Float64}(undef, 1, tierSize)
xmax = Array{Float64}(undef, 1, tierSize)
intervalSize = Array{Int64}(undef, 1, tierSize)
interval = Vector{Vector}(undef, tierSize)

readline(f)#SKIP "item[]" (line#8)

for i in 1:tierSize
    readline(f)#SKIP "item[i]" (line#9....)
    println("tier#",i)
    # class[i] = readline(f)[18:end-2]#println("class: " * readline(f)[18:end-2])
    # name[i] = readline(f)[17:end-2]#println("name: " * readline(f)[17:end-2])
    # xmin[i] = parse(Float64,readline(f)[15:end])#println("xmin: " * readline(f)[15:end])
    # xmax[i] = parse(Float64,readline(f)[15:end])#println("xmax: " * readline(f)[15:end])
    # intervalSize[i] = parse(Int64,readline(f)[26:end])#println("intervalSize: " * readline(f)[26:end])

    println("class: " * readline(f)[18:end-2])
    println("name: " * readline(f)[17:end-1])
    println("xmin: " * readline(f)[15:end])
    println("xmax: " * readline(f)[15:end])
    println("intervalSize: " * readline(f)[26:end])

    #intervalLines = Vector{Any}(undef,intervalSize[i])

    # for j in 1:intervalSize[i]
    #     #println(readline(f),readline(f),readline(f),readline(f),"\n")
    #     readline(f)#SKIP "intervals [j]:" (line#15)
    #     #println([readline(f)[19:end],readline(f)[19:end],readline(f)[21:end-2]])
    #     intervalLines[j] = [parse(Float64,readline(f)[19:end]),
    #                         parse(Float64,readline(f)[19:end]),
    #                         readline(f)[21:end-1]]
    #     #append!(interval,[readline(f)[19:end],readline(f)[19:end],readline(f)[21:end-2]])
    #     if(j == intervalSize[i])
    #         interval[i] = copy(intervalLines)
    #     end
    # end
end
close(f)

# for i in 1:length(interval[1])
#     test = 0
#     for n in 1:length(interval[7])
#         if ((interval[1][i][1]<interval[7][n][1]) & (interval[7][n][1]<interval[1][i][2]))
#             test += 1
#         end
#     end
#     val = test/(interval[1][i][2]-interval[1][i][1])*60
#     println("interval ",i," \"",interval[1][i][3],'"'," words: ",test," WPM: ",val)
# end

#
#
# function WPM(input)
#     f = open(input, "r")
#     #CHECK FILE TYPE
#     occursin("ooTextFile", readline(f)) ? print("Successuflly opened ooTextFile") : print("File type NOT: \"ooTextFile")#line#1
#     #CHECK OBJECT CLASS
#     if !(readline(f) == "Object class = \"TextGrid\"")#line#2
#         throw(UndefVarError(:NotTextGridFile))
#     end
#
#     readline(f)#SKIP BLANK LINE (line#3)
#     #READ START TIME
#     xmin = parse(Float64,readline(f)[8:end])#line#4
#     #READ END TIME
#     xmax = parse(Float64,readline(f)[8:end])#line#5
#
#     #READ TIERS INFO
#     tiers = occursin("tiers? <exists>", readline(f))#line#6
#     tierSize = parse(Int64,readline(f)[8:end])#line#7
#
#     #CREATE VECTOR VARIABLES
#     class = Array{String}(undef, 1, tierSize)
#     name = Array{String}(undef, 1, tierSize)
#     xmin = Array{Float64}(undef, 1, tierSize)
#     xmax = Array{Float64}(undef, 1, tierSize)
#     intervalSize = Array{Int64}(undef, 1, tierSize)
#     interval = Vector{Vector}(undef, tierSize)
#
#     readline(f)#SKIP "item[]" (line#8)
#
#     for i in 1:tierSize
#         readline(f)#SKIP "item[i]" (line#9....)
#         class[i] = readline(f)[18:end-2]#println("class: " * readline(f)[18:end-2])
#         name[i] = readline(f)[17:end-2]#println("name: " * readline(f)[17:end-2])
#         xmin[i] = parse(Float64,readline(f)[15:end])#println("xmin: " * readline(f)[15:end])
#         xmax[i] = parse(Float64,readline(f)[15:end])#println("xmax: " * readline(f)[15:end])
#         intervalSize[i] = parse(Int64,readline(f)[26:end])#println("intervalSize: " * readline(f)[26:end])
#
#         intervalLines = Vector{Any}(undef,intervalSize[i])
#         println("tier#",i)
#         for j in 1:intervalSize[i]
#             #println(readline(f),readline(f),readline(f),readline(f),"\n")
#             readline(f)#SKIP "intervals [j]:" (line#15)
#             #println([readline(f)[19:end],readline(f)[19:end],readline(f)[21:end-2]])
#             intervalLines[j] = [parse(Float64,readline(f)[19:end]),
#                                 parse(Float64,readline(f)[19:end]),
#                                 readline(f)[21:end-1]]
#             #append!(interval,[readline(f)[19:end],readline(f)[19:end],readline(f)[21:end-2]])
#             if(j == intervalSize[i])
#                 interval[i] = copy(intervalLines)
#             end
#         end
#     end
#     close(f)
# end
# input = raw"C:\Users\hemad\Desktop\Master\GoogleSpeech\Example1\file2\A032_Sarah_edited.TextGrid"
# WPM(input)
#
#
#
#
# #=
# for each data in item 1 (interval[1])
# find number of words in item 7
# =#
#
# for i in 1:length(interval[1])
#     test = 0
#     for n in 1:length(interval[7])
#         if ((interval[1][i][1]<interval[7][n][1]) & (interval[7][n][1]<interval[1][i][2]))
#             test += 1
#         end
#     end
#     val = test/(interval[1][i][2]-interval[1][i][1])*60
#     println("interval ",i," \"",interval[1][i][3],'"'," words: ",test," WPM: ",val)
# end


text:
# Main GaussianProcess type

mutable struct GPE{X<:MatF64,Y<:VecF64,M<:Mean,K<:Kernel,P<:AbstractPDMat,D<:KernelData} <: GPBase
    # Observation data
    "Input observations"
    x::X
    "Output observations"
    y::Y

    # Model
    "Mean object"
    mean::M
    "Kernel object"
    kernel::K
    "Log standard deviation of observation noise"
    logNoise::Float64

    # Auxiliary data
    "Dimension of inputs"
    dim::Int
    "Number of observations"
    nobs::Int
    "Auxiliary observation data (to speed up calculations)"
    data::D
    "`(k + exp(2*obsNoise))`"
    cK::P
    "`(k + exp(2*obsNoise))⁻¹y`"
    alpha::Vector{Float64}
    "Marginal log-likelihood"
    mll::Float64
    "Log target (marginal log-likelihood + log priors)"
    target::Float64
    "Gradient of marginal log-likelihood"
    dmll::Vector{Float64}
    "Gradient of log-target (gradient of marginal log-likelihood + gradient of log priors)"
    dtarget::Vector{Float64}

    function GPE{X,Y,M,K,P,D}(x::X, y::Y, mean::M, kernel::K, data::D, cK::P, logNoise::Float64) where {X,Y,M,K,P,D}
        dim, nobs = size(x)
        length(y) == nobs || throw(ArgumentError("Input and output observations must have consistent dimensions."))
        gp = new{X,Y,M,K,P,D}(x, y, mean, kernel, logNoise, dim, nobs, data, cK)
        initialise_target!(gp)
    end
end

"""
    GPE(x, y, mean, kernel[, logNoise])

Fit a Gaussian process to a set of training points. The Gaussian process is defined in
terms of its user-defined mean and covariance (kernel) functions. As a default it is
assumed that the observations are noise free.

# Arguments:
- `x::AbstractVecOrMat{Float64}`: Input observations
- `y::AbstractVector{Float64}`: Output observations
- `mean::Mean`: Mean function
- `kernel::Kernel`: Covariance function
- `logNoise::Float64`: Natural logarithm of the standard deviation for the observation
  noise. The default is -2.0, which is equivalent to assuming no observation noise.
"""
function GPE(x::MatF64, y::VecF64, mean::Mean, kernel::Kernel, data::KernelData, cK::AbstractPDMat, logNoise::Float64)
    GPE{typeof(x),typeof(y),typeof(mean),typeof(kernel),typeof(cK),typeof(data)}(x, y, mean, kernel, data, cK, logNoise)
end
function GPE(x::MatF64, y::VecF64, mean::Mean, kernel::Kernel, logNoise::Float64 = -2.0)
    nobs = length(y)
    kerneldata = KernelData(kernel, x)
    GPE(x, y, mean, kernel, kerneldata,
        PDMat(Σ_default(x, kernel, kerneldata, logNoise)),
        logNoise)
end
GPE(x::VecF64, y::VecF64, mean::Mean, kernel::Kernel, logNoise::Float64 = -2.0) =
    GPE(x', y, mean, kernel, logNoise)

"""
    GPE(; mean::Mean = MeanZero(), kernel::Kernel = SE(0.0, 0.0), logNoise::Float64 = -2.0)

Construct a [GPE](@ref) object without observations.
"""
function GPE(; mean::Mean = MeanZero(), kernel::Kernel = SE(0.0, 0.0), logNoise::Float64 = -2.0)
    x = Array{Float64}(undef, 1, 0) # ElasticArrays don't like length(x) = 0.
    y = Array{Float64}(undef, 0)
    GPE(x, y, mean, kernel, logNoise)
end

"""
    GP(x, y, mean::Mean, kernel::Kernel[, logNoise::Float64=-2.0])

Fit a Gaussian process that is defined by its `mean`, its `kernel`, and the logarithm
`logNoise` of the standard deviation of its observation noise to a set of training points
`x` and `y`.

See also: [`GPE`](@ref)
"""
GP(x::AbstractVecOrMat{Float64}, y::VecF64, mean::Mean, kernel::Kernel,
   logNoise::Float64 = -2.0) = GPE(x, y, mean, kernel, logNoise)

"""
    fit!(gp::GPE{X,Y}, x::X, y::Y)

Fit Gaussian process `GPE` to a training data set consisting of input observations `x` and
output observations `y`.
"""
function fit!(gp::GPE{X,Y}, x::X, y::Y) where {X,Y}
    length(y) == size(x,2) || throw(ArgumentError("Input and output observations must have consistent dimensions."))
    gp.x = x
    gp.y = y
    gp.data = KernelData(gp.kernel, x)
    gp.cK = PDMat(Σ_default(gp))
    gp.dim, gp.nobs = size(x)
    initialise_target!(gp)
end

fit!(gp::GPE, x::VecF64, y::VecF64) = fit!(gp, x', y)

#———————————————————————————————————————————————————————————
#Fast memory allocation function

"""
    get_ααinvcKI!(ααinvcKI::Matrix{Float64}, cK::AbstractPDMat, α::Vector)

Write `ααᵀ - cK * eye(nobs)` to `ααinvcKI` avoiding any memory allocation, where `cK` and
`α` are the covariance matrix and the alpha vector of a Gaussian process, respectively.
Hereby `α` is defined as `cK \\ (Y - μ)`.
"""
function get_ααinvcKI!(ααinvcKI::MatF64, cK::AbstractPDMat, α::Vector)
    nobs = length(α)
    size(ααinvcKI) == (nobs, nobs) || throw(ArgumentError(
                @sprintf("Buffer for ααinvcKI should be a %dx%d matrix, not %dx%d",
                         nobs, nobs,
                         size(ααinvcKI,1), size(ααinvcKI,2))))
    fill!(ααinvcKI, 0)
    @inbounds for i in 1:nobs
        ααinvcKI[i,i] = -1.0
    end
    ldiv!(cK.chol, ααinvcKI)
    BLAS.ger!(1.0, α, α, ααinvcKI)
end

#———————————————————————————————————————————————————————————————-
#Functions for calculating the log-target

Σ_default(gp) = Σ_default(gp.x, gp.kernel, gp.data, gp.logNoise)

Σ_default(x, kernel, data, logNoise) = make_posdef!(cov(kernel, x, data) + exp(2*logNoise)*I)[1]

"""
    update_cK!(gp::GPE)

Update the covariance matrix and its Cholesky decomposition of Gaussian process `gp`.
"""
function update_cK!(gp::GPE)
    old_cK = gp.cK
    Σbuffer = mat(old_cK)
    cov!(Σbuffer, gp.kernel, gp.x, gp.data)
    noise = (exp(2*gp.logNoise) + 1e-5)
    for i in 1:gp.nobs
        Σbuffer[i,i] += noise
    end
    chol_buffer = cholfactors(old_cK)
    copyto!(chol_buffer, Σbuffer)
    chol = cholesky!(Symmetric(chol_buffer))
    gp.cK = wrap_cK(gp.cK, Σbuffer, chol)
    gp.cK
end

wrap_cK(cK::PDMat, Σbuffer, chol) = PDMat(Σbuffer, chol)
mat(cK::PDMat) = cK.mat
cholfactors(cK::PDMat) = cK.chol.factors

# modification of initialise_target! that reuses existing matrices to avoid
# unnecessary memory allocations, which speeds things up significantly
function update_mll!(gp::GPE; noise::Bool=true, domean::Bool=true, kern::Bool=true)
    if kern | noise
        update_cK!(gp)
    end
    μ = mean(gp.mean, gp.x)
    y = gp.y - μ
    gp.alpha = gp.cK \ y
    # Marginal log-likelihood
    gp.mll = - (dot(y, gp.alpha) + logdet(gp.cK) + log2π * gp.nobs) / 2
    gp
end

function dmll_kern!(dmll::VecF64, k::Kernel, X::MatF64, data::KernelData, ααinvcKI::MatF64)
    dim, nobs = size(X)
    nparams = num_params(k)
    @assert nparams == length(dmll)
    dK_buffer = Vector{Float64}(undef, nparams)
    dmll[:] .= 0.0
    @inbounds for j in 1:nobs
        # off-diagonal
        for i in 1:j-1
            dKij_dθ!(dK_buffer, k, X, data, i, j, dim, nparams)
            @simd for iparam in 1:nparams
                dmll[iparam] += dK_buffer[iparam] * ααinvcKI[i, j]
            end
        end
        # diagonal
        dKij_dθ!(dK_buffer, k, X, data, j, j, dim, nparams)
        for iparam in 1:nparams
            dmll[iparam] += dK_buffer[iparam] * ααinvcKI[j, j] / 2.0
        end
    end
    return dmll
end
"""
     update_dmll!(gp::GPE, ...)

Update the gradient of the marginal log-likelihood of Gaussian process `gp`.
"""
function update_dmll!(gp::GPE, ααinvcKI::MatF64;
    noise::Bool=true, # include gradient component for the logNoise term
    domean::Bool=true, # include gradient components for the mean parameters
    kern::Bool=true, # include gradient components for the spatial kernel parameters
    )
    size(ααinvcKI) == (gp.nobs, gp.nobs) || throw(ArgumentError(
                @sprintf("Buffer for ααinvcKI should be a %dx%d matrix, not %dx%d",
                         gp.nobs, gp.nobs,
                         size(ααinvcKI,1), size(ααinvcKI,2))))
    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)
    gp.dmll = Array{Float64}(undef, noise + domean * n_mean_params + kern * n_kern_params)

    get_ααinvcKI!(ααinvcKI, gp.cK, gp.alpha)

    i=1
    if noise
        gp.dmll[i] = exp(2 * gp.logNoise) * tr(ααinvcKI)
        i += 1
    end

    if domean && n_mean_params>0
        Mgrads = grad_stack(gp.mean, gp.x)
        for j in 1:n_mean_params
            gp.dmll[i] = dot(Mgrads[:,j], gp.alpha)
            i += 1
        end
    end
    if kern
        dmll_k = @view(gp.dmll[i:end])
        dmll_kern!(dmll_k, gp.kernel, gp.x, gp.data, ααinvcKI)
    end
end

"""
    update_mll_and_dmll!(gp::GPE, ...)

Update the gradient of the marginal log-likelihood of a Gaussian
process `gp`.
"""
function update_mll_and_dmll!(gp::GPE,
        ααinvcKI::MatF64;
        kwargs...)
    update_mll!(gp; kwargs...)
    update_dmll!(gp, ααinvcKI; kwargs...)
end

"""
    initialise_target!(gp::GPE)

Initialise the log-posterior
```math
\\log p(θ | y) ∝ \\log p(y | θ) +  \\log p(θ)
```
of a Gaussian process `gp`.
"""
function initialise_target!(gp::GPE)
    update_mll!(gp)
        #HOW TO SET-UP A PRIOR FOR THE LOGNOISE?
    gp.target = gp.mll   + prior_logpdf(gp.mean) + prior_logpdf(gp.kernel)
    gp
end

"""
    update_target!(gp::GPE, ...)

Update the log-posterior
```math
\\log p(θ | y) ∝ \\log p(y | θ) +  \\log p(θ)
```
of a Gaussian process `gp`.
"""
function update_target!(gp::GPE; noise::Bool=true, domean::Bool=true, kern::Bool=true)
    update_mll!(gp; noise=noise, domean=domean, kern=kern)
    #HOW TO SET-UP A PRIOR FOR THE LOGNOISE?
    gp.target = gp.mll  + prior_logpdf(gp.mean) + prior_logpdf(gp.kernel)
    gp
end

function update_dtarget!(gp::GPE, Kgrad::MatF64, L_bar::MatF64; kwargs...)
    update_dmll!(gp, L_bar; kwargs...)
    gp.dtarget = gp.dmll + prior_gradlogpdf(gp; kwargs...)
    gp
end

"""
    update_target_and_dtarget!(gp::GPE, ...)

Update the log-posterior
```math
\\log p(θ | y) ∝ \\log p(y | θ) +  \\log p(θ)
```
of a Gaussian process `gp` and its derivative.
"""
function update_target_and_dtarget!(gp::GPE, Kgrad::MatF64, L_bar::MatF64; kwargs...)
    update_target!(gp; kwargs...)
    update_dtarget!(gp, Kgrad, L_bar; kwargs...)
end

function update_target_and_dtarget!(gp::GPE; kwargs...)
    ααinvcKI = Array{Float64}(undef, gp.nobs, gp.nobs)
    update_target!(gp; kwargs...)
    update_dmll!(gp, ααinvcKI; kwargs...)
    gp.dtarget = gp.dmll + prior_gradlogpdf(gp; kwargs...)
end

#——————————————————————#
# Predict observations #
#——————————————————————#

"""
    predict_y(gp::GPE, x::Union{Vector{Float64},Matrix{Float64}}[; full_cov::Bool=false])

Return the predictive mean and variance of Gaussian Process `gp` at specfic points which
are given as columns of matrix `x`. If `full_cov` is `true`, the full covariance matrix is
returned instead of only variances.
"""
function predict_y(gp::GPE, x::MatF64; full_cov::Bool=false)
    μ, σ2 = predict_f(gp, x; full_cov=full_cov)
    if full_cov
        return μ, σ2 + ScalMat(gp.nobs, exp(2*gp.logNoise))
    else
        return μ, σ2 .+ exp(2*gp.logNoise)
    end
end

# 1D Case for predictions
predict_y(gp::GPE, x::VecF64; full_cov::Bool=false) = predict_y(gp, x'; full_cov=full_cov)

@deprecate predict predict_y

## compute predictions
function _predict(gp::GPE, x::MatF64)
    cK = cov(gp.kernel, gp.x, x)
    Lck = whiten(gp.cK, cK)
    mu = mean(gp.mean, x) + cK'*gp.alpha        # Predictive mean
    Sigma_raw = cov(gp.kernel, x) - Lck'Lck # Predictive covariance
    # Add jitter to get stable covariance
    Sigma = tolerant_PDMat(Sigma_raw)
    return mu, Sigma
end

#———————————————————————————————————————————————————————————
# Sample from the GPE
function Random.rand!(gp::GPE, x::MatF64, A::DenseMatrix)
    nobs = size(x,2)
    n_sample = size(A,2)

    if gp.nobs == 0
        # Prior mean and covariance
        μ = mean(gp.mean, x);
        Σraw = cov(gp.kernel, x);
        Σ = tolerant_PDMat(Σraw)
    else
        # Posterior mean and covariance
        μ, Σ = predict_f(gp, x; full_cov=true)
    end
    return broadcast!(+, A, μ, unwhiten!(Σ,randn(nobs, n_sample)))
end

Random.rand(gp::GPE, x::MatF64, n::Int) = rand!(gp, x, Array{Float64}(undef, size(x, 2), n))

# Sample from 1D GPE
Random.rand(gp::GPE, x::VecF64, n::Int) = rand(gp, x', n)

# Generate only one sample from the GPE and returns a vector
Random.rand(gp::GPE, x::MatF64) = vec(rand(gp,x,1))
Random.rand(gp::GPE, x::VecF64) = vec(rand(gp,x',1))

#—————————————————————————————————————————————————————–
#Functions for setting and calling the parameters of the GP object

function get_params(gp::GPE; noise::Bool=true, domean::Bool=true, kern::Bool=true)
    params = Float64[]
    if noise; push!(params, gp.logNoise); end
    if domean
        append!(params, get_params(gp.mean))
    end
    if kern
        append!(params, get_params(gp.kernel))
    end
    return params
end

function num_params(gp::GPE; noise::Bool=true, domean::Bool=true, kern::Bool=true)
    n = 0
    noise && (n += 1)
    domean && (n += num_params(gp.mean))
    kern && (n += num_params(gp.kernelern))
    n
end

function set_params!(gp::GPE, hyp::VecF64;
                     noise::Bool=true, domean::Bool=true, kern::Bool=true)
    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)

    i = 1
    if noise
        gp.logNoise = hyp[1];
        i+=1
    end

    if domean && n_mean_params>0
        set_params!(gp.mean, hyp[i:i+n_mean_params-1])
        i += n_mean_params
    end
    if kern
        set_params!(gp.kernel, hyp[i:i+n_kern_params-1])
        i += n_kern_params
    end
end

function prior_gradlogpdf(gp::GPE; noise::Bool=true, domean::Bool=true, kern::Bool=true)
    grad = Float64[]
    if noise; push!(grad, 0.0); end # Noise does not have any priors
    if domean
        append!(grad, prior_gradlogpdf(gp.mean))
    end
    if kern
        append!(grad, prior_gradlogpdf(gp.kernel))
    end
    return grad
end

#———————————————————————————————————————————————————————————-
# Push function
function Base.push!(gp::GPE, x::MatF64, y::VecF64)
    @warn "push! method is currently inefficient as it refits all observations"
    if gp.nobs == 0
        GaussianProcesses.fit!(gp, x, y)
    elseif size(x,1) != size(gp.x,1)
        error("New input observations must have dimensions consistent with existing observations")
    else
        GaussianProcesses.fit!(gp, cat(2, gp.x, x), cat(1, gp.y, y))
    end
end

Base.push!(gp::GPE, x::VecF64, y::VecF64) = push!(gp, x', y)
Base.push!(gp::GPE, x::Float64, y::Float64) = push!(gp, [x], [y])
Base.push!(gp::GPE, x::VecF64, y::Float64) = push!(gp, reshape(x, length(x), 1), [y])

# —————————————————————————————————————————————————————————————
# Show function
function Base.show(io::IO, gp::GPE)
    println(io, "GP Exact object:")
    println(io, "  Dim = ", gp.dim)
    println(io, "  Number of observations = ", gp.nobs)
    println(io, "  Mean function:")
    show(io, gp.mean, 2)
    println(io, "\n  Kernel:")
    show(io, gp.kernel, 2)
    if gp.nobs == 0
        println("\n  No observation data")
    else
        println(io, "\n  Input observations = ")
        show(io, gp.x)
        print(io, "\n  Output observations = ")
        show(io, gp.y)
        print(io, "\n  Variance of observation noise = ", exp(2 * gp.logNoise))
        print(io, "\n  Marginal Log-Likelihood = ")
        show(io, round(gp.target; digits = 3))
    end
end


text:
module PointProcessGeneration

using Distributions
using Distances

# Write your package code here.
include("simple.jl")
include("metric.jl")
include("support.jl")
include("colored.jl")

# simple
export poisson, normal, circle_poisson, circle_uniform
# with support
export clusters, voronoi, circles
# colored
export normal_clusters, circle_clusters, anarchy, colored_clusters, target_circles

end


text:
#!/usr/bin/env julia

#struct session
#    ifvl::Bool
#    ifad::Bool
#    ifpr::Bool
#    ifps::Bool
#
#    nx1::Int
#    ny1::Int
#end

using SpectralElements

using FastGaussQuadrature
using Plots,LinearAlgebra,SparseArrays

using Zygote

#using NLopt,Optim
using DiffEqFlux
using Flux
using Statistics

#----------------------------------------------------------------------#
# size
#----------------------------------------------------------------------#
nx1 = 8; Ex = 5;
ny1 = 8; Ey = 5;

nx2 = nx1-2; nxd = Int(ceil(1.5*nx1));
ny2 = nx1-2; nyd = Int(ceil(1.5*ny1));

nxp = 5*nx1;
nyp = 5*ny1;

#----------------------------------------------------------------------#
# nodal operators
#----------------------------------------------------------------------#
zr1,wr1 = gausslobatto(nx1);  zs1,ws1 = gausslobatto(ny1);
zr2,wr2 = gausslobatto(nx2);  zs2,ws2 = gausslobatto(ny2);
zrd,wrd = gausslobatto(nxd);  zsd,wsd = gausslobatto(nyd);
zrp     = linspace(-1,1,nxp); zsp     = linspace(-1,1,nyp);

Jr1d = interpMat(zrd,zr1); Js1d = interpMat(zsd,zs1);
Jr2d = interpMat(zrd,zr2); Js2d = interpMat(zsd,zs2);
Jr21 = interpMat(zr1,zr2); Js21 = interpMat(zs1,zs2);
Jr1p = interpMat(zrp,zr1); Js1p = interpMat(zsp,zs1);

Dr1 = derivMat(zr1); Ds1 = derivMat(zs1);
Dr2 = derivMat(zr2); Ds2 = derivMat(zs2);
Drd = derivMat(zrd); Dsd = derivMat(zsd);

#----------------------------------------------------------------------#
# boundary conditions
#----------------------------------------------------------------------#
ifperiodicX = false
ifperiodicY = false

Ix1 = Matrix(I,Ex*nx1,Ex*nx1);
Iy1 = Matrix(I,Ey*ny1,Ey*ny1);

Rx1 = Ix1[2:end-1,:];
Ry1 = Iy1[2:end-1,:];

if(ifperiodicX) Rx1 = Ix1; end
if(ifperiodicY) Ry1 = Iy1; end

M1 = diag(Rx1'*Rx1) * diag(Ry1'*Ry1)'; # hom. dir BC

Ix1 = nothing
Iy1 = nothing
Rx1 = nothing
Ry1 = nothing
#----------------------------------------------------------------------#
# mapping
#----------------------------------------------------------------------#

# Q: global -> local op, Q': local -> global
Qx1 = semq(Ex,nx1,ifperiodicX); Qx2 = semq(Ex,nx2,ifperiodicX);
Qy1 = semq(Ey,ny1,ifperiodicY); Qy2 = semq(Ey,ny2,ifperiodicY);

# gather scatter op
QQtx1 = Qx1*Qx1'; QQtx2 = Qx2*Qx2';
QQty1 = Qy1*Qy1'; QQty2 = Qy2*Qy2';

# weight for inner products
mult1 = ones(nx1*Ex,ny1*Ey);
mult1 = gatherScatter(mult1,QQtx1,QQty1);
mult1 = @. 1 / mult1;

mult2 = ones(nx2*Ex,ny2*Ey);
mult2 = gatherScatter(mult2,QQtx2,QQty2);
mult2 = @. 1 / mult2;
#----------------------------------------------------------------------#
# geometry
#----------------------------------------------------------------------#
x1e,_ = semmesh(Ex,nx1); y1e,_ = semmesh(Ey,ny1);
x2e,_ = semmesh(Ex,nx2); y2e,_ = semmesh(Ey,ny2);
xde,_ = semmesh(Ex,nxd); yde,_ = semmesh(Ey,nyd);

x1,y1 = ndgrid(x1e,y1e);
x2,y2 = ndgrid(x2e,y2e);
xd,yd = ndgrid(xde,yde);

# deform grid with gordonhall
function usrgeom(r,s)
    x = @. 0.5*(r+1)
    y = @. 0.5*(s+1)
return x,y
end

x1,y1 = usrgeom(x1,y1)
x2,y2 = usrgeom(x2,y2)
xd,yd = usrgeom(xd,yd)

Jac1,Jaci1,rx1,ry1,sx1,sy1 = jac(x1,y1,Dr1,Ds1);
Jac2,Jaci2,rx2,ry2,sx2,sy2 = jac(x2,y2,Dr2,Ds2);
Jacd,Jacid,rxd,ryd,sxd,syd = jac(xd,yd,Drd,Dsd);

wx1 = kron(ones(Ex,1),wr1); wy1 = kron(ones(Ey,1),ws1);
wx2 = kron(ones(Ex,1),wr2); wy2 = kron(ones(Ey,1),ws2);
wxd = kron(ones(Ex,1),wrd); wyd = kron(ones(Ey,1),wsd);

# diagonal mass matrix
B1  = Jac1 .* (wx1*wy1');
B2  = Jac2 .* (wx2*wy2');
Bd  = Jacd .* (wxd*wyd');
Bi1 = 1 ./ B1;

#----------------------------------------------------------------------#
# debugging
#----------------------------------------------------------------------#

## don't use lumped mass matrices for now
#Jr1d = []
#Js1d = []
#rxd = rx1
#ryd = ry1
#sxd = sx1
#syd = sy1
#Bd = B1

#----------------------------------------------------------------------#
# case setup
#----------------------------------------------------------------------#
V  = 0.4
pa = 5
ε  = 1e-3
α  = 1e-8
f1 = @. 1e-2 + 0*x1
#M1[:,end] .= 1;
#M1[1,:]   .= 1;
af(p)   = @. 0.5*(tanh(p)+1)
kond(a) = @. ε + (1-ε)*a^pa

function problem(a) # topology parameter, forcing

    Ba = Bd .* ABu(Js1d,Jr1d,a);

    visc  = kond(a)
    viscd = ABu(Js1d,Jr1d,visc);
    G11 = @. viscd * Ba * (rxd * rxd + ryd * ryd);
    G12 = @. viscd * Ba * (rxd * sxd + ryd * syd);
    G22 = @. viscd * Ba * (sxd * sxd + syd * syd);

    #Ba = Zygote.hook(d->hmp(d),Ba);
    lhs(v) = lapl(v,M1,Jr1d,Js1d,QQtx1,QQty1,Dr1,Ds1,G11,G12,G22,mult1);

    f   = f1;

    rhs = mass(f,M1,Ba,Jr1d,Js1d,QQtx1,QQty1,mult1);

    return lhs, rhs
end

function opM(v) # preconditioner
    return v
end

function solver(opA,rhs,adj::Bool)
    if adj
        # zero out sensitivities at bounary points
        rhs = mask(rhs,M1);
        # ensure continuity and elem boundaries
        rhs = mult1 .* gatherScatter(rhs,QQty1,QQtx1);

        return pcg(rhs,opA,opM,mult1,false);
    else
        return pcg(rhs,opA,opM,mult1,false);
    end
end

function model(a)
    u = linsolve(a,problem,solver) # adjoint support thru Zygote
end

function loss(p)
    a = af(p)
    u = model(a) # du = 1 .* B1 .* a .* f
    f = f1
    adx,ady = grad(a,Dr1,Ds1,rx1,ry1,sx1,sy1);

    ad   = ABu(Js1d,Jr1d,a  );
    ud   = ABu(Js1d,Jr1d,u  );
    fd   = ABu(Js1d,Jr1d,f  );
    adxd = ABu(Js1d,Jr1d,adx);
    adyd = ABu(Js1d,Jr1d,ady);

    ll = @. fd*ud + α*(adxd^2+adyd^2);
    l  = sum(Bd.*ad.*ll);
    return l, u
end

#----------------------------------------------------------------------#
# DiffEqFlux.sciml_train
#----------------------------------------------------------------------#
p0 = @. 0.6 + 0.0*x1

global param = []
callback = function (p, l, pred; doplot = true)
  global param = p
  return false
end

result = DiffEqFlux.sciml_train(loss,p0
                                ,ADAM(1e-1),cb=Flux.throttle(callback,1)
                                ,maxiters=10)

p=param;dp=Zygote.gradient((p)->loss(p)[1], p)[1];
u= model(af(p));

Jp  = ABu(Js1p,Jr1p,p );
Jdp = ABu(Js1p,Jr1p,dp);
Ju  = ABu(Js1p,Jr1p,u );
pl=plot(heatmap(af(Jp)),heatmap(Jdp),heatmap(Ju));
display(pl)

display(loss(p0)[1])
display(loss(p )[1])
#----------------------------------------------------------------------#
# NLopt.Optimize
#----------------------------------------------------------------------#
#function fmin(a::Vector, grad::Vector)
#    a = reshape(a,nx1*Ex,ny1*Ey)
#    if length(grad)>0
#        grad[:] = Zygote.gradient((a)->loss(a)[1], a)[1][:]
#    end
#    return loss(a)[1]
#end

#function fc(a::Vector, grad::Vector)
#    a = reshape(a,nx1*Ex,ny1*Ey)
#    ineq(a) = sum(B1.*a) .- V
#    if length(grad)>0
#        grad[:] = Zygote.gradient((a)->ineq(a), a)[1][:]
#    end
#    return ineq(a)
#end

## set up and run optimizer
#a0 = rand(nx1*Ex,ny1*Ey)
#a0 = 0.4*ones(nx1*Ex,ny1*Ey)

#fig = mesh(x1,y1,model(a0)); display(fig);

#opt = NLopt.Opt(:LD_MMA, length(a0))
#opt.min_objective = fmin
#NLopt.inequality_constraint!(opt, fc)
#opt.lower_bounds = 0.
#opt.upper_bounds = 1.
#opt.maxeval = 20
#opt.xtol_abs = 1e-8
#opt.xtol_rel = 1e-8

#(optf,optx,ret) = NLopt.optimize(opt, a0[:])

#opta = reshape(optx,nx1*Ex,ny1*Ey)
#numevals = opt.numevals
#println("started with loss ",loss(a0)[1])
#println("got $optf after $numevals iterations (returned $ret)")
#fig = mesh(x1,y1,model(opta)); display(fig);
#----------------------------------------------------------------------#
nothing


text:
using LinearAlgebra: I, diagm
using KalmanFilters.UnscentedKalmanFilter
using KalmanFilters.SigmaPoints

function fx(
    ;x::AbstractArray{Float64, 1}, δt::Float64,
)::AbstractArray{Float64, 1}
    F = Float64[[1 δt 0 0];[0 1 0 0];[0 0 1 δt];[0 0 0 1]]
    F * x
end

hx(;x::AbstractArray{Float64, 1})::AbstractArray{Float64, 1} = x[[1, 3]]

function main()
    δt = 1.0
    σ_parameters = MerweScaled(n=4, α=0.1, β=2.0, κ=-1.0)
    ukf = UKFState(dim_x=4, dim_z=2, σ_parameters=σ_parameters)
    ukf.x .= [0, 1, 0, 1]
    ukf.P .*= 0.01
    Q = nothing
    R = diagm([0.01 ^ 2, 0.01 ^ 2])

    measurements = vcat([Float64[i i] for i = 1:10]...)
    for i = 1:size(measurements, 1)
        predict!(ukf=ukf, δt=δt, Q=Q, fx=fx)
        update!(ukf=ukf, z=measurements[i, :], R=R, hx=hx)
    end
    println(ukf.x)
end

main()

max_stars_repo_name: pxl-th/KalmanFilters

max_issues_repo_name: pxl-th/KalmanFilters

max_forks_repo_name: pxl-th/KalmanFilters


text:
"""
All stochastic iterated integrals are written in Stratonovich sense as indicated
by the J.
"""

abstract type AbstractJ end

# Iterated stochastic integral (for diagonal and commutative noise processes where the Levy area approximation is not required.)
abstract type AbstractJDiagonal <: AbstractJ end
abstract type AbstractJCommute <: AbstractJ end

struct JDiagonal_oop <: AbstractJDiagonal end

mutable struct JDiagonal_iip{JType} <: AbstractJDiagonal
  J::JType
  JDiagonal_iip(ΔW) = new{typeof(ΔW)}(false .* ΔW .* ΔW)
end

struct JCommute_oop <: AbstractJCommute end

mutable struct JCommute_iip{JType} <: AbstractJCommute
  J::JType
  function JCommute_iip(ΔW)
    J = false .* vec(ΔW) .* vec(ΔW)'
    new{typeof(J)}(J)
  end
end

function get_iterated_I(dt, dW, dZ, alg::JDiagonal_oop, p=nothing, c=1, γ=1//1)
  1//2 .* dW .* dW
end

function get_iterated_I!(dt, dW, dZ, alg::JDiagonal_iip, p=nothing, c=1, γ=1//1)
  @unpack J = alg
  if typeof(dW) <: Number
    alg.J = 1//2 .* dW .^ 2
  else
    @.. J = 1//2 * dW^2
  end
  return nothing
end

function get_iterated_I(dt, dW, dZ, alg::JCommute_oop, p=nothing, c=1, γ=1//1)
  J = 1//2 .* vec(dW) .* vec(dW)'
  J
end

function get_iterated_I!(dt, dW, dZ, alg::JCommute_iip, p=nothing, c=1, γ=1//1)
  @unpack J = alg
  mul!(J, vec(dW), vec(dW)')
  @.. J *= 1//2
  return nothing
end

# algs from LevyArea.jl # LevyArea.levyarea allocates random variables and then mutates these, see e.g.
# https://github.com/stochastics-uni-luebeck/LevyArea.jl/blob/68c5cb08ab103b4dcd3178651f7a5dd9ce8c666d/src/milstein.jl#L25
function get_iterated_I(dt, dW, dZ, alg::LevyArea.AbstractIteratedIntegralAlgorithm, p=nothing, c=1, γ=1//1)
  if isnothing(p)
    ε = c * dt^(γ + 1//2)
    p = terms_needed(length(dW), dt, ε, alg, MaxL2())
  end
  I = LevyArea.levyarea(dW / √dt, p, alg)
  I .= 1//2 * dW .* dW' .+ dt .* I
end

mutable struct IteratedIntegralAlgorithm_iip{JType,levyalgType} <: LevyArea.AbstractIteratedIntegralAlgorithm
  J::JType
  levyalg::levyalgType
  function IteratedIntegralAlgorithm_iip(ΔW, levyalg)
    J = false .* vec(ΔW) .* vec(ΔW)'
    new{typeof(J),typeof(levyalg)}(J, levyalg)
  end
end

function get_iterated_I!(dt, dW, dZ, alg::IteratedIntegralAlgorithm_iip, p=nothing, c=1, γ=1//1)
  @unpack J, levyalg = alg
  if isnothing(p)
    ε = c * dt^(γ + 1//2)
    p = terms_needed(length(dW), dt, ε, levyalg, MaxL2())
  end
  J .= LevyArea.levyarea(dW / √dt, p, levyalg)
  J .= 1//2 * dW .* dW' .+ dt .* J
  return nothing
end

# Default algorithms, keep KPWJ_oop() to have a non-mutating version
function get_Jalg(ΔW, dt, prob, alg)
  if alg.ii_approx isa IILevyArea
    if isinplace(prob)
      if typeof(ΔW) <: Number || is_diagonal_noise(prob)
        return JDiagonal_iip(ΔW)
      else
        # optimal_algorithm(dim, stepsize, eps=stepsize^(3/2), norm=MaxL2())
        return IteratedIntegralAlgorithm_iip(ΔW, LevyArea.optimal_algorithm(length(ΔW), dt))
      end
    else
      if typeof(ΔW) <: Number || is_diagonal_noise(prob)
        return JDiagonal_oop()
      else
        return LevyArea.optimal_algorithm(length(ΔW), dt)
      end
    end
  elseif alg.ii_approx isa IICommutative
    if isinplace(prob)
      if typeof(ΔW) <: Number || is_diagonal_noise(prob)
        return JDiagonal_iip(ΔW)
      else
        return JCommute_iip(ΔW)
      end
    else
      if typeof(ΔW) <: Number || is_diagonal_noise(prob)
        return JDiagonal_oop()
      else
        return JCommute_oop()
      end
    end
  else
    if isinplace(prob)
      IteratedIntegralAlgorithm_iip(ΔW, alg.ii_approx)
    else
      return alg.ii_approx
    end
  end
end

# Specific Levy area alg for an SDE solver
# function StochasticDiffEq.get_Jalg(ΔW,prob,alg::SOLVER)
#  return MronRoe()
# end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq

max_issues_repo_name: JuliaDiffEq/StochasticDiffEq

max_forks_repo_name: JuliaDiffEq/StochasticDiffEq


text:
using Test, GenericLinearAlgebra, LinearAlgebra, Quaternions

@testset "Singular value decomposition" begin
    @testset "Problem dimension ($m,$n)" for
        (m,n) in ((6,5)     , (6,6)     , (5,6),
                  (60, 50)  , (60, 60)  , (50, 60),
                  (200, 150), (200, 200), (150, 200))

        vals = reverse(collect(1:min(m,n)))
        U = qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:m, j = 1:min(m,n)]).Q
        V = qr(Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:n, j = 1:min(m,n)]).Q

        # FixMe! Using Array here shouldn't be necessary. Can be removed once
        # the bug in LinearAlgebra is fixed
        A = U*Array(Diagonal(vals))*V'

        @test size(A) == (m, n)
        @test vals ≈ svdvals(A)

        F = svd(A)
        @test vals ≈ F.S
        @show norm(F.U'*A*F.V - Diagonal(F.S), Inf)
        @test F.U'*A*F.V ≈ Diagonal(F.S)
    end

    @testset "The Ivan Slapničar Challenge" begin
        # This matrix used to hang (for n = 70). Thanks to Ivan Slapničar for reporting.
        n = 70
        J = Bidiagonal(0.5 * ones(n), ones(n-1), :U)
        @test GenericLinearAlgebra._svdvals!(copy(J)) ≈ svdvals(J)
        @test GenericLinearAlgebra._svdvals!(copy(J))[end] / svdvals(J)[end] - 1 < n*eps()
    end

    @testset "Compare to Base methods. Problem dimension ($m,$n)" for
        (m, n) in ((10,  9), # tall
                   (10, 10), # square
                   (9 , 10)) # wide

        A    = randn(m,n)
        Abig = big.(A)
        @test svdvals(A) ≈ Vector{Float64}(svdvals(Abig))
        @test cond(A)    ≈ Float64(cond(Abig))

        F    = svd(A)
        Fbig = svd(Abig)
        @test abs.(F.U'Float64.(Fbig.U)) ≈ I
        @test abs.(F.V'Float64.(Fbig.V)) ≈ I

        F    = svd(A, full=true)
        Fbig = svd(Abig, full=true)
        @test abs.(F.U'Float64.(Fbig.U)) ≈ I
        @test abs.(F.V'Float64.(Fbig.V)) ≈ I
    end

    @testset "Issue 54" begin
        U0, _, V0 = svd(big.(reshape(0:15, 4, 4)))
        A = U0[:, 1:3] * V0[:, 1:3]'

        U, S, V = svd(A)
        @test A ≈ U*Diagonal(S)*V'
    end

    @testset "Very small matrices. Issue 79" begin
        A = randn(1, 2)
        FA = svd(A)
        FAb = svd(big.(A))
        FAtb = svd(big.(A'))
        @test FA.S ≈ Float64.(FAb.S) ≈ Float64.(FAtb.S)
        @test abs.(FA.U'*Float64.(FAb.U))  ≈ I
        @test abs.(FA.U'*Float64.(FAtb.V)) ≈ I
        @test abs.(FA.V'*Float64.(FAb.V))  ≈ I
        @test abs.(FA.V'*Float64.(FAtb.U)) ≈ I
    end

    @testset "Issue 81" begin
        m = [1 0 0 0; 0 2 1 0; 0 1 2 0; 0 0 0 -1]
        @test Float64.(svdvals(big.(m))) ≈ svdvals(m)
    end
end


text:
using ESDL
using Test
using Dates
import Base.Iterators
using Distributed
using Statistics
addprocs(2)
@everywhere using ESDL, Statistics, NetCDF

@everywhere function sub_and_return_mean(xout1,xout2,xin)
    m=mean(skipmissing(xin))
    for i=1:length(xin)
        xout1[i]=xin[i]-m
    end
    xout2[1]=m
end
function sub_and_return_mean(c)
  mapCube(sub_and_return_mean,c,
  indims=InDims("Time"),
  outdims=(OutDims("Time"),OutDims()))
end

function doTests()
  # Test simple Stats first
  c=Cube()

  d = subsetcube(c,variable="air_temperature_2m",lon=(10,11),lat=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))

  dmem=readcubedata(d)

  @testset "Simple statistics using mapslices" begin
  # Basic statistics
  m=mapslices(mean∘skipmissing,d,dims = "Time")

  @test isapprox(readcubedata(m).data,[282.04703 281.84494 281.65927 282.2035;
    281.90726 281.59995 281.69135 282.38644;
    281.84158 281.35898 281.66083 282.49042;
    281.92462 281.12613 281.50513 282.4352])

  #Test Spatial meann along laitutde axis
  d1=subsetcube(c,variable="gross_primary_productivity",time=(Date("2002-01-01"),Date("2002-01-01")),lon=(30,30))

  dmem=readcubedata(d1)
  mtime=mapslices(mean∘skipmissing,dmem,dims = ("lon","lat"))

  end
  # Test Mean seasonal cycle retrieval
  @testset "Seasonal cycle statistics and anomalies" begin
  cdata=subsetcube(c,variable="land_surface_temperature",lon=10,lat=50.75)
  d=readcubedata(cdata)

  x2=getMSC(d)

  x3=getMedSC(d)

  a = d[:]
  a = a[3:46:end]
  @test isapprox(mean(skipmissing(a)),x2[3])
  @test isapprox(median(skipmissing(a)),x3[3])

  # Test gap filling
  cube_filled=readcubedata(gapFillMSC(d))
  imiss=findfirst(i->ismissing(i),d.data)
  @test !ismissing(cube_filled.data[imiss])
  its=mod(imiss-1,46)+1
  @test cube_filled.data[imiss]≈readcubedata(x2).data[its]
  @test !any(ismissing(cube_filled.data))

  # Test removal of MSC
  cube_anomalies=readcubedata(removeMSC(cube_filled))
  @test all(i->abs(i)<0.001, cube_anomalies.data[47:92].-(cube_filled.data[47:92].-readcubedata(x2).data[1:46]))

  # Test normalization
  anom_normalized=normalizeTS(cube_anomalies)[:]
  #@show cube_anomalies[:,:,:]
  @test mean(anom_normalized)<1e7
  @test 1.0-1e-6 <= std(anom_normalized) <= 1.0+1e-6

  #Test Polynomial fitting
  d = c[var = "land_surface_temperature"]
  dshort = d[time=2001:2003,lon=10.375,lat=51.125]
  dfill = gapFillMSC(dshort,complete_msc=true)
  @test all(!ismissing,dfill[:])
  @test dfill[1:10] == [262.4602f0,266.955f0,270.04968f0,279.886f0,268.8522f0,269.4936f0,269.0274f0,270.074f0,283.9174f0,284.86874f0]
  @test gapfillpoly(dshort)[80:90] == [283.9506f0, 285.94876f0, 284.3024f0, 282.85516f0, 282.19125f0, 276.15f0, 270.67f0, 273.86847f0, 278.861f0, 271.97f0, 271.06696f0]

  end

  d1=subsetcube(c,variable=["gross_primary_productivity","net_ecosystem_exchange"],lon=(10,10),lat=(50,50))
  d2=subsetcube(c,variable=["gross_primary_productivity","air_temperature_2m"],lon=(10,10),lat=(50,50))

  @testset "Multiple output cubes" begin
  #Test onvolving multiple output cubes
  c1=subsetcube(c,variable="gross_primary_productivity",lon=(10,11),lat=(50,51),time=Date(2001)..Date(2010))

  c2=readcubedata(c1)

  cube_wo_mean,cube_means=sub_and_return_mean(c2)

  @test isapprox(permutedims(c2[:,:,:].-mean(c2[:,:,:],dims=3),(3,1,2)),readcubedata(cube_wo_mean)[:,:,:])
  @test isapprox(mean(c2[:,:,:],dims=3)[:,:,1],cube_means[:,:])
  end
end

@testset "Parallel processing" begin
doTests()
end
rmprocs(workers())

@testset "Single proc processing" begin
doTests()
end


text:
using CSTParser: @cst_str, headof, valof

function test_iter_spans(x)
    n = 0
    for i = 1:length(x)
        a  = x[i]
        if !(a isa EXPR)
            @info i, headof(x), Expr(x)
        end
        @test a isa EXPR
        test_iter_spans(a)
        n += a.fullspan
    end
    length(x) > 0 && @test n == x.fullspan
end

@testset "iterators" begin
    @testset "const local global return" begin
        @testset "local" begin
            x = cst"local a = 1"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
        end

        @testset "global" begin
            x = cst"global a = 1"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
        end
        @testset "const" begin
            x = cst"const a = 1"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
        end
        @testset "simple" begin
            x = cst"global const a = 1"
            @test length(x) == 3
            @test x[1] === x.trivia[1]
            @test x[2] === x.trivia[2]
            @test x[3] === x.args[1]
            @test length(x[3]) == 1
            @test x[3][1] === x.args[1].args[1]
        end

        @testset "return" begin
            x = cst"return a"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
        end
    end

    @testset "datatype declarations" begin
        @testset "abstract" begin
            x = cst"abstract type T end"
            @test length(x) == 4
            @test x[1] === x.trivia[1]
            @test x[2] === x.trivia[2]
            @test x[3] === x.args[1]
            @test x[4] === x.trivia[3]
        end
        @testset "primitive" begin
            x = cst"primitive type T N end"
            @test length(x) == 5
            @test x[1] === x.trivia[1]
            @test x[2] === x.trivia[2]
            @test x[3] === x.args[1]
            @test x[4] === x.args[2]
            @test x[5] === x.trivia[3]
        end
        @testset "struct" begin
            x = cst"struct T body end"
            @test length(x) == 5
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.args[2]
            @test x[4] === x.args[3]
            @test x[5] === x.trivia[2]
        end
        @testset "mutable" begin
            x = cst"mutable struct T body end"
            @test length(x) == 6
            @test x[1] === x.trivia[1]
            @test x[2] === x.trivia[2]
            @test x[3] === x.args[1]
            @test x[4] === x.args[2]
            @test x[5] === x.args[3]
            @test x[6] === x.trivia[3]
        end
    end

    @testset ":quote" begin
        @testset "block" begin
            x = cst"""quote
                            ex1
                            ex2
                        end"""
            @test length(x) == 3
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.trivia[2]
        end
        @testset "op" begin
            x = cst""":(body + 1)"""
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
        end
    end

    @testset ":block" begin
        @testset "begin" begin
            x = cst"""begin
                            ex1
                            ex2
                        end"""
            @test length(x) == 4
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.args[2]
            @test x[4] === x.trivia[2]
        end
    end
    @testset ":for" begin
        x = cst"""for itr in itr
                        ex1
                        ex2
                    end"""
        @test length(x) == 4
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
    end

    @testset ":outer" begin
        x = cst"""for outer itr in itr
                        ex1
                        ex2
                    end""".args[1].args[1]
        @test length(x) == 2
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
    end

    @testset ":function" begin
        @testset "name only" begin
            x = cst"""function name end"""
            @test length(x) == 3
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.trivia[2]
        end
        @testset "full" begin
            x = cst"""function sig()
                            ex1
                            ex2
                        end"""
            @test length(x) == 4
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.args[2]
            @test x[4] === x.trivia[2]
        end
    end
    @testset ":braces" begin
        @testset "simple" begin
            x = cst"{a,b,c}"
            @test length(x) == 7
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.trivia[2]
            @test x[4] === x.args[2]
            @test x[5] === x.trivia[3]
            @test x[6] === x.args[3]
            @test x[7] === x.trivia[4]
        end

        @testset "params" begin
            x = cst"{a,b;c}"
            @test length(x) == 6
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[2]
            @test x[3] === x.trivia[2]
            @test x[4] === x.args[3]
            @test x[5] === x.args[1]
            @test x[6] === x.trivia[3]
        end
    end

    @testset ":curly" begin
        @testset "simple" begin
            x = cst"x{a,b}"
            @test length(x) == 6
            @test x[1] === x.args[1]
            @test x[2] === x.trivia[1]
            @test x[3] === x.args[2]
            @test x[4] === x.trivia[2]
            @test x[5] === x.args[3]
            @test x[6] === x.trivia[3]
        end
        @testset "params" begin
            x = cst"x{a,b;c}"
            @test length(x) == 7
            @test x[1] === x.args[1]
            @test x[2] === x.trivia[1]
            @test x[3] === x.args[3]
            @test x[4] === x.trivia[2]
            @test x[5] === x.args[4]
            @test x[6] === x.args[2]
            @test x[7] === x.trivia[3]
        end
    end

    @testset ":comparison" begin
        x = cst"a < b < c"
        @test length(x) == 5
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
        @test x[3] === x.args[3]
        @test x[4] === x.args[4]
        @test x[5] === x.args[5]
    end

    @testset ":using" begin
        @testset ":using" begin
            x = cst"using a"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]

            x = cst"using a, b, c"
            @test length(x) == 6
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
            @test x[3] === x.trivia[2]
            @test x[4] === x.args[2]
            @test x[5] === x.trivia[3]
            @test x[6] === x.args[3]

            x = cst"using .a"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]

            x = cst"using a: b, c"
            @test length(x) == 2
            @test x[1] === x.trivia[1]
            @test x[2] === x.args[1]
        end

        @testset ":" begin
            x = cst"using a: b, c".args[1]
            @test length(x) == 5
            @test x[1] === x.args[1]
            @test x[2] === x.head
            @test x[3] === x.args[2]
            @test x[3] === x.args[2]
        end

        @testset "." begin
            x = cst"using a".args[1]
            @test length(x) == 1
            @test x[1] === x.args[1]

            x = cst"using a.b".args[1]
            @test length(x) == 3
            @test x[1] === x.args[1]
            @test x[2] === x.trivia[1]
            @test x[3] === x.args[2]

            x = cst"using ..a.b".args[1]
            @test length(x) == 5
            @test x[1] === x.args[1]
            @test x[2] === x.args[2]
            @test x[3] === x.args[3]
            @test x[4] === x.trivia[1]
            @test x[5] === x.args[4]

            x = cst"using .a.b".args[1]
            @test length(x) == 4
            @test x[1] === x.args[1]
            @test x[2] === x.args[2]
            @test x[3] === x.trivia[1]
            @test x[4] === x.args[3]

            x = cst"using ...a.b".args[1]
            @test length(x) == 6
            @test x[1] === x.args[1]
            @test x[2] === x.args[2]
            @test x[3] === x.args[3]
            @test x[4] === x.args[4]
            @test x[5] === x.trivia[1]
            @test x[6] === x.args[5]
        end
    end

    @testset ":kw" begin
        x = cst"f(a=1)".args[2]
        @test length(x) == 3
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
    end

    @testset ":tuple" begin
        x = cst"a,b"
        @test length(x) == 3
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]

        x = cst"(a,b)"
        @test length(x) == 5
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
        @test x[3] === x.trivia[2]
        @test x[4] === x.args[2]
        @test x[5] === x.trivia[3]

        x = cst"(a,b,)"
        @test length(x) == 6
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
        @test x[3] === x.trivia[2]
        @test x[4] === x.args[2]
        @test x[5] === x.trivia[3]
        @test x[6] === x.trivia[4]
    end

    @testset ":call" begin
        x = cst"f()"
        @test length(x) == 3
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.trivia[2]

        x = cst"f(a)"
        @test length(x) == 4
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]

        x = cst"f(a,)"
        @test length(x) == 5
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
        @test x[5] === x.trivia[3]

        x = cst"f(;)"
        @test length(x) == 4
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]

        x = cst"f(a, b;c = 1)"
        @test length(x) == 7
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[3]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[4]
        @test x[6] === x.args[2]
        @test x[7] === x.trivia[3]

        x = cst"f(a, b,;c = 1)"
        @test length(x) == 8
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[3]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[4]
        @test x[6] === x.trivia[3]
        @test x[7] === x.args[2]
        @test x[8] === x.trivia[4]
    end

    @testset ":where" begin
        x = cst"a where {b,c;d}"
        @test length(x) == 8
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.trivia[2]
        @test x[4] === x.args[3]
        @test x[5] === x.trivia[3]
        @test x[6] === x.args[4]
        @test x[7] === x.args[2]
        @test x[8] === x.trivia[4]
    end

    @testset ":quotenode" begin
        x = cst"a.b".args[2]
        @test length(x) == 1
        @test x[1] === x.args[1]

        x = cst"a.:b".args[2]
        @test length(x) == 2
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
    end

    @testset ":if" begin
        x = cst"if cond end"
        @test length(x) == 4
        @test headof(x[1]) === :IF
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]
        @test headof(x[4]) === :END

        x = cst"if cond else end"
        @test length(x) == 6
        @test headof(x[1]) === :IF
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]
        @test headof(x[4]) === :ELSE
        @test x[5] === x.args[3]
        @test headof(x[6]) === :END

        x = cst"if cond args elseif a end"
        @test length(x) == 5
        @test headof(x[1]) === :IF
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]
        @test headof(x[4]) === :elseif
        @test headof(x[5]) === :END

        x = cst"a ? b : c"
        @test length(x) == 5
        @test valof(x[1]) === "a"
        @test CSTParser.isoperator(x[2])
        @test valof(x[3]) === "b"
        @test CSTParser.isoperator(x[4])
        @test valof(x[5]) === "c"
    end

    @testset ":elseif" begin
        x = cst"if cond elseif c args end".args[3]
        @test length(x) == 3
        @test headof(x[1]) === :ELSEIF
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]

        x = cst"if cond elseif c args else args end".args[3]
        @test length(x) == 5
        @test headof(x[1]) === :ELSEIF
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]
        @test headof(x[4]) === :ELSE
        @test x[5] === x.args[3]
    end

    @testset ":string" begin
        x = cst"\"txt$interp txt\""
        @test length(x) == 4
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.args[3]

        x = cst"\"txt1 $interp1 txt2 $interp2 txt3\""
        @test length(x) == 7
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.args[3]
        @test x[5] === x.trivia[2]
        @test x[6] === x.args[4]
        @test x[7] === x.args[5]

        x = cst"\"$interp\""
        @test length(x) == 4
        @test x[1] === x.trivia[1]
        @test x[2] === x.trivia[2]
        @test x[3] === x.args[1]
        @test x[4] === x.trivia[3]

        x = cst"\"$interp txt\""
        @test length(x) == 4
        @test x[1] === x.trivia[1]
        @test x[2] === x.trivia[2]
        @test x[3] === x.args[1]
        @test x[4] === x.args[2]

        x = cst"\"$(interp)\""
        @test length(x) == 6
        @test x[1] === x.trivia[1]
        @test x[2] === x.trivia[2]
        @test x[3] === x.trivia[3]
        @test x[4] === x.args[1]
        @test x[5] === x.trivia[4]
        @test x[6] === x.trivia[5]

        x = cst"\"a$b$c \""
        @test length(x) == 6
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[3]

        x = cst"\"$(a)$(b)$(c)d\""
        @test length(x) == 14
        @test x[1] === x.trivia[1]
        @test x[2] === x.trivia[2]
        @test x[3] === x.trivia[3]
        @test x[4] === x.args[1]
        @test x[5] === x.trivia[4]
        @test x[6] === x.trivia[5]
        @test x[7] === x.trivia[6]
        @test x[8] === x.args[2]
        @test x[9] === x.trivia[7]
        @test x[10] === x.trivia[8]
        @test x[11] === x.trivia[9]
        @test x[12] === x.args[3]
        @test x[13] === x.trivia[10]
        @test x[14] === x.args[4]

        x = cst"""
        "$(()$)"
        """
        @test x[6] === x.trivia[5]

        x = cst"\"$(\"\")\""
        @test length(x) == 6
        @test x[1] === x.trivia[1]
        @test x[2] === x.trivia[2]
        @test x[3] === x.trivia[3]
        @test x[4] === x.args[1]
        @test x[5] === x.trivia[4]
        @test x[6] === x.trivia[5]

        x = EXPR(:string, EXPR[cst"\" \"", EXPR(:errortoken, 0, 0), EXPR(:errortoken, 0, 0)], EXPR[cst"$"])
        @test x[4] == x.args[3]
    end

    @testset ":macrocall" begin
        x = cst"@mac a"
        @test length(x) == 3
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
        @test x[3] === x.args[3]

        x = cst"@mac(a)"
        @test length(x) == 5
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
        @test x[3] === x.trivia[1]
        @test x[4] === x.args[3]
        @test x[5] === x.trivia[2]

        x = cst"@mac(a, b)"
        @test length(x) == 7
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
        @test x[3] === x.trivia[1]
        @test x[4] === x.args[3]
        @test x[5] === x.trivia[2]
        @test x[6] === x.args[4]
        @test x[7] === x.trivia[3]

        x = cst"@mac(a; b = 1)"
        @test length(x) == 6
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
        @test x[3] === x.trivia[1]
        @test x[4] === x.args[4]
        @test x[5] === x.args[3]
        @test x[6] === x.trivia[2]

        x = cst"@mac(a, b; x)"
        @test length(x) == 8
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
        @test x[3] === x.trivia[1]
        @test x[4] === x.args[4]
        @test x[5] === x.trivia[2]
        @test x[6] === x.args[5]
        @test x[7] === x.args[3]
        @test x[8] === x.trivia[3]
    end

    @testset ":brackets" begin
        x = cst"(x)"
        @test length(x) == 3
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
        @test x[3] === x.trivia[2]
    end

    @testset ":ref" begin
        x = cst"x[i]"
        @test length(x) == 4
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]

        x = cst"x[i, j]"
        @test length(x) == 6
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[3]
        @test x[6] === x.trivia[3]

        x = cst"x[i, j; k = 1]"
        @test length(x) == 7
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[3]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[4]
        @test x[6] === x.args[2]
        @test x[7] === x.trivia[3]
    end
    @testset ":typed_vcat" begin
        x = cst"x[i;j]"
        @test length(x) == 5
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.args[3]
        @test x[5] === x.trivia[2]
    end

    @testset ":row" begin
        x = cst"[a b; c d ]".args[1]
        @test length(x) == 2
        @test x[1] === x.args[1]
        @test x[2] === x.args[2]
    end

    @testset ":module" begin
        x = cst"module a end"
        @test length(x) == 5
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
        @test x[3] === x.args[2]
        @test x[4] === x.args[3]
        @test x[5] === x.trivia[2]
    end

    @testset ":export" begin
        x = cst"export a, b, c"
        @test length(x) == 6
        @test x[1] === x.trivia[1]
        @test x[2] === x.args[1]
        @test x[3] === x.trivia[2]
        @test x[4] === x.args[2]
        @test x[5] === x.trivia[3]
        @test x[6] === x.args[3]
    end

    @testset ":parameters" begin
        x = cst"f(a; b=1, c=1, d=1)"[4]
        @test length(x) == 5
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[3]
    end

    @testset "lowered iterator" begin
        x = cst"for a in b end".args[1]
        @test length(x) == 3
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
    end

    @testset ":do" begin
        x = cst"f(x) do arg something end"
        @test length(x) == 4
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
    end

    @testset ":generator" begin
        x = cst"(a for a in A)".args[1]
        @test length(x) == 3
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]

        x = cst"(a for a in A, b in B)".args[1]
        @test length(x) == 5
        @test x[1] === x.args[1]
        @test x[2] === x.trivia[1]
        @test x[3] === x.args[2]
        @test x[4] === x.trivia[2]
        @test x[5] === x.args[3]
    end

    @testset ":flatten" begin
        function flatten(x)
            if length(x) == 0
                [x]
            else
                vcat([flatten(a) for a in x]...)
            end
        end
        function testflattenorder(s)
            x = CSTParser.parse(s)[2]
            issorted([Base.parse(Int, a.val) for a in flatten(x) if a.head === :INTEGER])
        end

        @test testflattenorder("(1 for 2 in 3)")
        @test testflattenorder("(1 for 2 in 3 for 4 in 5)")
        @test testflattenorder("(1 for 2 in 3, 4 in 5 for 6 in 7)")
        @test testflattenorder("(1 for 2 in 3 for 4 in 5, 6 in 7)")
        @test testflattenorder("(1 for 2 in 3 for 4 in 5, 6 in 7 if 8)")
    end

    @testset ":filter" begin
        x = cst"(a for a in A if a)".args[1].args[2]
        @test length(x) == 3
        @test valof(headof(x[1])) == "="
        @test headof(x[2]) === :IF
        @test valof(x[3]) == "a"
    end

    @testset ":try" begin
        x = cst"try expr catch e end"
        @test length(x) == 6
        @test headof(x[1]) === :TRY
        @test headof(x[2]) === :block
        @test headof(x[3]) === :CATCH
        @test valof(x[4]) == "e"
        @test headof(x[5]) === :block
        @test headof(x[6]) === :END

        x = cst"try expr finally expr2 end"
        @test length(x) == 8
        @test headof(x[1]) === :TRY
        @test headof(x[2]) === :block
        @test headof(x[3]) === :CATCH
        @test x[3].fullspan == 0
        @test headof(x[4]) === :FALSE
        @test headof(x[5]) === :FALSE
        @test headof(x[6]) === :FINALLY
        @test headof(x[7]) === :block
        @test headof(x[8]) === :END

        x = cst"try expr catch err expr2 finally expr3 end"
        @test length(x) == 8
        @test headof(x[1]) === :TRY
        @test headof(x[2]) === :block
        @test headof(x[3]) === :CATCH
        @test valof(x[4]) == "err"
        @test headof(x[5]) === :block
        @test headof(x[6]) === :FINALLY
        @test headof(x[7]) === :block
        @test headof(x[8]) === :END
    end

    @testset ":comprehension" begin
        x = cst"[a for a in A]"
        @test length(x) == 3
        @test headof(x[1]) === :LSQUARE
        @test headof(x[2]) === :generator
        @test headof(x[3]) === :RSQUARE
    end

    @testset ":typed_comprehension" begin
        x = cst"T[a for a in A]"
        @test length(x) == 4
        @test headof(x[1]) === :IDENTIFIER
        @test headof(x[2]) === :LSQUARE
        @test headof(x[3]) === :generator
        @test headof(x[4]) === :RSQUARE
    end

    @testset "unary syntax" begin
        x = cst"<:a"
        @test length(x) == 2
        @test headof(x[1]) === :OPERATOR
        @test headof(x[2]) === :IDENTIFIER

        x = cst">:a"
        @test length(x) == 2
        @test headof(x[1]) === :OPERATOR
        @test headof(x[2]) === :IDENTIFIER

        x = cst"::a"
        @test length(x) == 2
        @test headof(x[1]) === :OPERATOR
        @test headof(x[2]) === :IDENTIFIER

        x = cst"&a"
        @test length(x) == 2
        @test headof(x[1]) === :OPERATOR
        @test headof(x[2]) === :IDENTIFIER

        x = cst"a..."
        @test length(x) == 2
        @test headof(x[1]) === :IDENTIFIER
        @test headof(x[2]) === :OPERATOR

        x = cst"$a"
        @test length(x) == 2
        @test headof(x[1]) === :OPERATOR
        @test headof(x[2]) === :IDENTIFIER
    end
end

@testset "self test" begin
    test_iter_spans(CSTParser.parse(String(read("parser.jl")), true))
    [(@info f test_iter_spans(CSTParser.parse(String(read(f)), true))) for f in joinpath.(abspath("../src"), readdir("../src")) if endswith(f, ".jl")]
    [(@info f test_iter_spans(CSTParser.parse(String(read(f)), true))) for f in joinpath.(abspath("../src"), readdir("../src")) if endswith(f, ".jl")]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiDiscoveryV1beta1EndpointSlice <: SwaggerModel
    addressType::Any # spec type: Union{ Nothing, String } # spec name: addressType
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    endpoints::Any # spec type: Union{ Nothing, Vector{IoK8sApiDiscoveryV1beta1Endpoint} } # spec name: endpoints
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiDiscoveryV1beta1EndpointPort} } # spec name: ports

    function IoK8sApiDiscoveryV1beta1EndpointSlice(;addressType=nothing, apiVersion=nothing, endpoints=nothing, kind=nothing, metadata=nothing, ports=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("addressType"), addressType)
        setfield!(o, Symbol("addressType"), addressType)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("endpoints"), endpoints)
        setfield!(o, Symbol("endpoints"), endpoints)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiDiscoveryV1beta1EndpointSlice, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        o
    end
end # type IoK8sApiDiscoveryV1beta1EndpointSlice

const _property_map_IoK8sApiDiscoveryV1beta1EndpointSlice = Dict{Symbol,Symbol}(Symbol("addressType")=>Symbol("addressType"), Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("endpoints")=>Symbol("endpoints"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("ports")=>Symbol("ports"))
const _property_types_IoK8sApiDiscoveryV1beta1EndpointSlice = Dict{Symbol,String}(Symbol("addressType")=>"String", Symbol("apiVersion")=>"String", Symbol("endpoints")=>"Vector{IoK8sApiDiscoveryV1beta1Endpoint}", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("ports")=>"Vector{IoK8sApiDiscoveryV1beta1EndpointPort}")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1EndpointSlice))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1EndpointSlice[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1EndpointSlice[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1EndpointSlice)
    (getproperty(o, Symbol("addressType")) === nothing) && (return false)
    (getproperty(o, Symbol("endpoints")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiDiscoveryV1beta1EndpointSlice }, name::Symbol, val)
end


text:
using Base.Iterators, Base.Threads

function part1(filename, phases=6)
    function readinput()
        lines = readlines(filename)
        chars = collect.(lines)
        state = hcat(map(l -> Int.(l .== '#'), chars)...)
        plane = reshape(state, size(state)..., 1)
        L, M, N = size(plane)
        grid = zeros(Int, L + 2phases, M + 2phases, N + 2phases)
        grid[phases .+ (1:L), phases .+ (1:M), phases .+ (1:N)] = plane[:,:,:]

        grid
    end

    function update(grid)
        output = zeros(Int, size(grid)...)
        L, M, N = size(output)
        for (i, j, k) in product(map(i -> 1:i, size(output))...)
            nn = 0
            for (δi, δj, δk) in product(-1:1, -1:1, -1:1)
                if δi == δj == δk == 0
                    continue
                end
                try
                    nn += grid[i+δi, j+δj, k+δk]
                catch
                end
            end
            if grid[i, j, k] == 1
                output[i, j, k] = (nn == 2 || nn == 3)
            else
                output[i, j, k] = (nn == 3)
            end
        end
        output
    end

    grid = readinput()
    for _ in 1:phases
        grid = update(grid)
    end
    sum(grid)
end

function part2(filename, phases=6)
    function readinput()
        lines = readlines(filename)
        chars = collect.(lines)
        state = hcat(map(l -> Int.(l .== '#'), chars)...)
        plane = reshape(state, size(state)..., 1, 1)
        K, L, M, N = size(plane)
        grid = zeros(Int, K + 2phases, L + 2phases, M + 2phases, N + 2phases)
        grid[phases .+ (1:K), phases .+ (1:L), phases .+ (1:M), phases .+ (1:N)] = plane[:,:,:,:]

        grid
    end

    function update(grid)
        output = zeros(Int, size(grid)...)
        I, J, K, L = size(output)
        outer = collect(product(1:I, 1:J))
        Threads.@threads :static for (i, j) in outer
            for (k, l) in product(1:K, 1:L)
                nn = 0
                for (δi, δj, δk, δl) in product(-1:1, -1:1, -1:1, -1:1)
                    if δi == δj == δk == δl == 0
                        continue
                    end
                    try
                        nn += grid[i+δi, j+δj, k+δk, l+δl]
                    catch
                    end
                end
                if grid[i, j, k, l] == 1
                    output[i, j, k, l] = (nn == 2 || nn == 3)
                else
                    output[i, j, k, l] = (nn == 3)
                end
            end
        end
        output
    end

    grid = readinput()
    for _ in 1:phases
        grid = update(grid)
    end
    sum(grid)
end

@time begin
    @info "Solution" part1("input.txt") part2("input.txt")
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
export Minimizer, minimize!

include("./agent.jl")
include("./constraint.jl")
include("./optimizer.jl")

struct Minimizer
    xsol::Vector{Float64}
    xerr::Vector{Float64}
    buff::Vector{Float64}
    fork::Vector{Int64}
    pool::Vector{Agent}
    NP::Int
    NR::Int

    function Minimizer(ND::Int, NP::Int; NR::Int=0)
        NR   = NR < ND+1 ? ND+1 : NR
        xsol = Vector{Float64}(undef, ND)
        xerr = Vector{Float64}(undef, ND)
        buff = Vector{Float64}(undef, ND)
        fork = Vector{Int}(undef, NR)
        pool = return_agents(ND, NP)
        return new(xsol, xerr, buff, fork, pool, NP, NR)
    end
end

# @code_warntype ✓
function inits!(agents::VecO{Agent}, lb::NTuple, ub::NTuple)
    for agent in agents
        born!(agent.x, lb, ub)
    end
end

# @code_warntype ✓
function inits!(agents::VecO{Agent}, f::Function, cons::NTuple)
    fmax = -Inf
    for agent in agents
        violation = eval_violation(agent.x, cons)
        violation > 0.0 && (agent.c = violation; continue) # agent is infeasible
        agent.v = true
        agent.f = fcall(f, agent.x)
        fmax    = max(fmax, agent.f)
    end
    for agent in agents
        !agent.v && (agent.f = agent.c + fmax) # agent is infeasible
    end
end

# @code_warntype ✓
function group!(fork::VecO{Int}, agents::VecI{Agent}, NR::Int, NC::Int)
    diversity = 0.0
    @inbounds for i in eachindex(fork)
        diversity += agents[NR + 1].f - agents[i].f
    end
    if iszero(diversity) || isnan(diversity)
        fill!(fork, 1)
    else
        @inbounds for i in eachindex(fork)
            fork[i] = max(1, round(Int, NC * (agents[NR + 1].f - agents[i].f) / diversity))
        end
    end
    res = NC - sum(fork) # residue
    idx = 2
    while res > 0
        @inbounds fork[idx] += 1; res -= 1
        idx < NR ? idx += 1 : idx = 2
    end
    while res < 0
        @inbounds fork[idx] = max(1, fork[idx] - 1); res += 1
        idx < NR ? idx += 1 : idx = 2
    end
end

# @code_warntype ✓
function minimize!(obj::Minimizer, fn::Function, lb::NTuple{ND,T}, ub::NTuple{ND,T}, itmax::Int, dmax::Real, avgtimes::Int) where {ND,T<:Real}
    NP = obj.NP
    NR = obj.NR
    NC = NP - NR

    cons = boxbounds(lb, ub)
    xsol = obj.xsol
    xerr = obj.xerr
    buff = obj.buff
    fork = obj.fork

    agents = obj.pool
    elites = return_elites(agents, NR)
    throng = return_throng(agents, NR, NP)

    inits!(agents, lb, ub)
    inits!(agents, fn, cons)
    biinsertsort!(agents)

    generation = 0
    while generation < avgtimes
        generation += 1
        itcount     = 0
        @inbounds while itcount < itmax
            itcount += 1
            ss = logistic(itcount, 0.5 * itmax, -0.618, 20.0 / itmax, 2.0)
            group!(fork, agents, NR, NC)
            #### Moves throng/elites to elites/the best
            rx = 1
            fx = fork[rx]
            # move agents in throng toward elites
            for ix in eachindex(throng)
                sco_move!(buff, elites[rx].x, throng[ix].x, ss)
                check!(buff, agents, elites, throng, rx, ix, fn, cons)
                fx -= 1
                iszero(fx) && (rx += 1; fx = fork[rx])
            end
            # move agents in elites and find the best one
            for rx in 2:NR
                sco_move!(buff, elites[1].x, elites[rx].x, ss)
                check!(buff, elites, rx, fn, cons)
            end
            #### Random searching process
            for ix in 1:fork[1]
                if !(dmax < nrm2(agents[1].x, throng[ix].x, buff))
                    wca_move!(buff, agents[1].x)
                    check!(buff, agents, elites, throng, 1, ix, fn, cons)
                end
            end
            for rx in 2:NR
                if !(dmax < nrm2(agents[1].x, elites[rx].x, buff)) || !(0.1 < rand())
                    born!(buff, lb, ub)
                    check!(buff, elites, rx, fn, cons)
                end
            end
            #### Update the function-value of infeasible candidates
            fmax = -Inf
            for agent in agents
                agent.v && (fmax = max(fmax, agent.f))
            end
            for agent in agents
                !agent.v && (agent.f = agent.c + fmax)
            end

            biinsertsort!(agents)
            dmax -= dmax / itmax
        end
        @inbounds xnew = agents[1].x
        @inbounds for i in eachindex(xsol)
            xsol[i], xerr[i] = welford_step(xsol[i], xerr[i], xnew[i], generation)
        end
    end
end


text:
module Plotters

#using Interact
using Random
using RecipesBase
using DataFrames
using Dates

using ..AbsTypes
using ..Utils
import ..AbsTypes: fit!, transform!
using ..ValDateFilters

export fit!,transform!
export Plotter

## setup plotting for publication
#function setupplot(pdfoutput::Bool)
#  Plots.gr();
#  fntsm = Plots.font("sans-serif", 8);
#  fntlg = Plots.font("sans-serif", 8);
#  Plots.default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm);
#  if pdfoutput == true
#    Plots.default(size=(390,200)); #Plot canvas size
#  else
#    Plots.default(size=(500,300)); #Plot canvas size
#  end
#  return nothing
#end

struct DateVal{D<:DateTime,V<:Union{Missing,Real}}
   dtime::Vector{D}
   values::Vector{V}
end

@recipe function plot(dt::DateVal,sz=(390,200))
   marker := :none
   markertype --> :none
   seriestype := :path
   fontfamily := "sans-serif"
   titlefont := 8
   size := sz
   return (dt.dtime,dt.values)
end

"""
Plotter(
  Dict(
    :interactive => false,
    :pdfoutput => false
  )
)

Plots a TS by default but performs interactive plotting if specified during instance creation.
- `:interactive` => boolean to indicate whether to use interactive plotting with `false` as default
- `:pdfoutput` => boolean to indicate whether ouput will be saved as pdf with `false` as default

Example:

csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
pltr = Plotter(Dict(:interactive => false))

mpipeline = @pipeline csvfilter |> pltr
myplot = fit_transform!(mpipeline)

Implements: `fit!`, `transform!`
"""
mutable struct Plotter <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function Plotter(args=Dict())
      default_args = Dict(
         :name => "pltr",
         :interactive => false,
         :pdfoutput => false
      )
      #setupplot(margs[:pdfoutput])
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

"""
fit!(pltr::Plotter, features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}

Check validity of `features`: 2-column Date,Val data
"""
function fit!(pltr::Plotter, features::DataFrame, labels::Vector=[])
   ncol(features) == 2 || throw(ArgumentError("dataframe must have 2 columns: Date, Val"))
end

"""
transform!(pltr::Plotter, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Convert `missing` into `NaN` to allow plotting of discontinuities.
"""
function transform!(pltr::Plotter, features::DataFrame)
   features != DataFrame() || return DataFrame()
   ncol(features) == 2 || throw(ArgumentError("dataframe must have 2 columns: Date, Val"))
   sum(names(features) .== ("Date","Value"))  == 2 || throw(ArgumentError("wrong column names"))
   # covert missing to NaN
   df = deepcopy(features)
   df.Value = Array{Union{Missing, Float64,eltype(features.Value)},1}(missing,nrow(df))
   df.Value .= features.Value
   ndxmissing = findall(x->ismissing(x),df.Value)
   df.Value[ndxmissing] .= NaN
   #setupplot(pltr.model[:pdfoutput])
   #if pltr.model[:interactive] == true && pltr.model[:pdfoutput] == false
   #  # disable interactive plot due to Knockout.jl badly maintained (Interact.jl deps)
   #  #interactiveplot(df)
   #  pl=Plots.plot(df.Date,df.Value,xlabel="Date",ylabel="Value",legend=false,show=false);
   #  return pl
   #else
   #  pl=Plots.plot(df.Date,df.Value,xlabel="Date",ylabel="Value",legend=false,show=false);
   #  return pl
   #end
   sz = (500,300) # default
   if pltr.model[:pdfoutput] == true
      sz = (390,200)
   end
   dtval = DateVal(df.Date,df.Value)
   RecipesBase.plot(dtval,sz)
end

#function interactiveplot(df::Union{Vector,Matrix,DataFrame})
#  mlength = length(df.Value)
#  @manipulate for min in slider(1:mlength,label="min",value=1),max in slider(1:mlength,label="max",value=mlength)
#     Plots.plot(df.Date[min:max],df.Value[min:max],xlabel="Date",ylabel="Value",legend=false,show=false);
#  end
#end

end