text:
function add_mass_matrix!{T<:Number}(bas::GenBasis1d, el::Element1d, mass::AbstractMatrix{T}, λ=one(T))
    M = nmodes(bas)
    Q = nquad(bas)
    ϕ = qbasis(bas)
    wJ = jacweights(el)
    for k = 1:M
        for i = k:M
            m = 0.0
            for q = 1:Q
                m += ϕ[q,k] * ϕ[q, i] * wJ[q]
            end
            mass[k,i] += m*λ
            if k != i
                mass[i,k] += m*λ
            end
        end
    end
    return mass
end

function add_mass_matrix!{T<:Number}(bas::SEM1d{T}, el::Element1d, mass::AbstractMatrix{T}, λ=one(T))

    M = nmodes(bas)
    wJ = jacweights(el)

    for i = 1:M
        mass[i,i] += wJ[i]*λ
    end
end

mass_matrix{T<:Number}(bas::GenBasis1d, el::Element1d{T}, λ=one(T)) =
    add_mass_matrix!(bas, el, zeros(T, nmodes(bas), nmodes(bas)), λ)

function add_mass_matrix!{T<:Number}(bas::SEM1d{T}, el::Element1d, mass::AbstractMatrix{T}, λ::AbstractVector{T})

    M = nmodes(bas)
    wJ = jacweights(el)

    for i = 1:M
        mass[i,i] += wJ[i]*λ[i]
    end
end

function add_mass_matrix!{T<:Number}(bas::GenBasis1d, el::Element1d, mass::AbstractMatrix{T}, λ::AbstractVector{T})
    M = nmodes(bas)
    Q = nquad(bas)
    ϕ = qbasis(bas)
    wJ = jacweights(el)
    for k = 1:M
        for i = k:M
            m = zero(T)
            for q = 1:Q
                m += ϕ[q,k] * ϕ[q, i] * wJ[q] * λ[q]
            end
            mass[k,i] += m
            if k != i
                mass[i,k] += m
            end
        end
    end
    return mass
end

function add_stiff_matrix!{T<:Number}(bas::GenBasis1d, el::Element1d, mat::AbstractMatrix{T}, λ=one(T))
    M = nmodes(bas)
    Q = nquad(bas)
    dϕ = dqbasis(bas)
    D = diffmat(bas)
    wJ = jacweights(el)
    dξdx = deriv_ξ(el)
    for k = 1:M
        for i = k:M
            L = zero(T)
            for q = 1:Q
                L += dϕ[q,i] * dϕ[q,k] *wJ[q] * (dξdx[q]^2)*λ
            end
            mat[i,k] += L
            if k != i
                mat[k,i] += L
            end
        end
    end
    return mat
end

stiff_matrix{T<:Number}(bas::GenBasis1d, el::Element1d{T}, λ=one(T)) =
    add_stiff_matrix!(bas, el, zeros(T, nmodes(bas), nmodes(bas)), λ)

function add_stiff_matrix!{T<:Number}(bas::GenBasis1d, el::Element1d, mat::AbstractMatrix{T},
                                      λ::AbstractVector{T})
    M = nmodes(bas)
    Q = nquad(bas)
    dϕ = dqbasis(bas)
    D = diffmat(bas)
    wJ = jacweights(el)
    dξdx = deriv_ξ(el)

    for k = 1:M
        for i = k:M
            L = zero(T)
            for q = 1:Q
                L += dϕ[q,i] * dϕ[q,k] *wJ[q] * (dξdx[q]^2)*λ[q]
            end
            mat[i,k] += L
            if k != i
                mat[k,i] += L
            end
        end
    end

    return mat
end

# RHS

function add_rhs!{T<:Number}(bas::GenBasis1d{T}, el::Element1d{T}, f::AbstractVector{T},
                            Fe::AbstractVector{T})

    wJ = jacweights(el)
    ϕ = qbasis(bas)
    M = nmodes(bas)
    Q = nquad(bas)

    for k = 1:M
        F = zero(T)
        for q = 1:Q
            F += f[q] * ϕ[q,k] * wJ[q]
        end
        Fe[k] += F
    end
    return Fe
end

function add_rhs!{T<:Number}(bas::SEM1d{T}, el::Element1d{T}, f::AbstractVector{T},
                            Fe::AbstractVector{T})

    wJ = jacweights(el)
    M = nmodes(bas)

    for k = 1:M
        Fe[k] = f[k] * wJ[k]
    end
    return Fe
end


text:
"Used for testing run_inline_stability_checks"
macro enable_inline_checks(enabled::Bool)
    enable_inline_stability_checks(enabled)

    quote end
end

"""Copies the variable `x` to the global variable `x_inline_check_testing_copy`"""
macro copy(vars...)
    copy_names = @. Symbol(String(vars) * "_inline_check_testing_copy")
    esc(quote
        global $(copy_names...)
        $((:($copy_name = $var)
            for (copy_name, var) in zip(copy_names, vars))...)
    end)
end

"""For variable `x` `@test's` if the global variable
`x_inline_check_testing_copy` is a function"""
macro check(vars...)
    copy_names = @. Symbol(String(vars) * "_inline_check_testing_copy")
    esc(quote
        global $(copy_names...)
        $((quote
                @test $copy_name isa Function
                @test $(QuoteNode(var)) == nameof($copy_name)
            end
            for (copy_name, var) in zip(copy_names, vars))...)
    end)
end

macro create_func()
    esc(quote
        function created_func()
            return 0
        end
    end)
end

@enable_inline_checks true

@testset "Inline Checker" begin

    @test_nowarn begin
        @stable_function([()],
        function f()
            1
        end)

        @copy f
    end
    @check f

    @test_nowarn begin
        @stable_function([()], nothing,
        function fb()
            1
        end)

        @copy fb
    end
    @check fb

    @test_nowarn begin
        @stable_function [()] begin
            function f1()
                1
            end
            function f2()
                3
            end
        end
        @copy f1 f2
    end
    @check f1 f2

    @test_nowarn begin
        @stable_function [()] nothing begin
            function fb1()
                1
            end
            function fb2()
                3
            end
        end
        @copy fb1 fb2
    end
    @check fb1 fb2

    @test_nowarn begin
        @stable_function([(Float64,)],
        function g(x)
            if x > 0
                x
            else
                0.0
            end
        end)
        @copy g
    end
    @check g

    @test_nowarn begin
        @stable_function [(Complex{Float64},)] begin
            function g1(x)
                if x == 0
                    x
                else
                    0.0im
                end
            end

            function g2(x::Complex{T}) where T
                if x == 0
                    x
                else
                    Complex{T}(2.5)
                end
            end
        end
        @copy g1 g2
    end
    @check g1 g2

    @test_warn [r".*not stable.*", r".*eturn.*"] begin
        @stable_function([(UInt8, Float64), (UInt16, Float32)],
        function h(x, y)
            if x > 1
                1
            else
                1.0
            end
        end)
        @copy h
    end
    @check h

    @test_nowarn begin
        @stable_function([(UInt8, Float64), (UInt16, Float32)], Dict(:return => Number),
        function hb(x, y)
            if x > 1
                1
            else
                1.0
            end
        end)
        @copy hb
    end
    @check hb

    @test_warn [r".*h1.*", r".*h2.*", (result)->!occursin(result, "h3"), r".*not stable.*", r".*eturn.*"] begin
        @stable_function [(UInt8, Float64), (UInt16, Float32)] begin
            function h1(x, y)
                if x > 1
                    1
                else
                    1.0
                end
            end
            function h2(x, y)
                if x > 1
                    1
                else
                    2.0
                end
            end
            function h3(x, y)
                Float64(x)
            end
        end
        @copy h1 h2 h3
    end
    @check h1 h2 h3

    @test_nowarn @stable_function [()] i() = 1
    @test_nowarn @stable_function [()] begin
                    i1() = 2
                    i2() = 3
                end

    @test_warn [r".*not stable.*", r".*eturn.*"] begin
        @stable_function([(UInt8,), (UInt16,)],
        h(x) = if x > 1
                   1
               else
                   1.0
               end
        )
        @copy h
    end
    @check h

    @test_nowarn @stable_function [()] @create_func

    @test_nowarn @stable_function [()] begin
        @create_func
        @copy created_func
    end
    @check created_func

    function foo(x)
        if x > 0
            x
        else
            0.0
        end
    end

    @test_nowarn @stable_function [(Float64,)] foo
    @test_nowarn @stable_function [(Float64,)] nothing foo
    @test_warn [r".*not stable.*", r".*eturn.*"] @stable_function [(Int,)] foo
    @test_warn [r".*not stable.*", r".*eturn.*"] @stable_function [(Int,)] nothing foo
    @test_nowarn @stable_function [(Int,)] Dict(:return => Number) foo

    @enable_inline_checks false
    @test_nowarn @stable_function [(Int,)] foo

    @test_nowarn begin
        @stable_function([(Int,)],
        function bar(x)
            if x > 1
                1
            else
                1.0
            end
        end)
        @copy bar
    end
    @check bar

    @test_nowarn begin
        @stable_function [(Int,)] begin
            function bar1(x)
                if x > 1
                    1
                else
                    1.0
                end
            end
            function bar2(x)
                if x > 1
                    1
                else
                    12
                end
            end
        end
        @copy bar1 bar2
    end
    @check bar1 bar2
end

ext: jl


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

