text:
# the plan is to NOT define a Blotter as an array but rather as a Timestamp with defined type T <: FinancialInstrument
typealias Blotter{T<:FinancialAsset} Timestamp{T}

# designed to show one blotter at a time, and lined up so an array of these objects look pretty
function show{T<:FinancialAsset}(io::IO, b::Blotter{T})

    # window size for column alignment

    # quantity max out at 1_000_000, ergo = 9 (7 + 2 for min spacing)
    q_len = 9
    # ShortFuture, ergo = 13 (11 + 2)
    u_len = 13
    # price max out at  9999.99, ergo = 7
    # this is also a prepend, so the suffix is fixed at 2
    p_len = 7
    # currency max out at AU$ = 5 (3 + 2)
    c_len = 5

    # unit variable
    unit = split(string(typeof(b.value)), ".")[2]  # asset type as a string

    # currency variable
    cvar = show_currency(b)

    if b.value.quantity >= 0
        print_with_color(:blue, io, string(b.timestamp))
        print_with_color(:blue, io, " | ")
        print_with_color(:green, io, " " * string(b.value.quantity) * ^(" ", q_len - strwidth(string(b.value.quantity))))
        print_with_color(:blue, io, unit * ^(" ", u_len - strwidth(unit)))
        print_with_color(:blue, io, ^(" ", p_len - strwidth(string(b.value.basis)))  * string(b.value.basis) * "  ") # rpad fixed b.value two spaces
        print_with_color(:blue, io, cvar * ^(" ", c_len - strwidth(cvar)))
        print_with_color(:blue, io, string(b.value.ticker))
    else
        print_with_color(:blue, io, string(b.timestamp))
        print_with_color(:blue, io, " | ")
        print_with_color(:red, io, string(b.value.quantity) * ^(" ", q_len - strwidth(string(b.value.quantity)) + 1))  # the +1 for negative value padding
        #print_with_color(:blue, io, unit * ^(" ", u_len - strwidth(unit) - 1))  # adj for negative value
        print_with_color(:blue, io, unit * ^(" ", u_len - strwidth(unit) - 0))  # adj for negative value
        #print_with_color(:blue, io, ^(" ", p_len - strwidth(string(b.value.basis)))  * string(b.value.basis) * "   ") # rpad fixed b.value three spaces
        print_with_color(:blue, io, ^(" ", p_len - strwidth(string(b.value.basis)))  * string(b.value.basis) * "  ") # rpad fixed b.value two spaces
        print_with_color(:blue, io, cvar * ^(" ", c_len - strwidth(cvar)))
        print_with_color(:blue, io, string(b.value.ticker))
    end
end

# helper method for show
function show_currency(b::Blotter)

    if b.value.currency == USD
        cvar = "\$"
    elseif b.value.currency == EUR
        cvar = "€"
    elseif b.value.currency == GBP
        cvar = "£"
    elseif b.value.currency == JPY
        cvar = "¥"
    elseif b.value.currency == AUD
        cvar = "AU\$"
    elseif b.value.currency == NZD
        cvar = "NZ\$"
    else
        cvar = ""
    end
    cvar
end
# Not sure why these getindex methods only work on type FinancialAsset, but it works so hey

# ticker
function getindex(b::Vector{Timestamp}, s::Symbol)
    bval = [blot.value.ticker == s for blot in b]
    b[bval]
end

# asset type single value
function getindex(b::Vector{Timestamp}, d::DataType)
    tval = [typeof(blot.value) == d for blot in b]
    b[tval]
end

# array of asset types
function getindex(b::Vector{Timestamp}, dv::Vector{DataType})
    counter = falses(length(b))
    for i in 1:length(b)
        if typeof(b[i].value) in dv
            counter[i] = true
        else
            nothing
        end
    end
    b[counter]
end

# currency
function getindex(b::Vector{Timestamp}, c::Currency)
    bval = [blot.value.currency == c for blot in b]
    b[bval]
end

# date

function getindex(b::Vector{Timestamp}, d::Union(Date, DateTime))
    ds = falses(length(b))
    for i in 1:length(b)
        if b[i].timestamp == d
            ds[i] = true
        end
    end
    b[ds]
end

# vector of dates
function getindex(b::Vector{Timestamp}, ds::Union(Vector{Date}, Vector{DateTime}))
    counter = falses(length(b))
    for i in 1:length(b)
        if b[i].timestamp in ds
            counter[i] = true
        else
            nothing
        end
    end
    b[counter]
end

# range of dates -> converts to vector and calls above method
function getindex(b::Vector{Timestamp}, r::Union(StepRange{Date}, StepRange{DateTime}))
    b[[r;]]
end


text:
struct Newton{IL, L} <: SecondOrderOptimizer
    alphaguess!::IL
    linesearch!::L
end

"""
# Newton
## Constructor
```julia
Newton(; alphaguess = LineSearches.InitialStatic(),
linesearch = LineSearches.HagerZhang())
```

## Description
The `Newton` method implements Newton's method for optimizing a function. We use
a special factorization from the package `PositiveFactorizations.jl` to ensure
that each search direction is a direction of descent. See Wright and Nocedal and
Wright (ch. 6, 1999) for a discussion of Newton's method in practice.

## References
 - Nocedal, J. and S. J. Wright (1999), Numerical optimization. Springer Science 35.67-68: 7.
"""
function Newton(; alphaguess = LineSearches.InitialStatic(), # Good default for Newton
                linesearch = LineSearches.HagerZhang())    # Good default for Newton
    Newton(alphaguess, linesearch)
end

Base.summary(::Newton) = "Newton's Method"

mutable struct NewtonState{Tx, T, F<:Cholesky} <: AbstractOptimizerState
    x::Tx
    x_previous::Tx
    f_x_previous::T
    F::F
    s::Tx
    @add_linesearch_fields()
end

function initial_state(method::Newton, options, d, initial_x)
    T = eltype(initial_x)
    n = length(initial_x)
    # Maintain current gradient in gr
    s = similar(initial_x)

    value_gradient!!(d, initial_x)
    hessian!!(d, initial_x)

    NewtonState(copy(initial_x), # Maintain current state in state.x
                similar(initial_x), # Maintain previous state in state.x_previous
                T(NaN), # Store previous f in state.f_x_previous
                Cholesky(similar(d.H, T, 0, 0), :U, BLAS.BlasInt(0)),
                similar(initial_x), # Maintain current search direction in state.s
                @initial_linesearch()...)
end

function update_state!(d, state::NewtonState, method::Newton)
    # Search direction is always the negative gradient divided by
    # a matrix encoding the absolute values of the curvatures
    # represented by H. It deviates from the usual "add a scaled
    # identity matrix" version of the modified Newton method. More
    # information can be found in the discussion at issue #153.
    T = eltype(state.x)

    if typeof(NLSolversBase.hessian(d)) <: AbstractSparseMatrix
        state.s .= -NLSolversBase.hessian(d)\convert(Vector{T}, gradient(d))
    else
        state.F = cholesky!(Positive, NLSolversBase.hessian(d))
        if typeof(gradient(d)) <: Array
            # is this actually StridedArray?
            ldiv!(state.s, state.F, -gradient(d))
        else
            # not Array, we can't do inplace ldiv
            gv = Vector{T}(undef, length(gradient(d)))
            copyto!(gv, -gradient(d))
            copyto!(state.s, state.F\gv)
        end
    end
    # Determine the distance of movement along the search line
    lssuccess = perform_linesearch!(state, method, d)

    # Update current position # x = x + alpha * s
    @. state.x = state.x + state.alpha * state.s
    lssuccess == false # break on linesearch error
end

function assess_convergence(state::NewtonState, d, options)
  default_convergence_assessment(state, d, options)
end

function trace!(tr, d, state, iteration, method::Newton, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(gradient(d))
        dt["h(x)"] = copy(NLSolversBase.hessian(d))
        dt["Current step size"] = state.alpha
    end
    g_norm = norm(gradient(d), Inf)
    update!(tr,
            iteration,
            value(d),
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

ext: jl


text:
#how to run
#include("Single_Robust.jl")
#generate_robust
joshpath = "C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy"

PROJECT_ROOT = joshpath
currentpath = joshpath

using AmbulanceDeployment
using DataFrames, Winston, JLD, CSV, Gurobi, JuMP, GLPK, Dates
using JSON

function generate_robust(alpha::Float64 = .1, namb::Int64 = 40, cross_validation::Int64 = 1)
   adjacent_nbhd = CSV.File("C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy/test/austin_data/adjacent_nbhd.csv") |> DataFrame
   coverage = CSV.File("C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy/test/austin_data/coverage.csv") |> DataFrame
   hourly_calls = CSV.File("C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy/test/austin_data/Full_WeekdayCalls.csv") |> DataFrame

   print("\n coverage ",size(coverage))
   print("\n hourly_calls ",size(hourly_calls))
   print("\n adjacent_nbhd ",size(adjacent_nbhd))

   num_regions = 209
   regions = collect(1:209)

   #reformat data frames as matrix data
   locations = collect(1:size(coverage,2))
   coverage = Matrix(coverage)
   coverage = convert(Array{Bool, 2}, coverage[:, :])
   adjacent = Matrix(adjacent_nbhd[!,2:end])
   adjacent = adjacent[regions,regions] .> 0.5
   demand = Matrix(hourly_calls[!,5:end-1]); #we were trimming first 6 collumns but it seems to me it only the first 5 contain extraneous data

   print("\n coverage ",size(coverage))
   print("\n demand ",size(demand))
   print("\n adjacent ",size(adjacent))

    # We focus on emergency calls during the "peak period" (8AM - 8PM),
    # with the emergency calls from the first 3 month as our training set,
    # and the subsequent emergency calls from the remaining months as our test set
    peak_period = (hourly_calls[!,:hour] .>= 8) .* (hourly_calls[!,:hour] .<= 20)
    indices = 1:DataFrames.nrow(hourly_calls);
    if(cross_validation == 1)
          train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 3)
    elseif(cross_validation == 2)
          train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 6) .* (hourly_calls[!,:month] .>= 3)
    elseif(cross_validation == 3)
          train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 9) .* (hourly_calls[!,:month] .>= 6)
    end
    test_filter  = .~train_filter;
    train_indices = indices[train_filter]
    test_indices = indices[test_filter];

    # we distinguish between peak and offpeak hours
    train_peak = indices[peak_period .* train_filter]
    train_offpeak = indices[.~peak_period .* train_filter]

    test_peak = indices[peak_period .* test_filter]
    test_offpeak = indices[.~peak_period .* test_filter]

    p = DeploymentProblem(30, length(locations), length(regions), demand, train_indices,
          test_indices, coverage[regions,:], Array{Bool,2}(adjacent));
    #p.demand[p.train,:]

    #test_inc_peak = inc_peak_period .* inc_test_filter
    #test_inc_offpeak = .~inc_peak_period .* inc_test_filter;

    #data structures to record results
    scenarios = Dict{Symbol, Dict{Int, Vector{Vector{Int}}}}()
    generated_deployment = Dict{Symbol, Dict{Int, Vector{Vector{Int}}}}()
    upperbounds = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    lowerbounds = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    upptiming = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    lowtiming = Dict{Symbol, Dict{Int, Vector{Float64}}}()
    amb_deployment = Dict{Symbol, Dict{Int, Vector{Int}}}()

    (deployment_model, name) = (next_dp -> RobustDeployment(next_dp, α=alpha), :Robust01)
    amb_deployment[name] = Dict{Int, Vector{Int}}()
    scenarios[name] = Dict{Int, Vector{Vector{Int}}}()
    generated_deployment[name] = Dict{Int, Vector{Vector{Int}}}()
    upperbounds[name] = Dict{Int, Vector{Float64}}()
    lowerbounds[name] = Dict{Int, Vector{Float64}}()
    upptiming[name] = Dict{Int, Vector{Float64}}()
    lowtiming[name] = Dict{Int, Vector{Float64}}()
    println("$namb ")
    p.nambulances = namb

    model = deployment_model(p)
    set_optimizer(model.m, Gurobi.Optimizer)
    @time AmbulanceDeployment.optimize_robust!(model, p);

    amb_deployment[name][namb] = deployment(model)

    # for tracking purposes
    scenarios[name][namb] = model.scenarios
    generated_deployment[name][namb] = model.deployment
    upperbounds[name][namb] = model.upperbounds
    lowerbounds[name][namb] = model.lowerbounds
    upptiming[name][namb] = model.upptiming
    lowtiming[name][namb] = model.lowtiming

    solver_stats = Dict{String, Any}()
    push!(solver_stats, "amb_deployment" => amb_deployment)
    push!(solver_stats, "scenarios" => scenarios)
    push!(solver_stats, "generated_deployment" => generated_deployment)
    push!(solver_stats, "upperbounds" => upperbounds)
    push!(solver_stats, "lowerbounds" => lowerbounds)
    push!(solver_stats, "upptiming" => upptiming)
    push!(solver_stats, "lowtiming" => lowtiming)
    json_string = JSON.json(solver_stats)

    open(string(PROJECT_ROOT , "/src/outputs/solver_stats_8_7.json"),"w") do f
     write(f, json_string)
    end

    print(sum(model.scenarios[1]))
    model
