text:
export LoadParamsHook, FetchParamsHook

#####
# StopMsg
#####

(msg::StopMsg)(x::StopSignal) = x[] = true

(msg::StopMsg)(x) = nothing

function (msg::StopMsg)(x::ComposedStopCondition)
    for s in x.stop_conditions
        msg(s)
    end
end

#####
# LoadParamMsg
#####

(msg::LoadParamMsg)(x) = nothing

function (msg::LoadParamMsg)(x::ComposedHook)
    for h in x.hooks
        msg(h)
    end
end

struct LoadParamsHook <: AbstractHook
    buffer::Channel
end

LoadParamsHook() = LoadParamsHook(Channel(DEFAULT_MAILBOX_SIZE))

function (hook::LoadParamsHook)(::PostActStage, agent, env)
    ps = nothing
    while isready(hook.buffer)
        ps = take!(hook.buffer).data
    end
    isnothing(ps) || Flux.loadparams!(agent.policy, ps)
end

(msg::LoadParamMsg)(x::LoadParamsHook) = put!(x.buffer, msg)

#####
# FetchParamsHook
#####

Base.@kwdef mutable struct FetchParamsHook <: AbstractHook
    target::RemoteChannel
    buffer::RemoteChannel = RemoteChannel(() -> Channel(DEFAULT_MAILBOX_SIZE))
    n::Int = 0
    freq::Int = 1
    is_blocking::Bool = false
end

(msg::LoadParamMsg)(x::FetchParamsHook) = put!(x.buffer, msg)

function (hook::FetchParamsHook)(::PostActStage, agent, env)
    hook.n += 1
    if hook.n % hook.freq == 0
        if isready(hook.buffer)
            # some other workers have sent the request since our last request
            # so we just reuse it instead of creating new requests to reduce
            # message traffic
            while isready(hook.buffer)
                ps = take!(hook.buffer).data
            end
            Flux.loadparams!(agent.policy, ps)
            put!(target, FetchParamMsg(hook.buffer))
        else
            put!(target, FetchParamMsg(hook.buffer))
            if hook.is_blocking
                ps = take!(hook.buffer).data
                Flux.loadparams!(agent.policy, ps)
            end
        end
    end
end


text:
using Crayons
using Crayons.Box
using InteractiveUtils
using ColorTypes
using ImageInTerminal
using ImageCore

const SPACE = " " # one white space
BG_COLORS=[
    #BLACK_BG         ,
    RED_BG           ,
    GREEN_BG         ,
    YELLOW_BG        ,
    BLUE_BG          ,
    MAGENTA_BG       ,
    CYAN_BG          ,
    LIGHT_GRAY_BG    ,
    DEFAULT_BG       ,
    DARK_GRAY_BG     ,
    LIGHT_RED_BG     ,
    LIGHT_GREEN_BG   ,
    LIGHT_YELLOW_BG  ,
    LIGHT_BLUE_BG    ,
    LIGHT_MAGENTA_BG ,
    LIGHT_CYAN_BG    ,
    WHITE_BG         ,
]

function goma(c,H,W; init::Bool=false)
    H,W = 5,10
    buf = IOBuffer()
    if !init
        print(buf, "\x1b[999D\x1b[$(H)A") #rollback H-times
    end
    for i in 1:H
        msg = repeat("$(rand(1:9))",W)
        println(buf, c, msg, Crayon(reset=true))
    end
    print(buf |> take! |> String)
end

function clean(H)
    buf = IOBuffer()
    for i in 1:H
        print(buf, "\x1b[2K") # clear line
        print(buf, "\x1b[999D\x1b[$(1)A") # rollback
    end
    print(buf |> take! |> String)
end

function v()
    H,W=5,10
    goma(BLACK_BG,H,W,init=true)
    for c in BG_COLORS
        goma(c,H,W)
        sleep(0.1)
    end
    clean(H)
end

function iv()
    inH,inW=240,320
    imgA =colorview(RGB, zeros(3, inH,inW))
    imgB =colorview(RGB, ones(3, inH,inW))
    outH,outW=displaysize()
    init = true
    for img in repeat([imgA,imgB],3)
        buf=IOBuffer()
        if !init
            print(buf, "\x1b[999D\x1b[$(outH)A")
        else
            init = false
        end
        ImageInTerminal.imshow(buf,img)
        println(buf |> take! |> String)
        sleep(0.5)
    end
    clean(outH)
end

v()

iv()

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# The RQL type hierarchy is as follows:
# Top
type ReqlTop
  # t =  na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlTop})
  Union{
    ReqlDatum,
      ReqlNull, Nullable,
      ReqlBool, Bool,
      ReqlNumber, Number,
      ReqlString, ByteString,
      ReqlObject, Associative,
        ReqlSingleSelection,
      ReqlArray, Array,
    ReqlSequence,
      ReqlArray, Array,
      ReqlStream,
        ReqlStreamSelection,
          ReqlTable,
    ReqlDatabase,
    ReqlFunction,ReqlFunction1,
    ReqlOrdering,
    ReqlPathSpec}
end

# Top/DATUM
type ReqlDatum
  # t =  na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlDatum})
  Union{
    ReqlDatum,
      ReqlNull, Nullable,
      ReqlBool, Bool,
      ReqlNumber, Number,
      ReqlString, ByteString,
      ReqlObject, Associative,
        ReqlSingleSelection,
      ReqlArray, Array}
end

# Top/DATUM/NULL
type ReqlNull
  # t = 1
  # value::Nullable
  op_code::Number
  value
end

function reql_type(a::Type{ReqlNull})
  Union{ReqlNull, Void}
end

# Top/DATUM/BOOL
type ReqlBool
  # t = 2
  # value::Bool
  op_code::Number
  value
end

function reql_type(a::Type{ReqlBool})
  Bool
end

# Top/DATUM/NUMBER
type ReqlNumber
  # t = 3
  # value::Number
  op_code::Number
  value
end

function reql_type(a::Type{ReqlNumber})
  Number
end

function reql_arr_type(a::Type{ReqlNumber})
  Union{Array{Int64},Array{Float64}}
end

function reql_arr_type(a)
  Array
end

# Top/DATUM/STRING
type ReqlString
  # t = 4
  # value::ByteString
  op_code::Number
  value
end

function reql_type(a::Type{ReqlString})
  ByteString
end

# Top/DATUM/OBJECT
type ReqlObject
  # t = 6
  # value::Associative
  op_code::Number
  value
end

function reql_type(a::Type{ReqlObject})
  Union{
    ReqlObject, Associative,
      ReqlSingleSelection}
end

# Top/DATUM/OBJECT/SingleSelection
type ReqlSingleSelection
  # t = na
  # value::Associative
  op_code::Number
  value
end

function reql_type(a::Type{ReqlSingleSelection})
  ReqlSingleSelection
end

# Top/DATUM/ARRAY
type ReqlArray
  # t = 5
  # value::Array
  op_code::Number
  value
end

function reql_type(a::Type{ReqlArray})
  Union{ReqlArray,Array}
end

# Top/Sequence
type ReqlSequence
  # t = na
  # value::Array
  op_code::Number
  value
end

function reql_type(a::Type{ReqlSequence})
  Union{
    ReqlSequence,
      ReqlArray, Array,
      ReqlStream,
        ReqlStreamSelection,
          ReqlTable}
end

# Top/Sequence/ARRAY
# type ReqlArray
#   t = 5
#   value::Array
#   value
# end

# Top/SequenceStream
type ReqlStream
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlStream})
  Union{
    ReqlStream,ReqlStreamSelection,ReqlTable}
end

# Top/SequenceStream/StreamSelection
type ReqlStreamSelection
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlStreamSelection})
  ReqlStreamSelection
end

# Top/SequenceStream/StreamSelection/Table
type ReqlTable
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlTable})
  ReqlTable
end

# Top/Database
type ReqlDatabase
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlDatabase})
  ReqlDatabase
end

# Top/Function
type ReqlFunction
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlFunction})
  ReqlFunction
end

# Top/Function
type ReqlFunction1
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlFunction1})
  ReqlFunction1
end

# Top/Ordering - used only by ORDER_BY
type ReqlOrdering
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlOrdering})
  ReqlOrdering
end

# Top/Pathspec -- an object, string, or array that specifies a path
type ReqlPathSpec
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlPathSpec})
  ReqlPathSpec
end

# Error
type ReqlError
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlError})
  ReqlError
end

type ReqlTime
  # t = na
  op_code::Number
  value
end

function reql_type(a::Type{ReqlTime})
  ReqlTime
end

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb


text:
# ranges
@test size(10:1:0) == (0,)
@test length(1:.2:2) == 6
@test length(1.:.2:2.) == 6
@test length(2:-.2:1) == 6
@test length(2.:-.2:1.) == 6
@test length(2:.2:1) == 0
@test length(2.:.2:1.) == 0

@test length(1:0) == 0
@test length(0.0:-0.5) == 0
@test length(1:2:0) == 0
L32 = linspace(int32(1), int32(4), 4)
L64 = linspace(int64(1), int64(4), 4)
@test L32[1] == 1 && L64[1] == 1
@test L32[2] == 2 && L64[2] == 2
@test L32[3] == 3 && L64[3] == 3
@test L32[4] == 4 && L64[4] == 4

r = [5:-1:1]
@test r[1]==5
@test r[2]==4
@test r[3]==3
@test r[4]==2
@test r[5]==1

@test length(.1:.1:.3) == 3
@test length(1.1:1.1:3.3) == 3
@test length(1.1:1.3:3) == 2
@test length(1:1:1.8) == 1

let
    span = 5:20
    r = -7:3:42
    @test findin(r, span) == 5:10
    r = 15:-2:-38
    @test findin(r, span) == 1:6
end
#@test isempty(findin(5+0*(1:6), 2:4))
#@test findin(5+0*(1:6), 2:5) == 1:6
#@test findin(5+0*(1:6), 2:7) == 1:6
#@test findin(5+0*(1:6), 5:7) == 1:6
#@test isempty(findin(5+0*(1:6), 6:7))
#@test findin(5+0*(1:6), 5:5) == 1:6

@test intersect(1:5, 2:3) == 2:3
@test intersect(-3:5, 2:8) == 2:5
@test intersect(-8:-3, -8:-3) == -8:-3
@test intersect(1:5, 5:13) == 5:5
@test isempty(intersect(-8:-3, -2:2))
@test isempty(intersect(-3:7, 2:1))
@test intersect(1:11, -2:3:15) == 1:3:10
@test intersect(1:11, -2:2:15) == 2:2:10
@test intersect(1:11, -2:1:15) == 1:11
@test intersect(1:11, 15:-1:-2) == 1:11
@test intersect(1:11, 15:-4:-2) == 3:4:11
@test intersect(-20:-5, -10:3:-2) == -10:3:-7
@test isempty(intersect(-5:5, -6:13:20))
@test isempty(intersect(1:11, 15:4:-2))
@test isempty(intersect(11:1, 15:-4:-2))
#@test intersect(-5:5, 1+0*(1:3)) == 1:1
#@test isempty(intersect(-5:5, 6+0*(1:3)))
@test intersect(-15:4:7, -10:-2) == -7:4:-3
@test intersect(13:-2:1, -2:8) == 7:-2:1
@test isempty(intersect(13:2:1, -2:8))
@test isempty(intersect(13:-2:1, 8:-2))
#@test intersect(5+0*(1:4), 2:8) == 5+0*(1:4)
#@test isempty(intersect(5+0*(1:4), -7:3))
@test intersect(0:3:24, 0:4:24) == 0:12:24
@test intersect(0:4:24, 0:3:24) == 0:12:24
@test intersect(0:3:24, 24:-4:0) == 0:12:24
@test intersect(24:-3:0, 0:4:24) == 24:-12:0
@test intersect(24:-3:0, 24:-4:0) == 24:-12:0
@test intersect(1:3:24, 0:4:24) == 4:12:16
@test intersect(0:6:24, 0:4:24) == 0:12:24
@test isempty(intersect(1:6:2400, 0:4:2400))
@test intersect(-51:5:100, -33:7:125) == -26:35:79
@test intersect(-51:5:100, -32:7:125) == -11:35:94
#@test intersect(0:6:24, 6+0*(0:4:24)) == 6:6:6
#@test intersect(12+0*(0:6:24), 0:4:24) == Range(12, 0, 5)
#@test isempty(intersect(6+0*(0:6:24), 0:4:24))
@test intersect(-10:3:24, -10:3:24) == -10:3:23
@test isempty(intersect(-11:3:24, -10:3:24))

@test !(3.5 in 1:5)
@test (3 in 1:5)
@test (3 in 5:-1:1)
#@test (3 in 3+0*(1:5))
#@test !(4 in 3+0*(1:5))

r = 0.0:0.01:1.0
@test (r[30] in r)
r = (-4*int64(maxintfloat(is(Int,Int32) ? Float32 : Float64))):5
@test (3 in r)
@test (3.0 in r)

@test !(1 in 1:0)
@test !(1.0 in 1.0:0.0)

# indexing range with empty range (#4309)
@test (3:6)[5:4] == 7:6
@test_throws BoundsError (3:6)[5:5]
@test_throws BoundsError (3:6)[5]
@test (0:2:10)[7:6] == 12:2:10
@test_throws BoundsError (0:2:10)[7:7]

# avoiding intermediate overflow (#5065)
@test length(1:4:typemax(Int)) == div(typemax(Int),4) + 1

# overflow in length
@test_throws OverflowError length(0:typemax(Int))
@test_throws OverflowError length(typemin(Int):typemax(Int))
@test_throws OverflowError length(-1:typemax(Int)-1)

let s = 0
    # loops ending at typemax(Int)
    for i = (typemax(Int)-1):typemax(Int)
        s += 1
        @test s <= 2
    end
    @test s == 2

    s = 0
    for i = (typemax(Int)-2):(typemax(Int)-1)
        s += 1
        @test s <= 2
    end
    @test s == 2

    s = 0
    for i = typemin(Int):(typemin(Int)+1)
        s += 1
        @test s <= 2
    end
    @test s == 2

    # loops covering the full range of smaller integer types
    s = 0
    for i = typemin(Uint8):typemax(Uint8)
        s += 1
    end
    @test s == 256

    # loops past typemax(Int)
    n = 0
    s = int128(0)
    for i = typemax(Uint64)-2:typemax(Uint64)
        n += 1
        s += i
    end
    @test n == 3
    @test s == 3*int128(typemax(Uint64)) - 3

    # loops over empty ranges
    s = 0
    for i = 0xff:0x00
        s += 1
    end
    @test s == 0

    s = 0
    for i = int128(typemax(Int128)):int128(typemin(Int128))
        s += 1
    end
    @test s == 0
end

# sums (see #5798)
if WORD_SIZE == 64
    @test sum(int128(1:10^18)) == div(10^18 * (int128(10^18)+1), 2)
    @test sum(int128(1:10^18-1)) == div(10^18 * (int128(10^18)-1), 2)