struct UsersProfileApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_usersProfileGet(_api::UsersProfileApi; include_labels=nothing, user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.profile.get", ["slackAuth"])
    Swagger.set_param(_ctx.query, "include_labels", include_labels)  # type Bool
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Retrieves a user's profile information.
Param: include_labels::Bool
Param: user::String
Return: Dict{String, Any}
"""
function usersProfileGet(_api::UsersProfileApi; include_labels=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileGet(_api; include_labels=include_labels, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersProfileGet(_api::UsersProfileApi, response_stream::Channel; include_labels=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileGet(_api; include_labels=include_labels, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersProfileSet(_api::UsersProfileApi; name=nothing, profile=nothing, user=nothing, value=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/users.profile.set", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_param(_ctx.form, "profile", profile)  # type String
    Swagger.set_param(_ctx.form, "user", user)  # type String
    Swagger.set_param(_ctx.form, "value", value)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set the profile information for a user.
Param: name::String
Param: profile::String
Param: user::String
Param: value::String
Return: Dict{String, Any}
"""
function usersProfileSet(_api::UsersProfileApi; name=nothing, profile=nothing, user=nothing, value=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileSet(_api; name=name, profile=profile, user=user, value=value, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersProfileSet(_api::UsersProfileApi, response_stream::Channel; name=nothing, profile=nothing, user=nothing, value=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileSet(_api; name=name, profile=profile, user=user, value=value, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export usersProfileGet, usersProfileSet


text:
#TYPES
##java.util
JArray = @jimport java.util.Array
JArrayList = @jimport java.util.ArrayList
JCollection = @jimport java.util.Collection
JHashMap = @jimport java.util.HashMap
JIterator = @jimport java.util.Iterator
JList = @jimport java.util.List
JProperties = @jimport java.util.Properties
JSet = @jimport java.util.Set

##java.io
#JReader = @jimport java.io.Reader
#JStringReader = @jimport java.io.StringReader

##edu.stanford.nlp.ling
JCoreLabel = @jimport edu.stanford.nlp.ling.CoreLabel

##edu.stanford.nlp.pipeline
JAnnotation = @jimport edu.stanford.nlp.pipeline.Annotation
JStanfordCoreNLP = @jimport edu.stanford.nlp.pipeline.StanfordCoreNLP

##edu.stanford.nlp.trees
JTree = @jimport edu.stanford.nlp.trees.Tree
JLabeledScoredTreeNode = @jimport edu.stanford.nlp.trees.LabeledScoredTreeNode
JConstituent = @jimport edu.stanford.nlp.trees.Constituent
JGrammaticalStructure = @jimport edu.stanford.nlp.trees.GrammaticalStructure
JTypedDependency = @jimport edu.stanford.nlp.trees.TypedDependency
#JTreeAnnotation = @jimport edu.stanford.nlp.trees.TreeCoreAnnotations.TreeAnnotation

##edu.stanford.nlp.parser
JDependencyParser = @jimport edu.stanford.nlp.parser.nndep.DependencyParser
JLexicalizedParser = @jimport edu.stanford.nlp.parser.lexparser.LexicalizedParser

##edu.stanford.nlp.coref
JCorefCoreAnnotations = @jimport edu.stanford.nlp.coref.CorefCoreAnnotations
JCorefChain = @jimport edu.stanford.nlp.coref.data.CorefChain
JMention = @jimport edu.stanford.nlp.coref.data.Mention

##Misc.
JSemanticGraph = @jimport edu.stanford.nlp.semgraph.SemanticGraph
JRelationTriple = @jimport edu.stanford.nlp.ie.util.RelationTriple
#JDocumentPreprocessor = @jimport edu.stanford.nlp.process.DocumentPreprocessor
#JNaturalLogicAnnotations = @jimport edu.stanford.nlp.naturalli.NaturalLogicAnnotations
#JMaxentTagger = @jimport edu.stanford.nlp.tagger.maxent.MaxentTagger

#CLASSES
## For sentence, token, and tree annotations
JSentencesAnnotationClass = classforname("edu.stanford.nlp.ling.CoreAnnotations\$SentencesAnnotation")
JTokensAnnotationClass = classforname("edu.stanford.nlp.ling.CoreAnnotations\$TokensAnnotation")
JTreeAnnotationClass = classforname("edu.stanford.nlp.trees.TreeCoreAnnotations\$TreeAnnotation")

##For dependency annotation
JBasicDependenciesAnnotationClass = classforname("edu.stanford.nlp.semgraph.SemanticGraphCoreAnnotations\$BasicDependenciesAnnotation") #inner? outer?
JCollapsedDependenciesAnnotationClass = classforname("edu.stanford.nlp.semgraph.SemanticGraphCoreAnnotations\$CollapsedDependenciesAnnotation") #inner? outer?
JCollapsedCCProcessedDependenciesAnnotationClass = classforname("edu.stanford.nlp.semgraph.SemanticGraphCoreAnnotations\$CollapsedCCProcessedDependenciesAnnotation")

##For coreference annotation
JCorefMentionsAnnotationClass = classforname("edu.stanford.nlp.coref.CorefCoreAnnotations\$CorefMentionsAnnotation")
JCorefChainAnnotationClass = classforname("edu.stanford.nlp.coref.CorefCoreAnnotations\$CorefChainAnnotation")
JRelationTriplesAnnotationClass = classforname("edu.stanford.nlp.naturalli.NaturalLogicAnnotations\$RelationTriplesAnnotation")

##For sentiment analysis
JSentimentAnnotatedTreeClass = classforname("edu.stanford.nlp.sentiment.SentimentCoreAnnotations\$SentimentAnnotatedTree")


text:
# Symbolics.jl

# [[file:~/Documents/Julia/scrap.org::*Symbolics.jl][Symbolics.jl:1]]
module Symbolics

using LinearAlgebra
using Match
using Lazy: @>
using DataStructures
using DiffRules, SpecialFunctions

include("types.jl")
include("Utils.jl")
include("SymbolicAlgebra.jl")
include("Simplification.jl")
include("Calculus.jl")
include("FunctionAlgebra.jl")
include("Structure.jl")
include("NewNumbers.jl")
include("LaTeX.jl")

export Sym, SymExpr, AbstractSym, AbstractSymExpr, @sym, Symbolic, D, ∂, simplify, UpTuple, DownTuple
export up, down, square
export @new_number

end
# Symbolics.jl:1 ends here


text:
module TestInterface
using Test

# Checking EcoBase interface
using Diversity
using EcoBase

numspecies = 10;
numcommunities = 8;
manyweights = rand(numspecies, numcommunities);
manyweights /= sum(manyweights);

@testset "EcoBase interface" begin
    species = map(n -> "Species $n", 1:numspecies)
    communities = map(n -> "SC $n", 1:numcommunities)
    ut = UniqueTypes(species)
    sc = Subcommunities(communities)
    mc = Metacommunity(manyweights, ut, sc)

    @test all(thingnames(mc) .== species)
    @test all(placenames(mc) .== communities)
    @test all(occurrences(mc) .≈ manyweights)
    @test all(richness(mc) .== repeat([numspecies], inner=numcommunities))
    @test all(occupancy(mc) .== repeat([numcommunities], inner=numspecies))
    fewerweights = deepcopy(manyweights)
    fewerweights[1, 1] = 0
    fewerweights /= sum(fewerweights)
    fmc = Metacommunity(fewerweights, ut, sc)

    @test noccupied(fmc) == numcommunities
    @test noccurring(fmc) == numspecies
    @test noccupied(fmc, 1) == numcommunities - 1
    @test noccurring(fmc, 1) == numspecies - 1
    @test nthings(fmc) == numspecies
    @test nplaces(fmc) == numcommunities
end

end


text:
# ---
# title: 14. Longest Common Prefix
# id: problem14
# author: zhwang
# date: 2022-03-15
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/longest-common-prefix/description/>
# hidden: true
# ---
#
# Write a function to find the longest common prefix string amongst an array of
# strings.
#
# If there is no common prefix, return an empty string `""`.
#
#
#
# **Example 1:**
#
#
#
#     Input: strs = ["flower","flow","flight"]
#     Output: "fl"
#
#
# **Example 2:**
#
#
#
#     Input: strs = ["dog","racecar","car"]
#     Output: ""
#     Explanation: There is no common prefix among the input strings.
#
#
#
#
# **Constraints:**
#
#   * `0 <= strs.length <= 200`
#   * `0 <= strs[i].length <= 200`
#   * `strs[i]` consists of only lower-case English letters.
#
#
## @lc code=start
using LeetCode

function longest_common_prefix(strs::Vector{String})::String
    s1, s2 = minimum(strs), maximum(strs)
    pos = findfirst(i -> s1[i] != s2[i], 1:length(s1))
    return isnothing(pos) ? "" : s1[1:(pos - 1)]
end

## @lc code=end

## @lc test=start
@testset "14.longest-common-prefix.jl" begin
    @test longest_common_prefix(["flower","flow","flight"]) == "fl"
    @test longest_common_prefix(["dog","racecar","car"]) == ""
    @test longest_common_prefix(["reflower","flow","flight"]) == ""
end
## @lc test=end

## @lc info=start
# link: [solution to 14]()
## @lc info=end


text:
export ⊙

@concrete terse struct PointwiseProductMeasure{T} <: AbstractMeasure
    data :: T

    PointwiseProductMeasure(μs...) = new{typeof(μs)}(μs)
    PointwiseProductMeasure(μs) = new{typeof(μs)}(μs)
end

Base.size(μ::PointwiseProductMeasure) = size(μ.data)

function Base.show(io::IO, μ::PointwiseProductMeasure)
    io = IOContext(io, :compact => true)
    print(io, join(string.(μ.data), " ⊙ "))
end

function Base.show_unquoted(io::IO, μ::PointwiseProductMeasure, indent::Int, prec::Int)
    io = IOContext(io, :compact => true)
    if Base.operator_precedence(:*) ≤ prec
        print(io, "(")
        show(io, μ)
        print(io, ")")
    else
        show(io, μ)
    end
    return nothing
end

Base.length(m::PointwiseProductMeasure{T}) where {T} = length(m.data)

function ⊙(μ::PointwiseProductMeasure{X}, ν::PointwiseProductMeasure{Y}) where {X,Y}
    data = (μ.data..., ν.data...)
    PointwiseProductMeasure(data...)
end

function ⊙(μ::AbstractMeasure, ν::PointwiseProductMeasure)
    data = (μ, ν.data...)
    PointwiseProductMeasure(data...)
end

function ⊙(μ::PointwiseProductMeasure, ν::N) where {N <: AbstractMeasure}
    data = (μ.data..., ν)
    PointwiseProductMeasure(data...)
end

function ⊙(μ::M, ν::N) where {M <: AbstractMeasure, N <: AbstractMeasure}
    data = (μ, ν)
    PointwiseProductMeasure(data...)
end

function ⊙(μ::AbstractMeasure, ℓ::Likelihood)
    data = (μ, ℓ)
    PointwiseProductMeasure(data...)
end

function logdensity(d::PointwiseProductMeasure, x)
    sum((logdensity(dⱼ, x) for dⱼ in d.data))
end

function sampletype(d::PointwiseProductMeasure)
    @inbounds sampletype(first(d.data))
end

basemeasure(d::PointwiseProductMeasure) =  @inbounds basemeasure(first(d.data))


text:
abstract type AbstractPlayer end

mutable struct Player <: AbstractPlayer
    id::Int
    score::Int
    position::Int
end

mutable struct DeterministicDie
    sides::Int
    value::Int
    rolls::Int
    DeterministicDie(sides) = new(sides, 0, 0)
end

struct Game
    players::Vector{AbstractPlayer}
    winning_score::Int
end

function move!(player::AbstractPlayer, spaces::Int)
    player.position = mod1(player.position + spaces, 10)
    player.score += player.position
    player
end

function roll!(die::DeterministicDie)
    die.rolls += 1
    die.value = mod1(die.value + 1, die.sides)
    die
end

function play!(game::Game, die::DeterministicDie)
    while true
        for player in game.players
            move!(player, sum(roll!(die).value for _ in 1:3))
            player.score >= game.winning_score && return player.id
        end
    end
end

function Base.parse(::Type{Player}, str::String)
    Player(parse(Int, split(str, " ")[2]), 0, parse(Int, last(split(str, ":"))))
end

players = parse.(Player, readlines("day21/start.txt"))
game = Game(players, 1000)
die = DeterministicDie(100)
play!(game, die)

p1 = die.rolls * minimum(player.score for player in game.players)

@show p1

mutable struct MultiversePlayer <: AbstractPlayer
    id::Int
    score::Int
    position::Int
    universes::Int
end

mutable struct DiracDie
    sides::Int
end

function roll3(die::DiracDie)
    outcomes = Dict{Int, Int}()
    for sum_value in [i + j + k for i in 1:die.sides for j in 1:die.sides for k in 1:die.sides]
        outcomes[sum_value] = get(outcomes, sum_value, 0) + 1
    end
    outcomes
end

function play(game::Game, die::DiracDie)
    wins = Dict((player.id, 0) for player in game.players)
    die_outcomes = roll3(die)
    countwins!(wins, game, die_outcomes)
    wins
end

function countwins!(wins::Dict{Int, Int}, game::Game, die_outcomes::Dict{Int, Int})
    player = first(game.players)
    for (sum_value, n) in die_outcomes
        player_copy = move!(copy(player), sum_value)
        player_copy.universes *= n
        if player_copy.score >= game.winning_score
            wins[player_copy.id] += player_copy.universes
        else
            other_player_copy = copy(last(game.players))
            other_player_copy.universes *= n
            countwins!(wins, Game([other_player_copy, player_copy], game.winning_score), die_outcomes)
        end
    end
end

function Base.parse(::Type{MultiversePlayer}, str::String)
    MultiversePlayer(parse(Int, split(str, " ")[2]), 0, parse(Int, last(split(str, ":"))), 1)
end

Base.copy(player::MultiversePlayer) = MultiversePlayer(player.id, player.score, player.position, player.universes)

players = parse.(MultiversePlayer, readlines("day21/start.txt"))
game = Game(players, 21)
die = DiracDie(3)
wins = play(game, die)

p2 = maximum(values(wins))

@show p2


text:
# Function for Gibbs sampling multivariate density via slice sampling
function mv_slice_sampler(x::Array{Float64}, f::Function, opts::Options)
  P = length(x)
  gxp = 0
  for p = 1:length(x)
    function fp(w)  # want to fix all but the p'th value of this array
      x[p] = w
      f(x)
    end
    # TODO: initialize at last point
    # TODO: Pass on options to univariate slice sampler
    x[p], gxp = slice_sampler(x[p], fp)
  end
  x, gxp
end
mv_slice_sampler(x::Array{Float64}, d::Union(Density,DifferentiableDensity), opts::Options) = mv_slice_sampler(x, d.f, opts)


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

# Compute coboundary_1 in 2D via product FV * EV^t with fixing of redundancies.
# The cellular complex includes non-convex cells.
# The coboundary_1 matrix is generated coherently oriented, and with the outer cell on the first row.

V = [5. 0  6  9  5  3  0  5  8  10  10  10  0  0  3  9  8  6  5;
     9  9  5  5  5  6  6  0  0   9   5   0  0  3  3  8  3  8  3]

EV = [[1,2],[3,4],[3,5],[1,10],[16,18],[7,14],[6,15],[8,13],[13,14],[10,11],
[11,12],[8,19],[9,17],[1,5],[17,19],[6,7],[14,15],[4,11],[3,18],[4,16],[8,9],[9,12],[2,7]];

FV = [[1,2,7,6,15,14,13,8,19,17,9,12,11,4,3,5],
[3,4,16,18],[6,7,14,15],[1,3,4,5,10,11,16,18],
[8,9,17,19],[1,2,7,14,13,8,9,12,11,10]]

VV = [[k] for k in 1:size(V,2)];
GL.VIEW( GL.numbering(3)((V,[VV, EV, FV])) )

# testing Lar.coboundary_1  #TODO DEBUG
convex = false;
exterior = true;
copFE = Lar.coboundary_1( V::Lar.Points, FV::Lar.Cells, EV::Lar.Cells, convex,exterior );

[ collect(1:size(copFE,2))'; Matrix(copFE) ]


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -0.7481803354345802 <= q <= 3.1247048530266452)

                     add_NL_constraint(m, :(gelu(0.30283339323429637 + -0.6237495659853991*gelu(-0.21938578694289834 + -0.903262746252536*gelu(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + 0.16468819845671145*gelu(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4]))) + 0.9377177734653888*gelu(0.5819650565710264 + -0.47161154460717736*gelu(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + -0.23723822344605994*gelu(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4])))) + gelu(-0.5162811028857681 + 0.4410139812710443*gelu(-0.21938578694289834 + -0.903262746252536*gelu(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + 0.16468819845671145*gelu(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4]))) + 0.2178047619382908*gelu(0.5819650565710264 + -0.47161154460717736*gelu(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + -0.23723822344605994*gelu(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#gtk_base.jl: Base structures functionnality for gtk GUI
#-------------------------------------------------------------------------------

import Gtk: get_gtk_property, set_gtk_property!, signal_connect, @guarded
import Gtk: GConstants.GtkOrientation, GConstants.GtkSelectionMode
import Gtk: GConstants, GAccessor
import Gtk: GConstants.GtkAlign
import Gtk: GConstants.GtkButtonsType #NONE, OK, CLOSE, CANCEL, YES_NO, OK_CANCEL
import Gtk: GConstants.GtkResponseType #NONE, REJECT, ACCEPT, DELETE_EVENT, OK, CANCEL, CLOSE, YES, NO, APPLY, HELP
import Gtk: GConstants.GtkMessageType #INFO, WARNING, QUESTION, ERROR, OTHER
import Gtk: GConstants.GtkDialogFlags
import Gtk: GConstants.GtkDirectionType #TAB_FORWARD, TAB_BACKWARD, UP, DOWN, LEFT, RIGHT
import Gtk: GConstants.GtkPackType #START, END
import Gtk: GConstants.GtkWrapMode #NONE, CHAR, WORD, WORD_CHAR
import Gtk: GConstants.GtkScrollablePolicy #MINIMUM, NATURAL
import Gtk: GConstants.GtkPolicyType #ALWAYS, AUTOMATIC, NEVER, EXTERNAL
import Gtk: GConstants.GdkModifierType #CONTROL, LOCK, SHIFT, MOD1, ...
import Gtk: GConstants.GtkAccelFlags #VISIBLE, LOCKED, MASK

#==Types
===============================================================================#
abstract type AbstractDialog; end
struct NoDialog <: AbstractDialog; end

abstract type AbstractState end #Generic state object

abstract type AbstractWndState <: AbstractState end #Identifies current state of an ExploreWnd.
struct WSNormal <: AbstractWndState; end #Default state
struct WSUpdating <: AbstractWndState; end #Updating

"Exception: Invalid state"
struct InvalidState <: Exception
	curstate::AbstractState
end

mutable struct ExploreWnd
	state::AbstractWndState
	sel::ExploreSelection
	db::HDF5.HDF5File
	wnd::_Gtk.Window
	cb_subjects::_Gtk.ComboBoxText
	rb_grade_list::Vector{_Gtk.RadioButton}
	ent_sourcedoc::_Gtk.Entry
#	rg_grades::_Gtk.RadioButtonGroup
	ls_domains::_Gtk.ListStore
	tv_domains::_Gtk.TreeView
	ls_attentes::_Gtk.ListStore
	tv_attentes::_Gtk.TreeView
	ent_attentesdesc::_Gtk.Entry
	ls_content::_Gtk.ListStore
	tv_content::_Gtk.TreeView
	editdlg::AbstractDialog
end

#==Menu builders:
===============================================================================#
function Gtk_addmenu(parent::Union{_Gtk.Menu, _Gtk.MenuBar}, name::String)
	item = _Gtk.MenuItem(name)
	mnu = _Gtk.Menu(item)
	push!(parent, item)
	return mnu
end
function Gtk_addmenuitem(mnu::_Gtk.Menu, name::String)
	item = _Gtk.MenuItem(name)
	push!(mnu, item)
	return item
end
function Gtk_addsep(mnu::_Gtk.Menu)
	push!(mnu, _Gtk.SeparatorMenuItem())
	return
end

#==Copy functions for different field views
===============================================================================#
#=Actually, data will be copied from freshly read data.  This makes the code
cleaner/simpler - but it means what is copied could get out of sync with the GUI.
=#
function _build_domain_str(v::Tuple, opt::FieldViewOpts)
	abuf = String[]
	if opt.include_id; push!(abuf, v[2]); end
	if opt.include_shortdescr||opt.include_longdescr
		push!(abuf, v[3])
	end
	return join(abuf, " ")
end
function _build_attente_str(v::Tuple, opt::FieldViewOpts)
	abuf = String[]
	if opt.include_id; push!(abuf, v[2]); end
	if opt.include_shortdescr; push!(abuf, v[3]); end
	if opt.include_longdescr; push!(abuf, v[4]); end
	return join(abuf, " ")
end
function _build_content_str(v::Tuple, opt::FieldViewOpts)
	abuf = String[]
	if opt.include_id; push!(abuf, v[2]); end
	if opt.include_shortdescr||opt.include_longdescr
		push!(abuf, v[3])
	end
	return join(abuf, " ")
end
copyfv(::AbstractFieldView, explore::ExploreWnd, opt::FieldViewOpts) = nothing
copyfv(::FViewSrcDoc, explore::ExploreWnd, opt::FieldViewOpts) =
	read_sourcedoc(explore.db, explore.sel)
function copyfv(::FViewDomain, explore::ExploreWnd, opt::FieldViewOpts)
	sel = explore.sel
	dlist = read_domain_list(explore.db, sel)
	if isnothing(dlist); return nothing; end
	if sel.domain_idx in 1:length(dlist)
		return _build_domain_str(dlist[sel.domain_idx], opt)
	end
	return nothing
end
function copyfv(::FViewAttente, explore::ExploreWnd, opt::FieldViewOpts)
	sel = explore.sel
	dlist = read_attente_list(explore.db, sel)
	if isnothing(dlist); return nothing; end
	if sel.attente_idx in 1:length(dlist)
		return _build_attente_str(dlist[sel.attente_idx], opt)
	end
	return nothing
end
function copyfv(::FViewAttente_LongDescr, explore::ExploreWnd, opt::FieldViewOpts)
	sel = explore.sel
	opt = FieldViewOpts(opt.include_id, false, true) #Copy specifically long description
	dlist = read_attente_list(explore.db, sel)
	if isnothing(dlist); return nothing; end
	if sel.attente_idx in 1:length(dlist)
		return _build_attente_str(dlist[sel.attente_idx], opt)
	end
	return nothing
end
function copyfv(::FViewContent, explore::ExploreWnd, opt::FieldViewOpts)
	sel = explore.sel
	dlist = read_content_list(explore.db, sel)
	if isnothing(dlist); return nothing; end
	if sel.content_idx in 1:length(dlist)
		return _build_content_str(dlist[sel.content_idx], opt)
	end
	return nothing
end

copyallfv(fv::AbstractFieldView, explore::ExploreWnd, opt::FieldViewOpts) =
	copyfv(fv, explore, opt)
function copyallfv(::FViewDomain, explore::ExploreWnd, opt::FieldViewOpts)
	abuf = String[]
	dlist = read_domain_list(explore.db, explore.sel)
	if isnothing(dlist); return nothing; end
	for data in dlist
		push!(abuf, _build_domain_str(data, opt))
	end
	return join(abuf, "\n")
end
function copyallfv(::FViewAttente, explore::ExploreWnd, opt::FieldViewOpts)
	abuf = String[]
	dlist = read_attente_list(explore.db, explore.sel)
	if isnothing(dlist); return nothing; end
	for data in dlist
		push!(abuf, _build_attente_str(data, opt))
	end
	return join(abuf, "\n")
end
function copyallfv(::FViewContent, explore::ExploreWnd, opt::FieldViewOpts)
	abuf = String[]
	dlist = read_content_list(explore.db, explore.sel)
	if isnothing(dlist); return nothing; end
	for data in dlist
		push!(abuf, _build_content_str(data, opt))
	end
	return join(abuf, "\n")
end

#==Helper functions
===============================================================================#
function clipboard_set_text(explore::ExploreWnd, ::Nothing)
	#Do nothing; maybe clear clipboard??
end
function clipboard_set_text(explore::ExploreWnd, txtstr::String)
	clipboard = GAccessor.clipboard(explore.wnd, GTK_SELECTION_CLIPBOARD)
	clipboard_set_text(clipboard, txtstr, -1)
end

#Returns desired field view - depending on active widget
function _activefieldview(explore::ExploreWnd)
	widgetmap = [
		(explore.ent_sourcedoc, FViewSrcDoc()),
		(explore.tv_domains, FViewDomain()),
		(explore.tv_attentes, FViewAttente()),
		(explore.ent_attentesdesc, FViewAttente_LongDescr()),
		(explore.tv_content, FViewContent()),
	]
	for (w, fview) in widgetmap
		if Gtk.get_gtk_property(w, "is-focus", Bool)
			return fview
		end
	end
	return FViewNone()
end

#==Internal "public" interface
===============================================================================#
#ensure_updating: Make sure code doesn't run unless in a "WSUpdating" state.
#otherwise, an infinite recursion might happen.
function ensure_updating(explore)
	if !isa(explore.state, WSUpdating); throw(InvalidState(explore.state)); end
	return
end
function state_set!(explore::ExploreWnd, state::AbstractWndState)
	explore.state = state
	return explore
end

subjects_updatecb!(::WSUpdating, explore::ExploreWnd) = nothing #Inhibit infinite recursion
subjects_updatecb(explore::ExploreWnd) = subjects_updatecb(explore.state, explore::ExploreWnd)

#==Public interface
===============================================================================#

#Overwrite show (inhibit dumping Gtk info):
function Base.show(io::IO, ::MIME"text/plain", explore::ExploreWnd)
	print(io, ExploreWnd, "(\"", explore.db.filename, "\")")
end

refresh!(::WSUpdating, explore::ExploreWnd) = nothing #Inhibit infinite recursion
refresh!(explore::ExploreWnd) = refresh!(explore.state, explore)

Base.close(dlg::NoDialog) = nothing
Base.close(dlg::AbstractDialog) = window_close(dlg.wnd)
Base.close(explore::ExploreWnd) = window_close(explore.wnd)

#Last line


text:
"""
   Struct for snapshot which include four fields: vz, vx, pz, px
"""
struct Snapshot{Tv<:AbstractFloat}

    vz :: Vector{Tv}
    vx :: Vector{Tv}
    pz :: Vector{Tv}
    px :: Vector{Tv}
end

"""
   constructor of snapshot
"""
function Snapshot(params::TdParams)

    data_format = params.data_format
    N           = params.Nz * params.Nx

    return Snapshot(zeros(data_format,N), zeros(data_format,N),
                    zeros(data_format,N), zeros(data_format,N))
end

"""
   Struct for wavefiled which include three fields(vz, vx, p), PML boundary are truncated
"""
struct Wavefield{Tv<:AbstractFloat}
     vz :: Vector{Tv}
     vx :: Vector{Tv}
     p  :: Vector{Tv}
end

"""
   constructor of wavefield
"""
function Wavefield(params::TdParams)

    data_format = params.data_format
    N           = params.nz * params.nx

    return Wavefield(zeros(data_format,N), zeros(data_format,N), zeros(data_format,N))
end

"""
   sample snapshot to get wavefield
"""
function sample_spt2wfd(spt::Snapshot, params::TdParams)

    # initial a empty wavefield
    wfd = Wavefield(params)
    N   = length(params.spt2wfd)

    for i = 1 : N
        j = params.spt2wfd[i]
        wfd.vz[i] = spt.vz[j]
        wfd.vx[i] = spt.vx[j]
        wfd.p[i]  = spt.pz[j] + spt.px[j]
    end

    return wfd
end

"""
   sample snapshot to get pressure field
"""
function sample_adjoint2pre!(p::Vector{Tv}, spt::Snapshot,
         params::TdParams) where {Tv <: AbstractFloat}

    i = 0
    for j in params.spt2wfd
        i = i + 1
        p[i] = spt.px[j] # pz = px = p/2 inside the modeling area
    end
    return nothing
end

"""
   Create header of regularly sampled data for snapshot
"""
function snapshot_header(params::TdParams, interval::Int64)

    # the size of data
    n1 = params.Nz
    n2 = params.Nx
    n3 = 4       # include vz, vx, pz, px
    n4 = floor(Int64, (params.nt-1)/interval) + 1

    # label
    label1="z-axis"
    label2="x-axis"
    label3="component"
    label4="time step"

    return RegularSampleHeader(n1=n1, n2=n2, n3=n3, n4=n4,
           o1=0.0      , o2=0.0      , o3=1.0, o4=0.0    ,
           d1=params.dz, d2=params.dx, d3=1.0, d4=params.dt*interval,
           label1=label1, label2=label2, label3=label3, label4=label4,
           title="snapshots", data_format=params.data_format)
end

"""
   append one snapshot to the end of provided stream
"""
function append_one_snapshot(fid::IOStream, spt::Snapshot)

    write(fid, spt.vz); write(fid, spt.vx);
    write(fid, spt.pz); write(fid, spt.px);
    flush(fid)

    return nothing
end

"""
   read one snapshot from binary file in Regularly sampled format
"""
function read_one_snapshot(path::String, it::Int)

    hdr = read_RSheader(path)

    N   = hdr.n1 * hdr.n2
    esize = sizeof(hdr.data_format)

    location = field_location[:data] + (it-1)*N*4*esize
    fid = open(path, "r")
    seek(fid, location)

    vz = zeros(hdr.data_format, N); read!(fid, vz);
    vx = zeros(hdr.data_format, N); read!(fid, vx);
    pz = zeros(hdr.data_format, N); read!(fid, pz);
    px = zeros(hdr.data_format, N); read!(fid, px);
    close(fid)

    return Snapshot(vz, vx, pz, px)
end

"""
   create a header for wavefield in regularly sampled data format
"""
function wavefield_header(params::TdParams, interval::Int64)

    # the size of data
    n1 = params.nz
    n2 = params.nx
    n3 = 3      # include(vz, vx, p)
    n4 = floor(Int64, (params.nt-1)/interval) + 1

    # label
    label1="z-axis"
    label2="x-axis"
    label3="component"
    label4="time step"

    return RegularSampleHeader(n1=n1, n2=n2, n3=n3, n4=n4,
           o1=0.0, o2=0.0, o3=1.0, o4=0.0,
           d1=params.dz, d2=params.dx, d3=1.0, d4=params.dt*interval,
           label1=label1, label2=label2, label3=label3, label4=label4,
           title="wavefield", data_format=params.data_format)
end

"""
   append one wavefield to the end of provided stream
"""
function append_one_wavefield(fid::IOStream, spt::Snapshot, params::TdParams)

    # length of one wavefield
    N = params.nz * params.nx

    # write vz
    for i = 1 : N
        idx = params.spt2wfd[i]
        write(fid, spt.vz[idx])
    end

    # write vx
    for i = 1 : N
        idx = params.spt2wfd[i]
        write(fid, spt.vx[idx])
    end

    # write p (pz+px)
    for i = 1 : N
        idx = params.spt2wfd[i]
        write(fid, spt.pz[idx]+spt.px[idx])
    end
    flush(fid)

    return nothing
end

"""
   save the wavefield at the last time step in a independent file
"""
function write_wavefield(path::String, wfd::Wavefield, params::TdParams)

    # length of wavefield
    n1 = params.nz
    n2 = params.nx
    data_format = eltype(wfd.vz)

    hdr = RegularSampleHeader(n1=n1, n2=n2, n3=3,
                              o1=0.0, o2=0.0, o3=1.0,
                              d1=params.dz, d2=params.dx, d3=1.0,
                              label1="z-axis", label2="x-axis", label3="component",
                              title="wavefield at one time step", data_format=data_format)

    fid = write_RSheader(path, hdr)
    write(fid, wfd.vz)
    write(fid, wfd.vx)
    write(fid, wfd.p )
    close(fid)

    return nothing
end

"""
   read one wavefield at the time step it
"""
function read_one_wavefield(path::String, it::Int64)

    hdr = read_RSheader(path)
    N   = hdr.n1 * hdr.n2
    esize = sizeof(hdr.data_format)

    location = field_location[:data] + (it-1)*N*3*esize
    fid = open(path, "r")
    seek(fid, location)

    vz = zeros(hdr.data_format, N); read!(fid, vz);
    vx = zeros(hdr.data_format, N); read!(fid, vx);
    p  = zeros(hdr.data_format, N); read!(fid, p );
    close(fid)

    return Wavefield(vz, vx, p)
end

"""
   create a header of regularly sampled data for pressure
"""
function pressure_header(params::TdParams, interval::Int64)

    # the size of data
    n1 = params.nz
    n2 = params.nx
    n3 = floor(Int64, (params.nt-1)/interval) + 1

    # label
    label1="z-axis"
    label2="x-axis"
    label3="time step"

    return RegularSampleHeader(n1=n1, n2=n2, n3=n3,
           o1=0.0, o2=0.0, o3=1.0,
           d1=params.dz, d2=params.dx, d3=params.dt*interval,
           label1=label1, label2=label2, label3=label3,
           title="pressure", data_format=params.data_format)
end

"""
   append one pressure to the end of provided stream
"""
function append_one_pressure(fid::IOStream, spt::Snapshot, params::TdParams)

    # length of pressure vector
    N = params.nz * params.nx

    for i = 1 : N
        idx = params.spt2wfd[i]
        tmp = spt.pz[idx] + spt.px[idx]
        write(fid, tmp)
    end
    flush(fid)

    return nothing
end

"""
   read one pressure at a particular time step specified by it
"""
function read_one_pressure(path::String, it::Int64)

    hdr = read_RSheader(path)
    N   = hdr.n1 * hdr.n2
    esize = sizeof(hdr.data_format)

    location = field_location[:data] + (it-1)*N*esize
    fid = open(path, "r")
    seek(fid, location)

    p = zeros(hdr.data_format, N)
    read!(fid, p)

    return reshape(p, hdr.n1, hdr.n2)
end

"""
   Copy the content in spt2 to spt1
"""
function copy_snapshot!(spt1::Snapshot, spt2::Snapshot)

    copyto!(spt1.vz, spt2.vz)
    copyto!(spt1.vx, spt2.vx)
    copyto!(spt1.pz, spt2.pz)
    copyto!(spt1.px, spt2.px)
    return nothing
end

"""
   copy the content in wfd2 to wfd1
"""
function copy_wavefield!(wfd1::Wavefield, wfd2::Wavefield)

    copyto!(wfd1.vz, wfd2.vz)
    copyto!(wfd1.vx, wfd2.vx)
    copyto!(wfd1.p , wfd2.p )
    return nothing
end

"""
   structure save the boundary values of wavefield at all time step. The three
component vz, vx, p are matrix size of nb-by-(nt+1), the one extra time step is used
for born approximation.
"""
struct WavefieldBound{Tv<:AbstractFloat}

    vz :: Union{Matrix{Tv}, Vector{Tv}}
    vx :: Union{Matrix{Tv}, Vector{Tv}}
    p  :: Union{Matrix{Tv}, Vector{Tv}}
end

"""
   Constructor of WavefieldBound, it provide two options, 1: include all time steps, 2: just one time step.
by setting the keyword argument flag=1 or flag=2.
"""
function WavefieldBound(params::TdParams; step_flag=1)

    # length of boundary elements for one time step
    N = length(params.spt2bnd)

    # initialize empty struct
    if step_flag == 1
       return WavefieldBound(zeros(params.data_format, N),
                             zeros(params.data_format, N),
                             zeros(params.data_format, N))

    elseif step_flag == 2
       return WavefieldBound(zeros(params.data_format, N, params.nt),
                             zeros(params.data_format, N, params.nt),
                             zeros(params.data_format, N, params.nt))
    else
       error("non-supported flag, set flag=1 include all time steps, flag=2 just one time step")
    end
end

"""
   Create the  RS header for wavefield boundary
"""
function boundary_header(params::TdParams)

    # the size of data
    n1 = length(params.spt2bnd)
    n2 = 3                      # include(vz, vx, p)

    # the extra boundary value is used to compute time derivative of source-side wavefield
    # via central finite difference
    n3 = params.nt

    # label
    label1="elements"
    label2="component"
    label3="time step"

    return RegularSampleHeader(n1=n1, n2=3, n3=n3,
                               o1=0.0, o2=0.0, o3=1.0, o4=0.0,
                               d1=0.0, d2=0.0, d3=params.dt,
                               label1=label1, label2=label2, label3=label3,
                               title="wavefield boundaries", data_format=params.data_format)
end

"""
   write the boundary values for all time steps
"""
function write_boundary(path_bnd::Ts, bnd::WavefieldBound) where {Ts<:String}

    # create header for wavefield boundaries
    hdr = boundary_header(params)

    # check the size of bnd
    if (hdr.n1, hdr.n3) != size(bnd.vz) || (hdr.n1, hdr.n3) != size(bnd.vx) || (hdr.n1, hdr.n3) != size(bnd.p)
       error("size mismatch")
    end

    # write the header of wavefield boundary
    fid = write_RSheader(path_bnd, hdr)

    # write the data of wavefield boundary
    for i = 1 : hdr.n3
        write(fid, view(bnd.vz, :, i))
        write(fid, view(bnd.vx, :, i))
        write(fid, view(bnd.p , :, i))
    end
    close(fid)

    return nothing
end

"""
   append the wavefield boundary at one time step to the end of fid
"""
function append_one_boundary(fid::IOStream, spt::Snapshot, params::TdParams)

    # length of one wavefield
    N = length(params.spt2bnd)

    # write vz
    for i = 1 : N
        idx = params.spt2bnd[i]
        write(fid, spt.vz[idx])
    end

    # write vx
    for i = 1 : N
        idx = params.spt2bnd[i]
        write(fid, spt.vx[idx])
    end

    # write p (pz+px)
    for i = 1 : N
        idx = params.spt2bnd[i]
        write(fid, spt.pz[idx]+spt.px[idx])
    end
    flush(fid)

    return nothing
end

"""
   read all the boundary value at all time steps
"""
function read_boundary(path::String)

    (hdr, d) = read_RSdata(path)
    return bnd = WavefieldBound(d[:,1,:], d[:,2,:], d[:,3,:])
end

"""
  update the content in WavefieldBound, specifically
"""
function read_one_boundary!(bnd::WavefieldBound, fid::IOStream,
         it::Int64, params::TdParams)

    # byte length of one boundary
    byte_length = length(params.spt2bnd) * sizeof(params.data_format) * 3

    # location of file pointer
    location = field_location[:data] + byte_length * (it-1)
    seek(fid, location)

    # read the boundary value
    read!(fid, bnd.vz)
    read!(fid, bnd.vx)
    read!(fid, bnd.p )

    return nothing
end

"""
   ||snapshot||
"""
function l2norm_snapshot(spt::Snapshot)

    L = zero(eltype(spt.vz))

    L = L + dot(spt.vz, spt.vz)
    L = L + dot(spt.vx, spt.vx)
    L = L + dot(spt.pz, spt.pz)
    L = L + dot(spt.px, spt.px)

    return sqrt(L)
end

"""
   spt = spt1 - spt2
"""
function minus_snapshot(spt1::Snapshot, spt2::Snapshot)

    N = length(spt1.vz)
    data_format = eltype(spt1.vz)

    vz = zeros(data_format, N)
    vx = zeros(data_format, N)
    pz = zeros(data_format, N)
    px = zeros(data_format, N)

    for i = 1 : N
        vz[i] = spt1.vz[i] - spt2.vz[i]
        vx[i] = spt1.vx[i] - spt2.vx[i]
        pz[i] = spt1.pz[i] - spt2.pz[i]
        px[i] = spt1.px[i] - spt2.px[i]
    end

    return Snapshot(vz, vx, pz, px)
end

"""
   spt = spt1 + spt2
"""
function add_snapshot(spt1::Snapshot, spt2::Snapshot)

    N = length(spt1.vz)
    data_format = eltype(spt1.vz)

    vz = zeros(data_format, N)
    vx = zeros(data_format, N)
    pz = zeros(data_format, N)
    px = zeros(data_format, N)

    for i = 1 : N
        vz[i] = spt1.vz[i] + spt2.vz[i]
        vx[i] = spt1.vx[i] + spt2.vx[i]
        pz[i] = spt1.pz[i] + spt2.pz[i]
        px[i] = spt1.px[i] + spt2.px[i]
    end

    return Snapshot(vz, vx, pz, px)
end

"""
   reverse the order of snapshot or wavefield or pressure (path_tmp) and save them in the new
file(path)
"""
function reverse_order(path::String, path_tmp::String; save_flag="snapshot")

    # read header
    hdr = read_RSheader(path_tmp)

    # write header
    fid = write_RSheader(path, hdr)

    # the length of field at one time step
    if save_flag == "snapshot" || save_flag == "wavefield"
       N = hdr.n1 * hdr.n2 * hdr.n3
       nt= hdr.n4
    elseif save_flag == "pressure"
       N = hdr.n1 * hdr.n2
       nt= hdr.n3
    end

    elength = sizeof(hdr.data_format) * N
    fid_tmp = open(path_tmp, "r")
    tmp     = zeros(hdr.data_format, N)

    for it = 1 : nt

        position = field_location[:data] + elength * (nt-it)
        seek(fid_tmp, position)

        read!(fid_tmp, tmp);
        write(fid, tmp)
    end

    close(fid_tmp)
    close(fid)

    return nothing
end


text:
abstract type BehaviorGenerator end

struct DiscreteBehaviorSet <: BehaviorGenerator
    models::Vector{BehaviorModel}
    weights::Weights
end

rand(rng::AbstractRNG, s::DiscreteBehaviorSet) = sample(rng, s.models, s.weights)

function max_accel(gen::DiscreteBehaviorSet)
    m = gen.models
    w = gen.weights
    len = length(m)
    return maximum(max_accel(m[i]) for i in 1:len if w[i] > 0.0)
end

mutable struct UniformIDMMOBIL <: BehaviorGenerator
    min_idm::IDMParam
    max_idm::IDMParam
    min_mobil::MOBILParam
    max_mobil::MOBILParam
    next_idx::Int
end

function rand(rng::AbstractRNG, g::UniformIDMMOBIL)
    idm = IDMParam(
        g.min_idm.a + rand(rng)*(g.max_idm.a - g.min_idm.a),
        g.min_idm.b + rand(rng)*(g.max_idm.b - g.min_idm.b),
        g.min_idm.T + rand(rng)*(g.max_idm.T - g.min_idm.T),
        g.min_idm.v0 + rand(rng)*(g.max_idm.v0 - g.min_idm.v0),
        g.min_idm.s0 + rand(rng)*(g.max_idm.s0 - g.min_idm.s0),
        g.min_idm.del
    )
    mobil = MOBILParam(
        g.min_mobil.p + rand(rng)*(g.max_mobil.p - g.min_mobil.p),
        g.min_mobil.b_safe + rand(rng)*(g.max_mobil.b_safe - g.min_mobil.b_safe),
        g.min_mobil.a_thr + rand(rng)*(g.max_mobil.a_thr - g.min_mobil.a_thr),
    )
    g.next_idx += 1
    return IDMMOBILBehavior(idm, mobil, g.next_idx-1)
end

#=
N   T   A
120 100 140 v0  27.8 38.9   33.3    5.6
1.5 1.8 1.0 T   1.0  2.0     1.5     0.5
2.0 4.0 1.0 s0  0.0  4.0     2.0     2.0
1.4 1.0 2.0 a   0.8  2.0     1.4     0.6
2.0 1.0 3.0 b   1.0  3.0     2.0     1.0
=#

function standard_uniform(factor=1.0; correlation::Union{Bool,Float64}=false)
    ma = 1.4;    da = 0.6
    mb = 2.0;    db = 1.0
    #mT = 1.5;    dT = 0.5
    #mT = 0.6;    dT = 0.2 ####paper1
    mT = 0.4;    dT = 0.1
    #mv0 = 33.3;  dv0 = 5.63
    mv0 = 20;  dv0 = 5.63
    #ms0 = 2.0;   ds0 = 2.0
    ms0 = 0.3;   ds0 = 0.1
    del = 4.0
    #mp = 0.5;    dp = 0.5
    mp = 0.2;    dp = 0.1
    mbsafe = 2.0;dbsafe = 1.0
    mathr = 0.1; dathr = 0.1
    max_idm = IDMParam(
        ma + factor*da,
        mb + factor*db,
        mT + factor*dT,
        mv0 + factor*dv0,
        ms0 + factor*ds0,
        del)
    min_idm = IDMParam(
        ma - factor*da,
        mb - factor*db,
        mT - factor*dT,
        mv0 - factor*dv0,
        ms0 - factor*ds0,
        del)
    max_mobil = MOBILParam(
        mp + factor*dp,
        mbsafe + factor*dbsafe,
        mathr + factor*dathr
    )
    min_mobil = MOBILParam(
        mp - factor*dp,
        mbsafe - factor*dbsafe,
        mathr - factor*dathr
    )
    if correlation == true
        return CorrelatedIDMMOBIL(min_idm, max_idm, min_mobil, max_mobil, 2)
    elseif correlation == false
        return UniformIDMMOBIL(min_idm, max_idm, min_mobil, max_mobil, 2)
    else
        return CopulaIDMMOBIL(min_idm, max_idm, min_mobil, max_mobil, correlation)
    end
end

function clip(b::IDMMOBILBehavior, gen::UniformIDMMOBIL)
    bi = b.p_idm
    bm = b.p_mobil
    mini = gen.min_idm
    maxi = gen.max_idm
    minm = gen.min_mobil
    maxm = gen.max_mobil
    return IDMMOBILBehavior(
        IDMParam(
            max(min(bi.a, maxi.a), mini.a),
            max(min(bi.b, maxi.b), mini.b),
            max(min(bi.T, maxi.T), mini.T),
            max(min(bi.v0, maxi.v0), mini.v0),
            max(min(bi.s0, maxi.s0), mini.s0),
            max(min(bi.del, maxi.del), mini.del)
        ),
        MOBILParam(
            max(min(bm.p, maxm.p), minm.p),
            max(min(bm.b_safe, maxm.b_safe), minm.b_safe),
            max(min(bm.a_thr, maxm.a_thr), minm.a_thr)
        ),
        b.idx
    )
end

mutable struct CorrelatedIDMMOBIL <: BehaviorGenerator
    min_idm::IDMParam
    max_idm::IDMParam
    min_mobil::MOBILParam
    max_mobil::MOBILParam
    next_idx::Int
end

function rand(rng::AbstractRNG, g::CorrelatedIDMMOBIL)
    agg = rand(rng)
    return create_model(g, agg)
end

function create_model(g::CorrelatedIDMMOBIL, agg::Float64)
    idm = IDMParam(
        g.min_idm.a + agg*(g.max_idm.a - g.min_idm.a),
        g.min_idm.b + agg*(g.max_idm.b - g.min_idm.b),
        g.max_idm.T + agg*(g.min_idm.T - g.max_idm.T), # T is lower for more aggressive
        g.min_idm.v0 + agg*(g.max_idm.v0 - g.min_idm.v0),
        g.max_idm.s0 + agg*(g.min_idm.s0 - g.max_idm.s0), # s0 is lower for more aggressive
        g.min_idm.del
    )
    mobil = MOBILParam(
        g.max_mobil.p + agg*(g.min_mobil.p - g.max_mobil.p), # p is lower for more aggressive
        g.min_mobil.b_safe + agg*(g.max_mobil.b_safe - g.min_mobil.b_safe),
        g.max_mobil.a_thr + agg*(g.min_mobil.a_thr - g.max_mobil.a_thr), # a_thr is lower for more aggressive
    )
    g.next_idx += 1
    return IDMMOBILBehavior(idm, mobil, g.next_idx-1)
end

function aggressiveness(gen::CorrelatedIDMMOBIL, b::IDMMOBILBehavior)
    return (b.p_idm.v0 - gen.min_idm.v0)/(gen.max_idm.v0 - gen.min_idm.v0)
end

mutable struct CopulaIDMMOBIL <: BehaviorGenerator
    min_idm::IDMParam
    max_idm::IDMParam
    min_mobil::MOBILParam
    max_mobil::MOBILParam
    copula::GaussianCopula
    next_idx::Int
end

function CopulaIDMMOBIL(min_idm::IDMParam,
                        max_idm::IDMParam,
                        min_mobil::MOBILParam,
                        max_mobil::MOBILParam,
                        cor::Float64)
    return CopulaIDMMOBIL(min_idm, max_idm,
                          min_mobil, max_mobil,
                          GaussianCopula(8, cor), 2)
end

function rand(rng::AbstractRNG, g::CopulaIDMMOBIL)
    agg = rand(rng, g.copula)
    return create_model(g, agg)
end

function create_model(g::CopulaIDMMOBIL, agg::Vector{Float64})
    @assert length(agg) == 8
    idm = IDMParam(
        g.min_idm.a + agg[1]*(g.max_idm.a - g.min_idm.a),
        g.min_idm.b + agg[2]*(g.max_idm.b - g.min_idm.b),
        g.max_idm.T + agg[3]*(g.min_idm.T - g.max_idm.T), # T is lower for more aggressive
        g.min_idm.v0 + agg[4]*(g.max_idm.v0 - g.min_idm.v0),
        g.max_idm.s0 + agg[5]*(g.min_idm.s0 - g.max_idm.s0), # s0 is lower for more aggressive
        g.min_idm.del
    )
    mobil = MOBILParam(
        g.max_mobil.p + agg[6]*(g.min_mobil.p - g.max_mobil.p), # p is lower for more aggressive
        g.min_mobil.b_safe + agg[7]*(g.max_mobil.b_safe - g.min_mobil.b_safe),
        g.max_mobil.a_thr + agg[8]*(g.min_mobil.a_thr - g.max_mobil.a_thr), # a_thr is lower for more aggressive
    )
    g.next_idx += 1
    return IDMMOBILBehavior(idm, mobil, g.next_idx-1)
end

CorrelatedIDMMOBIL(gen::Union{CopulaIDMMOBIL, UniformIDMMOBIL}) = CorrelatedIDMMOBIL(
    gen.min_idm,
    gen.max_idm,
    gen.min_mobil,
    gen.max_mobil,
    gen.next_idx
)

max_accel(gen::BehaviorGenerator) = 1.5*gen.max_idm.a

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
using Test

using CommonLH.FilesLH

struct saveS
    x1
    x2
    x3
end

function file_path()
    return joinpath(testDir, "load_save_test.jld")
end

@testset "FilesLH" begin
    xS = saveS(1.23, "abc", [1.2 2.3; 3.4 4.5]);

    fPath = file_path();

    FilesLH.save(fPath, xS);

    yS = FilesLH.load(fPath);

    @test isequal(xS.x1, yS.x1)
    @test isequal(xS.x2, yS.x2)
end

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


text:
using ChristensenSet
using Test

include("test_PolynomialIterator.jl")
include("test_RootsPlot.jl")
include("test_Computations.jl")


text:
const ∞ =  Inf
@testset "Define, construct & flatten free parameters" begin

  @parameters struct FooParam
    a::Float64 =0.5 ∈(0,1)
    b::Float64 =5.0 ∈(-∞,+∞)
    c::Float64 =4.0
    d::Float64 ∈(2,+∞)
    e::Float64
  end

  foo = FooParam(d=4.0,e=5.0)

#   v = FreeParameters.flatten(Foo, foo) # Currently failing (Foo not define:)

#   FreeParameters.flatten(Foo, foo) # Currently failing (Foo not define:)

end


text:
using SoftRobots
using MultiPoly

function test_snake()
    world = SoftRobots.World2D()
    world_state = Dict{SoftRobots.Object, SoftRobots.ObjectState}()

    r, snake_state = SoftRobots.snake()
    push!(world.objects, r)
    world_state[r] = snake_state

    x, y = generators(MPoly{Float64}, :x, :y)
    terrain = SoftRobots.FixedObject()
    terrain_state = SoftRobots.FixedObjectState(-(0.1 + 0.5x^2 - 1y))
    push!(world.objects, terrain)
    world_state[terrain] = terrain_state

    vis = SoftRobots.PyPlotVisualizer(world)

    pygui(true)
    PyPlot.ion()

    SoftRobots.update!(world, world_state, 0.001)
    for j = 1:1000
        if mod(j, 100) == 0
            SoftRobots.draw(vis, world_state)
            xlim([0,1])
            ylim([0,1])
        end
        SoftRobots.update!(world, world_state, 0.001)
    end
end

test_snake()


text:
struct EmacsDisplay <: AbstractDisplay end

function Base.display(d::EmacsDisplay, x)
    display(d, "image/png", x)
end

function register_EmacsDisplay_backend()
    for d in Base.Multimedia.displays
        if typeof(d) <: EmacsDisplay
            return
        end
    end
    # register as backend
    pushdisplay(EmacsDisplay())
    # popdisplay()
    nothing
end

register_EmacsDisplay_backend()

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
function construct_reflection(m, i, r_v)
    params=[getfield(m,pn) for pn in propertynames(m)]
    params[2]=i
    params[6]=r_v
    return typeof(m)(params...)
end


text:
using Unitful

@testset "Dipoles" begin

    dips = Dipole[]

    @testset "Create" begin

        dip1 = Dipole("Talairach", 1u"m", 2u"m", 1u"m", 0, 0, 0, 1, 1, 1)
        dip2 = Dipole("Talairach", 1u"m", 2u"m", 3u"m", 0, 0, 0, 2, 2, 2)

        dips = push!(dips, dip1)
        dips = push!(dips, dip2)

        @test length(dips) == 2
    end

    @testset "Show" begin

        show(dips[1])
        show(dips)
    end

    @testset "Mean" begin

        b = mean(dips)
        @test b.x == 1.0u"m"
        @test b.y == 2.0u"m"
        @test b.z == 2.0u"m"
    end

    @testset "Std" begin

        b = std(dips)
        @test b.x == 0.0u"m"
        @test b.y == 0.0u"m"
        @test b.z == std([1, 3])u"m"
    end
end


text:
using NDTensors
using LinearAlgebra
using Test

@testset "random_orthog" begin
  n, m = 10, 4
  O1 = random_orthog(n, m)
  @test eltype(O1) == Float64
  @test norm(transpose(O1) * O1 - Diagonal(fill(1.0, m))) < 1E-14
  O2 = random_orthog(m, n)
  @test norm(O2 * transpose(O2) - Diagonal(fill(1.0, m))) < 1E-14
end

@testset "random_unitary" begin
  n, m = 10, 4
  U1 = random_unitary(n, m)
  @test eltype(U1) == ComplexF64
  @test norm(U1' * U1 - Diagonal(fill(1.0, m))) < 1E-14
  U2 = random_unitary(m, n)
  @test norm(U2 * U2' - Diagonal(fill(1.0, m))) < 1E-14
end

nothing


text:
@testset verbose=true "Base" begin
    @testset "Numbers $T" for T in Types
        # Note that, because `Num` from Symbolics is a weird type, we have to
        # be a little more explicit below than we normally would be.  Also,
        # because of signed zeros in the float types, we have to take the
        # absolute value of the difference before comparing to zero.

        # Define basis elements
        u = Quaternion{T}(1)
        i = Quaternion{T}(𝐢)
        j = Quaternion{T}(𝐣)
        k = Quaternion{T}(𝐤)
        basis = [u, i, j, k]

        # Check basis elements
        for (index, element) in enumerate(basis)
            for component in 1:4
                if component == index
                    @test element[component] == one(T)
                else
                    @test element[component] == zero(T)
                end
            end
        end
        @test u.w==one(T) && u.x==zero(T) && u.y==zero(T) && u.z==zero(T)
        @test i.w==zero(T) && i.x==one(T) && i.y==zero(T) && i.z==zero(T)
        @test j.w==zero(T) && j.x==zero(T) && j.y==one(T) && j.z==zero(T)
        @test k.w==zero(T) && k.x==zero(T) && k.y==zero(T) && k.z==one(T)

        # Check equality with constants; note that these are *equal*, but not *identical*
        @test u == one(T) + zero(T)*𝐢 == one(T)
        @test i == 𝐢 == imx
        @test j == 𝐣 == imy
        @test k == 𝐤 == imz

        @test Quaternion{T}(1, 2, 3, 4) == Quaternion{T}(SVector{4, T}(1, 2, 3, 4))
        @test Quaternion{T}(0, 2, 3, 4) == Quaternion{T}(2, 3, 4)
        @test Quaternion{T}(1, 0, 0, 0) == Quaternion{T}(1)
        @test Quaternion(T.([1, 2, 3, 4])...) == Quaternion(SVector{4, T}(1, 2, 3, 4))
        @test Quaternion(T.([0, 2, 3, 4])...) == Quaternion(T(2), T(3), T(4))
        @test Quaternion(T.([1, 0, 0, 0])...) == Quaternion(T(1))
        if !(T<:Integer)
            @test Rotor{T}(1, 2, 3, 4) == Rotor{T}(SVector{4, T}(1, 2, 3, 4)/√T(30))
            @test Rotor{T}(0, 2, 3, 4) == Rotor{T}(2, 3, 4)
            @test Rotor{T}(1, 0, 0, 0) == Rotor{T}(1)
            @test Rotor(T.([1, 2, 3, 4])...) == Rotor(SVector{4, T}(1, 2, 3, 4)/√T(30))
            @test Rotor(T.([0, 2, 3, 4])...) == Rotor(T(2), T(3), T(4))
            @test Rotor(T.([1, 0, 0, 0])...) == Rotor(T(1))
        end
        @test QuatVec{T}(1, 2, 3, 4) == QuatVec{T}(SVector{4, T}(0, 2, 3, 4))
        @test QuatVec{T}(0, 2, 3, 4) == QuatVec{T}(2, 3, 4)
        @test QuatVec{T}(1, 0, 0, 0) == QuatVec{T}(1)
        @test QuatVec(T.([0, 2, 3, 4])...) == QuatVec(SVector{4, T}(0, 2, 3, 4))
        @test QuatVec(T.([0, 2, 3, 4])...) == QuatVec(T(2), T(3), T(4))
        @test QuatVec(T.([0, 0, 0, 0])...) == QuatVec(T(0))

        # Test indexing
        q = Quaternion(T(1), T(2), T(3), T(4))
        @test q[[1, 2]] == [T(1), T(2)]
        @test q[[3, 4]] == [T(3), T(4)]
        @test q[[4, 2]] == [T(4), T(2)]
        @test q[[4, 2, 3]] == [T(4), T(2), T(3)]
        @test q[[4, 2, 3, 1]] == [T(4), T(2), T(3), T(1)]

        # Test copy constructor and self-equality
        @test Quaternion(u) == Quaternion{T}(u) == u
        @test Quaternion(i) == Quaternion{T}(i) == i
        @test Quaternion(j) == Quaternion{T}(j) == j
        @test Quaternion(k) == Quaternion{T}(k) == k
        @test Quaternionic.wrapper(Quaternion(u)) == Quaternion
        @test Quaternionic.wrapper(Rotor(u)) == Rotor
        @test Quaternionic.wrapper(QuatVec(u)) == QuatVec
        @test u == one(T)
        @test one(T) == u
        @test i != one(T)
        @test one(T) != i
        @test j != one(T)
        @test one(T) != j
        @test k != one(T)
        @test one(T) != k
        @test isequal(u, u)
        @test !isequal(u, i)
        @test !isequal(u, j)
        @test !isequal(u, k)
        @test !isequal(i, u)
        @test isequal(i, i)
        @test !isequal(i, j)
        @test !isequal(i, k)
        @test !isequal(j, u)
        @test !isequal(j, i)
        @test isequal(j, j)
        @test !isequal(j, k)
        @test !isequal(k, u)
        @test !isequal(k, i)
        @test !isequal(k, j)
        @test isequal(k, k)

        # Check isapprox
        @test u ≈ one(T)
        @test one(T) ≈ u
        if T !== Num  # Num(1) ≉ Num(2) doesn't work
            @test i ≉ one(T)
            @test one(T) ≉ i
            @test j ≉ one(T)
            @test one(T) ≉ j
            @test k ≉ one(T)
            @test one(T) ≉ k
        end

        # Check non-sensical basis elements are not allowed
        @test_throws DomainError one(QuatVec)
        @test_throws DomainError one(QuatVec{T})
        @test_throws DomainError one(QuatVec{T}(1, 2, 3))
        @test_throws DomainError zero(Rotor)
        @test_throws DomainError zero(Rotor{T})
        @test_throws DomainError zero(Rotor{T}(1))

        # Check "real" part
        @test real(u) == one(T)
        @test real(i) == zero(T)
        @test real(j) == zero(T)
        @test real(k) == zero(T)

        # Check "imag" part
        @test imag(u) == [zero(T), zero(T), zero(T)]
        @test imag(i) == [one(T), zero(T), zero(T)]
        @test imag(j) == [zero(T), one(T), zero(T)]
        @test imag(k) == [zero(T), zero(T), one(T)]

        # Check "isreal"
        @test isreal(u)
        @test !isreal(i)
        @test !isreal(j)
        @test !isreal(k)

        # Check "isinteger"
        if T != Num
            @test isinteger(u)
            @test !isinteger(1.2u)
        end
        @test !isinteger(i)
        @test !isinteger(j)
        @test !isinteger(k)

        if T<:AbstractFloat
            # Check "isnan"
            @test !isnan(u)
            @test !isnan(i)
            @test !isnan(j)
            @test !isnan(k)
            @test isnan(T(NaN) + 0imx)
            @test isnan(T(NaN)imx)
            @test isnan(T(NaN)imy)
            @test isnan(T(NaN)imz)

            # Check "isinf"
            @test !isinf(u)
            @test !isinf(i)
            @test !isinf(j)
            @test !isinf(k)
            @test isinf(T(Inf) + 0imx)
            @test isinf(T(Inf)imx)
            @test isinf(T(Inf)imy)
            @test isinf(T(Inf)imz)
        end

        # Check "isone"
        @test isone(u)
        @test !isone(2u)
        @test !isone(i)
        @test !isone(j)
        @test !isone(k)

        # Check "flipsign"
        @test flipsign(u, 1) == u
        @test flipsign(u, -1) == -u
        @test flipsign(i, 1) == i
        @test flipsign(i, -1) == -i
        @test flipsign(j, 1) == j
        @test flipsign(j, -1) == -j
        @test flipsign(k, 1) == k
        @test flipsign(k, -1) == -k

        if T != Num
            # Check "in"
            @test u ∈ 0:2
            @test u ∉ 2:4
            @test i ∉ 0:2
            @test j ∉ 0:2
            @test k ∉ 0:2
        end

    end

    @testset "bswap" begin
        @test bswap(Quaternion(1)) == Quaternion(2^56)
        @test bswap(1imx) == (2^56)imx
        @test bswap(1imy) == (2^56)imy
        @test bswap(1imz) == (2^56)imz
        @test bswap(bswap(Quaternion(1))) == Quaternion(1)
        @test bswap(bswap(1imx)) == 1imx
        @test bswap(bswap(1imy)) == 1imy
        @test bswap(bswap(1imz)) == 1imz
    end

    @testset "hash" begin
        for T1 in [FloatTypes...; IntTypes...]
            for T2 in [FloatTypes...; IntTypes...]
                q1 = Quaternion{T1}(1, 2, 3, 4)
                q2 = Quaternion{T2}(1, 2, 3, 4)
                @test isequal(q1, q2) && hash(q1)==hash(q2)
            end
        end

        for T1 in FloatTypes
            for T2 in FloatTypes
                q1 = Quaternion(T1(0.0))
                q2 = Quaternion(T2(-0.0))
                @test !isequal(q1, q2) && hash(q1)!=hash(q2)
                q1 = Quaternion(T1(NaN))
                q2 = Quaternion(T2(NaN))
                @test isequal(q1, q2) && hash(q1)==hash(q2)
            end
        end
    end

    @testset "io" begin
        io = IOBuffer()

        Base.show(io, MIME("text/plain"), Quaternion{Float64}(1, 2, 3, 4))
        @test String(take!(io)) == "1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤"
        Base.show(io, MIME("text/plain"), Quaternion{Int64}(1, 2, 3, 4))
        @test String(take!(io)) == "1 + 2𝐢 + 3𝐣 + 4𝐤"
        Base.show(io, MIME("text/plain"), Quaternion(a-b, b*c, c/d, d+e))
        @test String(take!(io)) == "a - b + b*c𝐢 + (c / d)𝐣 + (d + e)𝐤"
        Base.show(io, MIME("text/latex"), Quaternion{Float64}(1, 2, 3, 4))
        @test String(take!(io)) == "\$1.0 + 2.0\\,\\mathbf{i} + 3.0\\,\\mathbf{j} + 4.0\\,\\mathbf{k}\$"
        Base.show(io, MIME("text/latex"), Quaternion{Int64}(1, 2, 3, 4))
        @test String(take!(io)) == "\$1 + 2\\,\\mathbf{i} + 3\\,\\mathbf{j} + 4\\,\\mathbf{k}\$"
        Base.show(io, MIME("text/latex"), Quaternion(a-b, b*c, c/d, d+e))
        @test String(take!(io)) == "\$a - b + b c\\,\\mathbf{i} + \\frac{c}{d}\\,\\mathbf{j} + \\left(d + e\\right)\\,\\mathbf{k}\$"

        for T in PrimitiveTypes
            io = IOBuffer()
            q = Quaternion{T}(1, 2, 3, 4)
            write(io, q)
            seekstart(io)
            p = read(io, typeof(q))
            @test q == p
        end
    end

end


text:
module RefractiveIndexDatabase

using Pkg.Artifacts

using YAML
using HTTP: request
using DelimitedFiles: readdlm

using Memoize
using Unitful: @u_str, uparse, uconvert, ustrip, AbstractQuantity

using Dierckx: Spline1D

abstract type Formula end
include("formulas.jl")

export get_material, load_file, load_url, search, info
include("database.jl")

export RefractiveIndex

# Data strucutres for materials stored in the database

struct Metadata
    reference::String
    comment::String
    specs::Dict{Any,Any}
end

struct FormulaN{T<:Formula}
    metadata::Metadata
    n::T
end

struct FormulaNTabulatedK{T<:Formula}
    metadata::Metadata
    n::T
    λ_k::Vector{Float64}
    k::Vector{Float64}
end

struct TabulatedN
    metadata::Metadata
    λ::Vector{Float64}
    n::Vector{Float64}
end

struct TabulatedNK
    metadata::Metadata
    λ::Vector{Float64}
    n::Vector{Float64}
    k::Vector{Float64}
end

struct TabulatedNTabulatedK
    metadata::Metadata
    λ_n::Vector{Float64}
    n::Vector{Float64}
    λ_k::Vector{Float64}
    k::Vector{Float64}
end

# Interpolation of tabulated data

struct RefractiveIndex{Tn<:Union{Formula,Spline1D},Tk<:Union{Spline1D,Nothing}}
    n::Tn
    k::Tk
end

function RefractiveIndex(m::FormulaN)
    return RefractiveIndex(m.n, nothing)
end

function RefractiveIndex(m::FormulaNTabulatedK)
    k_itp = Spline1D(m.λ_k, m.k, bc = "error")
    return RefractiveIndex(m.n, k_itp)
end

function RefractiveIndex(m::TabulatedN)
    n_itp = Spline1D(m.λ, m.n, bc = "error")
    return RefractiveIndex(n_itp, nothing)
end

function RefractiveIndex(m::TabulatedNK)
    n_itp = Spline1D(m.λ, m.n, bc = "error")
    k_itp = Spline1D(m.λ, m.k, bc = "error")
    return RefractiveIndex(n_itp, k_itp)
end

function RefractiveIndex(m::TabulatedNTabulatedK)
    n_itp = Spline1D(m.λ_n, m.n, bc = "error")
    k_itp = Spline1D(m.λ_k, m.k, bc = "error")
    return RefractiveIndex(n_itp, k_itp)
end

(m::RefractiveIndex{<:Formula,Nothing})(λ) = m.n(λ)
(m::RefractiveIndex{<:Spline1D,Nothing})(λ) = m.n(λ)
(m::RefractiveIndex{<:Formula,<:Spline1D})(λ) = m.n(λ) + m.k(λ) * im
(m::RefractiveIndex{<:Spline1D,<:Spline1D})(λ) = m.n(λ) + m.k(λ) * im

# Calulate the refractive index of a material at a given wavelength

@memoize _dim_to_micron(dim) = ustrip(uconvert(u"μm", 1.0uparse(dim)))
(m::RefractiveIndex)(λ, dim::String) = m(λ * _dim_to_micron(dim))
(m::RefractiveIndex)(λ::AbstractQuantity) = m(Float64(ustrip(uconvert(u"μm", λ))))

function Base.show(io::IO, ::MIME"text/plain", g::Metadata)
    println(io, typeof(g))
    println(io, "  reference: ", g.reference)
    println(io, "  comment: ", g.comment)
    print(io, "  specs: ", g.specs)
end

function Base.show(io::IO, ::MIME"text/plain", g::FormulaN)
    println(io, typeof(g))
    println(io, "  n: ", g.n)
    print(io, "  metadata: ", g.metadata)
end

function Base.show(io::IO, ::MIME"text/plain", g::FormulaNTabulatedK)
    println(io, typeof(g))
    println(io, "  n: ", g.n)
    println(io, "  λ_k: ", g.λ_k)
    println(io, "  k: ", g.k)
    print(io, "  metadata: ", g.metadata)
end

function Base.show(io::IO, ::MIME"text/plain", g::TabulatedN)
    println(io, typeof(g))
    println(io, "  λ: ", g.λ)
    println(io, "  n: ", g.n)
    print(io, "  metadata: ", g.metadata)
end

function Base.show(io::IO, ::MIME"text/plain", g::TabulatedNK)
    println(io, typeof(g))
    println(io, "  λ: ", g.λ)
    println(io, "  n: ", g.n)
    println(io, "  k: ", g.k)
    print(io, "  metadata: ", g.metadata)
end

function Base.show(io::IO, ::MIME"text/plain", g::TabulatedNTabulatedK)
    println(io, typeof(g))
    println(io, "  λ_n: ", g.λ_n)
    println(io, "  n: ", g.n)
    println(io, "  λ_k: ", g.λ_k)
    println(io, "  k: ", g.k)
    print(io, "  metadata: ", g.metadata)
end

end # module


text:
module ccall_test

export test_small_routine, test_large_routine, test_array_routine, test_julia

const small_routine  = joinpath(dirname(@__FILE__), "..", "deps", "small_routine")

const large_routine  = joinpath(dirname(@__FILE__), "..", "deps", "large_routine")

const array_routine  = joinpath(dirname(@__FILE__), "..", "deps", "array_routine")

function test_small_routine(invar, outvar, n)

    for i = 1:n
        ccall((:small_routine_, small_routine), Void, (Ptr{Float64},Ptr{Float64}), invar, outvar)
    end

end

function test_large_routine(invar, outvar, n)

    for i = 1:n
        ccall((:large_routine_, large_routine), Void, (Ptr{Float64},Ptr{Float64}), invar, outvar)
    end

end

function test_array_routine(invar, outvar, n)

    ccall((:array_routine_, array_routine), Void, (Ptr{Float64},Ptr{Float64},Ptr{Int64}), invar, outvar, &n)

end

function test_julia(invar, outvar, n)

    for i = 1:n
        outvar[1] = 2.0 * invar[1]
    end

end

end


text:
"""
	all_trajectories(posh::PosEvo; fixed_thr = 0.95, lost_thr = 0.05, keep_unfinished=false)

Find all active trajectories for given position `posh`.
An active trajectory is defined by being above fixing threshold (resp. below absence threshold) for 2 timebins in a row, and below (resp. above) at the next timebin.
Those can then be filtered for given characteristics, such as being seen in some frequency bin, being absent or fixed in the past, etc...

## Technical note
`start_indices`: `Y[i-1,a]>=fixed_thr && Y[i,a]>=fixed_thr` (*i.e.* mutation was fixed) AND `Y[i+1,a]<fixed_thr, *i.e.* it is now not fixed. This and the reverse for `lost_thr`.
`end_indices`: `Y[i,a] >= fixed_thr && Y[i+1,a] >= fixed_thr` (*i.e.* mutation is now fixed) AND `Y[i-1,a] < fixed_thr` (*i.e.* it was not fixed before).
"""
function all_trajectories(posh::PosEvo{A}; fixed_thr = 0.95, lost_thr = 0.05, keep_unfinished=false) where A
	X,Y,N,αβ = frequency_series(posh)
	out = Array{FrequencyTraj{A}, 1}(undef, 0)
	for a in 1:size(Y,2)
		#
		start_indices = findall(collect(2:length(X)-1)) do i
			if Y[i-1,a] >= fixed_thr && Y[i,a] >= fixed_thr # Was fixed
				if Y[i+1,a] < fixed_thr # Now not fixed
					return true
				else
					return false
				end
			elseif Y[i-1,a] <= lost_thr && Y[i,a] <= lost_thr # Was lost
				if Y[i+1,a] > lost_thr  # Now exists
					return true
				else
					return false
				end
			else
				return false
			end
		end .+ 1
		#
		end_indices = findall(collect(2:length(X)-1)) do i
			if Y[i,a] >= fixed_thr && Y[i+1,a] >= fixed_thr # Now fixed
				if Y[i-1,a] < fixed_thr # Was not fixed
					return true
				else
					return false
				end
			elseif Y[i,a] <= lost_thr && Y[i+1,a] <= lost_thr # Now lost
				if Y[i-1,a] > lost_thr # Was not lost
					return true
				else
					return false
				end
			else
				return false
			end
		end .+ 1

		# Fixing problems with boundaries : trajectories that we did not see start
		while !isempty(end_indices) && !isempty(start_indices) && end_indices[1] < start_indices[1]
			splice!(end_indices,1)
		end
		# ... or that we will not see end.
		if !isempty(end_indices) && !isempty(start_indices) && end_indices[end] < start_indices[end]
			if keep_unfinished
				push!(end_indices, length(X)-1)
			elseif !isempty(end_indices) && !isempty(start_indices)
				splice!(start_indices, length(start_indices))
			end
		elseif isempty(end_indices) && !isempty(start_indices)
			if length(start_indices)>1
				error("Empty `end_indices` and more than one in `start_indices`")
			elseif keep_unfinished
				push!(end_indices, length(X)-1)
			else
				start_indices=Int64[]
			end
		end

		# If nothing found, break
		if isempty(start_indices) || isempty(end_indices)
			continue
		end
		# Storing found trajectory
		for (is,ie) in zip(start_indices, end_indices)
			val = αβ[a]
			date = X[is]
			t = X[is:ie] .- X[is]
			freq = Y[is:ie,a]
			pop = N[is:ie]
			strains = [Array{AbstractString, 1}(undef, x) for x in pop]
			index = Dict(:start=>1, :end=>ie-is+1	, :active=>missing)
			fixation = :poly
			if Y[ie,a] >= fixed_thr
				fixation = :fixed
			elseif Y[ie,a] <= lost_thr
				fixation = :lost
			end
			push!(out, FrequencyTraj(posh.i, val, date, t, freq, pop, strains, index, fixation, Dict()))
		end
	end
	return out
end
"""
	all_trajectories(posh::Array{<:PosEvo,1}; fixed_thr = 0.95, lost_thr = 0.05, keep_unfinished=false)
"""
function all_trajectories(posh::Array{<:PosEvo,1}; fixed_thr = 0.95, lost_thr = 0.05, keep_unfinished=false)
	return cat(all_trajectories.(posh, fixed_thr=fixed_thr, lost_thr=lost_thr, keep_unfinished=keep_unfinished)..., dims=1)
end

"""
	previous_state_condition(traj::Array{<:FrequencyTraj,1}, past_state::Symbol)
"""
function previous_state_condition(traj::Array{<:FrequencyTraj,1}, past_state::Symbol)
	idx = zeros(Int64, 0)
	for (id, ft) in enumerate(traj)
		if (ft.freq[1] < 0.5 && past_state==:lost) || (ft.freq[1] > 0.5 && past_state == :fixed)
			push!(idx,id)
		end
	end
	return traj[idx]
end

"""
	frequency_condition(traj::Array{FrequencyTraj{A},1}, α; dα=0.05, shift_time=true) where A
"""
function frequency_condition(traj::Array{FrequencyTraj{A},1}, α; dα=0.05, shift_time=true) where A
	traj_ = Array{FrequencyTraj{A},1}(undef,0)
	for (id, ft) in enumerate(traj)
		for (i,(t,f)) in enumerate(zip(ft.t, ft.freq))
			if f >= α - dα && f < α + dα
				ft_ = deepcopy(ft)
				if shift_time
					ft_.t = ft_.t .- t
					ft_.date = ft_.date + t
					ft_.index[:active] = i
				end
				push!(traj_, ft_)
				break
			end
		end
	end
	return traj_
end
"""
	min_frequency_condition(traj::Array{<:FrequencyTraj,1}, α; shift_time=true)

Return trajectories that reach at least frequency `α`.
"""
function min_frequency_condition(traj::Array{<:FrequencyTraj,1}, α; shift_time=true)
	dα = (1-α)/2
	return frequency_condition(traj, α+dα, dα=dα, shift_time=shift_time)
end

"""
	population_size_condition(traj::Array{FrequencyTraj,1}, minpop)

Select trajectories that are based on a total population of at least `minpop` for their whole duration.
Values of `mode`: `:overall`, `active` or `:end`.
"""
function population_size_condition(traj::Array{<:FrequencyTraj,1}, minpop; mode=:overall)
	idx = zeros(Int64, 0)
	for (id, ft) in enumerate(traj)
		flag = true
		if mode == :overall
			for p in ft.pop[2:end-1]
				if p < minpop
					flag = false
				end
			end
		elseif mode == :active
			flag = ft.pop[ft.index[:active]] >= minpop
		elseif mode == :end
			flag = ft.pop[end] >= minpop
		else
			@error "Unrecognized mode $mode"
		end
		if flag
			push!(idx, id)
		end
	end
	return traj[idx]
end

"""
	derivative_condition(traj::Array{<:FrequencyTraj,1}; direction=:positive)

Note: Use the `:active` index of trajectories for computing derivative. Zero derivative are ignored.
"""
function derivative_condition(traj::Array{<:FrequencyTraj,1}; direction=:positive)
	idx = Int64[]
	for (id, ft) in enumerate(traj)
		df = (ft.index[:active] < 2) ? 0. : (ft.freq[ft.index[:active]] - ft.freq[ft.index[:active] - 1])
		if direction == :positive && df > 0
			push!(idx, id)
		elseif direction == :negative && df < 0
			push!(idx, id)
		end
	end
	return traj[idx]
end

"""
	get_strains!(traj::FrequencyTraj, fp::FluPop)
	get_strains!(traj::Array{<:FrequencyTraj}, fp::FluPop)

Find strains corresponding to trajectory `traj`, filling the `traj.strains` field.
"""
function get_strains!(traj::FrequencyTraj, fp::FluPop)
	for (i,t) in enumerate(traj.t)
		db = find_datebin(traj.date + t, fp)
		traj.strains[i] = [x[:strain] for x in find_strains(fp.datebin[db], traj.i, traj.val)]
	end
end
function get_strains!(traj::Array{<:FrequencyTraj}, fp::FluPop)
	for t in traj
		get_strains!(t, fp)
	end
end

"""
	compute_tree_spread!(traj::FrequencyTraj, tree::Tree)

Find how `traj` is spread in the tree: does come from one clade, two, more? This is quantified using either entropy or the sum sum of squared frequencies for the distribution of `traj` into clades.

## Notes
- This only makes sense for mutations that are *rising*. Mutations falling from a fixed position to a finite frequency will not appear in the tree at the moment they start falling!
- This function assumes `traj.strains` is already known. If it hasn't been determined before, call compute_tree_spread!(traj::FrequencyTraj, tree::Tree, fp::FluPop)
"""
function compute_tree_spread!(traj::FrequencyTraj, tree::Tree)
	traj.data[:treespread] = [Int64[] for i in 1:length(traj)]
	for (i, strains) in enumerate(traj.strains)
		class = Dict()
		for label in strains
			if haskey(tree.lleaves, label)
				# Find branch upstream of given strain that bears mutation
				r = TreeTools.find_mut_root(tree.lleaves[label], traj.i, traj.val)
				class[r] = get(class, r, 0) + 1
			end
		end
		traj.data[:treespread][i] = collect(values(class))
	end
end
"""
	compute_tree_spread!(traj::FrequencyTraj, tree::Tree, fp::FluPop)
"""
function compute_tree_spread!(traj::FrequencyTraj, tree::Tree, fp::FluPop)
	get_strains!(traj, fp)
	compute_tree_spread!(traj, tree)
end

"""
	get_regions!(traj::FrequencyTraj, fp::FluPop)
	get_regions!(traj::Array{<:FrequencyTraj}, fp::FluPop)
"""
function get_regions!(traj::FrequencyTraj, fp::FluPop)
	get_strains!(traj, fp)
	traj.data[:regions] = [Dict{String,Int64}() for i in 1:length(traj)]
	traj.data[:countries] = [Dict{String,Int64}() for i in 1:length(traj)]
	for (i,strains) in enumerate(traj.strains)
		# traj.regionspread[i] = Dict(r=>0 for r in nextstrain_regions)
		for s in strains
			!in(fp.strains[s][:region], ["?",'?']) && (traj.data[:regions][i][fp.strains[s][:region]] = get(traj.data[:regions][i], fp.strains[s][:region], 0) + 1 )
			!in(fp.strains[s][:country], ["?",'?']) && (traj.data[:countries][i][fp.strains[s][:country]] = get(traj.data[:countries][i], fp.strains[s][:country], 0) + 1 )
		end
	end
end
function get_regions!(traj::Array{<:FrequencyTraj}, fp::FluPop)
	for t in traj
		get_regions!(t, fp)
	end
end

"""
	get_date_idx(traj::FrequencyTraj, date::Date)
	get_date_idx(traj::FrequencyTraj, date)

Index at which `date` appears in `traj`. Return `nothing` if nothing is found. `date` should be convertable to a `Date` type.
"""
function get_date_idx(traj::FrequencyTraj, date::Date)
	for i in 1:length(traj)
		if traj.date + traj.t[i] == date
			return i
		end
	end
	return nothing
end
get_date_idx(traj::FrequencyTraj, date) = get_date_idx(traj, Date(date))

"""
	find_common_dates(t1::FrequencyTraj, t2::FrequencyTraj)

Return two vectors of indices corresponding to the common dates of `t1` and `t2`.
"""
function find_common_dates(t1::FrequencyTraj, t2::FrequencyTraj)
	if t1.date > t2.date
		idx2, idx1 = find_common_dates(t2, t1)
		return idx1, idx2
	end
	idx1 = Int64[]; idx2 = Int64[]
	for i1 in 1:length(t1)
		i2 = get_date_idx(t2, t1.date + t1.t[i1])
		if !isnothing(i2)
			push!(idx1, i1)
			push!(idx2, i2)
		end
	end
	return idx1, idx2
end

"""
	traj_distance(t1::FrequencyTraj, t2::FrequencyTraj; dist_type=:l1, kwargs...)

Distance types are `(:l1, :linfit)`.
"""
function traj_distance(t1::FrequencyTraj, t2::FrequencyTraj; dist_type=:strains, kwargs...)
	!in(dist_type, traj_distance_types()) && @error("Unknown distance type $dist_type.")
	if dist_type == :l1
		return traj_distance_l1(t1, t2; kwargs...)
	elseif dist_type == :linfit
		return traj_distance_linfit(t1, t2; kwargs...)
	elseif dist_type == :strains
		return traj_distance_strains(t1,t2; kwargs...)
	end
end

function traj_distance_default_kwargs(dist_type::Symbol)
	!in(dist_type, traj_distance_types()) && @error("Unknown distance type $dist_type.")
	if dist_type == :l1
		return (normalize=:full_length, common_only=false)
	elseif dist_type == :linfit
		return (min_overlap=3,)
	elseif dist_type == :strains
		return (normalize=true, min_overlap=3)
	end
end
traj_distance_types() = (:l1, :linfit, :strains)

"""
	traj_distance_strains(t1::FrequencyTraj, t2::FrequencyTraj; min_overlap=3)
"""
function traj_distance_strains(t1::FrequencyTraj, t2::FrequencyTraj; normalize=true, min_overlap=3)
	if t1.date > t2.date
		return traj_distance_strains(t2, t1, min_overlap=min_overlap, normalize=normalize)
	end
	#
	idx1, idx2 = find_common_dates(t1,t2)
	if length(idx1) < min_overlap
		return 1.
	end
	#
	dif = 0.
	Z = 0.
	for (i1,i2) in zip(idx1, idx2)
		# if isempty(t1.strains[i1]) || isempty(t2.strains[i2])
		# 	@error("Trajectory without strain :$(t1.i)$(t1.val)")
		# end
		dif += 1 - jaccard(t1.strains[i1], t2.strains[i2])
		if !isempty(t1.strains[i1],) || !isempty(t2.strains[i2],)
			Z += 1
		end
	end
	if normalize
		return dif / Z
	else
		return dif
	end
end
function jaccard(A,B)
	try
		if isempty(A) && isempty(B)
			return 1.
		elseif !isempty(A) && !isempty(B)
			length(intersect(A,B)) / length(union(A,B))
		else
			return 0.
		end
	catch err
		println(A)
		println(B)
		@error(err)
	end
end

"""
	traj_distance_linfit(t1::FrequencyTraj, t2::FrequencyTraj; min_overlap=3)
"""
function traj_distance_linfit(t1::FrequencyTraj, t2::FrequencyTraj; min_overlap=3)
	if t1.date > t2.date
		return traj_distance_linfit(t2,t1,min_overlap=min_overlap)
	end
	#
	idx1, idx2 = find_common_dates(t1,t2)
	if length(idx1) < min_overlap
		return 1.
	end
	#
	f1 = t1.freq[idx1]
	f2 = t2.freq[idx2]
	line = Polynomials.fit(f1,f2,1)
	return round(abs(line[1] - 1.), digits=5)
end
"""
	traj_distance_l1(t1::FrequencyTraj, t2::FrequencyTraj; normalize=true, common_only=false)

Sum `abs(f1-f2)` for all the length of trajectories `t1` and `t2`. When only one trajectory is active, the value taken for the other is either 0. or 1.
"""
function traj_distance_l1(t1::FrequencyTraj, t2::FrequencyTraj; normalize=false, common_only=false)
	if t1.date > t2.date
		return traj_distance_l1(t2,t1,normalize=normalize, common_only=common_only)
	end
	#
	idx1, idx2 = find_common_dates(t1,t2)
	# If common support is null
	α = 0.
	if length(idx1) <= α * length(t1) || length(idx2) <= α * length(t2)
		if normalize
			return 1
		else
			return max(length(t1), length(t2))
		end
	end
	#
	f2 = t2.freq[1] < 0.5 ? 0. : 1. # Initial value of t2
	dif = 0.
	Z = 0 # Normalization
	# Before t2 starts
	if !common_only
		i1 = 1
		while i1 < idx1[1]
			dif += abs(t1.freq[i1] - f2)
			Z += 1
			i1 += 1
		end
	end

	# Common part of both trajectories
	for (i1,i2) in zip(idx1, idx2)
		dif += abs(t1.freq[i1] - t2.freq[i2])
		Z += 1
	end

	# Continue counting differences if one of the trajectories has ended and not the other
	if !common_only
		i1 = idx1[end]; i2 = idx2[end]
		while i1 < length(t1)-1 || i2 < length(t2)-1
			i1 += 1
			i2 += 1
			if i1 > length(t1)
				f1 = t1.fixation==:fixed ? 1. : 0.
			else
				f1 = t1.freq[i1]
			end
			if i2 > length(t2)
				f2 = t2.fixation==:fixed ? 1. : 0.
			else
				f2 = t2.freq[i2]
			end
			dif += abs(f1 - f2)
			Z += 1
		end
	end
	#
	if !normalize
		Z = 1
	end
	return sum(diff) / Z
end

"""
	compute_distance_matrix(f, trajectories)

Compute pairwise distance matrix of `trajectories` for distance function `f`.
"""
function compute_distance_matrix(f, trajectories)
	dist_mat = zeros(Float64, length(trajectories), length(trajectories))
	for (i,t1) in enumerate(trajectories)
	    for j in (i+1):length(trajectories)
	        dist_mat[i,j] = f(t1,trajectories[j])
	        dist_mat[j,i] = dist_mat[i,j]
	    end
	    dist_mat[i,i] = 0.
	end
	return dist_mat
end

"""
	cluster_trajectories_naive(f, trajectories; thr=0.)

Group together trajectories that are exactly the same. Distance function is `f`.
"""
function cluster_trajectories_naive(f, trajectories; thr=0.)
	dist_mat = compute_distance_matrix(f, trajectories)
	return cluster_trajectories_naive(dist_mat, trajectories, thr=thr)
end
"""
	cluster_trajectories_naive(dist_mat::Array{<:Real,2}, trajectories; thr=0.)
"""
function cluster_trajectories_naive(dist_mat::Array{<:Real,2}, trajectories; thr=0.)
	clusters = Dict{Int64,Array{Int64,1}}()
	clustered = zeros(Bool, length(trajectories)) # Is the trajectory in a cluster
	for i in 1:length(trajectories)
	    clustered[i] && continue
	    clusters[i] = Int64[i]
	    clustered[i] = true
	    increase_cluster_naive!(clusters[i], clustered, dist_mat, thr)
	end
	return clusters
end
function increase_cluster_naive!(cluster, clustered::Array{Bool,1}, dist_mat, thr::Real)
	cflag = true
	idx = 1
	while cflag
		cflag = false
		for i in idx:length(cluster)
			for j in 1:size(dist_mat,1)
				if !clustered[j] && dist_mat[cluster[i],j] <= thr # If j is not clustered already and is close to `cluster[i]`
					push!(cluster, j)
					clustered[j] = true
					cflag = true
				end
			end
		end
		idx += 1
	end
end

function traj_clustering_default_kwargs(clustering_type::Symbol)
	if clustering_type == :naive
		return (thr=0.,)
	else
		@error("Unknown clustering type $clustering_type")
	end
end

"""
	effective_trajectories(trajectories; clustering=:naive, distance=:l1, common_only=true)

Cluster trajectories based on similarity, and return one representative per cluster
"""
function effective_trajectories(trajectories;
	clustering=:naive, clust_kwargs=traj_clustering_default_kwargs(clustering),
	distance=:strains, dist_kwargs=traj_distance_default_kwargs(distance))

	# Could vary based on kwargs
	dfunc(t1, t2) = traj_distance(t1,t2; dist_type=distance, dist_kwargs...)
	# clusters -- could also vary based on kwargs
	clusters = cluster_trajectories_naive(dfunc, trajectories; clust_kwargs...)

	# Picking one representative per cluster
	efftraj = []
	for (i,idx) in clusters
		push!(efftraj, trajectories[idx[rand(1:length(idx))]])
	end
	return efftraj
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta2
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# [old version of findmax(f, domain)](
# https://github.com/JuliaLang/julia/pull/35316/files#diff-97ad8b900c30b39398c32d78754c404f4a6df9a153d45284b2c11ab54837deb9R797)
#
# ```julia
# findmax(f, domain) = mapfoldl(x -> (f(x), x), _rf_findmax, domain)
#  _rf_findmax((fm, m), (fx, x)) = isless(fm, fx) ? (fx, x) : (fm, m)
# ```
#
# [current version](https://github.com/JuliaLang/julia/pull/41076/files#diff-97ad8b900c30b39398c32d78754c404f4a6df9a153d45284b2c11ab54837deb9R803)
#
# ```julia
# findmax(f, domain) = mapfoldl( ((k, v),) -> (f(v), k), _rf_findmax, pairs(domain) )
#  _rf_findmax((fm, im), (fx, ix)) = isless(fm, fx) ? (fx, ix) : (fm, im)
# ```
#
# Ref. [`findmax` and friends: confusing behaviour to be introduced in 1.7](https://discourse.julialang.org/t/findmax-and-friends-confusing-behaviour-to-be-introduced-in-1-7/61904)

# %%
using OffsetArrays

# %%
VERSION

# %%
f(x, y) = cos(x)*sin(y) + 0.1(x - y)

X = range(-2, 2; length=401)
Y = range(-2, 2; length=401)
XtimesY = Iterators.product(X, Y)

# %%
val, idx = findmax(XtimesY) do (x, y); f(x, y) end

# %% [markdown]
# `pairs` function used in `findmax(f, domain)` function requires `keys` method.

# %% tags=[]
methods(keys)

# %%
keys(collect(XtimesY))

# %%
val, idx = findmax(collect(XtimesY)) do (x, y); f(x, y) end

# %% [markdown]
# But `collect` causes memory allocations and I do not feel like it.

# %%
axes(XtimesY)

# %%
CartesianIndices(axes(XtimesY))

# %%
Base.keys(pr::Iterators.ProductIterator) = CartesianIndices(axes(pr))
F((x, y)) = f(x, y)

@show val, idx = findmax(F, XtimesY)
@show f(X[idx[1]], Y[idx[2]])
@show argmax(F.(XtimesY));

# %%
#valargmax(f, X) = (x = argmax(f, X); (f(x), x))
struct ValArg{F} <:Function f::F end
(valarg::ValArg)(x) = (valarg.f(x), x)
valargmax(f, X) = mapfoldl(ValArg(f), Base._rf_findmax, X)
valargmax(X) = valargmax(Base.Fix1(getindex, X), keys(X))
F((x, y)) = f(x, y)

@show val, idx = findmax(F, XtimesY)
@show val, (X[idx[1]], Y[idx[2]])
@show argmax(F, XtimesY)
@show valargmax(F, XtimesY)
@show valargmax(F.(XtimesY))
@show findmax(F, XtimesY);

# %%
Base.keys(rv::Iterators.Reverse) = keys(reverse(rv.itr))

@show val, idx = findmax(sin, Iterators.Reverse(X))
@show sin(reverse(X)[idx]);

# %%
Base.keys(en::Iterators.Enumerate) = keys(en.itr)
G((i, x)) = sin(x)

@show val, idx = findmax(G, enumerate(X))
@show sin(X[idx]);

# %%
Base.keys(zp::Iterators.Zip) = Base.OneTo(length(zp))
H((x, y)) = cos(x) * sin(y)

@show val, idx = findmax(H, zip(X, reverse(Y)))
@show cos(X[idx]) * sin(reverse(Y)[idx]);

# %%
Base.keys(ac::Iterators.Accumulate) = keys(ac.itr)
A = OffsetArray(range(-2, 2; length=401), -200:200)

@show val, idx = findmax(sin, Iterators.accumulate(+, A; init=π))
@show sin(sum(A[begin:idx]) + π);

# %%
Base.keys(tk::Iterators.Take) = Base.OneTo(length(tk))

@show val, idx = findmax(sin, Iterators.take(reverse(X), 100))
@show sin(reverse(X)[idx]);

# %%
methods(keys, Main)

# %%
?keys

# %% [markdown]
# https://github.com/JuliaLang/julia/blob/16f433bb13cfc87eea21d26a797dac7b34a41d86/base/abstractdict.jl#L71

# %%
?Base.isgreater

# %%
?findmax

# %%
findmax([8, 0.1, -9, pi])

# %%
A = [1.0 2; 3 4]
findmax(A, dims=1)

# %%
methods(findmax)

# %%
findmax(-(-4:5).^2)

# %%
findmax(x -> -x^2, -4:5)

# %%
valindargmax(f, X) = valargmax(f∘last, pairs(X))

X = range(-2, 2; length=401)
@show m, (i, x) = valindargmax(sin, X)
m, (i, x)

# %%
F((x, y)) = f(x, y)

X = range(-2, 2; length=401)
Y = range(-2, 2; length=401)
XtimesY = Iterators.product(X, Y)
@show m, (i, x) = valindargmax(F, XtimesY)
m, (i, x)

# %%

max_stars_repo_path: 0008/new findmax(f, domain) spec and keys's methods for iterators.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module ScrollAreas

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template, register_normal_element

export scrollarea

register_normal_element("q__scroll__area", context = @__MODULE__)

function scrollarea(args...; wrap::Function = StippleUI.DEFAULT_WRAPPER, kwargs...)
  wrap() do
    q__scroll__area(args...; kwargs...)
  end
end

end


text:
abstract type AbstractBranchFormulation <: AbstractDeviceFormulation end

abstract type AbstractLineForm <: AbstractBranchFormulation end

abstract type AbstractDCLineForm <: AbstractBranchFormulation end

abstract type AbstractTransformerForm <: AbstractBranchFormulation end

abstract type PiLine <: AbstractLineForm end

abstract type SeriesLine <: AbstractLineForm end

abstract type SimpleHVDC <: AbstractDCLineForm end

function line_flow_limit(ps_m::CanonicalModel,
                         devices::Array{Br,1},
                         device_formulation::Type{D},
                         system_formulation::Type{S},
                         time_range::UnitRange{Int64}) where {Br <: PSY.MonitoredLine,
                                                               D <: AbstractBranchFormulation,
                                                               S <: PM.AbstractPowerFormulation}

    #rate_data = [(h.name, (min = -1*h.rate, max = h.rate) for h in devices]

    device_range(ps_m, range_data, time_range, :dc_rate_const, :Fbr)

    return

end

function line_rate_constraints(ps_m::CanonicalModel,
                               devices::Array{Br,1},
                               device_formulation::Type{D},
                                system_formulation::Type{PSI.StandardPTDFForm},
                          time_range::UnitRange{Int64}) where {Br <: PSY.Branch, D <: AbstractBranchFormulation}

    range_data = [(h.name, (min = -1*h.rate, max = h.rate)) for h in devices]

    device_range(ps_m, range_data, time_range, :line_rate_limit, :Fbr)

    return

end

function line_flow_limit(ps_m::CanonicalModel,
                         devices::Array{Br,1},
                         device_formulation::Type{D},
                         system_formulation::Type{S},
                         time_range::UnitRange{Int64}) where {Br <: PSY.MonitoredLine,
                                                               D <: AbstractBranchFormulation,
                                                               S <: PM.AbstractActivePowerFormulation}

    #rate_data = [(h.name, (min = -1*h.rate, max = h.rate) for h in devices]

    device_range(ps_m, range_data, time_range, :dc_rate_const, :Fbr)

    return

end

function flow_variables(ps_m::CanonicalModel,
                        system_formulation::Type{S},
                        devices::Array{B,1},
                        time_range::UnitRange{Int64}) where {B <: PSY.Branch,
                                                            S <: PM.DCPlosslessForm}

    add_variable(ps_m, devices, time_range, :Fbr, false)

end

function flow_variables(ps_m::CanonicalModel,
                        system_formulation::Type{S},
                        devices::Array{B,1},
                        time_range::UnitRange{Int64}) where {B <: PSY.Branch,
                                                             S <: PM.AbstractDCPLLForm}

    add_variable(ps_m, devices, time_range, :Fbr_to, false)
    add_variable(ps_m, devices, time_range, :Fbr_fr, false)

end

function flow_variables(ps_m::CanonicalModel,
                        system_formulation::Type{S},
                        devices::Array{B,1},
                        time_range::UnitRange{Int64}) where {B <: PSY.Branch,
                                                             S <:PM.AbstractPowerFormulation}

    add_variable(ps_m, devices, time_range, :PFbr_to, false)
    add_variable(ps_m, devices, time_range, :PFbr_fr, false)

    add_variable(ps_m, devices, time_range, :QFbr_to, false)
    add_variable(ps_m, devices, time_range, :QFbr_fr, false)

    return

end