end

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
##### Execute this file in the REPL from current Spotify.jl branch directory #####

# Note: lcov is required for generating HTML output. You can install it via:
# Debian/Ubuntu: sudo apt install lcov

# Load package environment from "test" folder
using Pkg
Pkg.activate("test")
Pkg.instantiate()
using LocalCoverage

# Switch back to Spotify.jl package environment
Pkg.activate(pwd())
Pkg.instantiate()

# Helper function
error_message(e, stage::String) = @info "Unable to $(stage) coverage data \n Check this error: $(e)"

# Generate coverage summary and html data
try
    generate_coverage("Spotify", genhtml=true, show_summary=true, genxml=false)
catch e
    error_message(e, "generate")
end

# Open coverage summary in a browser
try
    open_coverage("Spotify")
catch e
    error_message(e, "open")
end

# Clean up .cov files but keep coverage directory
try
    clean_coverage("Spotify", rm_directory=false)
catch e
    error_message(e, "clean up")
end


text:
@doc raw"""
    LogKernel([α [,γ]])

The Log Kernel is a negative definite kernel given by:

```math
\kappa(\mathbf{x},\mathbf{y}) =
\log \left(1 + \alpha\|\mathbf{x} - \mathbf{y} \|^{2\gamma}\right)
\qquad \alpha > 0, \; \gamma \in (0,1]
```
"""
struct LogKernel{T<:AbstractFloat} <: NegativeDefiniteKernel{T}
    alpha::HyperParameter{T}
    gamma::HyperParameter{T}
    LogKernel{T}(α::Real, γ::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
    )
end
function LogKernel(α::T1 = 1.0, γ::T2 = one(T1)) where {T1<:Real,T2<:Real}
    LogKernel{floattype(T1,T2)}(α, γ)
end

@inline logkernel(z::T, α::T, γ::T) where {T<:AbstractFloat} = log(α*z^γ+1)

@inline basefunction(::LogKernel) = SquaredEuclidean()
@inline function kappa(κ::LogKernel{T}, z::T) where {T}
    logkernel(z, getvalue(κ.alpha), getvalue(κ.gamma))
end


text:
using BlockBandedMatrices, BandedMatrices, BlockArrays, ArrayLayouts, LinearAlgebra, Test
import BlockBandedMatrices: MemoryLayout, TriangularLayout, BandedBlockBandedColumnMajor,
                        BandedColumnMajor, BlockSkylineSizes, blockrowstop, blockcolstop, ColumnMajor
import BlockArrays: BlockedUnitRange, blockisequal

@testset "triangular" begin
    @testset "triangular BandedBlockBandedMatrix mul" begin
        A = BandedBlockBandedMatrix{Float64}(undef, 1:10,1:10, (1,1), (1,1))
        A.data .= randn.()

        U = UpperTriangular(A)
        @test MemoryLayout(typeof(U)) == TriangularLayout{'U','N',BandedBlockBandedColumnMajor}()
        b = randn(size(U,1))
        @test U*b isa Vector{Float64}
        @test all(lmul(U,b) .=== U*b)
        @test U*b  ≈ Matrix(U)*b

        U = UnitUpperTriangular(A)
        @test MemoryLayout(typeof(U)) == TriangularLayout{'U','U',BandedBlockBandedColumnMajor}()
        b = randn(size(U,1))
        @test U*b isa Vector{Float64}
        @test all(lmul(U,b) .=== U*b)
        @test U*b  ≈ Matrix(U)*b

        L = LowerTriangular(A)
        @test MemoryLayout(typeof(L)) == TriangularLayout{'L','N',BandedBlockBandedColumnMajor}()
        b = randn(size(U,1))
        @test L*b isa Vector{Float64}
        @test all(lmul(L,b) .=== L*b)
        @test L*b  ≈ Matrix(L)*b

        L = UnitLowerTriangular(A)
        @test MemoryLayout(typeof(L)) == TriangularLayout{'L','U',BandedBlockBandedColumnMajor}()
        b = randn(size(L,1))
        @test L*b isa Vector{Float64}
        @test all(lmul(L,b) .=== L*b)
        @test L*b  ≈ Matrix(L)*b
    end

    @testset "Block by BlockIndex" begin
        A = BandedBlockBandedMatrix{Float64}(undef, 1:10,1:10, (1,1), (1,1))
        A.data .= randn.()
        V = view(A, Block(3), Block.(3:5))
        @test MemoryLayout(typeof(V)) == BandedBlockBandedColumnMajor()
        b = randn(size(V,2))
        c = similar(b, size(V,1))

        @test blockbandwidths(V) == (1,1)
        @test (c .= MulAdd(V,b)) ≈ (Matrix(V)*b)

        V = view(A, Block.(2:4), Block(3))

        @test MemoryLayout(typeof(V)) == BandedBlockBandedColumnMajor()
        @test blockbandwidths(V) == (2,0)
        @test blockisequal(axes(V), (blockedrange(2:4), Base.OneTo(3)))

        @test blocksize(V) == (3,1)
        V2 = view(V, Block(1), Block(1))
        @test MemoryLayout(typeof(V2)) == BandedColumnMajor()

        b = rand(size(V,2))
        @test (similar(b, size(V,1)) .= MulAdd(V, b)) ≈ V*b
    end

    @testset "triangular BandedBlockBandedMatrix ldiv" begin
        A = BandedBlockBandedMatrix{Float64}(undef, 1:10,1:10, (1,1), (1,1))
        A.data .= randn.()
        A = A+10I

        U = UpperTriangular(A)
        b = randn(size(U,1))
        @test copyto!(similar(b), Ldiv(U,b)) ≈ Matrix(U) \ b
        @test all((similar(b) .= Ldiv(U, b)) .=== copyto!(similar(b), Ldiv(U,b)))

        U = UnitUpperTriangular(A)
        b = randn(size(U,1))
        @test copyto!(similar(b), Ldiv(U,b)) ≈ (Matrix(U) \ b)
        @test all((similar(b) .= Ldiv(U, b)) .=== copyto!(similar(b), Ldiv(U,b)))

        L = LowerTriangular(A)
        b = randn(size(L,1))
        @test copyto!(similar(b), Ldiv(L,b)) ≈ (Matrix(L) \ b)
        @test all((similar(b) .= Ldiv(L, b)) .=== copyto!(similar(b), Ldiv(L,b)))
    end

    @testset "Rectangular blocks BlockBandedMatrix linear algebra" begin
        l , u = 0,1
        rows = [3,4,5]
        cols = [2,3,4,3]

        A = BlockBandedMatrix{Float64}(undef, rows,cols, (l,u))
        A.data .= randn(length(A.data))

        b = randn(size(A,1))

        V = view(A, Block.(1:3), Block.(1:4))
        @test blockisequal(axes(V), axes(A))
        @test all(Matrix(V) .=== Matrix(A))
        @test view(V, Block(2)[1:2], Block(3)) ≡ view(A, Block(2)[1:2], Block(3))

        @test similar(Ldiv(UpperTriangular(A), b)) isa PseudoBlockVector
        @test UpperTriangular(A) \ b ≈ UpperTriangular(V) \ b ≈ UpperTriangular(Matrix(A)) \ b

        V = view(A, 1:11, 1:11)
        b = randn(size(V,1))
        @test UpperTriangular(V) \ b ≈ UpperTriangular(Matrix(V)) \ b ≈ UpperTriangular(A[1:11,1:11]) \ b

        # bug from SingularIntegralEquations, fixed by k_old check in _squaredblocks_newbandwidth
        A = BlockBandedMatrix{Float64}(undef, [5; fill(4,470)], [7; fill(4,368)], (102,269))
        A.data .= randn(length(A.data))
        for k = 1:min(size(A,1), size(A,2))
            A[k,k] += 10
        end

        V = view(A, 1:400,1:400)
        b = randn(400)
        @test UpperTriangular(V) \ b ≈ UpperTriangular(Matrix(V)) \ b
    end

    @testset "SubBlockSkylineMatrix linear algebra" begin
        l , u = 1,1
        N = M = 5
        cols = rows = 1:N
        A = BlockBandedMatrix{Float64}(undef, rows,cols, (l,u))
        A.data .= randn(length(A.data))

        V = view(A, Block.(1:3), Block.(1:3))

        @test blockrowstop(V,1) == Block(2)
        @test blockcolstop(V,1) == Block(2)

        @test blockisequal(axes(V), blockedrange.((1:3, 1:3)))

        b = randn(size(V,1))
        r = UpperTriangular(Matrix(V)) \ b
        @test_broken BlockBandedMatrices.blockbanded_squareblocks_trtrs!(V, copy(b)) ≈ r

        @test_broken all(ldiv!(UpperTriangular(V), copy(b)) .=== BlockBandedMatrices.blockbanded_squareblocks_trtrs!(V, copy(b)))

        V = view(A, Block.(2:3), Block(3))
        @test unsafe_load(pointer(V)) == A[2,4]
        @test unsafe_load(pointer(V)+sizeof(Float64)*stride(V,2)) == A[2,5]
        @test MemoryLayout(typeof(V)) == ColumnMajor()

        @test size(V) == (5,3)
        b = randn(size(V,2))

        @test all((similar(b,size(V,1)) .= MulAdd(V,b)) .=== Matrix(V)*b .===
                    BLAS.gemv!('N', 1.0, V, b, 0.0, Vector{Float64}(undef, size(V,1))))

        V = view(A, Block.(1:3), Block(3)[2:3])
        @test_throws ArgumentError pointer(V)
        @test MemoryLayout(typeof(V)) == ColumnMajor()

        b = randn(size(A,1))
        @test UpperTriangular(A) \ b ≈ UpperTriangular(Matrix(A)) \ b

        V = view(A, Block.(2:3), Block(3)[2:3])
        @test unsafe_load(pointer(V)) == A[2,5]
        @test unsafe_load(pointer(V)+sizeof(Float64)*stride(V,2)) == A[2,6]
        @test MemoryLayout(typeof(V)) == ColumnMajor()

        @test size(V) == (5,2)
        b = randn(size(V,2))
        @test all((similar(b,size(V,1)) .= MulAdd(V,b)) .=== Matrix(V)*b .===
                    BLAS.gemv!('N', 1.0, V, b, 0.0, Vector{Float64}(undef, size(V,1))))

        V = view(A, Block.(1:3), Block(3)[2:3])
        @test_throws ArgumentError pointer(V)
        @test MemoryLayout(typeof(V)) == ColumnMajor()

        b = randn(size(A,1))
        @test UpperTriangular(A) \ b ≈ UpperTriangular(Matrix(A)) \ b

        V_22 = view(A, Block(N)[1:N],  Block(N)[1:N])
        @test MemoryLayout(typeof(V_22)) == ColumnMajor()

        V = view(A, Block(N),  Block(N))
        V_22 = view(A, Block(N)[1:N],  Block(N)[1:N])
        @test unsafe_load(pointer(V_22)) == V_22[1,1] == V[1,1]
        @test strides(V_22) == strides(V) == (1,9)
        b = randn(N)
        @test all(copyto!(similar(b) , MulAdd(V,b)) .=== copyto!(similar(b) , MulAdd(V_22,b)) .===
            Matrix(V)*b .===
            BLAS.gemv!('N', 1.0, V, b, 0.0, Vector{Float64}(undef, size(V,1))))

        @test_skip all(UpperTriangular(V_22) \ b .=== ldiv!(UpperTriangular(V_22) , copy(b)) .=== ldiv!(UpperTriangular(V) , copy(b)) .===
            ldiv!(UpperTriangular(Matrix(V_22)) , copy(b)))

        V = view(A, Block.(rows), Block.(cols))
        V2 = view(A, 1:size(A,1), 1:size(A,2))
        b = randn(size(V,1))

        @test all(Matrix(V) .=== Matrix(V2))
        @test UpperTriangular(V2) \ b ≈ UpperTriangular(V) \ b
    end
end


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

@doc raw"""FSGroupStrategyOptions defines the strategy type and options used to create the strategy. Deprecated: use FSGroupStrategyOptions from policy API Group instead.

    IoK8sApiExtensionsV1beta1FSGroupStrategyOptions(;
        ranges=nothing,
        rule=nothing,
    )

    - ranges::Vector{IoK8sApiExtensionsV1beta1IDRange} : ranges are the allowed ranges of fs groups.  If you would like to force a single fs group then supply a single range with the same start and end. Required for MustRunAs.
    - rule::String : rule is the strategy that will dictate what FSGroup is used in the SecurityContext.
"""
mutable struct IoK8sApiExtensionsV1beta1FSGroupStrategyOptions <: SwaggerModel
    ranges::Any # spec type: Union{ Nothing, Vector{IoK8sApiExtensionsV1beta1IDRange} } # spec name: ranges
    rule::Any # spec type: Union{ Nothing, String } # spec name: rule

    function IoK8sApiExtensionsV1beta1FSGroupStrategyOptions(;ranges=nothing, rule=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1FSGroupStrategyOptions, Symbol("ranges"), ranges)
        setfield!(o, Symbol("ranges"), ranges)
        validate_property(IoK8sApiExtensionsV1beta1FSGroupStrategyOptions, Symbol("rule"), rule)
        setfield!(o, Symbol("rule"), rule)
        o
    end