else
    @test sum(int64(1:10^9)) == div(10^9 * (int64(10^9)+1), 2)
    @test sum(int64(1:10^9-1)) == div(10^9 * (int64(10^9)-1), 2)
end

# operations with scalars
@test (1:3) - 2 == -1:1
@test (1:3) - 0.25 == 1-0.25:3-0.25
@test (1:3) + 2 == 3:5
@test (1:3) + 0.25 == 1+0.25:3+0.25
@test (1:2:6) + 1 == 2:2:6
@test (1:2:6) + 0.3 == 1+0.3:2:5+0.3
@test (1:2:6) - 1 == 0:2:4
@test (1:2:6) - 0.3 == 1-0.3:2:5-0.3

# operations between ranges and arrays
@test all(([1:5] + (5:-1:1)) .== 6)
@test all(((5:-1:1) + [1:5]) .== 6)
@test all(([1:5] - (1:5)) .== 0)
@test all(((1:5) - [1:5]) .== 0)

# tricky floating-point ranges

@test [0.1:0.1:0.3]   == [1:3]./10
@test [0.0:0.1:0.3]   == [0:3]./10
@test [0.3:-0.1:-0.1] == [3:-1:-1]./10
@test [0.1:-0.1:-0.3] == [1:-1:-3]./10
@test [0.0:0.1:1.0]   == [0:10]./10
@test [0.0:-0.1:1.0]  == []
@test [0.0:0.1:-1.0]  == []
@test [0.0:-0.1:-1.0] == [0:-1:-10]./10
@test [1.0:1/49:27.0] == [49:1323]./49
@test [0.0:0.7:2.1]   == [0:7:21]./10
@test [0.0:1.1:3.3]   == [0:11:33]./10
@test [0.1:1.1:3.4]   == [1:11:34]./10
@test [0.0:1.3:3.9]   == [0:13:39]./10
@test [0.1:1.3:4.0]   == [1:13:40]./10
@test [1.1:1.1:3.3]   == [11:11:33]./10
@test [0.3:0.1:1.1]   == [3:1:11]./10

@test [0.0:1.0:5.5]   == [0:10:55]./10
@test [0.0:-1.0:0.5]  == []
@test [0.0:1.0:0.5]   == [0.0]

@test [prevfloat(0.1):0.1:0.3] == [prevfloat(0.1), 0.2, 0.3]
@test [nextfloat(0.1):0.1:0.3] == [nextfloat(0.1), 0.2]
@test [prevfloat(0.0):0.1:0.3] == [prevfloat(0.0), 0.1, 0.2]
@test [nextfloat(0.0):0.1:0.3] == [nextfloat(0.0), 0.1, 0.2]
@test [0.1:0.1:prevfloat(0.3)] == [0.1, 0.2]
@test [0.1:0.1:nextfloat(0.3)] == [0.1, 0.2, nextfloat(0.3)]
@test [0.0:0.1:prevfloat(0.3)] == [0.0, 0.1, 0.2]
@test [0.0:0.1:nextfloat(0.3)] == [0.0, 0.1, 0.2, nextfloat(0.3)]
@test [0.1:prevfloat(0.1):0.3] == [0.1, 0.2, 0.3]
@test [0.1:nextfloat(0.1):0.3] == [0.1, 0.2]
@test [0.0:prevfloat(0.1):0.3] == [0.0, prevfloat(0.1), prevfloat(0.2), 0.3]
@test [0.0:nextfloat(0.1):0.3] == [0.0, nextfloat(0.1), nextfloat(0.2)]

for T = (Float32, Float64,),# BigFloat),
    a = -5:25, s = [-5:-1;1:25], d = 1:25, n = -1:15
    den   = convert(T,d)
    start = convert(T,a)/den
    step  = convert(T,s)/den
    stop  = convert(T,(a+(n-1)*s))/den
    @test [start:step:stop] == T[a:s:a+(n-1)*s]./den
end

# near-equal ranges
@test 0.0:0.1:1.0 != 0.0f0:0.1f0:1.0f0

# comparing and hashing ranges
let
    Rs = {1:2, int32(1:3:17), int64(1:3:17), 1:0, 17:-3:0,
          0.0:0.1:1.0, float32(0.0:0.1:1.0)}
    for r in Rs
        ar = collect(r)
        @test r != ar
        @test !isequal(r,ar)
        for s in Rs
            as = collect(s)
            @test !isequal(r,s) || hash(r)==hash(s)
            @test (r==s) == (ar==as)
        end
    end
end

# issue #2959
@test 1.0:1.5 == 1.0:1.0:1.5 == 1.0:1.0
#@test 1.0:(.3-.1)/.1 == 1.0:2.0

let r = typemin(Int64):2:typemax(Int64), s = typemax(Int64):-2:typemin(Int64)
    @test first(r) == typemin(Int64)
    @test last(r) == (typemax(Int64)-1)
    @test_throws OverflowError length(r)

    @test first(s) == typemax(Int64)
    @test last(s) == (typemin(Int64)+1)
    @test_throws OverflowError length(s)
end

@test length(typemin(Int64):3:typemax(Int64)) == 6148914691236517206
@test length(typemax(Int64):-3:typemin(Int64)) == 6148914691236517206

for s in 3:100
    @test length(typemin(Int):s:typemax(Int)) == length(big(typemin(Int)):big(s):big(typemax(Int)))
    @test length(typemax(Int):-s:typemin(Int)) == length(big(typemax(Int)):big(-s):big(typemin(Int)))
end

@test length(uint(1):uint(1):uint(0)) == 0
@test length(typemax(Uint):uint(1):(typemax(Uint)-1)) == 0
@test length(typemax(Uint):uint(2):(typemax(Uint)-1)) == 0
@test length((typemin(Int)+3):5:(typemin(Int)+1)) == 0

# issue #6364
@test length((1:64)*(pi/5)) == 64

max_stars_repo_name: Maple-Wang/julia

max_issues_repo_name: Maple-Wang/julia

max_forks_repo_name: Maple-Wang/julia


text:
"""model definition of capsule3d"""

const STEP_DIST_CAPSULE3D = 0.05

struct StateCapsule3D <: AbsState
    x::Float64
    y::Float64
    z::Float64
    ϕ::Float64  # roll
    θ::Float64  # pitch
    ψ::Float64  # yaw
end

function get_mid_status(p::StateCapsule3D, q::StateCapsule3D)::StateCapsule3D
    return StateCapsule3D(
        (p.x + q.x) / 2,
        (p.y + q.y) / 2,
        (p.z + q.z) / 2,
        diff_angles(q.ϕ, p.ϕ) / 2 + p.ϕ,
        diff_angles(q.θ, p.θ) / 2 + p.θ,
        diff_angles(q.ψ, p.ψ) / 2 + p.ψ,
    )
end

function dist(q1::StateCapsule3D, q2::StateCapsule3D)::Float64
    # normalize
    return norm([
        q1.x - q2.x,
        q1.y - q2.y,
        q1.z - q2.z,
        diff_angles(q1.ϕ, q2.ϕ) / π,
        diff_angles(q1.θ, q2.θ) / π,
        diff_angles(q1.ψ, q2.ψ) / π,
    ])
end

"""get rotation matrix for parameters roll, pitch, and yaw"""
function get_rotation_matrix(
    ϕ::Float64,  # roll
    θ::Float64,  # pitch
    ψ::Float64,   # yaw
)::Matrix{Float64}
    cϕ, sϕ = cos(ϕ), sin(ϕ)
    cθ, sθ = cos(θ), sin(θ)
    cψ, sψ = cos(ψ), sin(ψ)
    # Note: the last two columns are not use in Capsule3D model
    return [
        (cϕ*cθ) (-sϕ*cψ+cϕ*sθ*sψ) (sϕ*sψ+cϕ*sθ*cψ)
        (sϕ*cθ) (cϕ*cψ+sϕ*sθ*sψ) (-cϕ*sψ+sϕ*sθ*cψ)
        (-sθ) (cθ*sψ) (cθ*cψ)
    ]
end

function get_capsule3D_points(q::StateCapsule3D, axes::Float64)::Vector{Vector{Float64}}
    R = get_rotation_matrix(q.ϕ, q.θ, q.ψ)
    p1 = [q.x, q.y, q.z]
    p2 = R * [axes, 0, 0] + p1
    return [p1, p2]
end

function gen_connect(
    q::StateCapsule3D,  # to identify type
    obstacles::Vector{CircleObstacle3D},
    rads::Vector{Float64},
    axises::Vector{Float64};
    step_dist::Float64 = STEP_DIST_CAPSULE3D,
    max_dist::Union{Nothing,Float64} = nothing,
)::Function

    # check: q \in C_free
    f(q::StateCapsule3D, i::Int64)::Bool = begin
        a, b = get_capsule3D_points(q, axises[i])
        # outside
        any(x -> (x < rads[i] || 1 - rads[i] < x), vcat(a, b)) && return false
        # obstacles
        any(o -> dist(a, b, [o.x, o.y, o.z]) < o.r + rads[i], obstacles) && return false
        return true
    end

    f(q_from::StateCapsule3D, q_to::StateCapsule3D, i::Int64)::Bool = begin
        D = dist(q_from, q_to)
        !isnothing(max_dist) && D > max_dist && return false

        dϕ = diff_angles(q_to.ϕ, q_from.ϕ)
        dθ = diff_angles(q_to.θ, q_from.θ)
        dψ = diff_angles(q_to.ψ, q_from.ψ)

        for e in vcat(collect(0:step_dist:D) / D, 1.0)
            # intermediate: root
            a = (1 - e) * [q_from.x, q_from.y, q_from.z] + e * [q_to.x, q_to.y, q_to.z]
            # intermediate: rotation matrix
            R = get_rotation_matrix(q_from.ϕ + e * dϕ, q_from.θ + e * dθ, q_from.ψ + e * dψ)
            # intermediate: tip
            b = a + R * [axises[i], 0, 0]
            # outside
            any(x -> (x < rads[i] || 1 - rads[i] < x), vcat(a, b)) && return false
            # obstacles
            any(o -> dist(a, b, [o.x, o.y, o.z]) < o.r + rads[i], obstacles) &&
                return false
        end

        return true
    end

    return f
end

function gen_collide(
    q::StateCapsule3D,
    rads::Vector{Float64},
    axises::Vector{Float64};
    step_dist::Float64 = STEP_DIST_CAPSULE3D,
)::Function

    N = length(rads)

    f(
        q_i_from::StateCapsule3D,
        q_i_to::StateCapsule3D,
        q_j_from::StateCapsule3D,
        q_j_to::StateCapsule3D,
        i::Int64,
        j::Int64,
    ) = begin
        # check each pair of step
        D_i = dist(q_i_from, q_i_to)
        D_j = dist(q_j_from, q_j_to)

        dϕ_i = diff_angles(q_i_to.ϕ, q_i_from.ϕ)
        dθ_i = diff_angles(q_i_to.θ, q_i_from.θ)
        dψ_i = diff_angles(q_i_to.ψ, q_i_from.ψ)

        dϕ_j = diff_angles(q_j_to.ϕ, q_j_from.ϕ)
        dθ_j = diff_angles(q_j_to.θ, q_j_from.θ)
        dψ_j = diff_angles(q_j_to.ψ, q_j_from.ψ)

        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # root, rotation, tip
            a_i =
                (1 - e_i) * [q_i_from.x, q_i_from.y, q_i_from.z] +
                e_i * [q_i_to.x, q_i_to.y, q_i_to.z]
            R_i = get_rotation_matrix(
                q_i_from.ϕ + e_i * dϕ_i,
                q_i_from.θ + e_i * dθ_i,
                q_i_from.ψ + e_i * dψ_i,
            )
            b_i = a_i + R_i * [axises[i], 0, 0]

            for e_j in vcat(collect(0:step_dist:D_j) / D_j, 1.0)
                # root, rotation, tip
                a_j =
                    (1 - e_j) * [q_j_from.x, q_j_from.y, q_j_from.z] +
                    e_j * [q_j_to.x, q_j_to.y, q_j_to.z]
                R_j = get_rotation_matrix(
                    q_j_from.ϕ + e_j * dϕ_j,
                    q_j_from.θ + e_j * dθ_j,
                    q_j_from.ψ + e_j * dψ_j,
                )
                b_j = a_j + R_j * [axises[j], 0, 0]

                # check collision
                dist(a_i, b_i, a_j, b_j) < rads[i] + rads[j] && return true
            end
        end
        return false
    end

    f(q_i::StateCapsule3D, q_j::StateCapsule3D, i::Int64, j::Int64) = begin
        a_i, b_i = get_capsule3D_points(q_i, axises[i])
        a_j, b_j = get_capsule3D_points(q_j, axises[j])
        return dist(a_i, b_i, a_j, b_j) < rads[i] + rads[j]
    end

    f(Q_from::Vector{Node{StateCapsule3D}}, Q_to::Vector{Node{StateCapsule3D}}) = begin
        for i = 1:N, j = i+1:N
            f(Q_from[i].q, Q_to[i].q, Q_from[j].q, Q_to[j].q, i, j) && return true
        end
        return false
    end

    f(Q::Vector{Node{StateCapsule3D}}, q_i_to::StateCapsule3D, i::Int64) = begin
        q_i_from = Q[i].q
        D_i = dist(q_i_from, q_i_to)

        dϕ_i = diff_angles(q_i_to.ϕ, q_i_from.ϕ)
        dθ_i = diff_angles(q_i_to.θ, q_i_from.θ)
        dψ_i = diff_angles(q_i_to.ψ, q_i_from.ψ)

        positions = map(j -> get_capsule3D_points(Q[j].q, axises[j]), 1:N)
        for e_i in vcat(collect(0:step_dist:D_i) / D_i, 1.0)
            # root, rotation, tip
            a_i =
                (1 - e_i) * [q_i_from.x, q_i_from.y, q_i_from.z] +
                e_i * [q_i_to.x, q_i_to.y, q_i_to.z]
            R_i = get_rotation_matrix(
                q_i_from.ϕ + e_i * dϕ_i,
                q_i_from.θ + e_i * dθ_i,
                q_i_from.ψ + e_i * dψ_i,
            )
            b_i = a_i + R_i * [axises[i], 0, 0]

            for j in filter(k -> k != i, 1:N)
                # check collision
                dist(a_i, b_i, positions[j]...) < rads[i] + rads[j] && return true
            end
        end

        return false
    end

    return f
end

function gen_uniform_sampling(q::StateCapsule3D)::Function
    () -> begin
        return StateCapsule3D(rand(3)..., (rand(3) * 2π)...)
    end
end

function plot_motion!(
    q_from::StateCapsule3D,
    q_to::StateCapsule3D,
    rad::Float64,
    axes::Float64,
    params,
)
    plot!([q_from.x, q_to.x], [q_from.y, q_to.y], [q_from.z, q_to.z]; params...)
end

function plot_start_goal!(
    q_init::StateCapsule3D,
    q_goal::StateCapsule3D,
    rad::Float64,
    axes::Float64,
    params,
)
    plot!([q_init.x], [q_init.y], [q_init.z]; markershape = :hex, params...)
    plot!([q_goal.x], [q_goal.y], [q_goal.z]; markershape = :star, params...)