end # type IoK8sApiExtensionsV1beta1FSGroupStrategyOptions

const _property_map_IoK8sApiExtensionsV1beta1FSGroupStrategyOptions = Dict{Symbol,Symbol}(Symbol("ranges")=>Symbol("ranges"), Symbol("rule")=>Symbol("rule"))
const _property_types_IoK8sApiExtensionsV1beta1FSGroupStrategyOptions = Dict{Symbol,String}(Symbol("ranges")=>"Vector{IoK8sApiExtensionsV1beta1IDRange}", Symbol("rule")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1FSGroupStrategyOptions }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1FSGroupStrategyOptions))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1FSGroupStrategyOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1FSGroupStrategyOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1FSGroupStrategyOptions }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1FSGroupStrategyOptions[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1FSGroupStrategyOptions)
    true
end

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


text:
import StateSpaceReconstruction:
	cembed,
	assign_bin_labels

import PerronFrobenius:
	get_binvisits,
	estimate_transferoperator_from_binvisits,
	invariantmeasure

# Estimate transfer entropy from scratch from a random
# set of points `n_realizations` times.
ts_length = 100

# Estimate transfer entropy using the wrapper function,
# and manually (doing the steps performed inside the
# wrapper function) to make sure that they give the same
# result.
estimates_3D_wrapper = Vector{Float64}(undef, n_realizations)
estimates_3D_allsteps = Vector{Float64}(undef, n_realizations)

@testset "3D #$i" for i in 1:n_realizations
	E = cembed([diff(rand(ts_length)) for i = 1:3])
	ϵ = 3
	# Test by doing all the dirty work and providing the raw input to the estimator
	bins_visited_by_orbit = assign_bin_labels(E, ϵ)
	bininfo = get_binvisits(bins_visited_by_orbit)
	TO = estimate_transferoperator_from_binvisits(bininfo)
	iv = invariantmeasure(TO)
	v = TEVars([1], [2], [3], Int[])

	estimates_3D_wrapper[i] = tetogrid(E, ϵ, v)
	estimates_3D_allsteps[i] = tetogrid(bins_visited_by_orbit, iv, v)
	@test estimates_3D_wrapper[i] >= 0
	@test estimates_3D_allsteps[i] >= 0
end

ts_length = 200
estimates_4D_wrapper = Vector{Float64}(undef, n_realizations)
estimates_4D_allsteps = Vector{Float64}(undef, n_realizations)

@testset "4D #$i" for i in 1:n_realizations
	E = cembed([diff(rand(ts_length)) for i = 1:4])
	ϵ = 0.3
	bins_visited_by_orbit = assign_bin_labels(E, ϵ)
	bininfo = get_binvisits(bins_visited_by_orbit)
	TO = estimate_transferoperator_from_binvisits(bininfo)
	iv = invariantmeasure(TO)
	v = TEVars([1], [2], [3, 4], Int[])
	estimates_4D_wrapper[i] = tetogrid(E, ϵ, v)
	estimates_4D_allsteps[i] = tetogrid(bins_visited_by_orbit, iv, v)
	@test estimates_4D_wrapper[i] >= 0
	@test estimates_4D_allsteps[i] >= 0
end

ts_length = 300
estimates_5D_wrapper = Vector{Float64}(undef, n_realizations)
estimates_5D_allsteps = Vector{Float64}(undef, n_realizations)

@testset "5D #$i" for i in 1:n_realizations
	E = cembed([diff(rand(ts_length)) for i = 1:5])
	ϵ = [0.2, 0.2, 0.1, 0.2, 0.3]

	bins_visited_by_orbit = assign_bin_labels(E, ϵ)
	bininfo = get_binvisits(bins_visited_by_orbit)
	TO = estimate_transferoperator_from_binvisits(bininfo)
	iv = invariantmeasure(TO)
	v = TEVars([1], [2], [3, 4], [5])
	estimates_5D_wrapper[i] = tetogrid(E, ϵ, v)
	estimates_5D_allsteps[i] = tetogrid(bins_visited_by_orbit, iv, v)
	@test estimates_5D_wrapper[i] >= 0
	@test estimates_5D_allsteps[i] >= 0
end

# # If everything  works as expected, there should be no negative
# # transfer entropy values.
#@show estimates_3D_wrapper
#@show estimates_4D_wrapper
#@show estimates_5D_wrapper
#show estimates_3D_allsteps
#@show estimates_4D_allsteps
#@show estimates_5D_allsteps

#@test all(estimates_3D_wrapper .>= 0)
#@test all(estimates_4D_wrapper .>= 0)
#@test all(estimates_5D_wrapper .>= 0)
#@test all(estimates_3D_allsteps .>= 0)
#@test all(estimates_4D_allsteps .>= 0)
#@test all(estimates_5D_allsteps .>= 0)

# The invariant distribution is estimated independently in the calls
# to the different transfer entropy estimators. The distribution is
# estimated by repeated application of the transfer operator on a
# randomly initialised distribution until convergence is achieved.
# Because a different initial distribution is used for each function
# call, the estimated transfer entropy will be slightly different.
# However, we do not expect the differences to be very large.
#@test all(abs.(estimates_3D_wrapper .- estimates_3D_allsteps) .< 1e-2)
#@test all(abs.(estimates_4D_wrapper .- estimates_4D_allsteps) .< 1e-2)
#@test all(abs.(estimates_5D_wrapper .- estimates_5D_allsteps) .< 1e-2)


text:
#####
##### Table
#####

export Dealer, SmallBlind, BigBlind, FirstToAct
export Table
export move_buttons!

Base.@kwdef mutable struct Winners
    declared::Bool = false
    players::Union{Nothing,Tuple,Player} = nothing
end

function Base.show(io::IO, winners::Winners, include_type = true)
    include_type && println(io, typeof(winners))
    println(io, "Winners declared = $(winners.declared)")
    println(io, "Winners          = $(winners.players)")
end

struct Blinds{S,B}
    small::S
    big::B
end

Blinds() = Blinds(1,2) # default
default_dealer_id() = 1 # default

function Base.show(io::IO, blinds::Blinds, include_type = true)
    include_type && println(io, typeof(blinds))
    println(io, "Blinds           = (small=$(blinds.small),big=$(blinds.big))")
end

struct Buttons
    dealer::Int
    small_blind::Int
    big_blind::Int
    first_to_act::Int
end

dealer_id(b::Buttons) = b.dealer
small_blind_id(b::Buttons) = b.small_blind
big_blind_id(b::Buttons) = b.big_blind
first_to_act_id(b::Buttons) = b.first_to_act

buttons(b::Buttons) = (
    b.dealer,
    b.small_blind,
    b.big_blind,
    b.first_to_act,
)

mutable struct Table
    deck::PlayingCards.Deck
    players::Tuple
    cards::Union{Nothing,Tuple{<:Card,<:Card,<:Card,<:Card,<:Card}}
    blinds::Blinds
    pot::Float64
    state::AbstractGameState
    buttons::Buttons
    current_raise_amt::Float64
    initial_round_raise_amt::Float64
    transactions::TransactionManager
    winners::Winners
    play_out_game::Bool
    n_max_actions::Int
end

buttons(table::Table) = table.buttons
n_max_actions(table::Table) = table.n_max_actions

"""
    compute_n_max_actions(players, bb)

The total possible number of actions,
which happens when
 - pre-flop: all players check/call
 - flop: all players check
 - turn: all players check
 - river: 1) all but one player checks (n_players-1 actions)
          2) the last player raises (1 action)
          3) everyone, except the last player, calls, (n_players-1 actions)
          4) the last player re-raises (1 action)
          ...
          until all players are all-in.
"""
function compute_n_max_actions(players, bb)
    maxbr = maximum(bank_roll.(players))
    n_players = length(players)
    n_check_call_rounds = n_players*3 # PreFlop, Flop, Turn
    n_raise_rounds = findfirst(1:typemax(Int)) do i
        maxbr < bb^n_raises(i, n_players)
    end
    return n_check_call_rounds + n_raise_rounds
end
n_raises(i, n_players) = Int(floor(i/n_players))

function Table(;
    players::Tuple,
    deck = ordered_deck(),
    cards = nothing,
    blinds = Blinds(),
    pot = Float64(0),
    state = PreFlop(),
    dealer_id = default_dealer_id(),
    current_raise_amt = Float64(0),
    initial_round_raise_amt = blinds.small,
    transactions = nothing,
    winners = Winners(),
    play_out_game = false,
)
    if transactions == nothing
        transactions = TransactionManager(players)
    end
    buttons = Buttons(dealer_id, players)
    n_max_actions = compute_n_max_actions(players, blinds.big)
    @debug "n_max_actions = $n_max_actions"
    return Table(deck,
        players,
        cards,
        blinds,
        pot,
        state,
        buttons,
        current_raise_amt,
        initial_round_raise_amt,
        transactions,
        winners,
        play_out_game,
        n_max_actions)
end

function Buttons(dealer_id, players::Tuple)
    dealer_id       = this_or_next_valid_id(dealer_id       , players)
    small_blind_id  = this_or_next_valid_id(dealer_id     +1, players)
    big_blind_id    = this_or_next_valid_id(small_blind_id+1, players)
    first_to_act_id = this_or_next_valid_id(big_blind_id  +1, players)
    return Buttons(dealer_id, small_blind_id, big_blind_id, first_to_act_id)
end

valid_button(player::Player) = !zero_bank_roll(player) && still_playing(player)

"""
    this_or_next_valid_id(id, players::Tuple)

Given an index `id`, return this index, if valid,
otherwise find the next valid index.
"""
function this_or_next_valid_id(id, players::Tuple)
    n_players = length(players)
    id = circle_index(n_players, id) # so that we can pass in id+1
    if valid_button(players[id])
        return id
    else
        state = id
        id_need_defined = true
        while id_need_defined
            state > 2*n_players && error("Broken logic in this_or_next_valid_id")
            id = circle_index(n_players, state)
            if valid_button(players[id])
                return id
            end
            state+=1
        end
    end
end

function Base.show(io::IO, table::Table, include_type = true)
    include_type && println(io, typeof(table))
    show(io, blinds(table), false)
    show(io, table.winners, false)
    println(io, "Dealer           = $(dealer_id(table))")
    println(io, "Pot              = $(table.transactions)")
    println(io, "All cards        = $(table.cards)")
    println(io, "Observed cards   = $(observed_cards(table))")
end

get_table_cards!(deck::PlayingCards.Deck) =
    Iterators.flatten(ntuple(i->pop!(deck, 1), 5)) |> collect |> Tuple
cards(table::Table) = table.cards

observed_cards(table::Table) = observed_cards(table, table.state)
observed_cards(table::Table, ::PreFlop) = ()
observed_cards(table::Table, ::Flop) = table.cards[1:3]
observed_cards(table::Table, ::Turn) = table.cards[1:4]
observed_cards(table::Table, ::River) = table.cards
current_raise_amt(table::Table) = table.current_raise_amt
initial_round_raise_amt(table::Table) = table.initial_round_raise_amt
minimum_raise_amt(table::Table) = blinds(table).small

state(table::Table) = table.state

play_out_game(table::Table) = table.play_out_game

dealer_id(table::Table) = dealer_id(table.buttons)
small_blind_id(table::Table) = small_blind_id(table.buttons)
big_blind_id(table::Table) = big_blind_id(table.buttons)
first_to_act_id(table::Table) = first_to_act_id(table.buttons)

dealer(table::Table) = players_at_table(table)[dealer_id(table)]
small_blind(table::Table) = players_at_table(table)[small_blind_id(table)]
big_blind(table::Table) = players_at_table(table)[big_blind_id(table)]
first_to_act(table::Table) = players_at_table(table)[first_to_act_id(table)]

is_dealer(table::Table, player::Player) = seat_number(player) == seat_number(dealer(table))
is_small_blind(table::Table, player::Player) = seat_number(player) == seat_number(small_blind(table))
is_big_blind(table::Table, player::Player) = seat_number(player) == seat_number(big_blind(table))
is_first_to_act(table::Table, player::Player) = seat_number(player) == seat_number(first_to_act(table))

players_at_table(table::Table) = table.players
all_playing_checked(table::Table) = all(map(x -> not_playing(x) || checked(x), players_at_table(table)))
all_playing_all_in(table::Table) = all(map(x -> not_playing(x) || all_in(x), players_at_table(table)))
all_all_in_or_checked(table::Table) = all(map(x -> checked(x) || all_in(x), players_at_table(table)))

"""
    bank_roll_leader(table::Table)

Return the player who
 - Is still playing and
 - Has the highest bank roll
and a `Bool` indicating if there
are multiple players with the largest
bank roll.
"""
function bank_roll_leader(table::Table)
    max_rbr = 0
    players = players_at_table(table)
    br_leader = first(players)
    for player in players
        still_playing(player) || continue # only consider players still playing
        pbr = round_bank_roll(player)
        if pbr > max_rbr
            br_leader = player
            max_rbr = pbr
        end
    end
    multiple_leaders = count(map(players) do player
        round_bank_roll(player) ≈ max_rbr && still_playing(player)
    end) > 1
    return br_leader, multiple_leaders