end

function plot_agent!(
    q::StateCapsule3D,
    rad::Float64,
    axes::Float64,
    color::String;
    step::Float64 = 0.2,
)
    R = get_rotation_matrix(q.ϕ, q.θ, q.ψ)
    a = [q.x, q.y, q.z]
    b = a + R * [axes, 0, 0]

    plotP(P) = begin
        plot!(
            map(e -> e[1], P),
            map(e -> e[2], P),
            map(e -> e[3], P),
            label = nothing,
            color = color,
            linewidth = 3.0,
        )
    end

    P1 = map(t -> a + R * [0, cos(t) * rad, sin(t) * rad], 0:step:2π)
    P2 = map(t -> b + R * [0, cos(t) * rad, sin(t) * rad], 0:step:2π)
    plotP(P1)
    plotP(P2)

    P1 = map(t -> a + R * [cos(t) * rad, sin(t) * rad, 0], π/2:step:3π/2)
    plotP(P1)
    P2 = map(t -> b + R * [cos(t) * rad, sin(t) * rad, 0], -π/2:step:π/2)
    plotP(P2)
    plotP([P1[1], P2[end]])
    plotP([P1[end], P2[1]])

    P1 = map(t -> a + R * [sin(t) * rad, 0, cos(t) * rad], π:step:2π)
    P2 = map(t -> b + R * [sin(t) * rad, 0, cos(t) * rad], 0:step:π)
    plotP(P1)
    plotP(P2)
    plotP([P1[1], P2[end]])
    plotP([P1[end], P2[1]])
end

function gen_random_instance_StateCapsule3D(;
    N_min::Int64 = 2,
    N_max::Int64 = 8,
    N::Int64 = rand(N_min:N_max),
    num_obs_min::Int64 = 0,
    num_obs_max::Int64 = 10,
    num_obs::Int64 = rand(num_obs_min:num_obs_max),
    rad::Float64 = 0.025,
    rad_obs::Float64 = 0.05,
    rad_min::Float64 = rad,
    rad_max::Float64 = rad,
    rad_obs_min::Float64 = rad_obs,
    rad_obs_max::Float64 = rad_obs,
    TIME_LIMIT::Float64 = 0.5,
)::Tuple{
    Vector{StateCapsule3D},  # initial configuration
    Vector{StateCapsule3D},  # goal configuration
    Vector{CircleObstacle3D},
    Vector{Float64},  # radius
    Vector{Float64},  # axises
}

    _q = StateCapsule3D(0, 0, 0, 0, 0, 0)
    while true
        t_s = now()
        timeover() = elapsed_sec(t_s) > TIME_LIMIT

        # generate obstacles
        obstacles = gen_obstacles(3, num_obs, rad_obs_min, rad_obs_max)

        # determine rads
        rads = map(e -> rand() * (rad_max - rad_min) + rad_min, 1:N)

        # determine axises
        axises = map(e -> rand() * (rad_max - rad_min) + rad_min, 1:N)

        # generate starts & goals
        connect = gen_connect(_q, obstacles, rads, axises)
        collide = gen_collide(_q, rads, axises)
        config_init, config_goal = gen_config_init_goal(_q, N, connect, collide, timeover)

        timeover() && continue
        return (config_init, config_goal, obstacles, rads, axises)
    end
end


text:
@static if Sys.isunix()
    @info("building release 965...")
    run(`make libpicosat.so`)
end

@static if Sys.iswindows()
    @error("PicoSAT.jl does not currently work on Windows")
end


text:
# wrapper to convert Polyhedra solver into LPQP solver

# To enable LPQP support from a Polyhedra solver, define, e.g.,
# LinearQuadraticModel(s::CDDSolver) = PolyhedraToLPQPBridge(PolyhedraModel(s))

mutable struct PolyhedraToLPQPBridge{T<:Real} <: MPB.AbstractLinearQuadraticModel
    m::AbstractPolyhedraModel
    A::SparseMatrixCSC{T,Int}
    collb::Vector{T}
    colub::Vector{T}
    obj::Vector{T}
    rowlb::Vector{T}
    rowub::Vector{T}
    sense::Symbol
    rowoffset
    coloffset
end

PolyhedraToLPQPBridge(s::AbstractPolyhedraModel) = PolyhedraToLPQPBridge(s, sparse(Int[],Int[],Float64[]), Float64[], Float64[], Float64[], Float64[], Float64[], :Min, nothing, nothing)

export PolyhedraToLPQPBridge

# Loads the provided problem data to set up the linear programming problem:
# min c'x
# st  lb <= Ax <= ub
#      l <=  x <= u
# where sense = :Min or :Max
function MPBSI.loadproblem!(wrap::PolyhedraToLPQPBridge, A, collb, colub, obj, rowlb, rowub, sense)
    wrap.A = A
    wrap.collb = collb
    wrap.colub = colub
    wrap.obj = obj
    wrap.rowlb = rowlb
    wrap.rowub = rowub
    wrap.sense = sense
end