end

paid_current_raise_amount(table::Table, player::Player) =
    round_contribution(player) ≈ current_raise_amt(table)

# Can be true in exactly 2 cases:
#  1) Everyone (still playing) is all-in.
#  2) Everyone (still playing), except `player`, is all-in.
function all_oppononents_all_in(table::Table, player::Player)
    all_opp_all_in = true
    for opponent in players_at_table(table)
        seat_number(opponent) == seat_number(player) && continue
        not_playing(opponent) && continue
        if action_required(opponent)
            all_opp_all_in = false
        else
            all_opp_all_in = all_opp_all_in && all_in(opponent)
        end
    end
    return all_opp_all_in
end

# One case that we need to catch is when everyone, except
# the bank roll leader, has folded or gone all-in. In this
# case, while the bank roll leader may have required actions
# everyone else does not, so nobody can respond to their raise
# (if they chose to do so). Therefore, we must "play out" the
# entire game with no further actions.
# Note that this method is intended to be used _during_ a round.
function all_all_in_except_bank_roll_leader(table::Table)
    br_leader, multiple_leaders = bank_roll_leader(table)
    players = players_at_table(table)
    multiple_leaders && return false # the bank roll leader can go all-in

    @assert !multiple_leaders # We have a single bank roll leader

    return all(map(players) do player
        cond1 = not_playing(player)
        cond2 = all_in(player)
        cond3 = seat_number(player) == seat_number(br_leader) && !action_required(br_leader)
        any((cond1, cond2, cond3))
    end)
end

# Note that this method is only valid before or after a round has ended.
function set_play_out_game!(table::Table)
    br_leader, multiple_leaders = bank_roll_leader(table)
    players = players_at_table(table)
    if multiple_leaders
        return all(map(player -> not_playing(player) || all_in(player), players))
    end

    @assert !multiple_leaders # We have a single bank roll leader

    table.play_out_game = all(map(players) do player
        cond1 = not_playing(player)
        cond2 = all_in(player)
        cond3 = seat_number(player) == seat_number(br_leader)
        any((cond1, cond2, cond3))
    end)
end

blinds(table::Table) = table.blinds

function is_blind_call(table::Table, player::Player, amt = call_amount(table, player))
    pot_inv = pot_investment(player)
    @debug "amt = $amt, pot_investment(player) = $pot_inv"
    bb = blinds(table).big
    sb = blinds(table).small
    if is_small_blind(table, player)
        return amt ≈ sb && pot_inv ≈ sb
    elseif is_big_blind(table, player)
        return amt ≈ 0 && pot_inv ≈ bb
    else
        return amt ≈ bb && pot_inv ≈ 0
    end
end

function reset_round_bank_rolls!(table::Table)
    players = players_at_table(table)
    for player in players
        player.round_bank_roll = bank_roll(player)
    end
end

function reset_round!(table::Table)
    players = players_at_table(table)
    for player in players
        not_playing(player) && continue
        all_in(player) && continue
        player.checked = false
        player.action_required = true
        player.last_to_raise = false
        player.round_contribution = 0
    end
    table.initial_round_raise_amt = blinds(table).small
    table.current_raise_amt = 0
end

function set_state!(table::Table, state::AbstractGameState)
    table.state = state
end

# Check for winner, in case when only a single player remains
# playing
function check_for_and_declare_winner!(table::Table)
    players = players_at_table(table)
    n_players = length(players)
    table.winners.declared = count(not_playing.(players)) == n_players-1
    if table.winners.declared
        for player in players
            not_playing(player) && continue
            table.winners.players = player
        end
    end
end

#####
##### Circling the table
#####

"""
    move_buttons!(table::Table)

Move the dealer, small blind, big blind,
and first-to-act buttons to the next set
of players.
"""
function move_buttons!(table::Table)
    players = players_at_table(table)
    dealer = table.buttons.dealer
    dealer       = this_or_next_valid_id(dealer      + 1, players)
    small_blind  = this_or_next_valid_id(dealer      + 1, players)
    big_blind    = this_or_next_valid_id(small_blind + 1, players)
    first_to_act = this_or_next_valid_id(big_blind   + 1, players)
    table.buttons = Buttons(dealer, small_blind, big_blind, first_to_act)
end

"""
    circle_index(n_players, state)

A circular index for indexing into `players`.
`state = 1` corresponds to `player[1]`.
"""
circle_index(n_players, i) = mod(i-1, n_players)+1

any_actions_required(table::Table) = any(action_required.(players_at_table(table)))

abstract type TablePosition end
struct Dealer <: TablePosition end
struct SmallBlind <: TablePosition end
struct BigBlind <: TablePosition end
struct FirstToAct <: TablePosition end # (after BigBlind)

struct CircleTable{CircType,P, NITER}
    players::Tuple
    buttons::Buttons
    n_players::Int
    player::P
end

n_iterations(::CircleTable{CircType,P,NITER}) where {CircType,P,NITER} = NITER

function Base.length(ct::CircleTable)
    n_iter = n_iterations(ct)
    if !(n_iter === Inf)
        return n_iter
    else
        error("CircleTable has an infinite length")
    end
end

circle(table::Table, tp::TablePosition, n_iter = Inf) =
    CircleTable{typeof(tp),Nothing, n_iter}(table.players, buttons(table), length(table.players), nothing)

circle(table::Table, player::Player, n_iter = Inf) =
    CircleTable{typeof(player),typeof(player), n_iter}(table.players, buttons(table), length(table.players), player)

function Base.iterate(ct::CircleTable{Dealer}, state = dealer_id(ct.buttons))
    state > n_iterations(ct)+dealer_id(ct.buttons)-1 && return nothing
    i_circ = circle_index(ct.n_players, state)
    player = ct.players[i_circ]
    state == dealer_id(ct.buttons) && @assert active(player)
    return (player, state+1)
end

function Base.iterate(ct::CircleTable{SmallBlind}, state = small_blind_id(ct.buttons))
    state > n_iterations(ct)+small_blind_id(ct.buttons)-1 && return nothing
    i_circ = circle_index(ct.n_players, state)
    player = ct.players[i_circ]
    state == small_blind_id(ct.buttons) && @assert active(player)
    return (player, state+1)
end

function Base.iterate(ct::CircleTable{BigBlind}, state = big_blind_id(ct.buttons))
    state > n_iterations(ct)+big_blind_id(ct.buttons)-1 && return nothing
    i_circ = circle_index(ct.n_players, state)
    player = ct.players[i_circ]
    state == big_blind_id(ct.buttons) && @assert active(player)
    return (player, state+1)
end

function Base.iterate(ct::CircleTable{FirstToAct}, state = first_to_act_id(ct.buttons))
    state > n_iterations(ct)+first_to_act_id(ct.buttons)-1 && return nothing
    i_circ = circle_index(ct.n_players, state)
    player = ct.players[i_circ]
    state == first_to_act_id(ct.buttons) && @assert active(player)
    return (player, state+1)
end

function Base.iterate(ct::CircleTable{P},
        state = circle_index(ct.n_players, seat_number(ct.player))
    ) where {P <: Player}
    state > n_iterations(ct)+circle_index(ct.n_players, seat_number(ct.player))-1 && return nothing
    i_circ = circle_index(ct.n_players, state)
    player = ct.players[i_circ]
    state == circle_index(ct.n_players, seat_number(ct.player)) && @assert active(player)
    return (player, state+1)
end

show_cards(table::Table, player::Player{Human}) = player.cards
show_cards(table::Table, player::Player{LF}) where {LF <: AbstractAI} = "(??,??)"

#####
##### Deal
#####

function deal!(table::Table, blinds::Blinds)
    players = players_at_table(table)
    shuffle!(table.deck)
    call_blinds = true
    for (i, player) in enumerate(circle(table, SmallBlind()))

        i>length(players) && break # deal cards to each player once

        not_playing(player) && continue

        player.cards = pop!(table.deck, 2)

        if is_small_blind(table, player) && bank_roll(player) ≤ blinds.small
            contribute!(table, player, bank_roll(player), call_blinds)
            @info "$(name(player)) paid the small blind (all-in) and dealt cards: $(show_cards(table, player))"
        elseif is_big_blind(table, player) && bank_roll(player) ≤ blinds.big
            contribute!(table, player, bank_roll(player), call_blinds)
            @info "$(name(player)) paid the  big  blind (all-in) and dealt cards: $(show_cards(table, player))"
        else
            if is_small_blind(table, player)
                contribute!(table, player, blinds.small, call_blinds)
                @info "$(name(player)) paid the small blind and dealt cards: $(show_cards(table, player))"
            elseif is_big_blind(table, player)
                contribute!(table, player, blinds.big, call_blinds)
                @info "$(name(player)) paid the  big  blind and dealt cards: $(show_cards(table, player))"
            else
                @info "$(name(player)) dealt (free) cards:                   $(show_cards(table, player))"
            end
        end
    end

    table.cards = get_table_cards!(table.deck)
    @info "Table cards dealt (face-down)."
end


text:
cluster_block(isfirst::Val{true}) = chain(2, [repeat(2, H, 1:2), control(2, 1, 2=>Z)])
cluster_block(isfirst::Val{false}) = chain(2, [swap(2, 1, 2), put(2, 2=>H), control(2, 1, 2=>Z)])

function cluster_circuit(nrepeat::Int)
    sequence([cluster_block(Val(i==1)) for i=1:nrepeat])
end

function model(::Val{:cluster}; nbit, B=4096)
    nrepeat = nbit - 1
    c = cluster_circuit(nrepeat)
    chem = QuantumMPS(1, 1, 0, c, zero_state(2, B), zeros(Int, nbit))
    chem
end

max_stars_repo_name: Roger-luo/QuantumMPS

max_issues_repo_name: Roger-luo/QuantumMPS

max_forks_repo_name: Roger-luo/QuantumMPS


text:
using LinearAlgebra
include("kernel.jl")

function sample(K)
    n = size(K, 1)
    C = cholesky(K)
    z = C.L*(randn(n))
    return z
end

function getExample(dim, n, a, b, θ)
    if dim>1
        s = [3 6 1 3; 4 2 3 2; 2 1 1 5; 1 4 2 3;5 6 7 8]
        s0 = [1 2 3 2; 2 4 2 1; 3 1 2 6; 1 9 4 2; 2 3 8 6]
        X = [3 4 4; 9 3 5; 1 7 13; 4 1 2; 5 6 14]
        X0 = [1 1 2 ; 3 2 1; 4 -1 3; 5 5 4; 8 -3 5]
        z = [9; 11; 13; 6; 7]
        example = setting(s, s0, X, X0, z)
    else
        #s = [3 6 1 3; 4 2 3 2; 2 1 1 5; 1 4 2 3;5 6 7 8]
        #s0 = [1 2 3 2]
        s = rand(n, a)*100
        s0 = rand(1, a)*100
        X0 = rand(1, b)*100
        #X0 = [1 1 2]
        X = rand(n, b)*100
        kernelm = K(s, s, θ)
        z = sample(kernelm)
        z = reshape(z, length(z))
        #z = [9; 11; 13; 6; 7;]
        example = setting(s, s0, X, X0, z)
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module HGT

function testfunction()
    println("Testing, five, six, EIGHT")
end

end # Module HGT


text:
"""
`targetlist(targetexpr::Regex, dir2walk::AbstractString)` returns `targetpaths, allfiles`
"""
function targetlist(targetexpr::Regex, dir2walk::AbstractString)
    # targetexpr = r"(\.jl)$";
    fulllist = String[];
    for (root, folders, files) in walkdir(dir2walk)
        for file in files
            push!(fulllist, joinpath(root, file));
        end
    end
    allfiles = basename.(fulllist);
    targetid = occursin.(targetexpr, allfiles);

    targetpaths = fulllist[targetid];
    return targetpaths, allfiles
end


text:
"""
    ht_p2_p1_gen

The manufactured-solution colliding flow example from Elman et al 2014. The
Hood-Taylor formulation with quadratic triangles for the velocity and continuous
pressure on linear triangles.

The formulation is the general elasticity-like scheme with
strain-rate-displacement matrices. It can be manipulated into the one derived in
Reddy, Introduction to the finite element method, 1993. Page 486 ff.
"""
module ht_p2_p1_gen

using LinearAlgebra
using StaticArrays
using MeshCore.Exports
using MeshSteward.Exports
using Elfel.Exports
using UnicodePlots