function MPBSI.addquadconstr!(wrap::PolyhedraToLPQPBridge, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    error("For polyhedra solvers, no quadratic constraint is supported")
end

function computeoffsets(lp::LPHRepresentation)
    coloffset = Vector{Vector{Int}}(undef, size(lp.A, 2))
    rowoffset = Vector{Vector{Int}}(undef, size(lp.A, 1))
    # Assumes that hyperplanes are in first indices
    # and halfspaces in following indices
    curoffset = 0
    for i in 1:size(lp.A, 2)
        if i in lp.coleqs
            curoffset += 1
            coloffset[i] = [curoffset]
        end
    end
    for i in 1:size(lp.A, 1)
        if i in lp.roweqs
            curoffset += 1
            rowoffset[i] = [curoffset]
        end
    end
    for i in 1:size(lp.A, 2)
        if !(i in lp.coleqs)
            offsets = Int[]
            if i in lp.colleqs
                curoffset += 1
                push!(offsets, curoffset)
            end
            if i in lp.colgeqs
                curoffset += 1
                push!(offsets, -curoffset)
            end
            coloffset[i] = offsets
        end
    end
    for i in 1:size(lp.A, 1)
        if !(i in lp.roweqs)
            offsets = Int[]
            if i in lp.rowleqs
                curoffset += 1
                push!(offsets, curoffset)
            end
            if i in lp.rowgeqs
                curoffset += 1
                push!(offsets, -curoffset)
            end
            rowoffset[i] = offsets
        end
    end
    coloffset, rowoffset
end

function MPBSI.optimize!(wrap::PolyhedraToLPQPBridge)
    A = wrap.A
    collb = wrap.collb
    colub = wrap.colub
    obj = wrap.obj
    rowlb = wrap.rowlb
    rowub = wrap.rowub

    (nvar = length(collb)) == length(colub) || error("Unequal lengths for column bounds")
    (nrow = length(rowlb)) == length(rowub) || error("Unequal lengths for row bounds")

    lp = LPHRepresentation(A, collb, colub, rowlb, rowub)
    wrap.coloffset, wrap.rowoffset = computeoffsets(lp)
    MPBSI.loadproblem!(wrap.m, lp, obj, wrap.sense)
    MPBSI.optimize!(wrap.m)
end

MPBSI.getsolution(wrap::PolyhedraToLPQPBridge) = MPBSI.getsolution(wrap.m)
MPBSI.status(wrap::PolyhedraToLPQPBridge) = MPBSI.status(wrap.m)
MPBSI.getobjval(wrap::PolyhedraToLPQPBridge) = MPBSI.getobjval(wrap.m)

function MPBSI.getconstrsolution(wrap::PolyhedraToLPQPBridge)
    m = length(wrap.rowoffset)
    constrsol = MPBSI.getconstrsolution(wrap.m)
    map(offset -> sign(offset[1]) * constrsol[abs(offset[1])], wrap.rowoffset)
end

function getdualsaux(wrap::PolyhedraToLPQPBridge{T}, offsets) where T
    l = length(offsets)
    duals = zeros(T, l)
    constrduals = MPBSI.getconstrduals(wrap.m)
    for i in 1:l
        for offset in offsets[i]
            duals[i] = sign(offset) * constrduals[abs(offset)]
        end
    end
    duals
end

function MPBSI.getreducedcosts(wrap::PolyhedraToLPQPBridge)
    getdualsaux(wrap, wrap.coloffset)
end

function MPBSI.getconstrduals(wrap::PolyhedraToLPQPBridge)
    getdualsaux(wrap, wrap.rowoffset)
end

MPBSI.numconstr(wrap::PolyhedraToLPQPBridge) = size(wrap.A, 1)
MPBSI.numvar(wrap::PolyhedraToLPQPBridge) = size(wrap.A, 2)
MPBSI.getvarLB(wrap::PolyhedraToLPQPBridge) = wrap.collb
MPBSI.getvarUB(wrap::PolyhedraToLPQPBridge) = wrap.colub
MPBSI.getconstrLB(wrap::PolyhedraToLPQPBridge) = wrap.rowlb
MPBSI.getconstrUB(wrap::PolyhedraToLPQPBridge) = wrap.rowub
MPBSI.getobj(wrap::PolyhedraToLPQPBridge) = wrap.obj
MPBSI.getsense(wrap::PolyhedraToLPQPBridge) = wrap.sense
function MPBSI.setvarLB!(wrap::PolyhedraToLPQPBridge, l)
    wrap.collb = l
end
function MPBSI.setvarUB!(wrap::PolyhedraToLPQPBridge, u)
    wrap.colub = u
end
function MPBSI.setconstrLB!(wrap::PolyhedraToLPQPBridge, lb)
    wrap.rowlb = lb
end
function MPBSI.setconstrUB!(wrap::PolyhedraToLPQPBridge, ub)
    wrap.rowub = ub
end
function MPBSI.setobj!(wrap::PolyhedraToLPQPBridge, obj)
    wrap.obj = obj
end
function MPBSI.setsense!(wrap::PolyhedraToLPQPBridge, sense)
    wrap.sense = sense
end
function MPBSI.addvar!(wrap::PolyhedraToLPQPBridge, constridx::AbstractArray{T}, constrcoef, l, u, objcoef) where T<:Integer
    wrap.A = [wrap.A sparsevec(constridx, constrcoef, size(wrap.A, 1))]
    push!(wrap.collb, l)
    push!(wrap.colub, u)
    push!(wrap.obj, objcoef)
end
function MPBSI.addconstr!(wrap::PolyhedraToLPQPBridge, varidx::AbstractArray{T}, coef, lb, ub) where T<:Integer
    wrap.A = [wrap.A; sparsevec(varidx, coef, size(wrap.A, 2))']
    push!(wrap.rowlb, lb)
    push!(wrap.rowub, ub)
end


text:
export preallocate

"""
    PreallocatedMethod{F, Args <: Tuple, N, R}

Type for a preallocated method for method `F` with signature `Args`,
where `N` is number of threads, `R` is the record type.
"""
struct PreallocatedMethod{F, Args <: Tuple, N, R}
    f::F
    replay_ctxs::NTuple{N, R} # one replay context per thread

    function PreallocatedMethod{F, Args}(f::F, ctxs::NTuple{N, R}) where {F, Args, N, R}
        new{F, Args, N, R}(f, ctxs)
    end
end

function (f::PreallocatedMethod)(xs...)
    ctx = f.replay_ctxs[Threads.threadid()]
    reinitialize!(ctx.metadata)
    return Cassette.overdub(ctx, f.f, xs...)
end

"""
    preallocate(f, xs...)

Return the first run result and a preallocated version of the corresponding method of
callable object `f` given its argument. Note the behaviour of both the method generated record
and the method will be freezed, you shouldn't feed this function with a dynamic allocation
method.
"""
function preallocate(f::F, xs...) where F
    x, record = record_allocations(f, xs...)
    record = freeze(record)
    ctxs = ntuple(Threads.nthreads()) do k
        new_replay_ctx(k == 1 ? copy(record) : record)
    end
    return x, PreallocatedMethod{F, typeof(xs)}(f, ctxs)
end

function Base.show(io::IO, f::PreallocatedMethod{F, Args}) where {F, Args}
    print(io, "preallocate(", f.f)

    for each in Args.parameters
        print(io, ", ::", each)
    end

    print(io, ")")
end


text:
using FileManager
using Visualization
using Common

function view_octree(files, affineMatrix)
    mesh = []
    for file in files
        octree = FileManager.las2aabb(file)
        V,EV,FV = Common.getmodel(octree)
        PC = FileManager.las2pointcloud(file)
        push!(mesh,Visualization.points(Common.apply_matrix(affineMatrix,PC.coordinates), PC.rgbs))
        push!(mesh,Visualization.GLGrid(Common.apply_matrix(affineMatrix,V),EV,Visualization.COLORS[7]))
    end
    return mesh
end

potree = "D:/potreeDirectory/pointclouds/CASALETTO"
cloudMets = FileManager.CloudMetadata(potree)

vector_t = [cloudMets.boundingBox.x_min,cloudMets.boundingBox.y_min,cloudMets.boundingBox.z_min]

trie = FileManager.potree2trie(potree)
max_level = FileManager.max_depth(trie)
files = FileManager.get_files_in_potree_folder(potree,0, false)
leaves = FileManager.get_leaf(trie)
# PC = FileManager.las2pointcloud(files...)
affineMatrix = Common.t(-vector_t...)

mesh = view_octree(files,affineMatrix)
Visualization.VIEW(mesh)


text:
export FixedPatch

"""
    FixedPatch()

The `FixedPatch` is similar to the [`OrthogonalPatch`](@ref) but it doesn't change during
the tracking. Instead it only updates at the start of the tracking.
"""
struct FixedPatch <: AbstractAffinePatch end

struct FixedPatchState{T,N} <: AbstractAffinePatchState{N}
    v̄::PVector{T,N} # this is already conjugated
end

function state(::FixedPatch, x::PVector)
    v = copy(x)
    conj!(v.data)
    FixedPatchState(v)
end

function init!(state::FixedPatchState, x::AbstractVector)
    @boundscheck length(x) == length(state.v̄)
    LinearAlgebra.normalize!(x)
    @inbounds for i in eachindex(state.v̄)
        state.v̄[i] = conj(x[i])
    end
    state
end

onpatch!(x::AbstractVector, state::FixedPatchState) = onpatch!(x, state.v̄)
evaluate!(u, state::FixedPatchState, x::PVector) = evaluate_patch!(u, state.v̄, x)
jacobian!(U, state::FixedPatchState, x::PVector) = jacobian_patch!(U, state.v̄, x)


text:
@testset "tokenize" begin
    s = """
        A { B } C
        D } E { F
        """
    tokens = FP.find_tokens(s, FP.MD_TOKENS)
    names = [t.name for t in tokens]
    @test count(e -> e == :CU_BRACKET_OPEN, names) == 2
    @test count(e -> e == :CU_BRACKET_CLOSE, names) == 2
    @test count(e -> e == :LINE_RETURN, names) == 2
    @test names[end] == :EOS
    @test names[1] == :SOS
    @test length(tokens) == 2 + 2 + 2 + 1 + 1

    s = """
        A <!-- B --> C
        ---
        and +++
        """
    tokens = FP.find_tokens(s, FP.MD_TOKENS)
    names = [t.name for t in tokens]
    @test names == [
        :SOS,
        :COMMENT_OPEN, :COMMENT_CLOSE, :LINE_RETURN,
        :LINE_RETURN,  # hrule are processed as blocks
        :MD_DEF_BLOCK, :LINE_RETURN,
        :EOS
    ]
end

@testset "md-1" begin
    s = """
        { { } { } {{ }} }
        <!--
            hello
          bye
        -->"""
    tokens = FP.find_tokens(s, FP.MD_TOKENS)
    deleteat!(tokens, 1)

    check_tokens(tokens, [1, 2, 4, 6, 7], :CU_BRACKET_OPEN)
    check_tokens(tokens, [3, 5, 8, 9, 10], :CU_BRACKET_CLOSE)

    check_tokens(tokens, [11, 13, 14, 15], :LINE_RETURN)

    check_tokens(tokens, [12], :COMMENT_OPEN)
    check_tokens(tokens, [16], :COMMENT_CLOSE)

    check_tokens(tokens, [17], :EOS)

    @test length(tokens) == 17

    t = """
        abc
            def
        {{abc}}
        """ |> FP.default_md_tokenizer
    deleteat!(t, 1)
    @test t[1].name == t[2].name == :LINE_RETURN
end

@testset "md-base" begin
    s = """* *a _ b_ **a b**"""
    tokens = FP.find_tokens(s, FP.MD_TOKENS)
    deleteat!(tokens, 1)
    @test tokens[1].name == :EM_OPEN
    @test tokens[2].name == :EM_OPEN
    @test tokens[3].name == :EM_CLOSE
    @test tokens[4].name == :STRONG_OPEN
    @test tokens[5].name == :STRONG_CLOSE
    @test tokens[6].name == :EOS
    s = """--> ----"""
    tokens = FP.find_tokens(s, FP.MD_TOKENS)
    deleteat!(tokens, 1)
    @test tokens[1].name == :COMMENT_CLOSE
    @test tokens[2].name == :EOS
    s = """+++ +++
    """
    tokens = FP.find_tokens(s, FP.MD_TOKENS)
    deleteat!(tokens, 1)
    @test length(tokens) == 3 # 1 +++ 2 \n 3 EOS
    @test tokens[1].name == :MD_DEF_BLOCK
    tokens = """~~~ a""" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :RAW_HTML

    tokens = """:foobar: and: this: bar""" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test length(tokens) == 2
    @test tokens[1].name == :CAND_EMOJI

    # \\
    tokens = raw"\\, \ , \*" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test length(tokens) == 3 + 1
    @test tokens[1].name == :LINEBREAK
    @test tokens[2].name == :CHAR_92
    @test tokens[3].name == :CHAR_42
    tokens = raw"\_, \`, \@, \{, \}, \$, \#" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :CHAR_95
    @test tokens[2].name == :CHAR_96
    @test tokens[3].name == :CHAR_64
    @test tokens[4].name == :CHAR_123
    @test tokens[5].name == :CHAR_125
    @test tokens[6].name == :CHAR_36
    @test tokens[7].name == :CHAR_35
    tokens = raw"\[ \] \newenvironment{foo}{a}{b}" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :MATH_DISPL_B_OPEN
    @test tokens[2].name == :MATH_DISPL_B_CLOSE
    @test tokens[3].name == :LX_NEWENVIRONMENT
    @test tokens[4].name == :CU_BRACKET_OPEN
    tokens = raw"\newcommand{a}{b}\begin{a}\end{a}" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :LX_NEWCOMMAND
    @test tokens[2].name == :CU_BRACKET_OPEN
    @test tokens[4].name == :CU_BRACKET_OPEN
    @test tokens[6].name == :LX_BEGIN
    @test tokens[9].name == :LX_END
    tokens = raw"\foo\bar{a}" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :LX_COMMAND
    @test tokens[2].name == :LX_COMMAND
    @test tokens[1].ss == raw"\foo"
    @test tokens[2].ss == raw"\bar"
    @test tokens[4].name == :CU_BRACKET_CLOSE

    tokens = raw"@def @@d0 @@d1,d-2 @@" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :MD_DEF_OPEN
    @test tokens[2].name == :DIV_OPEN
    @test tokens[3].name == :DIV_OPEN
    @test tokens[3].ss == "@@d1,d-2"
    @test tokens[4].name == :DIV_CLOSE

    tokens = raw"& & foo" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :CHAR_HTML_ENTITY
    @test length(tokens) == 2

    tokens = raw"$ $$ $a$" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :MATH_INLINE
    @test tokens[2].name == :MATH_DISPL_A
    @test tokens[3].name == :MATH_INLINE
    @test tokens[4].name == :MATH_INLINE

    # tokens = raw"_$>_ _$<_ ___ **** ---" |> FP.default_md_tokenizer
    # deleteat!(tokens, 1)
    # @test tokens[1].name == :MATH_I_OPEN
    # @test tokens[2].name == :MATH_I_CLOSE
    # @test tokens[3].name == :EOS
end

@testset "MD-code" begin
    tokens = raw"` `` ``` ```` `````" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :CODE_SINGLE
    @test tokens[2].name == :CODE_DOUBLE
    @test tokens[3].name == :CODE_TRIPLE
    @test tokens[4].name == :CODE_QUAD
    @test tokens[5].name == :CODE_PENTA

    tokens = raw"``` ```! ```julia" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :CODE_TRIPLE
    @test tokens[2].name == :CODE_TRIPLE!
    @test tokens[3].name == :CODE_LANG3

    tokens = raw"````hello `````foo" |> FP.default_md_tokenizer
    deleteat!(tokens, 1)
    @test tokens[1].name == :CODE_LANG4
    @test tokens[2].name == :CODE_LANG5
end


text:
pattern_uppercase = r"[A-Z]+"
pattern_question = r"\?\s*$"
pattern_empty = r"^\s*$"

function bob(stimulus::AbstractString)
    if occursin(pattern_empty, stimulus)
        "Fine. Be that way!"
    elseif occursin(pattern_uppercase, stimulus) && uppercase(stimulus) == stimulus
        if occursin(pattern_question, stimulus)
            "Calm down, I know what I'm doing!"
        else
            "Whoa, chill out!"
        end
    elseif occursin(pattern_question, stimulus)
        "Sure."
    else
        "Whatever."
    end
end

max_stars_repo_name: cjtallman/exercism

max_issues_repo_name: cjtallman/exercism

max_forks_repo_name: cjtallman/exercism


text:
using AbstractPlotting, StatsMakie, StatsBase
using Test

using Random: seed!
using AbstractPlotting.GeometryTypes: HyperRectangle
using Distributions
using KernelDensity: kde

seed!(0)

@testset "boxplot" begin
    a = repeat(1:5, inner = 20)
    b = 1:100
    p = boxplot(a, b)
    plts = p[end].plots
    @test length(plts) == 4
    @test plts[1] isa Scatter
    @test isempty(plts[1][1][])

    # test categorical
    a = repeat(["a", "b", "c", "d", "e"], inner = 20)
    b = 1:100
    p = boxplot(a, b)
    plts = p[end].plots
    @test length(plts) == 4
    @test plts[1] isa Scatter
    @test isempty(plts[1][1][])

    @test plts[2] isa LineSegments
    pts = Point{2, Float32}[
        [1.0, 5.75], [1.0, 1.0], [0.6, 1.0], [1.4, 1.0], [1.0, 15.25],
        [1.0, 20.0], [1.4, 20.0], [0.6, 20.0], [2.0, 25.75], [2.0, 21.0],
        [1.6, 21.0], [2.4, 21.0], [2.0, 35.25], [2.0, 40.0], [2.4, 40.0],
        [1.6, 40.0], [3.0, 45.75], [3.0, 41.0], [2.6, 41.0], [3.4, 41.0],
        [3.0, 55.25], [3.0, 60.0], [3.4, 60.0], [2.6, 60.0], [4.0, 65.75],
        [4.0, 61.0], [3.6, 61.0], [4.4, 61.0], [4.0, 75.25], [4.0, 80.0],
        [4.4, 80.0], [3.6, 80.0], [5.0, 85.75], [5.0, 81.0], [4.6, 81.0],
        [5.4, 81.0], [5.0, 95.25], [5.0, 100.0], [5.4, 100.0], [4.6, 100.0]
    ]
    @test plts[2][1][] == pts

    @test plts[3] isa Poly

    poly = HyperRectangle{2,Float32}[
        HyperRectangle{2,Float32}(Float32[0.6, 5.75], Float32[0.8, 9.5]),
        HyperRectangle{2,Float32}(Float32[1.6, 25.75], Float32[0.8, 9.5]),
        HyperRectangle{2,Float32}(Float32[2.6, 45.75], Float32[0.8, 9.5]),
        HyperRectangle{2,Float32}(Float32[3.6, 65.75], Float32[0.8, 9.5]),
        HyperRectangle{2,Float32}(Float32[4.6, 85.75], Float32[0.8, 9.5]),
    ]

    @test plts[3][1][] == poly

    #notch
    p = boxplot(a, b, notch=true)
    plts = p[end].plots

    @test length(plts) == 4

    @test plts[2] isa LineSegments
    @test plts[2][1][] == pts

    @test plts[3] isa Poly

    notch_boxes = Vector{Point{2,Float32}}[map(Point2f0, [[0.6, 5.75], [1.4, 5.75], [1.4, 7.42817], [1.2, 10.5], [1.4, 13.5718], [1.4, 15.25], [0.6, 15.25], [0.6, 13.5718], [0.8, 10.5], [0.6, 7.42817], [0.6, 5.75]]),
                                           map(Point2f0, [[1.6, 25.75], [2.4, 25.75], [2.4, 27.4282], [2.2, 30.5], [2.4, 33.5718], [2.4, 35.25], [1.6, 35.25], [1.6, 33.5718], [1.8, 30.5], [1.6, 27.4282], [1.6, 25.75]]),
                                           map(Point2f0, [[2.6, 45.75], [3.4, 45.75], [3.4, 47.4282], [3.2, 50.5], [3.4, 53.5718], [3.4, 55.25], [2.6, 55.25], [2.6, 53.5718], [2.8, 50.5], [2.6, 47.4282], [2.6, 45.75]]),
                                           map(Point2f0, [[3.6, 65.75], [4.4, 65.75], [4.4, 67.4282], [4.2, 70.5], [4.4, 73.5718], [4.4, 75.25], [3.6, 75.25], [3.6, 73.5718], [3.8, 70.5], [3.6, 67.4282], [3.6, 65.75]]),
                                           map(Point2f0, [[4.6, 85.75], [5.4, 85.75], [5.4, 87.4282], [5.2, 90.5], [5.4, 93.5718], [5.4, 95.25], [4.6, 95.25], [4.6, 93.5718], [4.8, 90.5], [4.6, 87.4282], [4.6, 85.75]])]
    @test plts[3][1][] ≈ notch_boxes
end

@testset "density" begin
    v = randn(1000)
    d = density(v, bandwidth = 0.1)
    p1 = plot(d)
    p2 = lines(d.x, d.density)
    @test p1[end][1][] == p2[end][1][]
    p3 = plot(density(bandwidth = 0.1), v)
    @test p3[end] isa Lines
    @test p3[end][1][] == p1[end][1][]
    x = randn(1000)
    y = randn(1000)
    v = (x, y)
    d = density(v, bandwidth = (0.1, 0.1))
    p1 = heatmap(d)
    p2 = heatmap(d.x, d.y, d.density)
    @test p1[end][1][] == p2[end][1][]
    p3 = plot(density(bandwidth = (0.1, 0.1)), v)
    @test p3[end] isa Heatmap
    @test p3[end][1][] == p1[end][1][]
    p4 = surface(density(bandwidth = (0.1, 0.1)), v)
    @test p4[end] isa Surface
    @test p4[end][1][] == p1[end][1][]

    t = (x = x, y = y)
    p5 = surface(density(bandwidth = (0.1, 0.1)), Data(t), (:x, :y))
    plt = p5[end].plots[1]
    @test plt isa Surface
    @test plt[1][] == p1[end][1][]

    p6 = surface(density(bandwidth = (0.1, 0.1)), Data(t), [:x :y])
    plt = p6[end].plots[1]
    @test plt isa Surface
    @test plt[1][] == p1[end][1][]

    x = randn(1000)
    y = rand(1000)
    z = rand(1:3, 1000)
    t = (x=x, y=y, z=z)
    p7 = plot(density, Data(t), Group(:z), :x, (weights = :y,))

    m1 = findall(==(1), z)
    m2 = findall(==(2), z)
    m3 = findall(==(3), z)
    k1 = kde(x[m1], weights = y[m1])
    k2 = kde(x[m2], weights = y[m2])
    k3 = kde(x[m3], weights = y[m3])
    @test p7[end].plots[1][1][] ≈ Point2f0.(k1.x, k1.density)
    @test p7[end].plots[2][1][] ≈ Point2f0.(k2.x, k2.density)
    @test p7[end].plots[3][1][] ≈ Point2f0.(k3.x, k3.density)
end

@testset "distribution" begin
    d = Normal()
    rg = StatsMakie.support(d)
    @test minimum(rg) ≈ -3.7190164854556866
    @test maximum(rg) ≈ 3.719016485455714
    p = plot(d)
    plt = p[end]
    @test plt isa Lines
    @test !StatsMakie.isdiscrete(d)
    @test first(plt[1][][1]) ≈ minimum(rg) rtol = 1f-6
    @test first(plt[1][][end]) ≈ maximum(rg) rtol = 1f-6

    for (x, pd) in plt[1][]
        @test pd ≈ pdf(d, x) rtol = 1f-6
    end

    d = Poisson()
    rg = StatsMakie.support(d)
    @test rg == 0:6
    p = plot(d)
    @test p[end] isa ScatterLines
    plt = p[end].plots[1]
    @test StatsMakie.isdiscrete(d)

    @test first.(plt[1][]) == 0:6
    @test last.(plt[1][]) ≈ pdf.(d, first.(plt[1][]))
end

@testset "dodge" begin
    a = [1, 2, 3, 4]
    b = [1 11
         2 12
         3 13
         4 14]
    a_long = vcat(a, a)
    b_long = vec(b)
    c = [1, 1, 1, 1, 2, 2, 2, 2]

    p1 = barplot(Position.dodge, a, b)
    p2 = barplot(Position.dodge, b)
    p3 = barplot(Position.stack, a, b)
    p4 = barplot(Position.stack, Group(c), a_long, b_long)

    @test p1[end][1][] isa PlotList
    series = p1[end][1][][1]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] ≈ [0.8, 1.8, 2.8, 3.8]
    @test series[2] == [1, 2, 3, 4]
    @test series[:width] == 0.4

    series = p1[end][1][][2]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] ≈ [1.2, 2.2, 3.2, 4.2]
    @test series[2] == [11, 12, 13, 14]
    @test series[:width] == 0.4

    @test p2[end][1][] isa PlotList
    series = p2[end][1][][1]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] ≈ [0.8, 1.8, 2.8, 3.8]
    @test series[2] == [1, 2, 3, 4]
    @test series[:width] == 0.4

    series = p2[end][1][][2]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] ≈ [1.2, 2.2, 3.2, 4.2]
    @test series[2] == [11, 12, 13, 14]
    @test series[:width] == 0.4

    @test p3[end][1][] isa PlotList
    series = p3[end][1][][1]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] == [1, 2, 3, 4]
    @test series[2] == [11.0, 12.0, 13.0, 14.0]
    @test series[:fillto] == [12, 14, 16, 18]
    @test series[:width] == 0.8

    series = p3[end][1][][2]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] == [1, 2, 3, 4]
    @test series[2] == [0.0, 0.0, 0.0, 0.0]
    @test series[:fillto] == [11, 12, 13, 14]
    @test series[:width] == 0.8

    @test p4[end][1][] isa PlotList
    series = p4[end][1][][1]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] == [1, 2, 3, 4]
    @test series[2] == [11.0, 12.0, 13.0, 14.0]
    @test series[:fillto] == [12, 14, 16, 18]
    @test series[:width] == 0.8

    series = p4[end][1][][2]
    @test series isa PlotSpec
    @test AbstractPlotting.plottype(series) <: BarPlot
    @test series[1] == [1, 2, 3, 4]
    @test series[2] == [0.0, 0.0, 0.0, 0.0]
    @test series[:fillto] == [11, 12, 13, 14]
    @test series[:width] == 0.8