mu = 1.0 # dynamic viscosity
D = SMatrix{3, 3}(
    [2*mu 0 0
     0 2*mu 0
     0 0 mu])
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(2 * A, 2 * A, N, N), "velocity")
    ir = baseincrel(vmesh)
    transform(ir, x -> x .- A)
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(Uh, Ph, tndof, D)
    function integrateK!(ass, elits, qpits, D)
        B = (g, k) -> (k == 1 ?
            SVector{3}((g[1], 0, g[2])) :
            SVector{3}((0, g[2], g[1])))
        c = edofcompnt(Uh)
        unedof, pnedof = ndofsperel.((Uh, Ph))
        kuu = LocalMatrixAssembler(unedof, unedof, 0.0)
        kup = LocalMatrixAssembler(unedof, pnedof, 0.0)
        for el in zip(elits...)
            uel, pel = el
            init!(kuu, eldofs(uel), eldofs(uel))
            init!(kup, eldofs(uel), eldofs(pel))
            for qp in zip(qpits...)
                uqp, pqp = qp
                Jac, J = jacjac(uel, uqp)
                JxW = J * weight(uqp)
                gradNu = bfungrad(uqp, Jac)
                Np = bfun(pqp)
                for j in 1:unedof
                    DBj = D * B(gradNu[j], c[j])
                    for i in 1:unedof
                        Bi = B(gradNu[i], c[i])
                        kuu[i, j] += dot(Bi, DBj) * (JxW)
                    end
                end
                for j in 1:pnedof, i in 1:unedof
                    kup[i, j] += (-JxW * Np[j]) * gradNu[i][c[i]]
                end
            end
            assemble!(ass, kuu)
            assemble!(ass, kup)
            assemble!(ass, transpose(kup))
        end
        return ass
    end

    elits = (FEIterator(Uh), FEIterator(Ph))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(Uh, qargs), QPIterator(Ph, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrateK!(ass, elits, qpits, D)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(Ph)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(Ph)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(Ph, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(Uh)
    function integrate!(elit, qpit, trueux, trueuy)
        unedof = ndofsperel(elit)
        uedofcomp = edofcompnt(Uh)
        E = 0.0
        for el in elit
            udofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Nu = bfun(qp)
                uxt = trueux(location(el, qp)...)
                uyt = trueuy(location(el, qp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:unedof
                    (uedofcomp[j] == 1) && (uxa += (udofvals[j] * Nu[j]))
                    (uedofcomp[j] == 2) && (uya += (udofvals[j] * Nu[j]))
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(Uh)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(Uh, qargs)
    return integrate!(elit, qpit, trueux, trueuy)
end

function run(N)
    vmesh, pmesh = genmesh(N)
    # Velocity space: space with two components
    Uh = FESpace(Float64, vmesh, FEH1_T6(), 2)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(Uh, 0, i, 1, trueux(locs[i]...))
            setebc!(Uh, 0, i, 2, trueuy(locs[i]...))
        end
    end
    # Pressure space
    Ph = FESpace(Float64, pmesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(pmesh); nearestto = [0.0, 0.0])
    setebc!(Ph, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([Uh, Ph])
    tndof = ndofs(Uh) + ndofs(Ph)
    tnunk = nunknowns(Uh) + nunknowns(Ph)
    # Assemble the coefficient matrix
    K = assembleK(Uh, Ph, tndof, D)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [Uh, Ph])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([Uh, Ph], U)
    # Postprocessing
    makeattribute(Ph, "p", 1)
    makeattribute(Uh, "ux", 1)
    makeattribute(Uh, "uy", 2)
    ep = evaluate_pressure_error(Ph)
    ev = evaluate_velocity_error(Uh)
    vtkwrite("ht_p2_p1_gen-p", baseincrel(pmesh), [(name = "p",), ])
    vtkwrite("ht_p2_p1_gen-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    # geom = geometry(Ph.mesh)
    # ir = baseincrel(Ph.mesh)
    # pt = VecAttrib([truep(geom[i]...) for i in 1:length(geom)])
    # ir.right.attributes["pt"] = pt
    # vtkwrite("ht_p2_p1_gen-pt", baseincrel(pmesh), [(name = "pt",), ])
    return (ep, ev)
end

end

using .ht_p2_p1_gen
let
    N = 4
    for loop in 1:5
        ep, ev = ht_p2_p1_gen.run(N)
        @show N, ep, ev
        N = N * 2
    end
end


text:
"""
hasconverged(s::ConvergenceState)

Returns true if `s.converged` is true.
"""
hasconverged(s::ConvergenceState) = s.converged

"""
hasconverged(s::Solution)

Returns true if `hasconverged(s.convstate)` is true.
"""
hasconverged(s::Solution) = hasconverged(s.convstate)

"""
```
assess_convergence!(
    solution::Solution,
    model::AbstractModel,
    tol::Tolerance,
    criteria::ConvergenceCriteria,
)
```

Evaluates the convergence state `solution.convstate` given the current solution, `solution`, the tolerance, `tol`, and the convergence criteria `criteria`. `solution.convstate.converged` is then updated.

If `criteria` is an instance of `GenericCriteria`, `converged = (x_converged || f_converged) && infeas_converged`. `x_converged`, `f_converged` and `infeas_converged` are explained in [`Tolerance`](@ref). If `criteria` is an instance of `KKTCriteria` or `ScaledKKTCriteria`, `converged = kkt_converged && infeas_converged`. `kkt_converged` and `infeas_converged` are explained in [`Tolerance`](@ref). If `criteria` is an instance of `IpoptCriteria`, `converged = ipopt_converged && infeas_converged`. `ipopt_converged` and `infeas_converged` are explained in [`Tolerance`](@ref).
"""
function assess_convergence!(
    solution::Solution,
    model::AbstractModel,
    tol::Tolerance,
    criteria::ConvergenceCriteria,
)
    xtol, fabstol, freltol, kkttol, infeastol = tol.x, tol.fabs, tol.frel, tol.kkt, tol.infeas
    Δx, Δf, infeas = getresiduals(solution, model, GenericCriteria())
    relΔf = Δf / (abs(solution.f) + freltol)
    kkt_residual, infeas = getresiduals(solution, model, KKTCriteria())
    ipopt_residual, infeas = getresiduals(solution, model, IpoptCriteria())
    if show_residuals[]
        @show kkt_residual, ipopt_residual, kkttol
    end

    x_converged = Δx < xtol
    fabs_converged = Δf < fabstol
    frel_converged = relΔf < freltol
    if criteria isa ScaledKKTCriteria
        if debugging[]
            #@show get_objective_multiple(model)
        end
        m = get_objective_multiple(model)
        kkt_residual = kkt_residual / max(m, 1/m)
    end
    kkt_converged = kkt_residual < kkttol
    ipopt_converged = ipopt_residual < kkttol
    infeas_converged = infeas <= infeastol
    f_increased = solution.f > solution.prevf

    if criteria isa GenericCriteria
        converged = (x_converged || fabs_converged || frel_converged) && infeas_converged
    elseif criteria isa KKTCriteria || criteria isa ScaledKKTCriteria
        converged = kkt_converged && infeas_converged
    elseif criteria isa IpoptCriteria
        converged = ipopt_converged && infeas_converged
    else
        throw("Unsupported convergence criteria for MMA.")
    end
    @pack! solution.convstate = x_converged,
                                fabs_converged,
                                frel_converged,
                                kkt_converged,
                                ipopt_converged,
                                infeas_converged,
                                Δx,
                                Δf,
                                relΔf,
                                kkt_residual,
                                ipopt_residual,
                                infeas,
                                f_increased,
                                converged
    return solution
end

function getresiduals(solution::Solution, ::AbstractModel, ::GenericCriteria)
    @unpack prevx, x, prevf, f, g = solution
    Δx = maximum(abs(x[j] - prevx[j]) for j in 1:length(x))
    Δf = abs(f - prevf)
    infeas = length(g) == 0 ? zero(eltype(g)) : max(0, maximum(g))
    return Δx, Δf, infeas
end
function getresiduals(solution::Solution, model::AbstractModel, ::KKTCriteria)
    @unpack ∇f, g, ∇g, λ, x = solution
    xmin, xmax = getmin(model), getmax(model)
    T = eltype(x)
    res = maximum(1:length(x)) do j
        @views temp = ∇f[j] + dot(∇g[:,j], λ)
        if xmin[j] >= x[j]
            return abs(min(0, temp))
        elseif x[j] >= xmax[j]
            return max(0, temp)
        else
            return abs(temp)
        end
    end
    if debugging[]
        @show λ, g
    end
    res = length(g) == 0 ? res : max(res, maximum(abs.(λ .* g)))
    if debugging[]
        @show maximum(abs, g)
        @show maximum(abs, λ)
        @show maximum(x)
    end
    infeas = length(g) == 0 ? zero(eltype(g)) : max(maximum(g), 0)
    return res, infeas
end
function getresiduals(solution::Solution, model::AbstractModel, ::IpoptCriteria)
    @unpack ∇f, g, ∇g, λ, x = solution
    xmin, xmax = getmin(model), getmax(model)
    T = eltype(x)
    n, m, s = length(x), length(λ), zero(T)
    res = maximum(1:n) do j
        @views temp = ∇f[j] + dot(∇g[:,j], λ)
        if xmin[j] >= x[j]
            dj = temp
            s += max(dj, 0)
            return abs(min(0, dj))
        elseif x[j] >= xmax[j]
            yj = -temp
            s += max(yj, 0)
            return abs(min(0, yj))
        else
            return abs(temp)
        end
    end
    sd = max(100, (sum(abs, λ) + s) / (n + m)) / 100
    res = length(g) == 0 ? res : max(res, maximum(abs.(λ .* g)))
    res = res / sd
    infeas = length(g) == 0 ? zero(eltype(g)) : max(maximum(g), 0)
    if debugging[]
        println("Agg infeas = ", infeas)
    end
    return res, infeas
end


text:
module Yandex
using Gumbo
using Cascadia
using Franklin: globvar, locvar, pagevar, path
using JSON
using FranklinContent: walkfiles, islink, isldjson

const header = HTMLElement(:header)
const heading = HTMLText("")
const subheading = HTMLText("")
const figure = HTMLElement(:figure)
const image = HTMLElement(:image)
const menu = HTMLElement(:menu)
const crumbs_html = HTMLElement(:div)
const crumbs_links = Vector{Tuple{String, String}}()
const crumbs_io = IOBuffer()
const a_tag = HTMLElement(:a)
const page_url = Ref("")
const page_date = Ref("")
const page_id = Ref("")
const page_descr = Ref("")
const page_lang = Ref("")
const page_kws = []
const turbo_dir = Ref("")

skip_nodes = Set(HTMLElement{sym} for sym in [:style, :iframe, :applet, :audio, :canvas, :embed, :video, :img, :button, :form])

@doc "Convert a string to an HTMLText object."
macro ht_str(text)
    quote
        HTMLText($text)
    end
end

@doc "Set up yandex turbo page header elements."
function init_header()
    empty!(header.children)
    empty!(header.attributes)
    h1 = HTMLElement(:h1)
    push!(h1, heading)
    push!(header, h1)

    h2 = HTMLElement(:h2)
    push!(h2, subheading)
    push!(header, h2)

    push!(figure, image)
    push!(header, figure)

    setattr!(crumbs_html, "data-block", "breadcrumblist")
    push!(header, menu)
    push!(header, crumbs_html)
end

function set_header(;title::AbstractString,
                    subtitle::AbstractString="",
                    img_url::AbstractString="",
                    menu_links::AbstractVector=[])
    heading.text = title
    subheading.text = subtitle
    setattr!(image, "src", img_url)
    @debug @assert length(a_tag.attributes) === 0 && length(a_tag.children) === 0
    empty!(menu.children)
    for link in menu_links
        a = deepcopy(a_tag)
        setattr!(a, "src", link)
        push!(menu, a)
    end
    empty!(crumbs_html.children)
    for (_, link) in crumbs_links
        a = deepcopy(a_tag)
        setattr!(a, "src", link)
        push!(crumbs_html, a)
    end
end

@inline function is_script_id(el, tp, id="")
    tp === HTMLElement{:script} &&
        getattr(el, "id", "") === id
end

function process_head(in_head)
    canonical_unset = true
    title_unset = true
    title = ""
    subtitle = ""
    sub_unset = true
    empty!(crumbs_links)
    crumbs_unset = true
    date_unset = true
    empty!(page_kws)
    for el in in_head.children
        tp = typeof(el)
        if tp ∈ skip_nodes continue end
        if title_unset && tp === HTMLElement{:title}
            title = text(el)
            title_unset = false
            # use the meta "description" tag for the subtitle
        elseif canonical_unset &&
            tp === HTMLElement{:link} &&
            islink(el, "canonical")
            page_url[] = getattr(el, "href", "")
            canonical_unset = false
        elseif sub_unset && tp === HTMLElement{:meta} && hasattr(el, "description")
            subtitle = getattr(el, "description")
            sub_unset = false
            # get page published date
        elseif date_unset && is_script_id(el, tp, "ldj-webpage")
            data = el |> text |> JSON.parse
            page_date[] = data["datePublished"]
            date_unset= false
            page_id[] = data["mainEntityOfPage"]["@id"]
            append!(page_kws, data["keywords"])
            # parse the markup for breadcrumbs
        elseif crumbs_unset && is_script_id(el, tp, "ldj-breadcrumbs")
            data = el |> text |> JSON.parse
            for list_el in data["itemListElement"]
                push!(crumbs_links, (list_el["name"], list_el["item"]))
            end
            crumbs_unset = false
        end
        title_unset || sub_unset || crumbs_unset || break
    end
    menu_links = globvar(:menu; default=[])
    set_header(;title, subtitle, img_url=locvar(:image; default=globvar(:author_image)), menu_links)
end

function breadcrumbs_tags()
	for (name, link) in crumbs_links
        write(crumbs_io, "<breadcrumb url=\"")
        write(crumbs_io, link)
        write(crumbs_io, "\" text=\"")
        write(crumbs_io, name)
        write(crumbs_io, "\" />")
    end
    String(take!(crumbs_io))
end

# function yandex_related(file)

# end

const content_sel = Selector("body")
const head_sel = Selector("head")
const turbo_item_content = HTMLElement(:body)
const io_item = IOBuffer()
@doc "Generates a turbo page feed item from a html file."
function turboitem(file; cosel=content_sel)
    @debug file
    html = parsehtml(read(file, String))
    # content = Cascadia.matchFirst(cosel, html.root)
    head = Cascadia.matchFirst(head_sel, html.root)
    body = Cascadia.matchFirst(sel"body", html.root)

    empty!(turbo_item_content.children)
    empty!(turbo_item_content.attributes)

    page_lang[] = getattr(html.root, "lang", globvar(:lang_code; default="en"))
    process_head(head)
   """
<item turbo="true">
            <!-- Page information -->
            <turbo:extendedHtml>true</turbo:extendedHtml>
            <link>$(page_url[])</link>
            <language>$(page_lang[])</language>
            <!-- <turbo:source></turbo:source> -->
            <!-- <turbo:topic></turbo:topic> -->
            <pubDate>$(page_date[])</pubDate>
            <author>$(globvar(:author))</author>
            <metrics>
                <yandex schema_identifier="$(page_id[])">
                    <breadcrumblist>
                        $(breadcrumbs_tags())
                    </breadcrumblist>
                </yandex>
            </metrics>
            <!-- <yandex:related></yandex:related> -->
            <turbo:content>
                <![CDATA[
                    $(header)
                    $(body)
                ]]>
            </turbo:content>
</item>
"""
end

@inline get_file_feed(dir, n) = joinpath(dir, "turbo-$n.xml")

function init_feed(io)
	"""<?xml version = "1.0" encoding = "UTF-8"?>
    <rss xmlns:yandex="http://news.yandex.ru"
    xmlns:media="http://search.yahoo.com/mrss/"
    xmlns:turbo="http://turbo.yandex.ru"
    version="2.0">
    <channel>
        <!-- Information about the source site  -->
        <title>$(globvar(:website_title))</title>
        <link>$(globvar(:website_url))</link>
        <description>$(globvar(:website_description))</description>
        <language>$(globvar(:lang_code))</language>
        <!-- <turbo:analytics></turbo:analytics> <turbo:adNetwork></turbo:adNetwork> -->
    """ |> x -> write(io, x)
end

function write_feed(io, dir, file_feed, feed_counter)
end

@doc "Create yandex turbo pages for a directory, outputs an RSS feed.
These feeds need to be added to https://webmaster.yandex.com under 'Data Sources' tab."
function turbodir(target=nothing; extensions=[".html"], dirs=["posts", "tag", "reads", "_rss"], ex_dirs=["amp"])
    isnothing(target) && begin target = path(:site) end
    @assert isdir(target) "Path $target is not a valid directory"
    turbo_dir[] = dir = isdirpath(target) ? dirname(target) : target
    rpl = Regex("^$dir/") => ""
    cwd = pwd()
    cd(dir)

    item_counter = 0
    feed_counter = 1
    init_header()
    file_feed = get_file_feed(dir, feed_counter)
    io_feed = IOBuffer()
    init_feed(io_feed)
    try
        proc_dirs = Dict{String, Nothing}()
        for file in walkfiles(dir; exts=Set(extensions),
                              dirs=Set(dirs),
                              ex_dirs=Set(ex_dirs),
                              subdir=false)
            basename(file) === "404.html" && continue
            out_file = joinpath(dir, replace(file, rpl))
            out_dir = dirname(out_file)
            if out_dir ∉ keys(proc_dirs) && !isdir(out_dir)
                mkpath(out_dir)
                proc_dirs[out_dir] = nothing
            end
            item = turboitem(file)
            (item_counter > 999 ||
                io_feed.size + sizeof(item) > Int(15e6)) && begin
                item_counter = 0
                write(io_feed, "</channel></rss>")
                write(file_feed, String(take!(io_feed)))
                init_feed(io_feed)
                feed_counter += 1
                file_feed = get_file_feed(dir, feed_counter)
            end
            write(io_feed, item)
            item_counter += 1
        end
        write(io_feed, "</channel></rss>")
        write(file_feed, String(take!(io_feed)))
        if feed_counter > 100
            @warn "Yandex only support 100 feed lists."
        end
    finally
        cd(cwd)
        close(io_feed)
    end
end

end


text:
module SCS

using Libdl
import Requires
import SparseArrays

if haskey(ENV, "JULIA_SCS_LIBRARY_PATH") || VERSION < v"1.3"
    if !isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
        error(
            "SCS not properly installed. Please run `Pkg.build(\"SCS\")` and " *
            "restart julia",
        )
    end
    include(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    function __init__()
        version = VersionNumber(scs_version())
        if version < v"2.1.0"
            error(
                "Current SCS version installed is $version, but we require " *
                "version 2.1.*",
            )
        end
        return
    end
else
    import SCS_jll
    const indirect = SCS_jll.libscsindir
    const direct = SCS_jll.libscsdir
    function __init__()
        Requires.@require(
            CUDA_jll = "e9e359dc-d701-5aa8-82ae-09bbf812ea83",
            include("linear_solvers/gpu_indirect.jl"),
        )
        return
    end
end

include("c_wrapper.jl")
include("linear_solvers/direct.jl")
include("linear_solvers/indirect.jl")
include("MOI_wrapper/MOI_wrapper.jl")

const available_solvers = [DirectSolver, IndirectSolver]

export scs_solve

end


text:
#
println("Dr) Test of the MultiPhotonIonization module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-MultiPhotonIonization.sum")
wa = Atomic.Computation("xx",  Nuclear.Model(36.);
                        initialConfigs=[Configuration("1s^2 2s^2 2p^6")],
                        finalConfigs  =[Configuration("1s^2 2s^2 2p^5"), Configuration("1s^2 2s 2p^6") ],
                        process = MultiPI,
                        processSettings=MultiPhotonIonization.Settings() )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
##
## Types
##

MMI.@mlj_model mutable struct LSSVClassifier <: MMI.Deterministic
    kernel::Symbol = :rbf::(_ in (:rbf, :linear, :poly))
    γ::Float64 = 1.0::(_ > 0.0)
    σ::Float64 = 1.0::(_ > 0.0)
    degree::Int = 0::(_ >= 0)
end

MMI.@mlj_model mutable struct LSSVRegressor <: MMI.Deterministic
    kernel::Symbol = :rbf::(_ in (:rbf, :linear, :poly))
    γ::Float64 = 1.0::(_ > 0.0)
    σ::Float64 = 1.0::(_ > 0.0)
    degree::Int = 0::(_ >= 0)
end

##
## Fitting functions
##

function MMI.fit(model::LSSVClassifier, verbosity::Int, X, y)
    Xmatrix = MMI.matrix(X; transpose=true) # notice the transpose

    a_target_element = y[1]
    num_classes = length(MMI.classes(a_target_element))

    decode  = MMI.decoder(a_target_element) # for predict method
    y_plain = convert(Array{eltype(Xmatrix)}, MMI.int(y))

    cache = nothing

    if num_classes == 2 # binary classification
        new_y = broadcast(x -> x == 2.0 ? -1.0 : 1.0, y_plain)

        svm = LSSVC(;kernel=model.kernel, γ=model.γ, σ=model.σ)
        fitted = svmtrain(svm, Xmatrix, new_y)

        fitresult = (deepcopy(svm), fitted, decode)

        report = (kernel = model.kernel, γ = model.γ, σ = model.σ)
    else # multiclass classification
        svm = LSSVC(; kernel=model.kernel, γ=model.γ, σ=model.σ)
        fitted = svmtrain_mc(svm, Xmatrix, y_plain, num_classes)
        fitresult = (fitted, decode)
        report = (kernel = model.kernel, γ = model.γ, σ = model.σ)
    end

    return (fitresult, cache, report)
end

function MMI.fit(model::LSSVRegressor, verbosity::Int, X, y)
    Xmatrix = MMI.matrix(X; transpose=true) # notice the transpose

    cache = nothing

    svr = LSSVR(; kernel=model.kernel, γ=model.γ, σ=model.σ)
    fitted = svmtrain(svr, Xmatrix, y)

    fitresult = (deepcopy(svr), fitted)

    report = (kernel = model.kernel, γ = model.γ, σ = model.σ)

    return (fitresult, cache, report)
end

##
## Prediction functions
##

function MMI.predict(model::LSSVClassifier, fitresult, Xnew)
    Xmatrix = MMI.matrix(Xnew; transpose=true) # notice the transpose
    n_fits = length(fitresult) # number of elements from the fit step

    if n_fits == 3 # binary classification
        (svm, fitted, decode) = fitresult
        results = svmpredict(svm, fitted, Xmatrix)
        broadcast!(x -> x == -1.0 ? 2.0 : 1.0, results, results)
        y = convert(Array{UInt64}, results)
        predictions = decode(y)
    else
        (fitted, decode) = fitresult
        results = svmpredict_mc(fitted, Xmatrix)
        y = convert(Array{UInt64}, results)
        predictions = decode(y)
    end

    return predictions
end

function MMI.predict(model::LSSVRegressor, fitresult, Xnew)
    Xmatrix = MMI.matrix(Xnew; transpose=true) # notice the transpose
    (svr, fitted) = fitresult
    results = svmpredict(svr, fitted, Xmatrix)

    return results
end

##
## Metadata
##

# shared metadata
const LSSVM = (LSSVC, LSSVR)

MMI.metadata_pkg.(LSSVM,
    name="LeastSquaresSVM",
    uuid="6bfd0e71-701c-47cd-9c90-5bf8fe84640d",
    url="https://github.com/edwinb-ai/LeastSquaresSVM",
    julia=true,
    license="MIT",
    is_wrapper=false)

MMI.metadata_model(LSSVC,
    input=MMI.Table(MMI.Continuous),
    target=AbstractVector{MMI.Finite},
    weights=false,
    descr="A Least Squares Support Vector Classifier implementation.",
    path="LeastSquaresSVM.LSSVC",
)

MMI.metadata_model(LSSVR,
    input=MMI.Table(MMI.Continuous),
    target=AbstractVector{MMI.Continuous},
    weights=false,
    descr="A Least Squares Support Vector Regressor implementation.",
    path="LeastSquaresSVM.LSSVR",
)

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
#-------------------------------------------------------------------------------
# This function creates an instance of Mimi-FUND, given user specifications.
#-------------------------------------------------------------------------------

# Load packages.
using Mimi
using DataFrames
using CSVFiles
using MimiFUND

# Load necessary model and data files.
include(joinpath("new_components", "updated_welfare_fund.jl"))

# Create a function to construct an updated version of FUND.
function create_fund(ρ::Float64, η::Float64, welfare_year::Int)

    # ---------------------------------------------
    # Create MimiRICE2010 model and set parameters.
    # ---------------------------------------------

    # Load base version of Mimi-FUND and set time horizon to 1950-2300.
    m = MimiFUND.get_model()
    set_dimension!(m, :time, collect(1950:2300))

    # Add in new welfare component.
    add_comp!(m, fund_welfare, after = :impactaggregation)

    # Set parameters.
    set_param!(m, :fund_welfare, :ρ, ρ)
    set_param!(m, :fund_welfare, :η, η)
    set_param!(m, :fund_welfare, :welfare_year, welfare_year)

    # Create component connections.
    connect_param!(m, :fund_welfare, :pop, :population, :population)
    connect_param!(m, :fund_welfare, :c, :socioeconomic, :consumption)

    # Return user-specified model.
    return m
end

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism


text:
using KernelDensity

@testset "pdf and Gaussian state data generator" begin
    # StateMatrix
    state = SqueezedThermalState(ξ(1., π/4), 0.5)
    θs = LinRange(0, 2π, 10)
    xs = LinRange(-10, 10, 10)

    ground_truth_pdf = q_pdf(state, θs, xs)

    single_point_pdf = (θ, x) -> q_pdf(state, θ, x)
    @test single_point_pdf.(θs, xs') ≈ ground_truth_pdf

    n = 100000
    data = gaussian_state_sampler(state, n)
    sampled_pdf = pdf(kde((data[1, :], data[2, :])), θs, xs)
    @test sum(abs.(sampled_pdf .- ground_truth_pdf)) / n < 5e-5

    # StateVector
    state = VacuumState()
    θs = LinRange(0, 2π, 10)
    xs = LinRange(-10, 10, 10)

    ground_truth_pdf = q_pdf(state, θs, xs)

    single_point_pdf = (θ, x) -> q_pdf(state, θ, x)
    @test single_point_pdf.(θs, xs') ≈ ground_truth_pdf

    n = 100000
    data = gaussian_state_sampler(state, n)
    sampled_pdf = pdf(kde((data[1, :], data[2, :])), θs, xs)
    @test sum(abs.(sampled_pdf .- ground_truth_pdf)) / n < 5e-5
end


text:
module SuperLearn

##############################
## by Qin Yu, Nov 2018
## using Julia 1.0.1
## for Supervised Learning
##############################

using LinearAlgebra # Julia 1.0.1 only
using Plots
using Printf
using Statistics

export TRANS_BASIS,
       POLY_OR_SINE,
       sort_matrix_by_jth_col,
       trained_regression_line,
       trained_regression_line_M,
       test_error_k_dim_basis,
       training_error_k_dim_basis,
       get_se

POLY_OR_SINE = "poly"
TRANS_BASIS = true

# Formulae - Find 𝝎:
function phik(x, k)  # function phik(x, k; basis="poly")  # don't know how to implement this efficiently
    if POLY_OR_SINE == "poly"
        [xi^(k-1) for xi in x]  # ϕₖ(𝒙) = ..., for vector 𝒙, all inputs
    elseif POLY_OR_SINE == "sine"
        [sin(k * π * xi) for xi in x]  # ϕₖ(𝒙) = ..., for vector 𝒙, all inputs
    end
end

phi1tok(x, k) = [phik(x, i) for i in 1:k]  # Φ(𝒙) = 𝒙 ⋅ 𝝋 = 𝒙 ⋅ (ϕ₁(), ϕ₂(), ..., ϕₖ())

function transformed_x_kk(x, k)  # Φ(𝒙) as a matrix
    if TRANS_BASIS
        hcat(phi1tok(x, k)...)
    else
        size(x) == () ? [x, 1] : hcat(x, ones(size(x)[1]))
    end
end

𝒘(x, y, k) = transformed_x_kk(x, k) \ y  # 𝝎 = Φ\𝒚

# Formulae - Equation of Fitted Regression Line:
𝒘Φ_M(X_test, X_train, 𝒚_train, nothing) = transformed_x_kk(X_test, nothing) * 𝒘(X_train, 𝒚_train, nothing)
𝒘Φ(x_test, x_train, y, k) = dot(transformed_x_kk(x_test, k), 𝒘(x_train, y, k))  # ̂𝑦 = Φ(𝑥) ⋅ 𝝎
trained_regression_line = 𝒘Φ
trained_regression_line_M = 𝒘Φ_M

# MSE for Testing:
function test_error_k_dim_basis(𝒙_test, 𝒚_test, 𝒙, 𝒚, k)
    # SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)² = 𝚺ᵢ(𝑦ᵢ - Φ(𝑥ᵢ) ⋅ 𝝎)²
    # MSE = SSE/N, where N = number_of_rows(input_data_set), here is S
    sse = sum((𝒚_test - transformed_x_kk(𝒙_test, k) * 𝒘(𝒙, 𝒚, k)).^2)
    mse = sse / first(size(𝒚_test))
    return mse
end

# MSE for Training:
training_error_k_dim_basis(𝒙, 𝒚, k) = test_error_k_dim_basis(𝒙, 𝒚, 𝒙, 𝒚, k)

# Sort a Matrix by the jth Column:
function sort_matrix_by_jth_col(A, j)
    disassembled_A = [A[i,:] for i in 1:size(A, 1)]
    sort!(disassembled_A, by = x -> x[j])
    vcat(disassembled_A'...)
end

# My function of get standard_error from observations:
function get_se(observations)
    n = size(observations, 1)
    standard_deviation = std(observations)
    standard_error = standard_deviation / sqrt(n)
    return standard_error
end

end

max_stars_repo_name: qin-yu/julia-regression-boston-housing

max_issues_repo_name: qin-yu/ml-julia-boston-housing

max_forks_repo_name: qin-yu/ml-julia-boston-housing


text:
struct SmartMeter
    id::String
    dat::Array{Float32,2}
    tVec::Array{DateTime,1}
    tRes::Second
    nDays::Day
end

mutable struct Project_SMSPS
    path::Dict
    fullDataList::Array{String,1}
    dataIndices
    smartMeters::Array{SmartMeter,1}
    summaryData::DataFrame
    htgSchedules::DataFrame
    verbose::Array{DataFrame,1}
    Project_SMSPS() = new()
end


text:
using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, "..")))

using DataFrames
using Configurations
using HDF5
using Glob
using ArgParse
using Logging
using Printf
using ATools
using NReco

logger = SimpleLogger(stdout, Logging.Warn)
old_logger = global_logger(logger)

function makezoo(args)
	dr      = args["dir"]
	outd    = args["odir"]
	outf    = args["ofile"]
	file_i  = args["filei"]
	file_l  = args["filel"]
	detconf = args["config"]

	if isdir(outd) == false
		mkdir(outd)
	end

	output_path = joinpath(outd, outf)
	filenames   = glob("*.h5",dr)

	if detconf != "default"
		dconf = from_toml(NReco.DetConf, detconf)
	else
		dconf = NReco.DetConf()
	end

	println("number of files in data dir = $(length(filenames))")
	println("reading = $(file_l - file_i + 1) files")
	println("output file  = $output_path")

	h5open(output_path, "w") do h5out
		grp        = create_group(h5out, "EventCounts")
		dtype      = ATools.generate_hdf5_datatype(ATools.EventTypes)

		evt_counts = NReco.event_classifier(filenames, dconf, file_i, file_l)

		dspace     = dataspace([evt_counts])
		dset       = create_dataset(grp, "counts", dtype, dspace)
		write_dataset(dset, dtype, [evt_counts])
	end
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--dir", "-d"
            help = "directory with nema simulations"
            arg_type = String
		"--odir", "-o"
            help = "output directory"
            arg_type = String
		"--ofile", "-x"
            help = "output file"
            arg_type = String
            default = "evtdf.csv"
		"--filei", "-i"
	        help = "number of initial file in glob list"
	        default  = 1
			arg_type = Int
		"--filel", "-l"
		    help = "number of last file in glob list"
		    default  = 1
			arg_type = Int
		"--config", "-c"
			help = "detector configuration"
			default = "default"
			arg_type = String
    end

    return parse_args(s)
end

function main()
	parsed_args = parse_commandline()
	println("Running makezoo with arguments", parsed_args)
	makezoo(parsed_args)
end

@time main()

max_stars_repo_name: andLaing/NReco

max_issues_repo_name: andLaing/NReco

max_forks_repo_name: andLaing/NReco


text:
function r3frac!(epslat::Float64, v)
    for i in 1:3
        v[i] = v[i] - round(Int64, v[i])
        if v[i] < 0.0
            v[i] = v[i] + 1.0
        end
        if (1.0 - v[i]) < epslat
            v[i] = 0.0
        end
        if v[i] < epslat
            v[i] = 0.0
        end
    end
    return
end


text:
module CloudMicrophysics

include("Common.jl")
include("Microphysics_0M.jl")
include("Microphysics_1M.jl")
include("AerosolModel.jl")
include("AerosolActivation.jl")

end # module


text:
using LinearAlgebra, DataFrames, Statistics

n = 1000
x = [ones(n) rand(1:0.1:10, n)]
b,sigma = [2.5, -1.2], 1
y = x*b + sigma .* randn(n);

# Write types
abstract type StandardErrorEstimator end
struct Bootstrap <: StandardErrorEstimator
    bootstrap_samples::Int
    subsample_size::Int
end
struct Spherical <: StandardErrorEstimator end

# Estimator code
ols(x, y) = ols(x, y, Spherical())
function ols(x, y, se_estimator::Spherical)
    V_inv = inv(x'x)
    b_hat = V_inv*x'y
    y_hat = x*b_hat
    e = y - y_hat

    n = length(y)
    k = size(x, 2)
    s2  = e'e / (n-k)

    se=sqrt.(s2 .* diag(V_inv))

    return DataFrame(
        b_hat=b_hat,
        se=se,
    )
end

function ols(x, y, se_estimator::Bootstrap)
    n = length(y)
    n_samples = se_estimator.bootstrap_samples
    subsample_size = se_estimator.subsample_size
    b_hats = zeros(n_samples, size(x, 2))

    for k in 1:n_samples
        inds = rand(1:n, subsample_size)
        x_sub = x[inds,:]
        y_sub = y[inds]

        b_hats[k, :] = inv(x_sub'x_sub)*x_sub'y_sub
    end

    b_hat = vec(mean(b_hats, dims=1))
    se_hat = vec(std(b_hats, dims=1))

    return DataFrame(
        b_hat=b_hat,
        se_hat=se_hat,
    )
end

fit1 = ols(x, y, Spherical())
fit2 = ols(x, y, Bootstrap(100_000, 10))

println("Spherical")
display(fit1)

println("\nBootstrapped")
display(fit2)


text:
using DataFrames

# TYPE DEFINITIONS
type Generator
   index::Any
   bus_idx::Int
   g_P_max::Float64
   g_Q_max::Float64
   cost::Float64
   function Generator(index, bus_idx, g_P_max, g_Q_max)
      g = new()
      g.index  = index
      g.bus_idx = bus_idx
      g.cost = 1.00
      g.g_P_max = g_P_max
      g.g_Q_max = g_Q_max
      return g
   end
end

type Bus
   index::Any
   is_root::Bool
   d_P::Float64
   d_Q::Float64
   cosphi::Float64
   tanphi::Float64
   v_max::Float64
   v_min::Float64
   children::Vector{Int}
   ancestor::Vector{Int}
   generator::Generator
   function Bus(index, d_P, d_Q, v_max, v_min)
      b = new()
      b.index = index
      b.is_root = false
      b.d_P = d_P
      b.d_Q = d_Q
      b.v_max = v_max
      b.v_min = v_min
      b.children = Int[]
      b.ancestor = Int[]
      cosphi = d_P/(sqrt(d_P^2 + d_Q^2))
      tanphi = tan(acos(b.cosphi))
      if isnan(cosphi)
        b.cosphi = 0
        b.tanphi = 0
      else
        b.cosphi = cosphi
        b.tanphi = tan(acos(cosphi))
      end
      return b
   end
end

type Line
   index::Any
   to_node::Int # the "to" node
   from_node::Int # the "from" node
   r::Float64 # the resistance value
   x::Float64 # the reactance value
   b::Float64 # the susceptance value
   s_max::Float64 # the capacity of the line
   function Line(index, to_node, from_node, r, x, s_max)
      l = new()
      l.index = index
      l.to_node = to_node
      l.from_node = from_node
      l.r = r
      l.x = x
      l.b = (x/(r^2 + x^2))
      l.s_max = s_max
      return l
   end
end

function load_case_data(;datafile="")
# READ RAW DATA
info("Reading Data")

if datafile == ""
  data_dir = DATA_DIR
else
  data_dir = datafile
end

nodes_raw = readtable("data/$data_dir/nodes.csv")
sum(nonunique(nodes_raw, :index)) != 0 ? warn("Ambiguous Node Indices") : nothing

lines_raw = readtable("data/$data_dir/lines.csv")
sum(nonunique(lines_raw, :index)) != 0  ? warn("Ambiguous Line Indices") : nothing

generators_raw = readtable("data/$data_dir/generators.csv")
sum(nonunique(generators_raw, :index)) != 0 ? warn("Ambiguous Generator Indices") : nothing

# PREPARING MODEL DATA

buses = Dict()
for n in 1:nrow(nodes_raw)
    index = nodes_raw[n, :index]
    d_P = nodes_raw[n, :d_P]
    d_Q = nodes_raw[n, :d_Q]
    v_max = nodes_raw[n, :v_max]
    v_min = nodes_raw[n, :v_min]
    newb = Bus(index, d_P, d_Q, v_max, v_min)
    buses[newb.index] = newb
end

lines = Dict()
for l in 1:nrow(lines_raw)
    index = lines_raw[l, :index]
    from_node = lines_raw[l, :from_node]
    to_node = lines_raw[l, :to_node]
    r = lines_raw[l, :r]
    x = lines_raw[l, :x]
    b = lines_raw[l, :b]
    s_max = lines_raw[l, :s_max]
    newl = Line(index, to_node, from_node, r, x, s_max)

    push!(buses[newl.from_node].children, newl.to_node)
    push!(buses[newl.to_node].ancestor, newl.from_node)

    lines[newl.index] = newl
end

# Check topology
r = 0
root_bus = 0
for b in keys(buses)
    l = length(buses[b].ancestor)
    if l > 1
        warn("Network not Radial (Bus $(buses[b].index))")
    elseif l == 0
        buses[b].is_root = true
        root_bus = b
        r += 1
    end
end
if r == 0
    warn("No root detected")
    root_bus = 0
elseif r > 1
    warn("More than one root detected")
end

generators = Dict()
for g in 1:nrow(generators_raw)
    index = generators_raw[g, :index]
    bus_idx = generators_raw[g, :node]
    g_P_max = generators_raw[g, :p_max]
    g_Q_max = generators_raw[g, :q_max]
    cost = generators_raw[g, :cost]
    newg = Generator(index, bus_idx, g_P_max, g_Q_max)
    newg.cost = cost

    buses[newg.bus_idx].generator = newg

    generators[newg.index] = newg
end

# Radial PTDF
A = zeros(length(lines), length(buses))
for b in collect(keys(buses))
    a = b
    while a != root_bus
        A[a-1, b] = 1
        a = buses[a].ancestor[1]
    end
end

info("Done preparing Data")
return buses, lines, generators, A
end

function change_load_same_pf(α)
# always uses the GLOBAL BUSES as reference

  buses = deepcopy(BUSES)
  for i in keys(buses)
      P = buses[i].d_P
      Q = buses[i].d_Q
      PF = P/sqrt(P^2 + Q^2)
      if !(isnan(PF))
         buses[i].d_P = P*α
         buses[i].d_Q = P*α*sqrt((1-PF^2)/PF^2)
      end
   end

   return buses
end

function add_generator(index, bus, g_P_max, g_Q_max, cost)
    newg = Generator(index, bus, g_P_max, g_Q_max)
    newg.cost = cost
    BUSES[newg.bus_idx].generator = newg
    GENERATORS[newg.index] = newg
end


text:
using DecFP
using Base.Test

@test unsafe_load(DecFP.flags) == 0

import DecFP.isnanstr
@test isnanstr("nan") && isnanstr("  +NAN") && isnanstr("-NaN") && !isnanstr("nano")

for T in (Dec32, Dec64, Dec128)
    println(STDERR, "TESTING $T ...")

    if T == Dec32
        @test d32"3.2" * d32"4.5" == d32"14.4"
        @test eps(Dec32) == parse(Dec32, "1e-6")
        @test realmax(Dec32) == parse(Dec32, "+9999999E+90")
        @test realmin(Dec32) == parse(Dec32, "+1E-101")
        @test bswap(Dec32(1.5)) == reinterpret(Dec32, bswap(reinterpret(UInt32, Dec32(1.5))))
    elseif T == Dec64
        @test d"3.2" * d"4.5" == d"14.4"
        @test d64"3.2" * d64"4.5" == d64"14.4"
        @test eps(Dec64) == parse(Dec64, "1e-15")
        @test realmax(Dec64) == parse(Dec64, "+9999999999999999E+369")
        @test realmin(Dec64) == parse(Dec64, "+1E-398")
        @test bswap(Dec64(1.5)) == reinterpret(Dec64, bswap(reinterpret(UInt64, Dec64(1.5))))
    else
        @test d128"3.2" * d128"4.5" == d128"14.4"
        @test eps(Dec128) == parse(Dec128, "1e-33")
        @test realmax(Dec128) == parse(Dec128, "+9999999999999999999999999999999999E+6111")
        @test realmin(Dec128) == parse(Dec128, "+1E-6176")
        @test bswap(Dec128(1.5)) == reinterpret(Dec128, bswap(reinterpret(UInt128, Dec128(1.5))))
    end

    @test parse(T, "1.5")::T == 1.5
    @test isnan(parse(T, "NaN")::T) && isnan(T(NaN))
    @test isinf(parse(T, "inf")::T) && !isfinite(parse(T, "inf")::T)
    @test parse(T,"inf")::T == typemax(T)::T
    @test parse(T,"-inf")::T == typemin(T)::T

    x,y,z = 1.5, -3.25, 0.0625 # exactly represented in binary
    xd = T(x); yd = T(y); zd = T(z)

    @test fma(xd,yd,zd)::T == muladd(xd,yd,zd)::T == xd*yd+zd == fma(x,y,z)

    @test one(T)::T == 1
    @test zero(T)::T == 0
    @test signbit(xd) == signbit(x) && signbit(yd) == signbit(y) && signbit(parse(T,"-Inf"))
    @test sign(xd)::T == sign(x) && sign(yd)::T == sign(y)

    @test nextfloat(xd)::T == xd + eps(T)
    @test prevfloat(xd)::T == xd - eps(T)
    @test nextfloat(T(1.5e10)) == 1.5e10 + eps(T(1.5e10))
    @test prevfloat(T(1.5e10)) == 1.5e10 - eps(T(1.5e10))

    for f in (isnan,isinf,isfinite,issubnormal,abs)
        @test f(xd) == f(x)
    end

    for f in (+,-,*,copysign,==,>,>=,<,<=)
        @test f(xd,yd) == f(x,y)
    end

    for f in (/,hypot,atan2,^)
        @test_approx_eq f(xd,yd) f(x,y)
        if f != ^
            @test_approx_eq f(yd,xd) f(y,x)
            @test_approx_eq f(yd,yd) f(y,y)
        end
    end

    for f in (exp,log,sin,cos,tan,asin,acos,atan,sinh,cosh,tanh,asinh,acosh,atanh,log1p,expm1,log10,log2,exp2,exp10,erf,erfc,lgamma,sqrt,cbrt)
        v = f == acosh ? x : z
        @test_approx_eq f(T(v)) f(v)
    end

    for c in (π, e, γ, catalan, φ)
        @test_approx_eq T(c) Float64(c)
    end

    for Tf in (Float64, Float32, Float16)
        @test xd == Tf(x) == T(Tf(x)) == Tf(xd)
    end
    for Ti in (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64)
        @test parse(T, "17") == T(Ti(17)) == Ti(17) == Ti(T(17))
        @test floor(Ti, T(4.5)) == 4 == ceil(Ti, T(4.5)) - 1
        @test_throws InexactError convert(Ti, xd)
        @test_throws InexactError floor(Ti, realmax(T))
        @test_throws InexactError ceil(Ti, realmax(T))
        if Ti <: Signed
            @test parse(T, "-17") == T(Ti(-17)) == Ti(-17) == Ti(T(-17))
            @test floor(Ti, T(-4.5)) == -5 == ceil(Ti, T(-4.5)) - 1
            @test_throws InexactError convert(Ti, yd)
        end
    end

    TI = eval(Symbol(string("UInt", sizeof(T)*8)))
    @test bswap(xd) == reinterpret(T, bswap(reinterpret(TI, xd)))

    @test_throws InexactError parse(T, "1e10000")
    @test_throws DomainError asin(xd)
    @test_throws DomainError sqrt(yd)
    @test_throws DomainError acosh(zd)

    # @test ldexp(parse(T, "1"), 3) == 1000
    # @test exponent(parse(T, "1000")) == 3
    @test_approx_eq sqrt(complex(yd)) sqrt(complex(y))

    @test typeof(xd * pi) == T
    @test typeof((xd+yd*im)*pi) == Complex{T}
end

@test unsafe_load(DecFP.flags) == 0


text:
@testset "Spatial data" begin
  # basic checks
  @test coordinates(data3D) == Dict(var => Float64 for var in [:x,:y,:z])
  @test variables(data3D) == Dict(:value => Float64)
  @test npoints(data3D) == 100
  X, z = valid(data3D, :value)
  @test size(X,2) == 100
  @test length(z) == 100
  @test_throws AssertionError readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t', coordnames=[:a])

  # missing data and NaN
  X, z = valid(missdata, :value)
  @test size(X) == (2,1)
  @test length(z) == 1

  # show methods
  rawdata = DataFrame(x=[1,2,3],y=[4,5,6])
  geodata = GeoDataFrame(rawdata, [:x,:y])
  @test sprint(show, geodata) == "3×2 GeoDataFrame (x and y)"
  @test sprint(show, MIME"text/plain"(), geodata) == "3×2 GeoDataFrame (x and y)\n\n│ Row │ x     │ y     │\n│     │ Int64 │ Int64 │\n├─────┼───────┼───────┤\n│ 1   │ 1     │ 4     │\n│ 2   │ 2     │ 5     │\n│ 3   │ 3     │ 6     │"
  @test sprint(show, MIME"text/html"(), geodata) == "3×2 GeoDataFrame (x and y)\n

	x	y

	Int64	Int64

3 rows × 2 columns

1	1	4

2	2	5

3	3	6

"
end


text:
## 11.jl
# タブをスペースに変換

f = open("./dat/hightemp.txt")
o = open("./dat/11.txt", "w")
print(o, replace(readall(f), "\t", " "))
close(f)
close(o)


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic

"""
    ShowExampleAppMainMenuBar(p_open::Ref{Bool})
Create a fullscreen menu bar and populating it.
"""
function ShowExampleAppMainMenuBar(p_open::Ref{Bool})
    if CImGui.BeginMainMenuBar()
        if CImGui.BeginMenu("File")
            ShowExampleMenuFile()
            CImGui.EndMenu()
        end
        if CImGui.BeginMenu("Edit")
            if CImGui.MenuItem("Undo", "CTRL+Z")
                @info "Trigger Undo | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            if CImGui.MenuItem("Redo", "CTRL+Y", false, false)  # disabled item
                @info "Trigger Redo | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            CImGui.Separator()
            if CImGui.MenuItem("Cut", "CTRL+X")
                @info "Trigger Cut | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            if CImGui.MenuItem("Copy", "CTRL+C")
                @info "Trigger Copy | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            if CImGui.MenuItem("Paste", "CTRL+V")
                @info "Trigger Paste | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            CImGui.EndMenu()
        end
        CImGui.EndMainMenuBar()
    end
end

function ShowExampleMenuFile()
    CImGui.MenuItem("(dummy menu)", C_NULL, false, false)
    if CImGui.MenuItem("New")
        @info "Trigger New | find me here: $(@__FILE__) at line $(@__LINE__)"
    end
    if CImGui.MenuItem("Open", "Ctrl+O")
        @info "Trigger Open | find me here: $(@__FILE__) at line $(@__LINE__)"
    end
    if CImGui.BeginMenu("Open Recent")
        CImGui.MenuItem("fish_hat.c")
        CImGui.MenuItem("fish_hat.inl")
        CImGui.MenuItem("fish_hat.h")
        if CImGui.BeginMenu("More..")
            CImGui.MenuItem("Hello")
            CImGui.MenuItem("Sailor")
                if CImGui.BeginMenu("Recurse..")
                    ShowExampleMenuFile()
                    CImGui.EndMenu()
                end
            CImGui.EndMenu()
        end
        CImGui.EndMenu()
    end
    if CImGui.MenuItem("Save", "Ctrl+S")
        @info "Trigger Save | find me here: $(@__FILE__) at line $(@__LINE__)"
    end
    if CImGui.MenuItem("Save As..")
        @info "Trigger Save As.. | find me here: $(@__FILE__) at line $(@__LINE__)"
    end
    CImGui.Separator()
    if CImGui.BeginMenu("Options")
        @cstatic enabled=true begin
            @c CImGui.MenuItem("Enabled", "", &enabled)
        end
        CImGui.BeginChild("child", ImVec2(0, 60), true)
            foreach(i->CImGui.Text("Scrolling Text $i"), 0:9)
        CImGui.EndChild()

        @cstatic f=Cfloat(0.5) n=Cint(0) b=true begin
            @c CImGui.SliderFloat("Value", &f, 0.0, 1.0)
            @c CImGui.InputFloat("Input", &f, 0.1)
            @c CImGui.Combo("Combo", &n, "Yes\0No\0Maybe\0\0")
            @c CImGui.Checkbox("Check", &b)
        end
        CImGui.EndMenu()
    end
    if CImGui.BeginMenu("Colors")
        sz = CImGui.GetTextLineHeight()
        for i = 0:CImGui.ImGuiCol_COUNT-1
            name = CImGui.GetStyleColorName(i)
            p = CImGui.GetCursorScreenPos()
            CImGui.AddRectFilled(CImGui.GetWindowDrawList(), p, (p.x+sz,p.y+sz), CImGui.GetColorU32(i))
            CImGui.Dummy(sz, sz)
            CImGui.SameLine()
            CImGui.MenuItem(name)
        end
        CImGui.EndMenu()
    end
    if CImGui.BeginMenu("Disabled", false)  # disabled
        throw(AssertionError("unreachable reached."))
    end
    CImGui.MenuItem("Checked", "", true) && @info "Trigger Checked | find me here: $(@__FILE__) at line $(@__LINE__)"
    CImGui.MenuItem("Quit", "Alt+F4") && @info "Trigger Quit | find me here: $(@__FILE__) at line $(@__LINE__)"
end


text:
using SparseRelaxation
using Base.Test

function conv_test(smoother)
    # test zero input
    @test all(x -> x == 0, smoother(speye(100), zeros(100), zeros(100)))

    # test convergence
    function converges()
        A = spdiagm(ones(99), 1, 100, 100)
        A = A + transpose(A)
        L = spdiagm(vec(sum(A, 1))) - A
        x = zeros(size(L, 1))
        b = rand(size(L, 1))
        b = b - mean(b)
        x1 = smoother(L,x,b)
        norm(b - L*x1) < norm(b - L*x)
    end

    # check convergence
    @test converges()
end

@testset "gauss_seidel forwards" begin
    conv_test((A,x,b)->gauss_seidel(A,x,b,backwards=false))
end
@testset "gauss_seidel backwards" begin
    conv_test((A,x,b)->gauss_seidel(A,x,b,backwards=true))
end
@testset "weighted_jacobi" begin
    conv_test(weighted_jacobi)
end