end

@testset "group" begin
    c = repeat(1:2, inner = 50)
    m = repeat(1:2, outer = 50)
    p = scatter(
        Group(
            color = c,
            marker = m,
        ),
        1:100,
        1:100,
        color = [:blue, :red],
        marker = [:cross, :circle]
    )
    @test length(p[end].plots) == 4
    @test p[end].plots[1].color[] == :blue
    @test p[end].plots[2].color[] == :blue
    @test p[end].plots[3].color[] == :red
    @test p[end].plots[4].color[] == :red
    @test p[end].plots[1].marker[] == :cross
    @test p[end].plots[2].marker[] == :circle
    @test p[end].plots[3].marker[] == :cross
    @test p[end].plots[4].marker[] == :circle

    @test p[end].plots[1][1][] == Point{2, Float32}.(1:2:49, 1:2:49)
    @test p[end].plots[2][1][] == Point{2, Float32}.(2:2:50, 2:2:50)
    @test p[end].plots[3][1][] == Point{2, Float32}.(51:2:99, 51:2:99)
    @test p[end].plots[4][1][] == Point{2, Float32}.(52:2:100, 52:2:100)

    t = (x = 1:100, y = 1:100, z = 2:2:200, m = m, c = c)
    q = scatter(
        Data(t),
        Group(color = :c, marker = :m),
        :x, :y,
        color = [:blue, :red],
        marker = [:cross, :circle]
    )

    @test length(q[end].plots) == 4
    @test q[end].plots[1].color[] == :blue
    @test q[end].plots[2].color[] == :blue
    @test q[end].plots[3].color[] == :red
    @test q[end].plots[4].color[] == :red
    @test q[end].plots[1].marker[] == :cross
    @test q[end].plots[2].marker[] == :circle
    @test q[end].plots[3].marker[] == :cross
    @test q[end].plots[4].marker[] == :circle

    @test q[end].plots[1][1][] == Point{2, Float32}.(1:2:49, 1:2:49)
    @test q[end].plots[2][1][] == Point{2, Float32}.(2:2:50, 2:2:50)
    @test q[end].plots[3][1][] == Point{2, Float32}.(51:2:99, 51:2:99)
    @test q[end].plots[4][1][] == Point{2, Float32}.(52:2:100, 52:2:100)

    colors = AbstractPlotting.wong_colors
    r = scatter(
        Data(t),
        Group(color = :c, marker = bycolumn),
        :x, (:y, :z), color = colors, marker = [:cross, :circle]
    )

    @test r[end].plots[1].color[] == colors[1]
    @test r[end].plots[2].color[] == colors[1]
    @test r[end].plots[3].color[] == colors[2]
    @test r[end].plots[4].color[] == colors[2]
    @test r[end].plots[1].marker[] == :cross
    @test r[end].plots[2].marker[] == :circle
    @test r[end].plots[3].marker[] == :cross
    @test r[end].plots[4].marker[] == :circle

    @test r[end].plots[1][1][] == Point{2, Float32}.(1:50, 1:50)
    @test r[end].plots[2][1][] == Point{2, Float32}.(1:50, 2:2:100)
    @test r[end].plots[3][1][] == Point{2, Float32}.(51:100, 51:100)
    @test r[end].plots[4][1][] == Point{2, Float32}.(51:100, 102:2:200)
end

@testset "histogram" begin
    v = randn(1000)
    h = fit(Histogram, v)
    p = plot(h)

    plt = p[end]
    @test plt isa BarPlot
    x = h.edges[1]
    @test plt[1][] ≈ Point{2, Float32}.(x[1:end-1] .+ step(x)/2, h.weights)

    v = (randn(1000), randn(1000))
    h = fit(Histogram, v, nbins = 30)
    p = plot(h)
    plt = p[end]
    @test plt isa Heatmap
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    p = surface(h)
    plt = p[end]
    @test plt isa Surface
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    p = surface(histogram(nbins = 30), v)
    plt = p[end]
    @test plt isa Surface
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    p = surface(histogram(nbins = 30), v...)
    plt = p[end]
    @test plt isa Surface
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    w = rand(1000)
    h = fit(Histogram, v, weights(w), nbins = 30)
    p = surface(histogram(nbins = 30), v..., (weights = w,))
    plt = p[end]
    @test plt isa Surface
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    v = (randn(1000), randn(1000), randn(1000))
    h = fit(Histogram, v)
    p = plot(h)
    plt = p[end]
    @test plt isa Volume
    x = h.edges[1]
    y = h.edges[2]
    z = h.edges[3]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ z[1:end-1] .+ step(z)/2
    @test plt[4][] == h.weights
end

@testset "frequency" begin
    v = rand(1:3, 1000)
    p = plot(frequency, v)
    n1 = count(==(1), v)
    n2 = count(==(2), v)
    n3 = count(==(3), v)
    @test p[end] isa BarPlot
    @test p[end][1][] == Point{2, Float32}.([1, 2, 3], [n1, n2, n3])
end

@testset "qqplot" begin
    v = randn(1000)
    q = qqbuild(fit(Normal, v), v)
    p = qqnorm(v)

    @test length(p[end].plots) == 2
    plt = p[end].plots[1]
    @test plt isa Scatter
    @test first.(plt[1][]) ≈ q.qx rtol = 1e-6
    @test last.(plt[1][]) ≈ q.qy rtol = 1e-6

    plt = p[end].plots[2]
    @test plt isa LineSegments
    @test first.(plt[1][]) ≈ [extrema(q.qx)...] rtol = 1e-6
    @test last.(plt[1][]) ≈ [extrema(q.qx)...] rtol = 1e-6

    p = qqnorm(v, qqline = nothing)
    @test length(p[end].plots) == 1
    plt = p[end].plots[1]
    @test plt isa Scatter
    @test first.(plt[1][]) ≈ q.qx rtol = 1e-6
    @test last.(plt[1][]) ≈ q.qy rtol = 1e-6

    p = qqnorm(v, qqline = :fit)
    plt = p[end].plots[2]
    itc, slp = hcat(fill!(similar(q.qx), 1), q.qx) \ q.qy
    xs = [extrema(q.qx)...]
    ys = slp .* xs .+ itc
    @test first.(plt[1][]) ≈ xs rtol = 1e-6
    @test last.(plt[1][]) ≈ ys rtol = 1e-6

    p = qqnorm(v, qqline = :quantile)
    plt = p[end].plots[2]
    xs = [extrema(q.qx)...]
    quantx, quanty = quantile(q.qx, [0.25, 0.75]), quantile(q.qy, [0.25, 0.75])
    slp = diff(quanty) ./ diff(quantx)
    ys = quanty .+ slp .* (xs .- quantx)
    @test first.(plt[1][]) ≈ xs rtol = 1e-6
    @test last.(plt[1][]) ≈ ys rtol = 1e-6

end

@testset "violin" begin
    x = repeat(1:4, 250)
    y = x .+ randn.()
    p = violin(x, y, side = :left, color = :blue)
    @test p[end] isa Violin
    @test p[end].plots[1] isa Mesh
    @test p[end].plots[1][:color][] == :blue
    @test p[end].plots[2] isa LineSegments
    @test p[end].plots[2][:color][] == :white
    @test p[end].plots[2][:visible][] == :false

    # test categorical
    x = repeat(["a", "b", "c", "d"], 250)
    p = violin(x, y, side = :left, color = :blue)
    @test p[end] isa Violin
    @test p[end].plots[1] isa Mesh
    @test p[end].plots[1][:color][] == :blue
    @test p[end].plots[2] isa LineSegments
    @test p[end].plots[2][:color][] == :white
    @test p[end].plots[2][:visible][] == :false
end

@testset "errorbar" begin
    x = [1:4;]
    y =  [1:4;]
    Δx = fill(0.25, 4)
    Δy = fill(0.25, 4)
    p = errorbar(x,y,Δx,Δy,xcolor=:green,ycolor=:red)
    @test p[end] isa ErrorBar
    @test p[end][:ycolor][] == :red
    @test p[end][:xcolor][] == :green

    @test p[end].plots[1] isa ErrorBarX
    @test p[end].plots[1].plots[1] isa LineSegments
    @test p[end].plots[1][:color][] == :green

    @test p[end].plots[2] isa ErrorBarY
    @test p[end].plots[2].plots[1] isa LineSegments
    @test p[end].plots[2][:color][] == :red
end


text:
# Types for quadrature algorithms.
"""
Abstract type for quadrature algorithms without user-defined nodes (e.g., Gaussian quadrature.)
"""
abstract type QuadratureAlgorithm end

"""
Abstract type for quadrature algorithms with user-defined nodes (e.g., trapezoidal integration).
"""
abstract type ExplicitQuadratureAlgorithm end

# Concrete types for quadrature algorithms.
"""
Gaussian quadrature. See specific methods for what precise algorithm is used (e.g., Gauss-Legendre, Gauss-Hermite, etc.)
"""
struct Gaussian <: QuadratureAlgorithm end # Distribution-family specific quadrature.

"""
A custom quadrature scheme written by Spencer Lyon as part of the QuantEcon.jl library. Used with permission.

For detailed information, see: https://github.com/QuantEcon/QuantEcon.jl/blob/be0a32ec17d1f5b04ed8f2e52604c70c69f416b2/src/quad.jl#L918.
"""
struct QuantileRange <: QuadratureAlgorithm end

"""
A dot product of a (finite) PDF vector and a finite set of transformed nodes.
"""
struct FiniteDiscrete <: ExplicitQuadratureAlgorithm end # Dot-product basically.

"""
Trapezoidal integration.
"""
struct Trapezoidal <: ExplicitQuadratureAlgorithm end # For iterable expectations.

# Abstract types for expectations.
"""
Abstract type for all expectations.
"""
abstract type Expectation end # Supports E(f)

# Concrete types for expectations.

#= For an example of using abstract types named in this way, see: https://github.com/JuliaStats/Distributions.jl/blob/2d98eb6f31e9a92cce416e7391a84cff9bba7292/src/truncate.jl#L1-L10. We define a family of Truncated{blahblahblah} types parametrically, but use the abstract Truncated as a supertype for all Truncated distributions.
=#
"""
Expectations which are paramterized by a vector of nodes (e.g., a discretized support) and corresponding quadrature weights.
"""
struct IterableExpectation{NT, WT} <: Expectation # Supports E(f), nodes, weights, *
    nodes::NT
    weights::WT
end


text:
include("BayesTest.jl")
import Gadfly.plot, Mamba.hpd

type BinomTest <: BayesTest
    prior::Distribution
    data::Dict{AbstractString, Real}
    BF::Float64 # BF_{01}
    posterior::Any # either distribution, or NaN
    params::AbstractString

    function BinomTest(y::Integer, N::Integer; a = 1, b = 1)
        self = new()
        self.params = "θ"
        self.data = Dict("y" => y, "N" => N)
        self.prior = Distributions.Beta(a, b)
        self.posterior = Distributions.Beta(a + y, b + N - y)
        self.BF = pdf(self.prior, .5) / pdf(self.posterior, .5) # Savage-Dickey Trick
        return self
    end
end

function update(Model::BinomTest, y, N)
    BinomTest(Model.data["y"] + y, Model.data["N"] + N;
              a = Model.posterior.α, b = Model.posterior.β)
end

function posterior(Model::BinomTest; iter = 10000)
    return rand(Model.posterior, iter)
end

function hpd(samples::Array; low::Real = .025, up::Real = .975)
    ss = sort(samples); n = length(ss)
    return ss[Int(n * low)], ss[Int(n * up)]
end

function plot(Model::BinomTest)
    fun(x) = Distributions.pdf(Model.posterior, x)
    sup = Distributions.support(Model.posterior)
    Gadfly.plot(fun, sup.lb, sup.ub)
end

function plot(samples::Array)
    density = KernelDensity.kde(samples)
    f(x) = KernelDensity.pdf(density, x)
    Gadfly.plot(f, minimum(dist), maximum(dist))
end


text:
module TestStats

using Test: @test, @test, @testset
using Ensemble
using Ensemble.Stats: logsumexp

function testlogsumexp_random()
    for i in 1:100
        x = randn()
        y = randn()

        @test isapprox(logsumexp(x,y), log(exp(x)+exp(y)))
    end
end

function testlogsumexp_inf()
    @test logsumexp(-Inf, 0.0)==0.0
    @test logsumexp(0.0, -Inf)==0.0
    @test logsumexp(-Inf,-Inf)==-Inf
end

function testall()
    testlogsumexp_random()
    testlogsumexp_inf()
end

end


text:
# tests based on AMLCR 2nd ed,
@testset "AMLCR Math" begin

    # §5.7 payment frequency and certain (Tables 5.1 and 5.2)
    function f(age)
        life = SingleLife(mortality = AMLCR_std, issue_age = age)
        int = Yields.Constant(0.05)
        lc = [
            AnnuityImmediate(life, int)
            AnnuityImmediate(life, int, frequency = 4)
            AnnuityDue(life, int)
            AnnuityDue(life, int, frequency = 4)
            AnnuityImmediate(life, int, 10)
            AnnuityImmediate(life, int, 10, frequency = 4)
            AnnuityDue(life, int, 10, frequency = 4)
        ]

        return present_value.(lc)
    end

    # atol only to two decimals b/c AMLCR uses rounded rates?
    @test all(isapprox.(f(20), [18.966, 19.338, 19.966, 19.588, 7.711, 7.855, 7.952], atol = 0.01))
    @test all(isapprox.(f(80), [7.548, 7.917, 8.548, 8.167, 6.128, 6.373, 6.539], atol = 0.01))

    life = SingleLife(mortality = AMLCR_std, issue_age = 20)
    int = Yields.Constant(0.05)

    # property tests

    @test present_value(AnnuityDue(life, int)) > present_value(AnnuityDue(life, int, start_time = 5))

    @test present_value(AnnuityDue(life, int, 5, certain = 5, start_time = 5)) == sum([1 for t = 5:9] .* [1 / 1.05^t for t = 5:9])

    # relation on pg 124
    @test issorted(present_value.([
        AnnuityImmediate(life, int)
        AnnuityImmediate(life, int, frequency = 4)
        AnnuityDue(life, int, frequency = 4)
        AnnuityDue(life, int)
    ]))

end#


text:
using Pkg
Pkg.activate(joinpath(@__DIR__,".."))
using PdbTool
using ArgParse
using HDF5
using LinearAlgebra: norm
using Statistics
using UnicodePlots
using JLD
using ArDCA

function zero_sum_gauge!(couplings)
    q = size(couplings, 2)
    @assert q == size(couplings, 3)
    for l in 1:size(couplings, 1)
        cmat = copy(couplings[l, :, :])
        cmat -= repeat(mean(cmat, dims=1), q, 1)
        cmat -= repeat(mean(cmat, dims=2), 1, q)
        cmat .+= mean(cmat)
        couplings[l, :, :] = cmat[:, :]
    end
end

function create_hmm(fasta_path)
    hmm_path = tempname()
    run(pipeline(fasta_path, `sed '/^>/!s/[a-z,\.]//g'`, `hmmbuild -n hmm --symfrac=0.0 --informat=a2m $hmm_path - `))
    return hmm_path
end

function get_predictions(model_path)

    couplings = h5read(model_path, "couplings")
    fields = h5read(model_path, "fields")

    zero_sum_gauge!(couplings)

    N = size(fields, 1)

    F = zeros(N,N)
    l = 1
    for i in 1:N
        for j in (i+1):N
            F[i,j] = norm(couplings[l ,1:end-1, 1:end-1])
            F[j,i] = F[i,j]
            l += 1
        end
    end
    @assert l == size(couplings,1)+1

    Fi = sum(F, dims=1)
    Fj = sum(F, dims=2)
    Fa = sum(F) * (1-1/N)
    F -= (Fj * Fi) / Fa

    preds = Vector{Tuple{Int64, Int64, Float64}}()
    for i in 1:N
        for j in (i+1):N
            push!(preds, (i,j,F[i,j]))
        end
    end

    print(heatmap(F))

    sort!(preds, by=k->k[3], rev=true)
    println(preds[1:10])

    return preds
end

function get_roc(fasta_path, model_path, pdb_path, chain)

    pdb = PdbTool.parsePdb(pdb_path)
    hmm_path = create_hmm(fasta_path)
    PdbTool.mapChainToHmm(pdb.chain[chain], hmm_path)
    if occursin("ardca", model_path) && !occursin("extracted", model_path)
        preds = load(model_path)["episcore"]
    else
        preds = get_predictions(model_path)
    end
    roc, intraAlignDist = PdbTool.makeIntraRoc(preds, pdb.chain[chain]), PdbTool.intraAlignDist(pdb.chain[chain])

    # replace pdb ids with align ids
    align_roc = []
    for i in 1:length(roc)
        a1 = pdb.chain[chain].residue[roc[i][1]].alignmentPos
        a2 = pdb.chain[chain].residue[roc[i][2]].alignmentPos
        r = roc[i][3]
        push!(align_roc, (a1, a2, r))
    end

    return align_roc, intraAlignDist

end

s = ArgParseSettings()

@add_arg_table s begin
    "--fasta_path"
    required = true
    "--model_path"
    required = true
    "--pdb_path"
    required = true
    "--chain"
    default="A"
    "--roc_out_file"
    required = true
    "--cmap_out_file"
end

parsed_args = parse_args(ARGS, s)

roc, intraAlignDist = get_roc(parsed_args["fasta_path"], parsed_args["model_path"], parsed_args["pdb_path"], parsed_args["chain"])

open(parsed_args["roc_out_file"], "w") do fid
    for i in 1:length(roc)
        println(fid, roc[i][1], ' ', roc[i][2], ' ', roc[i][3])
    end
end

if "cmap_out_file" in keys(parsed_args)
        open(parsed_args["cmap_out_file"], "w") do fid
        N = size(intraAlignDist, 1)
        for i in 1:N
            for j in (i+1):N
                println(fid, i, ' ', j, ' ', intraAlignDist[i,j])
            end
        end
    end
end

max_stars_repo_name: christophfeinauer/PairwiseDistillations

max_issues_repo_name: christophfeinauer/PairwiseDistillations

max_forks_repo_name: christophfeinauer/PairwiseDistillations


text:
using NearestNeighbors
using StaticArrays
using Random
using Plots
using LightGraphs
using SimpleWeightedGraphs
using Distances
using SparseArrays

# some naive samples
function sampler(f)
    function(n)
        res = Vector{SVector{3, Float64}}(undef, n)
        for i in eachindex(res)
            θ, φ = rand(2)
            res[i] = f(θ, φ)
        end
        res
    end
end

torus(R=3, r=1) = sampler((θ, φ) -> @SVector[(R+r*cospi(2θ))cospi(2φ),
                                             (R+r*cospi(2θ))sinpi(2φ),
                                             r*sinpi(2θ)])

sphere(r=1) = sampler((θ, φ) -> @SVector[r*sinpi(θ)cospi(2φ),
                                         r*sinpi(θ)sinpi(2φ),
                                         r*cospi(θ)])

points3d(X; markersize=0.5, kwargs...) =
    scatter3d(getindex.(X, 1), getindex.(X, 2), getindex.(X, 3);
              markersize=markersize, kwargs...)

points2d(X; markersize=1.0, kwargs...) =
    scatter(getindex.(X, 1), getindex.(X, 2);
            markersize=markersize, kwargs...)

points(X; kwargs...) =
    if length(eltype(X)) == 2
        points2d(X; kwargs...)
    else
        points3d(X; kwargs...)
    end

function graphplot(P, G)
    plt = points3d(P)
    E = SVector{3, Float64}[]

    for e in edges(G)
        s = src(e)
        d = dst(e)
        push!(E, P[s])
        push!(E, P[d])
        push!(E, @SVector[NaN,NaN,NaN])
    end

    plot3d!(plt, getindex.(E, 1), getindex.(E, 2), getindex.(E, 3))
end

function uniform_downsample(P, r, kdt = KDTree(X))
    n = length(P)
    kdt = KDTree(P)
    free = trues(n)
    idxs = shuffle(1:n)
    chosen = falses(n)

    for i in idxs
        free[i] || continue
        free[i] = false
        chosen[i] = true
        free[inrange(kdt, P[i], r)] .= false
    end
    (1:n)[chosen]
end

struct GeodesicComplex{T, D, P<:SVector{D, T}, M<:Metric, K<:NNTree{P, M}} <:
                       AbstractSimpleWeightedGraph{Int, T}
    points    ::Vector{P}
    landmarks ::Vector{Int}
    graph     ::SimpleWeightedGraph{Int, T}
    tree      ::K
    cover     ::Vector{Set{Int}}
    radius    ::T
    metric    ::M
end

function GeodesicComplex(pts::AbstractVector{SVector{D, T}}, r;
                         metric=Euclidean(), tree=KDTree) where {D, T}
    kdt = tree(pts, metric)
    landmarks = uniform_downsample(pts, r, kdt)
    cover = [Set(inrange(kdt, pts[landmarks[i]], r)) for i in eachindex(landmarks)]

    I = Int[]
    J = Int[]
    V = T[]
    for i in 1:length(landmarks), j in 1:(i - 1)
        if i ≠ j && !isempty(intersect(cover[i], cover[j]))
            d = evaluate(metric, pts[landmarks[i]], pts[landmarks[j]])
            append!(I, (i, j))
            append!(J, (j, i))
            append!(V, (d, d))
        end
    end
    graph = SimpleWeightedGraph(sparse(I, J, V))

    GeodesicComplex(pts, landmarks, graph, kdt, cover, r, metric)
end

function geodesic_rips(P::Vector{SVector{N, T}}, r) where {N, T}
    S = P[uniform_downsample(P, r)]
    kdt = KDTree(S)
    n = length(S)
    I = Int[]
    J = Int[]
    V = T[]

    for i in 1:n
        for j in inrange(kdt, S[i], 2r)
            v = evaluate(Euclidean(), S[i], S[j])
            append!(I, (i, j))
            append!(J, (j, i))
            append!(V, (v, v))
        end
    end

    S, SimpleWeightedGraph(sparse(I, J, V))
end

function geodesic_witness(P::Vector{SVector{N, T}}, r) where {N, T}
    kdt = KDTree(P)
    chosen = uniform_downsample(P, r)
    n = length(chosen)
    I = Int[]
    J = Int[]
    V = T[]

    for i in chosen
        for j in chosen
            if !isempty(intersect(inrange(kdt, P[i], r), inrange(kdt, P[j], r)))
                v = evaluate(Euclidean(), P[i], P[j])
                append!(I, (i, j))
                append!(J, (j, i))
                append!(V, (v, v))
            end
        end
    end

    P[chosen], SimpleWeightedGraph(sparse(I, J, V))
end

# generate some plots
function generate(n; r = 0.2, out="img", ext=".pdf")
    S2 = sphere()(n)
    T2 = torus()(n)
    G  = [@SVector[randn(), randn()] for _ in 1:n]

    for (X, name) in zip([S2, T2, G], ["S^2", "T^2", "gauss"])
        xs = getindex.(X, 1)
        ys = getindex.(X, 2)
        if length(eltype(X)) > 2
            zs = getindex.(X, 3)
        else
            zs = 0
        end
        xlim = extrema(xs)
        ylim = extrema(xs)
        zlim = extrema(zs)

        plt = points(X, #title="$name n=$(length(X))",
                     xlim=xlim, ylim=ylim, zlim=zlim,
                     aspect_ratio = 1,
                     legend = false)
        savefig(plt, joinpath(out, "$(name)_full$ext"))

        Y = uniform_downsample(X, r)
        plt = points(Y, #title="$name uniform r=$r n=$(length(Y))",
                     xlim=xlim, ylim=ylim, zlim=zlim,
                     aspect_ratio = 1,
                     legend = false)
        savefig(plt, joinpath(out, "$(name)_uniform$ext"))

        Y = weighted_downsample(X, r)
        plt = points(Y, #title="$name weighted r=$r n=$(length(Y))",
                     xlim=xlim, ylim=ylim, zlim=zlim,
                     aspect_ratio = 1,
                     legend = false)
        savefig(plt, joinpath(out, "$(name)_weighted$ext"))
    end
end

normalize.(@SVector rand(3) for _ in 1:1000) |> points

function datasets()

    T3 = torus()(1_000_000)
    g3 = GeodesicComplex(T3, 1.)

    T4 = vcat(T3[1:500_000], map(x -> x.+2.5, T3[500_0001:end]))
    g4w = GeodesicComplex(T4, .5)
    g4r = GeodesicComplex(T4, .5, witness = false)
end


text:
abstract type BatchPQs <: PoolQs end
abstract type ExtendingBatchQs <: BatchPQs end
abstract type MultiObjectiveBatchQs <: BatchPQs end

"""
    select_batch(batch_query_strategy, query_data, label_map, candidate_indices)
Applies batch_query_strategy to select the most useful unlabeled observations
from query_data.

pools is used to find out which observations in query_data are still unlabeled,
candidate_indices is used to find out which observations in query_data are valid
candidates for the query.

Returns list of indices of most useful observations.
"""
function select_batch end

function get_query_objects(qs::BatchPQs,
                           query_data::Array{T, 2} where T <: Real,
                           pools::Vector{Symbol},
                           global_indices::Vector{Int},
                           history::Vector{Vector{Int}})::Vector{Int}
    pool_map = MLLabelUtils.labelmap(pools)
    haskey(pool_map, :U) || throw(ArgumentError("No more points that are unlabeled."))
    all_history_values = collect(Iterators.flatten(history))
    candidate_indices = [i for i in pool_map[:U] if global_indices[i] ∉ all_history_values]
    debug(getlogger(@__MODULE__), "[QS] Selecting best batch of $(qs.k) from $(length(candidate_indices)) candidates.")
    local_query_indices = select_batch(qs, query_data, pool_map, candidate_indices)
    return global_indices[local_query_indices]
end

"""
rep_measure computes representativeness of observations.
Only KDE is supported so far.
"""
function get_rep_measure(name::Symbol)::Function
    if (name == :KDE)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, labels::Dict{Symbol, Vector{Int}}, batch::Vector{Int}) -> begin
            γ = MLKernels.getvalue(model.kernel_fct.alpha)
            return multi_kde(data, γ)(data[:, batch])
        end
    else
        return throw(ArgumentError("Invalid representativeness measure $(name) specified."))
    end
end

"""
div_measure computes diversity of all candidates incrementally
Incremental computation: only diversity of the candidates to the last sample that
has been added to the batch (new_in_batch_idx) must be computed,
old_scores saves aggregated result, i.e., the diversity to the rest of the batch

currently two measures are implemented:
:AngleDiversity - Batch diversity is minimal angle between two batch samples in kernel space
:EuclideanDistance - Batch diversity is minimal euclidean distance between two batch samples in data space
"""
function get_incremental_div_measure(name::Symbol)::Function
    if (name == :AngleDiversity)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, new_in_batch_idx::Int, candidates::Vector{Int}, old_scores::Vector{Float64}) -> begin
            if model.data == data
                # reuse model kernel matrix
                K = SVDD.is_K_adjusted(model) ? model.K_adjusted : model.K
            else
                # compute a new kernel matrix because query data differs from model data
                K = MLKernels.kernelmatrix(Val(:col), SVDD.get_kernel(model), data)
            end
            div_scores = [-abs(K[i, new_in_batch_idx]) / (sqrt(K[i, i]) * sqrt(K[new_in_batch_idx, new_in_batch_idx])) for i in candidates]
            if !isempty(old_scores)
                div_scores = min.(div_scores, old_scores)
            end
            return div_scores
        end
    elseif (name == :EuclideanDistance)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, new_in_batch_idx::Int, candidates::Vector{Int}, old_scores::Vector{Float64}) -> begin
            div_scores = Distances.colwise(Distances.Euclidean(), data[:, candidates], data[:, new_in_batch_idx])
            if !isempty(old_scores)
                div_scores = min.(div_scores, old_scores)
            end
            return div_scores
        end
    else
        return throw(ArgumentError("Invalid diversity measure $(name) specified."))
    end
end

function pairwise_angle_distance(kernel::MLKernels.Kernel, data::Array{T, 2}) where T <: Real
    K = MLKernels.kernelmatrix(Val(:col), kernel, data)
    return pairwise_angle_distance(K)
end

"""
Calculate the pairwise angle between all observations in the kernel space.
Kernel matrix is copied.
"""
function pairwise_angle_distance(K::Array{T, 2}) where T <: Real
    angle_dist = deepcopy(K)
    min_div = Inf
    for i in 1:size(K, 1)
        for j in i+1:size(K, 2)
            @inbounds angle_dist[i, j] = -abs(K[i, j]) / (sqrt(K[i, i]) * sqrt(K[j, j]))
        end
    end
    return angle_dist
end

function angle_batch_diversity(kernel::MLKernels.Kernel, data::Array{T, 2}) where T <: Real
    K = MLKernels.kernelmatrix(Val(:col), kernel, data)
    return angle_batch_diversity(K)
end

"""
The angle batch diversity is minimal angle between all pairs in the batch and
is calculated in the kernel space.
No kernel matrix allocation required.
"""
function angle_batch_diversity(K::Array{T, 2}) where T <: Real
    min_div = Inf
    for i in 1:size(K, 1)
        for j in i+1:size(K, 2)
            @inbounds div = -abs(K[i, j]) / (sqrt(K[i, i]) * sqrt(K[j, j]))
            min_div = min(div, min_div)
        end
    end
    return min_div
end

"""
Pairwise euclidean distance with diagonal set to Inf.
"""
function pairwise_euclidean_distance(data::Array{T, 2}) where T <: Real
    distances = Distances.pairwise(Distances.Euclidean(), data, dims=2)
    distances[LinearAlgebra.diagind(distances)] .= Inf
    return distances
end

"""
The euclidean batch diversity is the minimal euclidean distance between all pairs
in the batch and is calculated in the data space
"""
function euclidean_batch_diversity(data::Array{T, 2}) where T <: Real
    distances = pairwise_euclidean_distance(data)
    return minimum(distances)
end

"""
div_measure computes diversity of the samples in a batch
enumerative computation: compute diversity for the whole batch from scratch

currently two measures are implemented:
:AngleDiversity - Batch diversity is minimal angle between two batch samples in kernel space
:EuclideanDistance - Batch diversity is minimal euclidean distance between two batch samples in data space
"""
function get_div_measure(name::Symbol)::Function
    if (name == :AngleDiversity)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, batch::Vector{Int}) -> begin
            if model.data == data
                # reuse model kernel matrix
                K = SVDD.is_K_adjusted(model) ? model.K_adjusted[batch, batch] : model.K[batch, batch]
                return angle_batch_diversity(K)
            else
                # compute a new kernel matrix because query data differs from model data
                return angle_batch_diversity(SVDD.get_kernel(model), data[:, batch])
            end
        end
    elseif (name == :EuclideanDistance)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, batch::Vector{Int}) -> begin
            return euclidean_batch_diversity(data[:, batch])
        end
    else
        throw(ArgumentError("Invalid diversity measure $(name) specified."))
    end
end

function get_pairwise_distance(name::Symbol)::Function
    if (name == :AngleDiversity)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, candidates::Vector{Int}) -> begin
            if model.data == data
                # reuse model kernel matrix
                K = SVDD.is_K_adjusted(model) ? model.K_adjusted[candidates, candidates] : model.K[candidates, candidates]
                return pairwise_angle_distance(K)
            else
                # compute a new kernel matrix because query data differs from model data
                return pairwise_angle_distance(SVDD.get_kernel(model), data[:, candidates])
            end
        end
    elseif (name == :EuclideanDistance)
        return (model::SVDD.OCClassifier, data::Array{T, 2} where T <: Real, candidates::Vector{Int}) -> begin
            return pairwise_euclidean_distance(data[:, candidates])
        end
    else
        throw(ArgumentError("Invalid distance measure $(name) specified."))
    end
end

function min_max_normalization(x::Vector{Float64})::Vector{Float64}
    return (x .- minimum(x)) ./ (maximum(x) - minimum(x))
end


text:
using RecipesBase
using LaTeXStrings
using ProgressMeter
using PushVectors
include("TerminalCache.jl")

const SimpleIIHist = AbstractVector{Int}
const SimpleIIInfoState = Vector{Vector{Int}}

struct NullVec <: AbstractVector{Int} end
Base.size(::NullVec) = (0,)
Base.length(::NullVec) = 0

"""
Simple Imperfect Information Game
- Formulation of a matrix game as in imperfect information extensive form game
- Solved with CFR
- Takes reward matrix as input

    `SimpleIIGame(R::Matrix{NTuple{2,T}})`

## Example
```
RPS = SimpleIIGame([
    (0,0) (-1,1) (1,-1);
    (1,-1) (0,0) (-1,1);
    (-1,1) (1,-1) (0,0)
])
```
"""
struct SimpleIIGame{T}
    R::Matrix{NTuple{2,T}}
    terminals::Matrix{Vector{Int}}
    _terminal_cache::TerminalCache
    _util_cache::PushVector{Int, Vector{Int}}
end

function SimpleInfoState(g::SimpleIIGame, i::Int)
    if i === 1
        return [Int[]]
    else
        return [Int[i] for i in 1:size(g.R,1)]
    end
end

function SimpleIIGame(R::Matrix{NTuple{2,T}}) where T
    s = size(R)
    terminals = [[i,j] for i in 1:s[1], j in 1:s[2]]
    return SimpleIIGame(R,terminals, TerminalCache(terminals), PushVector{Int}(2))
end

"""
Simple Imperfect Information Game Player

Instantiate with `SimpleIIPlayer(game, id[, initial_strategy])`
Default to uniform strategy

## Example
```
init_strategy = [0.1,0.3,0.6]
player1 = SimpleIIPlayer(game, 1, init_strategy)
```
"""
struct SimpleIIPlayer{T}
    id::Int
    game::SimpleIIGame{T}
    strategy::Vector{Float64}
    hist::Vector{Vector{Float64}}
    regret_avg::Vector{Float64}
end

player(g::SimpleIIGame, id::Int) = SimpleIIPlayer(g,id)
player(g::SimpleIIGame, id::Int, s::Vector{Float64}) = SimpleIIPlayer(g, id, s)

function SimpleIIPlayer(game::SimpleIIGame, id::Int)
    n_actions = size(game.R, id)
    strategy = fill(1/n_actions, n_actions)
    SimpleIIPlayer(
        id,
        game,
        strategy,
        [copy(strategy)],
        zeros(n_actions),
    )
end

function SimpleIIPlayer(game::SimpleIIGame, id::Int, strategy::Vector{Float64})
    n_actions = size(game.R, id)
    @assert size(game.R, id) == length(strategy)
    SimpleIIPlayer(
        id,
        game,
        copy(strategy),
        [copy(strategy)],
        zeros(n_actions),
    )
end

function terminals(game::SimpleIIGame, h::AbstractVector{Int})
    return game._terminal_cache[h]
end

player(::SimpleIIGame, h::AbstractVector{Int}) = length(h) < 1 ? 1 : 2
player(h) = length(h) < 1 ? 1 : 2

path_prob(σ::NTuple{2,Vector{Float64}}, i::Int, h′::AbstractVector{Int}) = path_prob(σ,i,NullVec(), h′)
path_prob(σ::NTuple{2,Vector{Float64}}, h′::AbstractVector{Int}) = path_prob(σ,NullVec(), h′)

function path_prob(σ::NTuple{2,Vector{Float64}}, i::Int, h::AbstractVector{Int}, h′::AbstractVector{Int})
    if i < 0
        return neg_path_prob(σ, -i, h, h′)
    else
        return path_prob(σ[i], i, h, h′)
    end
end

function neg_path_prob(σ::NTuple{2,Vector{Float64}}, i::Int, h::AbstractVector{Int}, h′::AbstractVector{Int})
    if h′ == h
        return 1.0
    else
        act = last(h′)
        h′ = @view h′[1:end-1]
        player_turn = player(h′)
        prob = (player_turn != i) ? σ[player_turn][act] : 1.0
        return prob*neg_path_prob(σ, i, h, h′)
    end
end

function path_prob(σ_i::Vector{Float64}, i::Int, h::AbstractVector{Int}, h′::AbstractVector{Int})
    if h′ == h
        return 1.0
    else
        act = last(h′)
        h′ = @view h′[1:end-1]
        prob = (player(h′) == i) ? σ_i[act] : 1.0
        return prob*path_prob(σ_i, i, h, h′)
    end
end

function path_prob(σ::NTuple{2,Vector{Float64}}, h::AbstractVector{Int}, h′::AbstractVector{Int})
    if h′ == h
        return 1.0
    else
        act = last(h′)
        h′ = @view h′[1:end-1]
        prob = σ[player(h′)][act]
        return prob*path_prob(σ, h, h′)
    end
end

function path_prob(σ::NTuple{2,Vector{Float64}}, I::SimpleIIInfoState, i::Int)
    s = 0.0
    for h in I
        s += path_prob(σ, i, h)
    end
    return s
end

function u(game::SimpleIIGame, i::Int, h′::AbstractVector{Int})
    game.R[h′[1], h′[2]][i]
end

function u(game::SimpleIIGame, i::Int, I::SimpleIIInfoState, σ::NTuple{2,Vector{Float64}})
    num = 0.0
    den = 0.0
    for h in I
        pi_I = path_prob(σ, -i, h)
        den += pi_I
        for h′ in terminals(game,h)
            num += pi_I*path_prob(σ,h,h′)*u(game,i,h′)
        end
    end
    return num/den
end

function u(game::SimpleIIGame, i::Int, I::SimpleIIInfoState, a::Int, σ::NTuple{2,Vector{Float64}})
    num = 0.0
    den = 0.0
    for hk in I
        h = append!(game._util_cache,  hk)
        push!(h, a)
        pi_I = path_prob(σ, -i, h)
        den += pi_I
        for h′ in terminals(game,h)
            num += pi_I*path_prob(σ,h,h′)*u(game,i,h′)
        end
        empty!(h)
    end
    return num/den
end

function sub_regret(game::SimpleIIGame,i::Int,I::SimpleIIInfoState,σ::NTuple{2, Vector{Float64}},a::Int)
    path_prob(σ, I, -i)*(u(game,i,I,a,σ) - u(game,i, I, σ))
end

function update_strategy!(game::SimpleIIGame, I::SimpleIIInfoState, p1::SimpleIIPlayer, p2::SimpleIIPlayer; pushp2::Bool=true)
    σ = p1.strategy
    ra = p1.regret_avg
    T = length(p1.hist)
    norm_sum = 0.0
    for a in 1:length(σ)
        sr = sub_regret(game,1,I,(last(p1.hist),last(p2.hist)),a)
        ra[a] += (sr-ra[a])/T
        if ra[a] > 0.0
            norm_sum += ra[a]
            σ[a] = ra[a]
        else
            σ[a] = 0.0
        end
    end
    norm_sum === 0.0 ? fill!(σ,1/length(σ)) : σ ./= norm_sum

    σ1′ = copy(σ)
    push!(p1.hist, σ1′)
    if pushp2
        σ2′ = copy(p2.strategy)
        push!(p2.hist, σ2′)
    end
    return σ
end

function update_strategies!(game::SimpleIIGame, Is::NTuple{2,SimpleIIInfoState}, p1::SimpleIIPlayer, p2::SimpleIIPlayer)
    σ1 = p1.strategy
    ra1 = p1.regret_avg
    T1 = length(p1.hist)
    norm_sum = 0.0
    for a in 1:length(σ1)
        sr = sub_regret(game,1,Is[1],(last(p1.hist),last(p2.hist)),a)
        ra1[a] += (sr-ra1[a])/T1
        if ra1[a] > 0.0
            norm_sum += ra1[a]
            σ1[a] = ra1[a]
        else
            σ1[a] = 0.0
        end
    end

    norm_sum === 0.0 ? fill!(σ1,1/length(σ1)) : σ1 ./= norm_sum

    σ2 = p2.strategy
    ra2 = p2.regret_avg
    T2 = length(p2.hist)
    norm_sum = 0.0
    for a in 1:length(σ2)
        sr = sub_regret(game,2,Is[2],(last(p1.hist),last(p2.hist)),a)
        ra2[a] += (sr-ra2[a])/T2
        if ra2[a] > 0.0
            norm_sum += ra2[a]
            σ2[a] = ra2[a]
        else
            σ2[a] = 0.0
        end
    end
    norm_sum === 0.0 ? fill!(σ2,1/length(σ2)) : σ2 ./= norm_sum

    push!(p1.hist, copy(σ1))
    push!(p2.hist, copy(σ2))
    return σ1, σ2
end

function train_both!(p1::SimpleIIPlayer, p2::SimpleIIPlayer, N::Int; show_progress::Bool=false)
    game = p1.game
    I1 = SimpleInfoState(game, 1)
    I2 = SimpleInfoState(game, 2)
    L1 = length(p1.hist)
    L2 = length(p2.hist)

    prog = Progress(N; enabled=show_progress)
    for i in 1:N
        update_strategies!(game, (I1,I2), p1, p2)
        next!(prog)
    end

    finalize_strategy!(p1)
    finalize_strategy!(p2)
    return p1, p2
end

function train_one!(p1::SimpleIIPlayer, p2::SimpleIIPlayer, N::Int; show_progress::Bool=false)
    I = SimpleInfoState(p1.game, 1)

    prog = Progress(N; enabled=show_progress)
    for i in 1:N
        update_strategy!(p1.game, I, p1, p2)
        next!(prog)
    end

    finalize_strategy!(p1)
    return p1
end

function cumulative_strategies(p::SimpleIIPlayer)
    mat = Matrix{Float64}(undef, length(p.hist), length(p.strategy))
    σ = zeros(Float64, length(p.strategy))

    for (i,σ_i) in enumerate(p.hist)
        σ = σ + (σ_i - σ)/i
        mat[i,:] .= σ
    end
    return mat
end

function avg_strat(p::SimpleIIPlayer)
    return sum(p.hist)/length(p.hist)
end

function finalize_strategy!(p::SimpleIIPlayer) # `sum` causes gc
    σ = p.strategy .= 0.0
    for σ_i in p.hist
        σ .+= σ_i
    end
    σ ./= sum(σ)
end

function evaluate(p1::SimpleIIPlayer, p2::SimpleIIPlayer)
    # strategies assumed already finalized
    game = p1.game
    σ1 = p1.strategy
    σ2 = p2.strategy
    R = game.R
    s1,s2 = size(R)
    p1_eval = 0.0
    p2_eval = 0.0
    for i in 1:s1, j in 1:s2
        prob = σ1[i]*σ2[j]
        p1_eval += prob*R[i,j][1]
        p2_eval += prob*R[i,j][2]
    end
    return p1_eval, p2_eval
end

## extras

@recipe function f(p1::SimpleIIPlayer, p2::SimpleIIPlayer)
    layout --> 2
    link := :both
    framestyle := [:axes :axes]

    xlabel := "Training Steps"

    L1 = length(p1.strategy)
    labels1 = Matrix{String}(undef, 1, L1)
    for i in eachindex(labels1); labels1[i] = L"a_{%$(i)}"; end

    @series begin
        subplot := 1
        ylabel := "Strategy"
        title := "Player 1"
        labels := labels1
        cumulative_strategies(p1)
    end

    L2 = length(p2.strategy)
    labels2 = Matrix{String}(undef, 1, L2)
    for i in eachindex(labels2); labels2[i] = L"a_{%$(i)}"; end

    @series begin
        subplot := 2
        title := "Player 2"
        labels := labels2
        cumulative_strategies(p2)
    end
end

@recipe function f(p::SimpleIIPlayer)

    xlabel := "Training Steps"

    L = length(p.strategy)
    labels = Matrix{String}(undef, 1, L)
    for i in eachindex(labels); labels[i] = L"a_{%$(i)}"; end

    @series begin
        subplot := 1
        ylabel := "Strategy"
        title := "Player 1"
        labels := labels
        cumulative_strategies(p)
    end
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
using ModelParameters
using Base.Test

struct TestStruct{T<:Number}
    x::T
end

struct TestStruct2{T<:Number}
    x::T
    y::T
end

@testset "ModelParameters" begin
    @testset "Parameters" begin
        @test_throws ArgumentError Parameter(0.0, 1.0, -1.0)
        @test_throws ArgumentError Parameter(-2.0, -1.0, 1.0)
        p = Parameter(0.0)
        @test lower(p) == -Inf
        @test upper(p) == Inf
        p = Parameter(0)
        @test lower(p) == typemin(Int)
        @test upper(p) == typemax(Int)
        p = Parameter(-1.0, 1.0)
        @test initial(p) == 0.0
        @test isparameter(p)
        c = constant(1.0)
        @test !isparameter(c)
        @test isfixed(c)
        @test_throws ArgumentError push!(p, 3.0)
        @test_throws ArgumentError push!(p, -3.0)
        push!(p, 0.5)
        @test p == 0.5
        @test p ≈ 0.5
        reset!(p)
        @test 0 == p
        @test 0 ≈ p
        arr = [p, 1.0]
        @test isparameter.(arr) == [true, false]
        @test value.(arr) == [0.0, 1.0]
        @test getparameters(p) == [p]
        @test getparameters(c) == Parameter[]
        @test getparameters(arr) == [p]
        @test getparameters([Nullable(p), Nullable(c)]) == [p]
        @test arr ≈ arr

        p = one(p)
        @test norm(p) == value(p)
        @test imag(p) == 0im
        @test isequal(1.0, p)
        @test isequal(p, 1.0)
        @test isapprox(p, p)
        @test isapprox(1.0, p)
        @test zero(p) == 0.0

        p = Parameter(0.0)
        @test string(p) == "-Inf ≤ 0.0 ≤ Inf"
        @test string(c) == "1.0"

        @test c + c == 2
        @test c - c == 0
        @test c * c == 1
        @test c / c == 1

        @test c + 1 == 2
        @test c - 1 == 0
        @test c * 1 == 1
        @test c / 1 == 1
        @test 1 + c == 2
        @test 1 - c == 0
        @test 1 * c == 1
        @test 1 / c == 1

        a = [1.0, 2.0]
        @test c + a == [2.0, 3.0]
        @test c - a == [0.0, -1.0]
        @test c * a == [1.0, 2.0]
        @test a + c == [2.0, 3.0]
        @test a - c == [0.0, 1.0]
        @test a * c == [1.0, 2.0]
        @test a / c == [1.0, 2.0]
        @test isequal(p,p)
        @test !isequal(p,c)

        @test arr+3 == [3.0, 4.0]
        @test 3+arr == [3.0, 4.0]
        @test arr-3 == [-3.0, -2.0]
        @test 3-arr == [3.0, 2.0]
        @test arr*3 == [0.0, 3.0]
        @test 3*arr == [0.0, 3.0]
        @test arr/3 == [0.0, 1/3.0]
        @test p < c
        @test c > p
        @test -1 < c
        @test c < 3
        @test c > -1
        @test 3 > c

        @test [c, c] + [c, c] == [1, 1] + [1, 1]
        @test [1, 1] + [c, c] == [1, 1] + [1, 1]
        @test [c, c] + [1, 1] == [1, 1] + [1, 1]
        @test [c, c] - [c, c] == [1, 1] - [1, 1]
        @test [1, 1] - [c, c] == [1, 1] - [1, 1]
        @test [c, c] - [1, 1] == [1, 1] - [1, 1]

        @test [c, c] .* [c, c] == [1, 1] .* [1, 1]
        @test [1, 1] .* [c, c] == [1, 1] .* [1, 1]
        @test [c, c] .* [1, 1] == [1, 1] .* [1, 1]
        @test [c, c] ./ [c, c] == [1, 1] ./ [1, 1]
        @test [1, 1] ./ [c, c] == [1, 1] ./ [1, 1]
        @test [c, c] ./ [1, 1] == [1, 1] ./ [1, 1]

        p = Parameter(0.0)
        const fun = (
            :abs2, :acosh, :acot, :acotd, :acoth, :acsc, :acscd, :acsch, :asec,
            :asecd, :asinh, :atan, :atand, :cbrt, :cos, :cosd, :cosh, :cot, :cotd,
            :coth, :csc, :cscd, :csch, :exp, :exp2, :expm1, :gamma, :inv, :lgamma,
            :log, :log10, :log1p, :log2, :sec, :secd, :sech, :sin, :sind, :sinh,
            :sqrt, :tan, :tand, :tanh,
        )
        const fun2 = (:acos, :acosd, :asech, :asin, :asind, :atanh)
        val = pi/2
        push!(p, val)
        for f in fun
            @eval begin
                @test $f($p) == $f(float($val))
            end
        end
        val = 0.0
        push!(p, val)
        for f in fun2
            @eval begin
                @test $f($p) == $f(float($val))
            end
        end
    end
    @testset "Promotion" begin
        @test [Parameter(1.0), 1.0, 1] ==
            [Parameter(1.0), constant(1.0), constant(1.0)]
        @test promote(Parameter(1.0), 1.0, 1) ==
            (Parameter(1.0), constant(1.0), constant(1.0))
    end
    @testset "vary" begin
        exp = TestStruct(Parameter(0.5, 0.0, 1.0))
        act = @vary(x >= 0.0, x <= 1.0, x = 0.5, TestStruct(x))
        @test exp.x == act.x
        exp = TestStruct2(Parameter(0.5, 0.0, 1.0), Parameter(0.5, 0.0, 1.0))
        act = @vary(
            x >= 0.0,
            x <= 1.0,
            x = 0.5,
            y >= 0.0,
            y <= 1.0,
            y = 0.5,
            TestStruct2(x, y)
        )
        @test exp.x == act.x
        @test exp.y == act.y
        exp = TestStruct2(Parameter(0.5), Parameter(0.5))
        act = @vary(
            x = 0.5,
            y = 0.5,
            TestStruct2(x, y)
        )
        @test exp.x == act.x
        @test exp.y == act.y
    end
end


text:
@deprecate state_diff_jacobian!(args...) errstate_jacobian!(args...)
@deprecate ∇²differential!(args...) ∇errstate_jacobian!(args...)
@deprecate SampledTrajectory(X,U,dt::Vector{<:Real}) SampledTrajectory(X,U,dt=dt)
@deprecate SampledTrajectory(x::AbstractVector{<:Real}, u::AbstractVector{<:Real}, dt::AbstractVector{<:Real}) Traj([copy(x) for k = 1:length(dt)], [copy(u) for k = 1:length(dt)-1], dt=dt)
@deprecate get_data(x::SampledTrajectory) getdata(x)
Base.Base.@deprecate_binding Traj SampledTrajectory true

import Base.size
@deprecate size(fun::AbstractFunction) dims(fun)
@deprecate traj_size(Z) dims(Z)


text:
using Entropia
using Test

tol = 1e-3

@testset "bandt and pompe"  begin
    x = [4, 7, 9, 10, 6, 11, 3]

    @testset "test 00" begin
        e = bandt_and_pompe(x, 2)
        @test abs(e - 0.918) <= tol
    end

    @testset "test 01" begin
        e = bandt_and_pompe(x, 3)
        @test abs(e - 1.522) <= tol
    end
end

@testset "bandt and pompe normal" begin
    n = 2048
    t = [i/n for i in 0:n]

    @testset "test 00" begin
        x = sin.(2π * 15 * t)
        e = bandt_and_pompe_normal(x, 6)
        @test abs(e - 0.1579) <= tol
    end

    @testset "test 01" begin
        x = (0.5*sin.(2π*5*t).+1) .* (sin.(2π*50*t.^2+2π*10*t))
        e = bandt_and_pompe_normal(x, 6)
        @test abs(e - 0.2625) <= tol
    end
end

@testset "weight entropy" begin
    @testset "test 00" begin
        x = [4, 7, 9, 10, 6, 11, 3]
        e = weight_entropy(x)
        @test abs(e - 0.9801) <= tol
    end
end

@testset "complexity_entropy" begin
    @testset "test 00" begin
        x = [4, 7, 9, 10, 6, 11, 3]
        e, jsd, cjs = complexity_entropy(x)
        @test abs(e - 0.5888) <= tol
        @test abs(jsd - 0.2235) <= tol
        @test abs(cjs - 0.2900) <= tol
    end
end


text:
"""
"""
mutable struct StochasticGradDescent <: AbstractOptimizer

end

"""
"""
mutable struct ADAMOptimizer <: AbstractOptimizer

end


text:
using Weave
using Base.Test

weave("documents/winston_formats.txt", plotlib="Winston", doctype="tex")
result = readall(open("documents/winston_formats.tex"))
ref = readall(open("documents/winston_formats_ref.tex"))
@test result == ref

weave("documents/winston_formats.txt", plotlib="Winston", doctype="github")
result = readall(open("documents/winston_formats.md"))
ref = readall(open("documents/winston_formats_ref.md"))
@test result == ref

weave("documents/winston_formats.txt", plotlib="Winston", doctype="pandoc", fig_ext=".svg")
result = readall(open("documents/winston_formats.md"))
ref = readall(open("documents/winston_formats_svg_ref.md"))
@test result == ref

weave("documents/winston_formats.txt", plotlib="Winston", doctype="rst")
result = readall(open("documents/winston_formats.rst"))
ref = readall(open("documents/winston_formats_ref.rst"))
@test result == ref


text:
using ITensors
using ITensorUnicodePlots
using PastaQ
using Zygote

import ITensors: op

function op(
  ::OpName"CRy",
  ::SiteType"S=1/2"; θ)
  c = cos(θ/2)
  s = sin(θ/2)
  return [
  1 0 0 0
  0 1 0 0
  0 0 c -s
  0 0 s c]
end

i1 = Index(2, "S=1/2,i1")
i2 = Index(2, "S=1/2,i2")

Xp1 = state("Xp", i1)
Xm1 = state("Xm", i1)
Xp2 = state("Xp", i2)
Xm2 = state("Xm", i2)
Zp1 = state("Zp", i1)
Zm1 = state("Zm", i1)
Zp2 = state("Zp", i2)
Zm2 = state("Zm", i2)

Id1 = op("Id", i1)
Z1 = op("Z", i1)
X1 = op("X", i1)
Id2 = op("Id", i2)
Z2 = op("Z", i2)
X2 = op("X", i2)

i1 != i2

ZpZm = ITensor(i1, i2)
ZpZm[i1=>1, i2=>2] = 1

Zp1 = state("Zp", i1)
Zp2 = state("Zp", i2)
Zm1 = state("Zm", i1)
Zm2 = state("Zm", i2)

ZpZm = Zp1 * Zm2

ψ = ITensor(i1, i2)
ψ[i1=>1, i2=>2] = 1/√2
ψ[i1=>2, i2=>1] = 1/√2

ψ = (Zp1 * Zm2 + Zm1 * Zp2)/√2

inner(ZpZm, ψ) == 1/√2

_, S, _ = svd(ZpZm, i1);
diag(S) == [1, 0]

_, S, _ = svd(ψ, i1);
diag(S) == [1/√2, 1/√2]

H = ITensor(i1', i2',
i1, i2)
H[i1'=>2, i2'=>1,
i1=>2, i2=>1] = -1
# …

apply(X1, Zp1 * Zm2) ≈ Zm1 * Zm2

@visualize X1 * ZpZm

ZZ = Z1 * Z2
XI = X1 * Id2
IX = Id1 * X2
h = 0.5
H = -ZZ + h * (XI + IX)

inner(Zp1' * Zp2', H, Zp1 * Zp2)

D, _ = eigen(H);
diag(D) ≈ [-√2, -1, 1, √2]

CRy = op("CRy", i1, i2; θ=π/2)

apply(CRy, Zm1 * Zp2) ≈ Zm1 * Xp2

# Optimization
function minimize(f, ∂f, x; nsteps, γ)
  for n in 1:nsteps
    x = x - γ * ∂f(x)
    println("Step = $(n), f(x) = $(f(x))")
  end
  return x
end

# State optimization
function E(ψ)
  ψHψ = inner(ψ', H, ψ)
  ψψ = inner(ψ, ψ)
  return ψHψ / ψψ
end
∂E(ψ) = gradient(E, ψ)[1]

ψ0 = (Zp1 * Zp2 + Zm1 * Zm2)/√2
E(ψ0) == -1
norm(∂E(ψ0)) == 2

ψ = minimize(E, ∂E, ψ0; nsteps=10, γ=0.1)
E(ψ0) == -1
E(ψ) # ≈ -1.4142131 ≈ -√2

# Circuit optimization
U(θ,i1,i2) = [
  op("Ry",i1;θ=θ[1]),
  op("Ry",i2;θ=θ[2]),
  op("CX",i1,i2),
  op("Ry",i1;θ=θ[3]),
  op("Ry",i2;θ=θ[4])]

u(θ) = [
  ("Ry",1,(;θ=θ[1])),
  ("Ry",2,(;θ=θ[2])),
  ("CNOT",1,2),
  ("Ry",1,(;θ=θ[3])),
  ("Ry",2,(;θ=θ[4]))]

U(θ,i1,i2) = buildcircuit([i1,i2],u(θ))

ψ0 = Zp1 * Zp2
function E(θ)
  ψθ = apply(U(θ,i1,i2),ψ0)
  return inner(ψθ',H,ψθ)
end
∂E(ψ) = gradient(E, ψ)[1]

θ0 = [0, 0, 0, 0]
θ = minimize(E, ∂E, θ0; nsteps=40, γ=0.5)
E(θ0) == -1
E(θ) # ≈ -1.4142077 ≈ -√2

# Fidelity optimization
ψ0 = Zp1 * Zp2
ψ = (Zp1 * Zp2 + Zm1 * Zm2) / √2

function F(θ)
  ψθ = apply(U(θ,i1,i2),ψ0)
  return -abs2(inner(ψ,ψθ))
end
∂F(θ) = gradient(F, θ)[1]

θ0 = [0, 0, 0, 0]
θ = minimize(F, ∂F, θ0; nsteps=50, γ=0.1)
F(θ0) == -0.5
F(θ) # ≈ -0.9938992 ≈ -1