text:
# Unit tests to check that the implemented merge strategy from CHOMPACK works as expected
using COSMO, SparseArrays, LinearAlgebra, Test

function get_example_tree(strategy)

  # define example tree with known merge order
  snd1 = Set([15; 16; 17]); sep1 = Set{Int64}();
  snd2 = Set([5; 9]); sep2 = Set([15; 16]);
  snd3 = Set([3; 4]); sep3 = Set([5; 15]);
  snd4 = Set([1]); sep4 = Set([3]);
  snd5 = Set([2]); sep5 = Set([3; 4]);
  snd6 = Set([6]); sep6 = Set([9; 16]);
  snd7 = Set([7; 8]); sep7 = Set([9; 15]);
  snd8 = Set([12; 13; 14]); sep8 = Set([16; 17]);
  snd9 = Set([10; 11]); sep9 = Set([13; 14; 17]);

  supernodes = [snd1, snd2, snd3, snd4, snd5, snd6, snd7, snd8, snd9]
  seperators = [sep1, sep2, sep3, sep4, sep5, sep6, sep7, sep8, sep9]
  parents = [0; 1; 2; 3; 3; 2; 2; 1; 8]
  snd_post = collect(9:-1:1)
  post = collect(1:1:17)
  return COSMO.SuperNodeTree(supernodes, parents, snd_post, seperators, strategy, post = post)
end

function get_example_graph(strategy)
  t = get_example_tree(strategy)
  @. t.snd_par = -1
  @. empty!(t.snd_child)

  # add separators to snds
  for (i, snd) in enumerate(t.snd)
    union!(snd, t.sep[i])
    empty!(t.sep[i])
  end
  return t
end

function complexity_savings(t, c1, c2)
  dim_c1 = length(t.snd[c1])
  dim_c2 = length(t.snd[c2])
  dim_u = length(union(t.snd[c1], t.snd[c2]))
  return dim_c1^3 + dim_c2^3 - dim_u^3
end

# Correct solution for adjacency matrix with edge scores
function get_adjacency_matrix(t)
  rows = [2;  8; #1
  3; 6; 7; #2
  4; 5; #3
   5; #4
  9; #8
  ]
  cols = [ones(Int64, 2); 2*ones(Int64, 3); 3*ones(Int64, 2);  4; 8]
  N = length(rows)
  edge_weight = zeros(N)
  for iii = 1:N
    edge_weight[iii] = complexity_savings(t, rows[iii], cols[iii])
  end
  A = sparse(rows, cols, edge_weight, 9, 9)
  return A
end

# correct solution for clique graph with intersection weights
function get_intersection_matrix()
  rows = [2; 3; 6; 7; 8; 9; 3; 6; 7; 8; 4; 5; 7; 5; 7; 8; 9]
  cols = [ones(Int64, 6); 2 * ones(Int64, 4); 3 * ones(Int64, 3); 4; 6*ones(Int64, 2); 8]
  vals = [2.0; 1; 1; 1; 2; 1; 2; 2; 2; 1; 1; 2; 1; 1; 1; 1; 3]
  A = sparse(rows, cols, vals, 9, 9)
end
@testset "Clique merging Example" begin

  @testset "Clique merging (Parent-child merge strategy)" begin

    strategy = COSMO.ParentChildMerge()
    t = get_example_tree(strategy)
    COSMO.initialise!(t, strategy)

    # Check first  merge_cand: (1, 2)
    cand = COSMO.traverse(t, strategy)
    dim_par_snd, dim_par_sep = COSMO.clique_dim(t, cand[1])
    dim_clique_snd, dim_clique_sep = COSMO.clique_dim(t, cand[2])
    @test cand == [1; 2] && dim_clique_snd == 2 && dim_clique_sep == 2 && dim_par_snd == 3 && dim_par_sep == 0
    f_i = COSMO.fill_in(dim_clique_snd, dim_clique_sep, dim_par_snd, dim_par_sep)
    s = COSMO.max_snd_size(dim_clique_snd, dim_par_snd)
    @test f_i == 2 && s == 3
    @test COSMO.evaluate(t, strategy, cand)
    ordered_cand = COSMO.merge_two_cliques!(t, cand, strategy)
    @test isempty(t.snd[2]) && t.snd[1] == Set([15; 16; 17; 5; 9])
    @test isempty(t.sep[2]) && isempty(t.sep[1])
    @test t.num == 8 && t.snd_par[2] == -1 && t.snd_par[3] == t.snd_par[6]  == t.snd_par[7] == 1
    @test t.snd_child[1] == Set([8; 3; 6; 7]) && isempty(t.snd_child[2])

    # Let's go throught the whole tree and check if the merge log is as expected
    strategy = COSMO.ParentChildMerge()
    t = get_example_tree(strategy)
    COSMO.merge_cliques!(t, strategy)
    # considered clique pairs
    @test t.merge_log.clique_pairs == [1 2; 1 3; 1 4; 1 5; 1 6; 1 7; 1 8; 8 9]
    # decision if merge should be done
    @test t.merge_log.decisions == [true; true; true; true; true; false; false; true]
    # number of merges
    @test t.merge_log.num == 6

  end

  @testset "Clique merging (Clique graph based merge strategy)" begin

    strategy = COSMO.CliqueGraphMerge(edge_weight = COSMO.ComplexityWeight())
    t = get_example_tree(strategy)
    t.snd = union.(t.snd, t.sep)

    COSMO.initialise!(t, t.strategy)

    # since all values are < 0, this merge strategy wouldn't merge at all
    @test strategy.edges == get_adjacency_matrix(t)

    t = get_example_tree(strategy)
    t.snd = union.(t.snd, t.sep)
    COSMO.merge_cliques!(t, strategy)
    # number of merges
    @test t.merge_log.num == 0
    @test t.snd_par == [0; 1; 2; 3; 3; 2; 2; 1; 8]
  end
end


text:
using Test

include("sieve.jl")

@testset "no primes under two" begin
    @test sieve(1) == []
end

@testset "find first prime" begin
    @test sieve(2) == [2]
end

@testset "find primes up to 10" begin
    @test sieve(10) == [2, 3, 5, 7]
end

@testset "limit is prime" begin
    @test sieve(13) == [2, 3, 5, 7, 11, 13]
end

@testset "find primes up to 1000" begin
    @test sieve(1000) ==  [
        2,   3,   5,   7,   11,  13,  17,  19,  23,  29,  31,  37,  41,
        43,  47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97,  101,
        103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
        173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239,
        241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313,
        317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
        401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
        479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569,
        571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643,
        647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
        739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823,
        827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
        919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997
      ]
end

ext: jl


text:
using BinaryBuilder, Pkg

# GCC version for building the whole system
gcc_version = v"6"

# Versions of various Coin-OR libraries
Cbc_version = v"2.10.3"
Cbc_gitsha = "6fe3addaa76436d479d4431add67b371e11d3e83"

Cgl_version = v"0.60.2"
Cgl_gitsha = "6377b88754fafacf24baac28bb27c0623cc14457"
Cgl_packagespec = PackageSpec(; name = "Cgl_jll",
                              uuid = "3830e938-1dd0-5f3e-8b8e-b3ee43226782",
                              version = Cgl_version)

Clp_version = v"1.17.6"
Clp_gitsha = "756ddd3ed813eb1fa8b2d1b4fe813e6a4d7aa1eb"
Clp_packagespec = PackageSpec(; name = "Clp_jll",
                              uuid = "06985876-5285-5a41-9fcb-8948a742cc53",
                              version = Clp_version)

Osi_version = v"0.108.5"
Osi_gitsha = "2bd34ae6b8c93d342d54fd19d4d773f07194583c"
Osi_packagespec = PackageSpec(; name = "Osi_jll",
                              uuid = "7da25872-d9ce-5375-a4d3-7a845f58efdd",
                              version = Osi_version)

CoinUtils_version = v"2.11.3"
CoinUtils_gitsha = "ea66474879246f299e977802c94a0e45334e7afb"
CoinUtils_packagespec = PackageSpec(; name = "CoinUtils_jll",
                                    uuid = "be027038-0da8-5614-b30d-e42594cb92df",
                                    version = CoinUtils_version)

MUMPS_seq_version = v"4.10.0"
MUMPS_seq_packagespec = PackageSpec(; name = "MUMPS_seq_jll",
                                uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d",
                                version = MUMPS_seq_version)

METIS_version = v"4.0.3"
METIS_packagespec = PackageSpec(; name = "METIS_jll",
                                uuid = "d00139f3-1899-568f-a2f0-47f597d42d70",
                                version = METIS_version)

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())
platforms = [p for p in platforms if !(typeof(p) <: FreeBSD)]


text:
#=
Model: projectile.jl
=#
using Reachability, MathematicalSystems

function compute(input_options::Pair{Symbol,<:Any}...)
    # =====================
    # Problem specification
    # =====================
    # dynamics
    A = [0. 0.5 0. 0. ; 0. 0. 0. 0. ; 0. 0. 0. 0.7 ; 0. 0. 0. 0.]
    b = [0., 0., 0., -9.81]

    # initial set
    X0 = Singleton([0.,5.,100.,0])

    # instantiate continuous LTI system
    S = IVP(AffineContinuousSystem(A, b), X0)

    # ===============
    # Problem solving
    # ===============

    # define solver-specific options
    options = merge(Options(
        :mode => "reach",
        :blocks => [1, 2],
        :plot_vars => [1, 3]
        ), Options(input_options...))

    result = solve(S, options)

    # ========
    # Plotting
    # ========
    if options[:mode] == "reach"
        println("Plotting...")
        tic()
        plot(result)
        @eval(savefig(@filename_to_png))
        toc()
    end
end # function

compute(:N => 10, :T => 20.0); # warm-up
compute(:δ => 0.5, :T => 20.0); # benchmark settings (long)

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
function part1(input)
    return maximum(decode.(readlines(input)))
end

function part2(input)
    seats = decode.(readlines(input))
    return only(setdiff(minimum(seats):maximum(seats), seats))
end

function decode(s)
    for pair in ["B" => "1", "R" => "1", "F" => "0", "L" => "0"]
        s = replace(s, pair)
    end
    return parse(Int, s, base = 2)
end


text:
export DirectedPath, DirectedCycle, DirectedComplete
export RandomDigraph, RandomTournament, TorusDigraph

"""
`DirectedPath(n)` creates a directed cycles with vertices `1:n`.
"""
function DirectedPath(n::Int)
    if n < 1
        error("n must be positive")
    end
    G = IntDigraph(n)
    for u=1:(n-1)
        add!(G,u,u+1)
    end
    return G
end

"""
`DirectedCycle(n)` creates a directed cycles with vertices `1:n`.
"""
function DirectedCycle(n::Int)
    G = DirectedPath(n)
    add!(G,n,1)
    return G
end

# Create a complete digraph (all possible edges)
"""
`DirectedComplete(n)` creates a directed complete graph with
all possible edges (including a loop at each vertex). Use
`DirectedComplete(n,false)` to supress the creation of loops.
"""
function DirectedComplete(n::Int, with_loops::Bool=true)
    G = IntDigraph(n)
    if !with_loops
        forbid_loops!(G)
    end
    for u=1:n
        for v=1:n
            add!(G,u,v)
        end
    end
    return G
end

# Create a random digraph (Erdos-Renyi style)
"""
`RandomDigraph(n,p)` creates an Erdos-Renyi style random directed
graph with vertices `1:n` and edge probability `p` (equal to 0.5 by
default). The possible edges `(u,v)` and `(v,u)` are independent. No
loops are created. To also create loops (each with probability `p`)
use `RandomDigraph(n,p,true)`.
"""
function RandomDigraph(n::Int, p::Real=0.5, with_loops=false)
    G = IntDigraph(n)
    if !with_loops
        forbid_loops!(G)
    end
    for u=1:n
        for v=1:n
            if rand() < p
                add!(G,u,v)
            end
        end
    end
    return G
end

# Create a random tournament (no loops!)

"""
`RandomTournament(n)` creates a random tournament with vertex set
`1:n`.  This is equivalent to randomly assigning a direction to every
edge of a simple complete graph.
"""
function RandomTournament(n::Int)
    G = IntDigraph()
    for u=1:n-1
        for v=u+1:n
            if rand() < 0.5
                add!(G,u,v)
            else
                add!(G,v,u)
            end
        end
    end
    return G
end

"""
`all_tuples(alphabet, n)` creates an iterator that produces all
length-`n` tuples of distinct elements in `alphabet`.
"""
function all_tuples(alphabet, n::Int)
    elts = collect(distinct(alphabet))
    src  = [ elts for _=1:n ]
    its  = Base.Iterators.product(src...)
    return its
end

export ShiftDigraph

"""
`ShiftDigraph(alphabet,n)` creates a `SimpleDigraph` whose vertices
are all length-`n` tuples of the elements in `alphabet` (which can be
an array such as `[0,1]` or a string such as `"abc"`). An edge from
`u` to `v` corresponds to an element dropped from the first position
in `u` and another element added to the end yielding `v`. For example,
in `ShiftDigraph([0,1],5)` there are two edges leaving vertex
`(0,1,0,1,1)`; one goes to `(1,0,1,1,0)` and the other to
`(1,0,1,1,1)`.
"""
function ShiftDigraph(alphabet=[0,1], n::Int=3)
    elts = collect(distinct(alphabet))
    vertex_iter = all_tuples(alphabet, n)
    vlist = collect(vertex_iter)
    T = typeof(vlist[1])
    G = SimpleDigraph{T}()
    for v in vlist
        add!(G,v)
    end

    # create edges here

    for v in vlist
        head = collect(IterTools.drop(v,1))
        for c in elts
            w = tuple([head;c]...)
            add!(G,v,w)
        end
    end

    return G

end

"""
function for creating a Torus Graph
"""
function TorusDigraph(n::Int=4, m::Int=3)

  G = SimpleDigraph();

  #create vertices
  vlist = Tuple{Int,Int}[]
  for i = 1:m
    for j = 1:n
      push!(vlist,(i,j))
    end
  end

  #create edges
  for v in vlist
    if v[1] + 1 <= m
      w = (v[1]+1,v[2])
      add!(G,v,w)
    end
    if v[2] + 1 <= n
      w = (v[1],v[2]+1)
      add!(G,v,w)
    end
    if v[1] == m
      w = (1,v[2])
      add!(G,v,w)
    end
    if v[2] == n
      w = (v[1],1)
      add!(G,v,w)
    end
  end

  return G
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EmpiricalHistogram(sdata, var, [s]; kwargs...)

Spatial histogram of variable `var` in spatial data `sdata`.
Optionally, specify the block side `s` and the keyword
arguments `kwargs` for the `fit(Histogram, ...)` call.
"""
struct EmpiricalHistogram{H}
  hist::H
end

EmpiricalHistogram(d, v::Symbol, w::WeightingMethod; kwargs...) =
  EmpiricalHistogram(fit(Histogram, d[v], weight(d, w); kwargs...))
EmpiricalHistogram(d, v::Symbol, s::Number; kwargs...) =
  EmpiricalHistogram(d, v, BlockWeighting(ntuple(i->s,embeddim(d))); kwargs...)
EmpiricalHistogram(d, v::Symbol; kwargs...) =
  EmpiricalHistogram(d, v, median_heuristic(d); kwargs...)
EmpiricalHistogram(d, w::WeightingMethod; kwargs...) =
  Dict(v => EmpiricalHistogram(d, v, w; kwargs...) for v in name.(variables(d)))
EmpiricalHistogram(d, s::Number; kwargs...) =
  EmpiricalHistogram(d, BlockWeighting(ntuple(i->s,embeddim(d))); kwargs...)

"""
    values(histogram)

Return the abscissa and ordinates of empirical `histogram`.
"""
Base.values(h::EmpiricalHistogram) = midpoints(first(h.hist.edges)), h.hist.weights


text:
using EconPDEs

Base.@kwdef  struct CampbellCochraneModel
    # consumption process parameters
    μ::Float64 = 0.0189
    σ::Float64 = 0.015

    # utility
    γ::Float64 = 2.0
    ρ::Float64 = 0.116

    # habit
    κs::Float64 = 0.138
    b::Float64 = 0.0
end
# I choose persistence so that monthly simulation of the model matches processes in CC (1999)
# ρ = 12 * (1 - 0.89^(1/12))
# κs = 12 * (1 - 0.87^(1/12))

function initialize_stategrid(m::CampbellCochraneModel; smin = -300.0, n = 1000)
    μ = m.μ ; σ = m.σ ; γ = m.γ ; ρ = m.ρ ; κs = m.κs ; b = m.b
    Sbar = σ * sqrt(γ / (κs - b / γ))
    sbar = log.(Sbar)
    smax =  sbar + 0.5 * (1 - Sbar^2)
    # corresponds to Grid 3 in Wachter (2005)
    shigh = log.(range(0.0, stop = exp(smax), length = div(n, 10)))
    slow = range(smin, stop = shigh[2], length = n - div(n, 10))
    OrderedDict(:s => vcat(slow[1:(end-1)], shigh[2:end]))
end

function initialize_y(m::CampbellCochraneModel, stategrid)
    OrderedDict(:p => ones(length(stategrid[:s])))
end

function (m::CampbellCochraneModel)(state::NamedTuple, y::NamedTuple)
    μ = m.μ ; σ = m.σ ; γ = m.γ ; ρ = m.ρ ; κs = m.κs ; b = m.b
    s = state.s
    p, ps, pss = y.p, y.ps, y.pss

    # drift and volatility of  s and p
    Sbar = σ * sqrt(γ / (κs - b / γ))
    sbar = log(Sbar)
    λ = 1 / Sbar * sqrt(1 - 2 * (s - sbar)) - 1
    μs = - κs * (s - sbar)
    σs = λ * σ
    σp = ps / p * σs
    μp = ps / p * μs + 0.5 * pss / p * σs^2

    # market price of risk κ
    κ = γ * (σ + σs)

    # risk free rate  r
    r = ρ + γ * μ - (γ * κs - b) / 2 + b * (sbar - s)

    # PDE
    pt = p * (1 / p + μ + μp + σp * σ - r - κ * (σ + σp))
    return (pt,), (μs, ), (p = p, κ = κ, λ = λ, r = r, σp = σp, μs = μs, σs = σs)
end

# Campbell Cochrane (1999)
m = CampbellCochraneModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0)

# Wachter (2005) calibration
# m = CampbellCochraneModel(μ = 0.022, σ = 0.0086, γ = 2.0, ρ = 0.073, κs = 0.116, b = 0.011)
# stategrid = initialize_stategrid(m)
# y0 = initialize_y(m, stategrid)
# y, result, distance = pdesolve(m, stategrid, y0)


text:
module Surrogates
using LinearAlgebra
using Distributions
using Requires

abstract type AbstractSurrogate <: Function end
include("utils.jl")
include("Radials.jl")
include("Kriging.jl")
include("Sampling.jl")
include("Optimization.jl")
include("Lobachevsky.jl")
include("LinearSurrogate.jl")
include("InverseDistanceSurrogate.jl")
include("SecondOrderPolynomialSurrogate.jl")

function __init__()
    @require Stheno="8188c328-b5d6-583d-959b-9690869a5511" begin
        include("SthenoKriging.jl")
    end
end

include("RandomForestSurrogate.jl")
include("NeuralSurrogate.jl")
include("Wendland.jl")
include("MOE.jl") #rewrite gaussian mixture with own algorithm to fix deps issue
include("VariableFidelity.jl")
include("PolynomialChaos.jl")
include("Earth.jl")
include("GEK.jl")

current_surrogates = ["Kriging","LinearSurrogate","LobachevskySurrogate","NeuralSurrogate",
                      "RadialBasis","RandomForestSurrogate","SecondOrderPolynomialSurrogate",
                      "Wendland","GEK","PolynomialChaosSurrogate"]

#Radial structure:
function RadialBasisStructure(;radial_function,scale_factor,sparse)
    return (name = "RadialBasis", radial_function = radial_function, scale_factor = scale_factor, sparse = sparse)
end

#Kriging structure:
function KrigingStructure(;p,theta)
    return (name = "Kriging", p = p, theta = theta)
end

function GEKStructure(;p,theta)
    return (name = "GEK", p = p, theta = theta)
end

#Linear structure
function LinearStructure()
    return (name = "LinearSurrogate")
end

#InverseDistance structure
function InverseDistanceStructure(;p)
    return (name = "InverseDistanceSurrogate", p = p)
end

#Lobachevsky structure
function LobachevskyStructure(;alpha,n,sparse)
    return (name = "LobachevskySurrogate", alpha = alpha, n = n, sparse = sparse)
end

#Neural structure
function NeuralStructure(;model,loss,opt,n_echos)
    return (name ="NeuralSurrogate", model = model ,loss = loss,opt = opt,n_echos = n_echos)
end

#Random forest structure
function RandomForestStructure(;num_round)
    return (name = "RandomForestSurrogate", num_round = num_round)
end

#Second order poly structure
function SecondOrderPolynomialStructure()
    return (name = "SecondOrderPolynomialSurrogate")
end

#Wendland structure
function WendlandStructure(; eps, maxiters, tol)
    return (name = "Wendland", eps = eps, maxiters = maxiters, tol = tol)
end

#Polychaos structure
function PolyChaosStructure(; op)
    return (name = "PolynomialChaosSurrogate", op = op)
end

export current_surrogates
export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure
export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure
export WendlandStructure
export AbstractSurrogate, SamplingAlgorithm
export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial
export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, LowDiscrepancySample
export RandomSample, KroneckerSample, GoldenSample
export SRBF,LCBS,EI,DYCORS,SOP,EGO,RTEA,SMB,surrogate_optimize
export LobachevskySurrogate, lobachevsky_integral, lobachevsky_integrate_dimension
export LinearSurrogate
export RandomForestSurrogate
export SVMSurrogate
export NeuralSurrogate
export InverseDistanceSurrogate
export SecondOrderPolynomialSurrogate
export SthenoKriging
export Wendland
export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure
export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure
export WendlandStructure
#export MOE
export VariableFidelitySurrogate
export PolynomialChaosSurrogate
export EarthSurrogate
export GEK
end


text:
#####
##### Smolyak implementation details
#####

####
#### Block sizes and shuffling
####

"""
$(SIGNATURES)

Length of each block `b`.

!!! note
    Smolyak grids use “blocks” of polynomials, each indexed by ``b == 0, …, B`, with an
    increasing number of points in each.
"""
@inline __block_length(b::Int) = b ≤ 1 ? b + 1 : 1 << (b - 1)

"""
$(SIGNATURES)

Cumulative block length at block `b`.
"""
@inline __cumulative_block_length(b::Int) = b == 0 ? 1 : ((1 << b) + 1)

"""
$(TYPEDEF)

An iterator of indices for picking elements from a grid of length `len`, which should be a
valid cumulative block length.
"""
struct SmolyakGridShuffle
    len::Int
end

Base.length(ι::SmolyakGridShuffle) = ι.len

Base.eltype(::Type{SmolyakGridShuffle}) = Int

function Base.iterate(ι::SmolyakGridShuffle)
    @unpack len = ι
    i = (len + 1) ÷ 2
    i, (0, 0)                   # step = 0 is special-cased
end

function Base.iterate(ι::SmolyakGridShuffle, (i, step))
    @unpack len = ι
    i == 0 && return len > 1 ? (1, (1, len - 1)) : nothing
    i′ = i + step
    if i′ ≤ len
        i′, (i′, step)
    else
        step′ = step ÷ 2
        if step′ ≥ 2
            i′ = step′ ÷ 2 + 1
            i′, (i′, step′)
        else
            nothing
        end
    end
end

####
#### index traversal
####

"""
$(SIGNATURES)

Internal implementation of the Smolyak indexing iterator.

# Arguments

- `M`: upper bound on each block index, remains constant during the iteration.

- `slack`: `B - sum(blocks)`, cached

- `indices`: current indices

- `blocks`: block indexes

- `limits`: `__cumulative_block_length.(blocks)`, cached

# Return values

- `valid::Bool`: `false` iff there is no next element, in which case the following values
  should be ignored

- `Δ::Int`: change in `slack`

- `indices′`, `blocks′, `limits′`: next values for corresponding arguments above, each an
  `::NTuple{N,Int}`
"""
@inline function __inc(cumulative_block_lengths::NTuple{M,Int}, slack::Int,
                       indices::NTuple{N,Int}, blocks::NTuple{N,Int},
                       limits::NTuple{N,Int}) where {M,N}
    i1, iτ... = indices
    b1, bτ... = blocks
    l1, lτ... = limits
    if i1 < l1                  # increment i1, same block
        true, 0, (i1 + 1, iτ...), blocks, limits
    elseif b1 < M && slack > 0  # increment i1, next block
        b1′ = b1 + 1
        true, -1, (i1 + 1, iτ...), (b1′, bτ...), (cumulative_block_lengths[b1′], lτ...)
    else
        if N == 1               # end of iteration, arbitrary value since !valid
            false, 0, indices, blocks, limits
        else                    # i1 = 1, increment tail if applicable
            Δ1 = b1
            valid, Δτ, iτ′, bτ′, lτ′ = __inc(cumulative_block_lengths, slack + Δ1, iτ, bτ, lτ)
            valid, Δ1 + Δτ, (1, iτ′...), (0, bτ′...), (__cumulative_block_length(0), lτ′...)
        end
    end
end

function __inc_init(::Val{N}, ::Val{B}) where {N,B}
    indices = ntuple(_ -> 1, Val(N))
    blocks = ntuple(_ -> 0, Val(N))
    l = __cumulative_block_length(0)
    limits = ntuple(_ -> l, Val(N))
    slack = B
    slack, indices, blocks, limits
end

"""
$(SIGNATURES)

Calculate the length of a [`SmolyakIndices`](@ref) iterator. Argument as in the latter.
"""
function __smolyak_length(::Val{N}, ::Val{B}, M::Int) where {N,B}
    # implicit assumption: M ≤ B
    c = zeros(MVector{B+1,Int}) # indexed as 0, …, B
    for b in 0:M
        c[b + 1] = __block_length(b)
    end
    for n in 2:N
        for b in B:(-1):0            # blocks with indices that sum to b
            s = 0
            for a in 0:min(b, M)
                s += __block_length(a) * c[b - a + 1]
            end
            # can safely overwrite since they will not be used again for n + 1
            c[b + 1] = s
        end
    end
    sum(c)
end


text:
# See the example file for a more commented version of this.

using StableRNGs
const data = rand(StableRNG(123), TDist(1.0), 10_000)
const (mle, observed_information_matrix) = fit_mle(BulkAndTailsDist, data, print_level=0)

@testset "Fitting" begin
  for _x in (0.1, 1.0, 10.0, 100.0, 1000.0)
    @test isapprox(pdf(TDist(1.0), _x), pdf(mle, _x), atol=1e-2)
  end
end


text:
# Go Text Protocol
# Standard go interface

# Based on implementation from Michi
function gtp(policy::Policy)
    known_commands = ["boardsize", "clear_board", "komi", "play", "genmove",
                      "final_score", "quit", "name", "version", "known_command",
                      "list_commands", "protocol_version", "tsdebug"]
    board = Board()
    while true
        println(STDERR, "start====")
        flush(STDERR)
        if eof(STDIN)
            return
        end
        line = readline()
        println(STDERR, line)
        println(STDERR, "========")
        flush(STDERR)
        if line == ""
            continue
        end
        command = split(lowercase(line))
        if match(r"\d+", command[1]) != nothing
            cmdid = command[1]
            command = command[2:end]
        else
            cmdid = ""
        end
        ret = ""
        if command[1] == "boardsize"
            if parse(Int, command[2]) != N
                warn("Warning: Trying to play on different board $(command[1]) != $(N)")
                ret = nothing
            end
        elseif command[1] == "clear_board"
            board = Board()
        elseif command[1] == "komi"
            # Set komi
            board.komi = parse(Float64, command[2])
        elseif command[1] == "play"
            if length(command) == 2
                move = parse_coord(command[2])
            else
                move = parse_coord(command[3])
            end
            println(STDERR, command)
            println(STDERR, move)
            play_move(board, move)
            println(STDERR, "played move")
        elseif command[1] == "genmove"
            move = choose_move(board, policy)
            play_move(board, move)
            ret = str_coord(move)
            println(STDERR, ret)
            println(STDERR, move)
            # TODO: Resigning, passing
        elseif command[1] == "final_score"
            score = calculate_score(board)
            if score == 0
                ret = "0"
            elseif score > 0
                ret = @sprintf("B+%f", score)
            elseif score < 0
                ret = @sprintf("W+%f", -score)
            end
        elseif command[1] == "name"
            ret = "go.jl"
        elseif command[1] == "version"
            ret = "Go bot written in Julia"
        elseif command[1] == "tsdebug"
            print_pos(board, output=STDOUT)
        elseif command[1] == "ldebug"
            # Print liberties of each stone
            print_liberties(board, output=STDOUT)
        elseif command[1] == "list_commands"
            ret = join(known_commands, "\n")
        elseif command[1] == "known_command"
            ret = command[2] in known_commands ? "true" : "false"
        elseif command[1] == "protocol_version"
            ret = "2"
        elseif command[1] == "quit"
            @printf("=%s \n\n", cmdid)
            break
        else
            warn("Warning: Ignoring unknown command - $(line)")
            ret = nothing
        end
        print_pos(board)
        if ret != nothing
            out = @sprintf("=%s %s\n\n", cmdid, ret)
        else
            out = @sprintf("?%s ???\n\n", cmdid)
        end
        @printf("%s", out)
        println(STDERR, out)
        flush(STDOUT)
        flush(STDERR)
    end
end


text:
function get_pu_bases(MVAbase, kVbase)
    eurobase = 1 #
    hourbase = 1 #

    Sbase = MVAbase * 1e6
    Vbase = kVbase  * 1e3
    #Zbase = (Vbase)^2 / (Sbase) # TODO
    kAbase = MVAbase / (sqrt(3) * kVbase)
    Zbase = 1/(kAbase^2 / MVAbase)
    Ibase = (Sbase)   / (Vbase)
    timebase = hourbase*3600
    Ebase = Sbase * timebase

    bases = Dict{String, Any}(
    "Z" => Zbase,         # Impedance (Ω)
    "I" => Ibase,         # Current (A)
    "€" => eurobase,      # Currency (€)
    "t" => timebase,      # Time (s)
    "S" => Sbase,         # Power (W)
    "V" => Vbase,         # Voltage (V)
    "E" => Ebase          # Energy (J)
    )
    return bases
end

function process_additional_data!(data)
    to_pu!(data)
    fix_data!(data)
    convert_matpowerdcline_to_branchdc!(data)
end

function is_single_network(data)
    return !haskey(data, "multinetwork") || data["multinetwork"] == false
end

function to_pu!(data)
    if is_single_network(data)
        to_pu_single_network!(data)
    else
        to_pu_multinetwork!(data)
    end
end

function to_pu_single_network!(data)
    MVAbase = data["baseMVA"]
    if haskey(data, "convdc")
        for (i, conv) in data["convdc"]
            dcbus = conv["busdc_i"]
            kVbase = conv["basekVac"]
            Zbase = get_pu_bases(MVAbase, kVbase)["Z"]
            Ibase = get_pu_bases(MVAbase, kVbase)["I"]

            set_conv_pu_power(conv, MVAbase)
            set_conv_pu_volt(conv, kVbase*sqrt(3))
            set_conv_pu_ohm(conv, Zbase)
        end
    end
    if haskey(data, "branchdc")
        for (i, branchdc) in data["branchdc"]
            set_branchdc_pu(branchdc, MVAbase)
        end
    end
    if haskey(data, "busdc")
        for (i, busdc) in data["busdc"]
            set_busdc_pu(busdc, MVAbase)
        end
    end
    if haskey(data, "convdc_ne")
        for (i, conv) in data["convdc_ne"]
            dcbus = conv["busdc_i"]
            kVbase = conv["basekVac"]
            Zbase = get_pu_bases(MVAbase, kVbase)["Z"]
            Ibase = get_pu_bases(MVAbase, kVbase)["I"]

            set_conv_pu_power(conv, MVAbase)
            set_conv_pu_volt(conv, kVbase*sqrt(3))
            set_conv_pu_ohm(conv, Zbase)
        end
    end
    if haskey(data, "branchdc_ne")
        for (i, branchdc) in data["branchdc_ne"]
            set_branchdc_pu(branchdc, MVAbase)
        end
    end
    if haskey(data, "busdc_ne")
        new_busdc_ne = Dict{String, Any}()
        for (i, busdc) in data["busdc_ne"]
            set_busdc_pu(busdc, MVAbase)
            new_bus = busdc["busdc_i"] # assigning new bus numbers: continous numbers from dc bus numbers
            if new_bus == i
                display("candidate dc buses and existing dc buses should have different bus numbers")
            end
            new_busdc_ne[string(new_bus)] = busdc # assigning new bus numbers: continous numbers from dc bus numbers
            busdc["index"] = new_bus
        end
        data["busdc_ne"] = new_busdc_ne # assigning new bus numbers: continous numbers from dc bus numbers
    end
end

function to_pu_multinetwork!(data)
    for (n, network) in data["nw"]
        MVAbase = network["baseMVA"]
        if haskey(data["nw"][n], "convdc")
            for (i, conv) in data["nw"][n]["convdc"]
                dcbus = conv["busdc_i"]
                kVbase = conv["basekVac"]
                Zbase = get_pu_bases(MVAbase, kVbase)["Z"]
                Ibase = get_pu_bases(MVAbase, kVbase)["I"]

                set_conv_pu_power(conv, MVAbase)
                set_conv_pu_volt(conv, kVbase*sqrt(3))
                set_conv_pu_ohm(conv, Zbase)
            end
        end
        if haskey(data["nw"][n], "branchdc")
            for (i, branchdc) in data["nw"][n]["branchdc"]
                set_branchdc_pu(branchdc, MVAbase)
            end
        end
        if haskey(data["nw"][n], "busdc")
            for (i, busdc) in data["nw"][n]["busdc"]
                set_busdc_pu(busdc, MVAbase)
            end
        end
        if haskey(data["nw"][n], "convdc_ne")
            for (i, conv) in data["nw"][n]["convdc_ne"]
                dcbus = conv["busdc_i"]
                kVbase = conv["basekVac"]
                Zbase = get_pu_bases(MVAbase, kVbase)["Z"]
                Ibase = get_pu_bases(MVAbase, kVbase)["I"]

                set_conv_pu_power(conv, MVAbase)
                set_conv_pu_volt(conv, kVbase*sqrt(3))
                set_conv_pu_ohm(conv, Zbase)
                # conv["cost"] = conv["cost"]/length(data["nw"])
            end
        end
        if haskey(data["nw"][n], "branchdc_ne")
            for (i, branchdc) in data["nw"][n]["branchdc_ne"]
                set_branchdc_pu(branchdc, MVAbase)
                # branchdc["cost"] = branchdc["cost"]/length(data["nw"])
            end
        end
        # if haskey(data["nw"][n], "ne_branch")
        #     for (i, ne_branch) in data["nw"][n]["ne_branch"]
        #         ne_branch["construction_cost"] = ne_branch["construction_cost"]/length(data["nw"])
        #     end
        # end
        if haskey(data["nw"][n], "busdc_ne")
            new_busdc_ne = Dict{String, Any}()
            for (i, busdc) in data["nw"][n]["busdc_ne"]
                set_busdc_pu(busdc, MVAbase)
                new_bus = busdc["busdc_i"] # assigning new bus numbers: continous numbers from exisiting dc bus numbers
                if new_bus == i
                    display("candidate dc buses and existing dc buses should have different bus numbers")
                end
                new_busdc_ne[string(new_bus)] = busdc # assigning new bus numbers: continous numbers from dc bus numbers
                busdc["index"] = new_bus
            end
            data["nw"][string(n)]["busdc_ne"] = new_busdc_ne # assigning new bus numbers: continous numbers from exisiting dc bus numbers
        end
    end
end

function convert_matpowerdcline_to_branchdc!(data)
    if is_single_network(data)
        convert_matpowerdcline_to_branchdc_single_network!(data)
    else
        convert_matpowerdcline_to_branchdc_multinetwork!(data)
    end
end

function convert_matpowerdcline_to_branchdc_single_network!(data)
    MVAbase = data["baseMVA"]
    if haskey(data, "dcline") && haskey(data["dcline"], "1")
        if !haskey(data, "convdc")
            data["convdc"] = Dict{String, Any}()
        end
        if !haskey(data, "branchdc")
            data["branchdc"] = Dict{String, Any}()
        end
        if !haskey(data, "busdc")
            data["busdc"] = Dict{String, Any}()
        end
        conv_i = length(data["convdc"])
        branch_i = length(data["branchdc"])
        bus_i = length(data["busdc"])

        for (i, dcline) in data["dcline"]
            # make DC bus from side
            bus_i = bus_i + 1
            data["busdc"]["$bus_i"] = get_busdc(bus_i)

            prev_bus = bus_i - 1
            if haskey(data["busdc"],["$prev_bus"])
                data["busdc"]["$bus_i"]["basekVdc"] = data["busdc"]["$prev_bus"]["basekVdc"]
            else
                data["busdc"]["$bus_i"]["basekVdc"] = 100 # arbitrary choice
            end

            # DC bus to
            bus_i = bus_i + 1
            data["busdc"]["$bus_i"] = get_busdc(bus_i)
            prev_bus = bus_i - 1
            if haskey(data["busdc"],["$prev_bus"])
                data["busdc"]["$bus_i"]["basekVdc"] = data["busdc"]["$prev_bus"]["basekVdc"]
            else
                data["busdc"]["$bus_i"]["basekVdc"] = 100 # arbitrary choice
            end

            branch_i = branch_i + 1
            conv_i = conv_i + 1
            converter1, converter2, branchdc = convert_to_dcbranch_and_converters(data, dcline, branch_i, conv_i, bus_i-1, bus_i)
            data["branchdc"]["$branch_i"] = branchdc
            data["convdc"]["$conv_i"] = converter1
            conv_i = conv_i + 1
            data["convdc"]["$conv_i"] = converter2
        end
    end
    data["dcline"] = []
end

function convert_matpowerdcline_to_branchdc_multinetwork!(data)
    for (n, network) in data["nw"]
        MVAbase = network["baseMVA"]
        if haskey(data["nw"][n], "dcline") && haskey(data["nw"][n]["dcline"], "1")
            if !haskey(data["nw"][n], "convdc")
                data["nw"][n]["convdc"] = Dict{String, Any}()
            end
            if !haskey(data["nw"][n], "branchdc")
                data["nw"][n]["branchdc"] = Dict{String, Any}()
            end
            if !haskey(data["nw"][n], "busdc")
                data["nw"][n]["busdc"] = Dict{String, Any}()
            end
            conv_i = length(data["nw"][n]["convdc"])
            branch_i = length(data["nw"][n]["branchdc"])
            bus_i = length(data["nw"][n]["busdc"])

            for (i, dcline) in data["nw"][n]["dcline"]
                # make DC bus from side
                bus_i = bus_i + 1
                data["nw"][n]["busdc"]["$bus_i"] = get_busdc(bus_i)

                prev_bus = bus_i - 1
                if haskey(data["nw"][n]["busdc"],["$prev_bus"])
                    data["nw"][n]["busdc"]["$bus_i"]["basekVdc"] = data["nw"][n]["busdc"]["$prev_bus"]["basekVdc"]
                else
                    data["nw"][n]["busdc"]["$bus_i"]["basekVdc"] = 100 # arbitrary choice
                end

                # DC bus to
                bus_i = bus_i + 1
                data["nw"][n]["busdc"]["$bus_i"] = get_busdc(bus_i)
                prev_bus = bus_i - 1
                if haskey(data["nw"][n]["busdc"],["$prev_bus"])
                    data["nw"][n]["busdc"]["$bus_i"]["basekVdc"] = data["nw"][n]["busdc"]["$prev_bus"]["basekVdc"]
                else
                    data["nw"][n]["busdc"]["$bus_i"]["basekVdc"] = 100 # arbitrary choice
                end

                branch_i = branch_i + 1
                conv_i = conv_i + 1
                converter1, converter2, branchdc = convert_to_dcbranch_and_converters(data["nw"][n], dcline, branch_i, conv_i, bus_i-1, bus_i)
                data["nw"][n]["branchdc"]["$branch_i"] = branchdc
                data["nw"][n]["convdc"]["$conv_i"] = converter1
                conv_i = conv_i + 1
                data["nw"][n]["convdc"]["$conv_i"] = converter2
            end
        end
        data["nw"][n]["dcline"] = []
    end
end

function fix_data!(data)

    rescale_energy_cost = x -> (MWhbase/dollarbase)*x

    if is_single_network(data)
        fix_data_single_network!(data)
    else
        fix_data_multinetwork!(data)
    end
end

function fix_data_single_network!(data)
    MVAbase = data["baseMVA"]
    @assert(MVAbase>0)
    if haskey(data, "convdc")
        for (i, conv) in data["convdc"]
            check_conv_parameters(conv)
        end
    end
    if haskey(data, "branchdc")
        for (i, branchdc) in data["branchdc"]
            check_branchdc_parameters(branchdc)
        end
    end
    if haskey(data, "busdc")
        new_busdc = Dict{String, Any}()
        for (i, busdc) in data["busdc"]
            new_bus = busdc["busdc_i"] # assigning new bus numbers: continous numbers from dc bus numbers
            new_busdc[string(new_bus)] = busdc # assigning new bus numbers: continous numbers from dc bus numbers
        end
        data["busdc"] = new_busdc # assigning new bus numbers: continous numbers from dc bus numbers
    end
    if !haskey(data, "dcpol")
        data["dcpol"] = 2
    end
    if haskey(data, "convdc_ne")
        for (i, conv) in data["convdc_ne"]
            check_conv_parameters(conv)
        end
    end
    if haskey(data, "branchdc_ne")
        for (i, branchdc) in data["branchdc_ne"]
            check_branchdc_parameters(branchdc)
        end
    end
end

function fix_data_multinetwork!(data)
    for (n, network) in data["nw"]
        MVAbase = network["baseMVA"]
        @assert(MVAbase>0)
        if haskey(data["nw"][n], "convdc")
            for (i, conv) in data["nw"][n]["convdc"]
                check_conv_parameters(conv)
            end
        end
        if haskey(data["nw"][n], "branchdc")
            for (i, branchdc) in data["nw"][n]["branchdc"]
                check_branchdc_parameters(branchdc)
            end
        end
        if haskey(data["nw"][n], "busdc")
            new_busdc = Dict{String, Any}()
            for (i, busdc) in data["nw"][n]["busdc"]
                new_bus = busdc["busdc_i"] # assigning new bus numbers: continous numbers from dc bus numbers
                new_busdc[string(new_bus)] = busdc # assigning new bus numbers: continous numbers from dc bus numbers
            end
            data["nw"][n]["busdc"] = new_busdc # assigning new bus numbers: continous numbers from dc bus numbers
        end
        if !haskey(data["nw"][n], "dcpol")
            data["nw"][n]["dcpol"] = 2
        end
        if haskey(data["nw"][n], "convdc_ne")
            for (i, conv) in data["nw"][n]["convdc_ne"]
                check_conv_parameters(conv)
            end
        end
        if haskey(data["nw"][n], "branchdc_ne")
            for (i, branchdc) in data["nw"][n]["branchdc_ne"]
                check_branchdc_parameters(branchdc)
            end
        end
    end
end

function check_branchdc_parameters(branchdc)
    @assert(branchdc["rateA"]>=0)
    @assert(branchdc["rateB"]>=0)
    @assert(branchdc["rateC"]>=0)
end

function set_branchdc_pu(branchdc, MVAbase)
    rescale_power = x -> x/MVAbase
    _PM._apply_func!(branchdc, "rateA", rescale_power)
    _PM._apply_func!(branchdc, "rateB", rescale_power)
    _PM._apply_func!(branchdc, "rateC", rescale_power)
end

function set_busdc_pu(busdc, MVAbase)
    rescale_power = x -> x/MVAbase
    _PM._apply_func!(busdc, "Pdc", rescale_power)
end

function set_conv_pu_power(conv, MVAbase)
    rescale_power = x -> x/MVAbase
    _PM._apply_func!(conv, "P_g", rescale_power)
    _PM._apply_func!(conv, "Q_g", rescale_power)
    _PM._apply_func!(conv, "Pdcset", rescale_power)
    _PM._apply_func!(conv, "LossA", rescale_power)
    _PM._apply_func!(conv, "Pacmax", rescale_power)
    _PM._apply_func!(conv, "Pacmin", rescale_power)
    _PM._apply_func!(conv, "Qacmax", rescale_power)
    _PM._apply_func!(conv, "Qacmin", rescale_power)
    _PM._apply_func!(conv, "Pacrated", rescale_power)
    _PM._apply_func!(conv, "Qacrated", rescale_power)
end

function set_conv_pu_volt(conv, kVbase)
    rescale_volt = x -> x  / (kVbase)
    _PM._apply_func!(conv, "LossB", rescale_volt)
end

function set_conv_pu_ohm(conv, Zbase)
    rescale_ohm = x -> x / Zbase
    _PM._apply_func!(conv, "LossCrec", rescale_ohm)
    _PM._apply_func!(conv, "LossCinv", rescale_ohm)
end

function check_conv_parameters(conv)
    @assert(conv["LossA"]>=0)
    @assert(conv["LossB"]>=0)
    @assert(conv["LossCrec"]>=0)
    @assert(conv["LossCinv"]>=0)
    conv_id = conv["index"]
    conv["Pacrated"] = max(abs(conv["Pacmax"]),abs(conv["Pacmin"]))
    conv["Qacrated"] = max(abs(conv["Qacmax"]),abs(conv["Qacmin"]))
    if conv["Imax"] < sqrt(conv["Pacrated"]^2 + conv["Qacrated"]^2)
        Memento.warn(_PM._LOGGER, "Inconsistent current limit for converter $conv_id, it will be updated.")
        conv["Imax"] = sqrt(conv["Pacrated"]^2 + conv["Qacrated"]^2)
    end
    if conv["LossCrec"] != conv["LossCinv"]
        Memento.warn(_PM._LOGGER, "The losses of converter $conv_id are different in inverter and rectifier mode, inverter losses are used.")
    end
    if conv["islcc"] == 1
        Memento.warn(_PM._LOGGER, "Converter $conv_id is an LCC, reactive power limits might be updated.")
        if abs(conv["Pacmax"]) >= abs(conv["Pacmin"])
            conv["phimin"] = 0
            conv["phimax"] = acos(conv["Pacmin"] / conv["Pacmax"])
        else
            conv["phimin"] = pi - acos(conv["Pacmax"] / conv["Pacmin"])
            conv["phimax"] = pi
        end
        conv["Qacmax"] = conv["Pacrated"]
        conv["Qacrated"] = conv["Pacrated"]
        conv["Qacmin"] =  0
    end
    @assert(conv["Pacmax"]>=conv["Pacmin"])
    @assert(conv["Qacmax"]>=conv["Qacmin"])
    @assert(conv["Pacrated"]>=0)
    @assert(conv["Qacrated"]>=0)
end

function get_branchdc(matpowerdcline, branch_i, fbusdc, tbusdc)
    branchdc = Dict{String, Any}()
    branchdc["index"] = branch_i
    branchdc["fbusdc"] = fbusdc
    branchdc["tbusdc"] = tbusdc
    branchdc["r"] = 0
    branchdc["l"] = 0
    branchdc["c"] = 0
    branchdc["rateA"] = max(matpowerdcline["pmaxf"], matpowerdcline["pmaxt"])
    branchdc["rateB"] = max(matpowerdcline["pmaxf"], matpowerdcline["pmaxt"])
    branchdc["rateC"] = max(matpowerdcline["pmaxf"], matpowerdcline["pmaxt"])
    branchdc["status"] = matpowerdcline["br_status"]
    return branchdc
end

function get_busdc(bus_i)
    busdc = Dict{String, Any}()
    busdc["index"] = bus_i
    busdc["busdc_i"] = bus_i
    busdc["grid"] = 1
    busdc["Pdc"] = 0
    busdc["Vdc"] = 1
    busdc["Vdcmax"] = 1.1
    busdc["Vdcmin"] = 0.9
    busdc["Cdc"] = 0
    return busdc
end

function get_converter(conv_i, dcbus, acbus, kVbaseAC, vmax, vmin, status, pac, qac, qmaxac, qminac, vmac, Imax, lossA, lossB, pmaxac, pminac)
    conv = Dict{String, Any}()
    conv["index"] = conv_i
    conv["busdc_i"] = dcbus
    conv["busac_i"] = acbus
    conv["type_dc"] = 1
    conv["type_ac"] = 2
    conv["P_g"] = pac
    conv["Q_g"] = qac
    conv["Vtar"] = vmac
    conv["rtf"] = 0
    conv["xtf"] = 0
    conv["transformer"] = 0
    conv["bf"] = 0
    conv["filter"] = 0
    conv["rc"] = 0
    conv["xc"] = 0
    conv["reactor"] = 0
    conv["basekVac"] = kVbaseAC
    conv["Vmmax"] =  vmax
    conv["Vmmin"] =  vmin
    conv["Imax"] =  Imax # assuming 1pu
    conv["status"] = status
    conv["LossA"] = lossA
    conv["LossB"] = lossB
    conv["LossCrec"] = 0
    conv["LossCinv"] = 0
    conv["droop"] = 0
    conv["islcc"] = 0
    conv["Pdcset"] = 0
    conv["Vdcset"] = 1
    conv["tm"] = 1
    conv["dVdcSet"] = 0
    conv["Qacmax"] = qmaxac
    conv["Qacmin"] = qminac
    conv["Pacmax"] = pmaxac
    conv["Pacmin"] = pminac
    check_conv_parameters(conv)
    return conv
end

function convert_to_dcbranch_and_converters(data, dcline, branchdc_id, conv_i, fbusdc, tbusdc)
    # make one more DC branch
    branchdc = get_branchdc(dcline, branchdc_id, fbusdc, tbusdc)
    vmaxf = data["bus"]["$fbusdc"]["vmax"]
    vminf = data["bus"]["$fbusdc"]["vmin"]
    vmaxt =  data["bus"]["$tbusdc"]["vmax"]
    vmint =  data["bus"]["$tbusdc"]["vmin"]

    pac = dcline["pf"]
    qac = dcline["qf"]
    vmac = dcline["vf"]
    acbus = dcline["f_bus"]
    kVbaseAC = data["bus"]["$fbusdc"]["base_kv"]
    Imax =  sqrt(dcline["pmaxf"]^2 + dcline["qmaxf"]^2) / vminf
    status = dcline["br_status"]
    qmaxac = dcline["qmaxf"]
    qminac = dcline["qminf"]
    lossA = dcline["loss0"] /2
    lossB = dcline["loss1"]
    vac = 1
    pmaxac = dcline["pmaxf"]
    pminac = dcline["pminf"]
    converter1 =  get_converter(conv_i, fbusdc, acbus, kVbaseAC, vmaxf, vminf, status, pac, qac, qmaxac, qminac, vac, Imax, lossA, lossB, pmaxac, pminac)

    # converter 2
    conv_i = conv_i + 1
    acbus = dcline["t_bus"]
    pac = dcline["pt"]
    qac = dcline["qt"]
    vmac = dcline["vt"]
    kVbaseAC = data["bus"]["$tbusdc"]["base_kv"]
    Imax =  sqrt(dcline["pmaxt"]^2 + dcline["qmaxt"]^2) / vmint# assuming 1pu
    status = dcline["br_status"]
    lossA = dcline["loss0"] / 2
    lossB = 0
    pmaxac =  dcline["pmaxt"]
    pminac =  dcline["pmint"]
    qmaxac = dcline["qmaxt"]
    qminac = dcline["qmint"]
    converter2 =  get_converter(conv_i, tbusdc, acbus, kVbaseAC, vmaxt, vmint, status, pac, qac, qmaxac, qminac, vac, Imax, lossA, lossB, pmaxac, pminac)

    return converter1, converter2, branchdc
end

function converter_bounds(pmin, pmax, loss0, loss1)
    if pmin >= 0 && pmax >=0
        pminf = pmin
        pmaxf = pmax
        pmint = loss0 - pmaxf * (1 - loss1)
        pmaxt = loss0 - pminf * (1 - loss1)
    end
    if pmin >= 0 && pmax < 0
        pminf = pmin
        pmint = pmax
        pmaxf = (-pmint + loss0) / (1-loss1)
        pmaxt = loss0 - pminf * (1 - loss1)
    end
    if pmin < 0 && pmax >= 0
        pmaxt = -pmin
        pmaxf = pmax
        pminf = (-pmaxt + loss0) / (1-loss1)
        pmint = loss0 - pmaxf * (1 - loss1)
    end
    if pmin < 0 && pmax < 0
        pmaxt = -pmin
        pmint = pmax
        pmaxf = (-pmint + loss0) / (1-loss1)
        pminf = (-pmaxt + loss0) / (1-loss1)
    end
    return pminf, pmaxf, pmint, pmaxt
end


text:
using CSV, DataFrames, StatsPlots, StanSample, Statistics

ProjDir = @__DIR__

old_df = CSV.read(joinpath(ProjDir, "data", "data_final.csv"), delim=' ')

df = DataFrame()
for row in eachrow(old_df)
	vals = split(row[1], "\t")
	yr = Meta.parse(row[1][1:4])
	wk = Meta.parse(row[1][5:6])
	s = Meta.parse(vals[2])
	p = Meta.parse(vals[3])
	q = Meta.parse(vals[4])
	append!(df, DataFrame(year=yr, week=wk, store=s, quantity=q, price=p))
end

df[!, :ppu] = df[:, :price] ./ df[:, :quantity]
df[!, :ppu_l] = log.(df[:, :ppu])
df[!, :quantity_l] = log.(df[:, :quantity])

# Clean up

df = df[df.ppu_l .> 0.04, :]
df = df[df.quantity_l .> 0.0, :]

CSV.write(joinpath(ProjDir, "data", "data.csv"), df)

p1 = scatter(df[:, :ppu], df[:, :quantity],
	xlabel="ppu", ylabel="quantity sold")
p2 = scatter(df[:, :ppu_l], df[:, :quantity_l],
	xlabel="log(ppu)", ylabel="log(quantity sold)")
plot(p1, p2, layout=(2, 1))
savefig(joinpath(ProjDir, "plots", "data.png"))

bglm = "
data {
  int<lower=0> N;
  vector[N] q; // Outcome
  vector[N] p; // Predictor
}
parameters {
  real<lower=0> a;
  real b;
  real<lower=0> sigma;
}
model {
  vector[N] mu;            		// mu is a vector
  a ~ normal(2.0, 10.0);		// intercept positive
  b ~ normal(-2.0, 10.0);		// demand drops as price increases
  sigma ~ uniform(0, 15);
  mu = a + b * p;
  q ~ normal(mu, sigma);
}
";

data = Dict(
  :N => size(df, 1),
  :q => df[:, :quantity_l],
  :p => df[:, :ppu_l]
)

sm = SampleModel("bglm", bglm)

rc = stan_sample(sm, data=data)

if success(rc)
	chns = read_samples(sm)
	plot(chns)
	savefig(joinpath(ProjDir, "plots", "chains.png"))

	dfs = DataFrame(chns)
	show(chns)

	p = 0.6:0.01:2.0
	p1 = scatter(df[:, :ppu_l], df[:, :quantity_l],
		xlabel="ppu_l", ylabel="quantity_l sold")
	a_mean = mean(dfs[:, :a])
	b_mean = mean(dfs[:, :b])
	mu_q = a_mean .+ b_mean .* p
	plot!(p1, p, mu_q)
	savefig(joinpath(ProjDir, "plots", "regression.png"))

	price_points_l = [1.05, 1.35, 1.5, 1,6, 1.78]
	scatter(price_points_l, a_mean .+ b_mean .* price_points_l,
		xlims = (0.6, 2.0), ylims= (0, 7), leg=false,
		xlabel="log(ppu)", ylabel="log(quantity sold)")
	for row in eachrow(dfs)
		plot!(p, row[:a] .+ row[:b] .* p, color=:gray)
	end
	scatter!(price_points_l, a_mean .+ b_mean .* price_points_l)
	savefig(joinpath(ProjDir, "plots", "sampling.png"))
end


text:
module _CondEn
export CondEn
using Statistics: std, mean
using StatsBase: Histogram, fit
    """
        Cond, SEw, SEz = CondEn(Sig)

    Returns the corrected conditional entropy estimates (`Cond`) and the
    corresponding Shannon entropies (m: `SEw`, m+1: `SEz`) for m = [1,2]
    estimated from the data sequence (`Sig`) using the default  parameters:
    embedding dimension = 2, time delay = 1, symbols = 6, logarithm = natural,
    normalisation = false
    *Note: CondEn(m=1) returns the Shannon entropy of `Sig`.*

        Cond, SEw, SEz = CondEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, c::Int=6, Logx::Real=exp(1), Norm::Bool=false)

    Returns the corrected conditional entropy estimates (`Cond`) from the data
    sequence (`Sig`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, an integer > 1  \n
    `tau`   - Time Delay, a positive integer  \n
    `c`     - # of symbols, an integer > 1  \n
    `Logx`  - Logarithm base, a positive scalar  \n
    `Norm`  - Normalisation of CondEn value:  \n
              [false]  no normalisation - default
              [true]   normalises w.r.t Shannon entropy of data sequence `Sig`

    # See also `XCondEn`, `MSEn`, `PermEn`, `DistEn`, `XPermEn`

    # References:
        [1] Alberto Porta, et al.,
            "Measuring regularity by means of a corrected conditional
            entropy in sympathetic outflow."
            Biological cybernetics
            78.1 (1998): 71-78.

    """
    function CondEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1,
        c::Int=6, Logx::Real=exp(1), Norm::Bool=false)

    (size(Sig)[1] > 10) ? nothing :   error("Sig:   must be a numeric vector")
    (m > 1) ? nothing :  error("m:     must be an integer > 1")
    (tau > 0) ? nothing :  error("tau:   must be an integer > 0")
    (Logx > 0) ? nothing : error("Logx:     must be a positive number > 0")
    (c > 1) ? nothing :  error("c:     must be an integer > 1")

    Edges = range(minimum(Sig),maximum(Sig),length=c+1)
    Sx = map(x -> sum(Edges[1:c].<=x), Sig)

    N = size(Sx)[1]
    SEw = zeros(m-1)
    SEz = zeros(m-1)
    Prcm = zeros(m-1)
    Xi = zeros(N,m)
    Xi[:,m] = Sx
    for k = 1:m-1
        Nx = N-(k*tau)
        Xi[1:Nx,end-k] = Sx[(k*tau)+1:N]
        Wi = Xi[1:Nx,m-k+1:m] * (c.^collect(k-1:-1:0)) # Maybe dot notation here???
        Zi = Xi[1:Nx,m-k:m] * (c.^collect(k:-1:0))

        Pw = fit(Histogram, Wi, minimum(Wi)-.5:maximum(Wi)+.5).weights
        Pz = fit(Histogram, Zi, minimum(Zi)-.5:maximum(Zi)+.5).weights
        Prcm[k] = sum(Pw.==1)/Nx

        if sum(Pw)!= Nx || sum(Pz)!= Nx
            @warn("Potential error estimating probabilities.")
        end

        Pw = Pw[Pw.!=0];    Pw /= N;
        Pz = Pz[Pz.!=0];    Pz /= N;
        SEw[k] = -transpose(Pw)*log.(Logx, Pw)
        SEz[k] = -transpose(Pz)*log.(Logx, Pz)

    end

    Temp = fit(Histogram,Sx,.5:c+.5).weights/N;
    Temp = Temp[Temp.!=0]
    S1 = -transpose(Temp)*log.(Logx, Temp)
    Cond = SEz - SEw + Prcm*S1;
    Cond = vcat(S1, Cond);
    if Norm
        Cond = Cond/S1;
    end

    return Cond, SEw, SEz
    end

end

    #= Sig = (Sig.-mean(Sig))./std(Sig,corrected=false);
    Edges = range(minimum(Sig),maximum(Sig),length=c+1)
    #Edges[1] -= .1; Edges[end] += .1
    Sx = map(x -> searchsortedfirst(Edges,x), Sig) .- 1
    Sx[Sx.==0] .= 1 =#

    """
    Copyright 2021 Matthew W. Flood, EntropyHub

    Licensed under the Apache License, Version 2.0 (the "License");
    you may not use this file except in compliance with the License.
    You may obtain a copy of the License at

         http://www.apache.org/licenses/LICENSE-2.0

    Unless required by applicable law or agreed to in writing, software
    distributed under the License is distributed on an "AS IS" BASIS,
    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    See the License for the specific language governing permissions and
    limitations under the License.

    For Terms of Use see https://github.com/MattWillFlood/EntropyHub
    """


text:
using ApproxFun

#The following solves the Airy ODE with dirichlet boundary conditions

x=Fun(identity,-1000..15)   # Fun corresponding to multiplication by x, on [-100,15]
d=domain(x)
D=Derivative()             # The derivative operator
B=Dirichlet()              # Dirichlet boundary conditions, [u(-100),u(15)]

#Construct operator

A=[B;D^2-x]                # This is dirichlet conditions and u'' - x u
b=[airyai(leftendpoint(d)),0]     # We want it to equal airyai(-100) at -100, and 0 at
                           # 10, with 0 rhs
#Solve ODE

u=A\[b,0]                  # u satisfies A*u = b, or in other words,
                           # B*u = [airyai(-100.),0.] and (D2 - x)*u = 0.

# Check the accuracy
norm(u - Fun(airyai,d))

## We now solve with Neumann conditions

B=Neumann()
A=[B;D^2-x]
b=[[airyaiprime(leftendpoint(d)),0.],0.]

u=A\b

# Check the accuracy
norm(u - Fun(airyai,d))

max_stars_repo_path: ODEs/Airy equation.jl

max_issues_repo_path: ODEs/Airy equation.jl

max_forks_repo_path: ODEs/Airy equation.jl


text:
module BinaryTraits

using MacroTools

export @trait, @assign, @implement, @check
export check, istrait, traits, required_contracts
export inittraits

include("types.jl")
include("misc.jl")
include("utils.jl")
include("trait.jl")
include("assignment.jl")
include("interface.jl")

end # module


text:
"""
    DeterministicEquivalent

Deterministic equivalent memory structure. Stochastic program is stored as one large optimization problem. Supported by any standard `AbstractOptimizer`.

"""
struct DeterministicEquivalent{N, M, S <: NTuple{M, Scenarios}} <: AbstractStochasticStructure{N}
    decisions::Decisions
    decision_variables::NTuple{M, Vector{DecisionRef}}
    scenarios::S
    sub_objectives::NTuple{M, Vector{MOI.AbstractScalarFunction}}
    model::JuMP.Model

    function DeterministicEquivalent(scenarios::NTuple{M, Scenarios}) where M
        N = M + 1
        decisions = Decisions()
        decision_variables = ntuple(Val(M)) do i
            Vector{DecisionRef}()
        end
        sub_objectives = ntuple(Val(M)) do i
            Vector{MOI.AbstractScalarFunction}()
        end
        S = typeof(scenarios)
        return new{N,M,S}(decisions, decision_variables, scenarios, sub_objectives, Model())
    end
end

function StochasticStructure(scenario_types::ScenarioTypes{M}, ::Deterministic) where M
    scenarios = ntuple(Val(M)) do i
        Vector{scenario_types[i]}()
    end
    return DeterministicEquivalent(scenarios)
end

function StochasticStructure(scenarios::NTuple{M, Vector{<:AbstractScenario}}, ::Deterministic) where M
    return DeterministicEquivalent(scenarios)
end

# Base overloads #
# ========================== #
function Base.print(io::IO, structure::DeterministicEquivalent)
    print(io, "Deterministic equivalent problem\n")
    print(io, structure.model)
end
# ========================== #

# MOI #
# ========================== #
function MOI.get(structure::DeterministicEquivalent, attr::MOI.AbstractModelAttribute)
    return MOI.get(backend(structure.model), attr)
end
function MOI.get(structure::DeterministicEquivalent, attr::MOI.AbstractVariableAttribute, index::MOI.VariableIndex)
    return MOI.get(backend(structure.model), attr, index)
end
function MOI.get(structure::DeterministicEquivalent, attr::Type{MOI.VariableIndex}, name::String)
    return MOI.get(backend(structure.model), attr, name)
end
function MOI.get(structure::DeterministicEquivalent, attr::MOI.AbstractConstraintAttribute, cindex::MOI.ConstraintIndex)
    return MOI.get(backend(structure.model), attr, cindex)
end

function MOI.set(structure::DeterministicEquivalent, attr::MOI.Silent, flag)
    MOI.set(backend(structure.model), attr, flag)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::MOI.AbstractModelAttribute, value)
    MOI.set(backend(structure.model), attr, value)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    MOI.set(backend(structure.model), attr, index, value)
    return nothing
end
function MOI.set(structure::DeterministicEquivalent, attr::MOI.AbstractConstraintAttribute,
                 cindex::MOI.ConstraintIndex, value)
    MOI.set(backend(structure.model), attr, cindex, value)
    return nothing
end

function MOI.is_valid(structure::DeterministicEquivalent, index::MOI.VariableIndex)
    return MOI.is_valid(backend(structure.model), index)
end

function MOI.add_constraint(structure::DeterministicEquivalent, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    return MOI.add_constraint(backend(structure.model), f, s)
end

function MOI.delete(structure::DeterministicEquivalent, index::MOI.Index)
    # TODO: more to do if index is decision
    MOI.delete(backend(structure.model), index)
    return nothing
end

# Getters #
# ========================== #
function structure_name(structure::DeterministicEquivalent)
    return "Deterministic equivalent"
end
function decision(structure::DeterministicEquivalent, index::MOI.VariableIndex)
    return decision(structure.decisions, index)
end
function decisions(structure::DeterministicEquivalent)
    return structure.decisions
end
function all_decisions(structure::DeterministicEquivalent)
    return structure.decisions.undecided
end
function num_decisions(structure::DeterministicEquivalent{N}, s::Integer = 2) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    s == N && error("The final stage does not have decisions.")
    return length(structure.decision_variables[s])
end
function scenario(structure::DeterministicEquivalent{N}, i::Integer, s::Integer = 2) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    s == 1 && error("The first stage does not have scenarios.")
    return structure.scenarios[s-1][i]
end
function scenarios(structure::DeterministicEquivalent{N}, s::Integer = 2) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    s == 1 && error("The first stage does not have scenarios.")
    return structure.scenarios[s-1]
end
function subproblem(structure::DeterministicEquivalent, i::Integer, s::Integer = 2)
    return subproblem(scenarioproblems(structure, s), i)
end
function subproblems(structure::DeterministicEquivalent, s::Integer = 2)
    return subproblems(scenarioproblems(structure, s))
end
function num_subproblems(structure::DeterministicEquivalent, s::Integer = 2)
    return 0
end
function deferred(structure::DeterministicEquivalent)
    return num_variables(structure.model) == 0
end
# ========================== #

# Setters
# ========================== #
function update_decisions!(structure::DeterministicEquivalent, change::DecisionModification)
    update_decisions!(structure.model, change)
end

function add_scenario!(structure::DeterministicEquivalent, scenario::AbstractScenario, stage::Integer = 2)
    push!(scenarios(structure, stage), scenario)
    return nothing
end
function add_worker_scenario!(structure::DeterministicEquivalent, scenario::AbstractScenario, w::Integer, stage::Integer = 2)
    add_scenario!(structure, scenario, stage)
    return nothing
end
function add_scenario!(scenariogenerator::Function, structure::DeterministicEquivalent, stage::Integer = 2)
    add_scenario!(structure, scenariogenerator(), stage)
    return nothing
end
function add_worker_scenario!(scenariogenerator::Function, structure::DeterministicEquivalent, w::Integer, stage::Integer = 2)
    add_scenario!(scenariogenerator, structure, stage)
    return nothing
end
function add_scenarios!(structure::DeterministicEquivalent, _scenarios::Vector{<:AbstractScenario}, stage::Integer = 2)
    append!(scenarios(structure, stage), _scenarios)
    return nothing
end
function add_worker_scenarios!(structure::DeterministicEquivalent, scenarios::Vector{<:AbstractScenario}, w::Integer, stage::Integer = 2)
    add_scenarios!(structure, scenarios, stasge)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, structure::DeterministicEquivalent, n::Integer, stage::Integer = 2)
    for i = 1:n
        add_scenario!(structure, stage) do
            return scenariogenerator()
        end
    end
    return nothing
end
function add_worker_scenarios!(scenariogenerator::Function, structure::DeterministicEquivalent, n::Integer, w::Integer, stage::Integer = 2)
    add_scenarios!(scenariogenerator, structure, n, stage)
    return nothing
end
function sample!(structure::DeterministicEquivalent, sampler::AbstractSampler, n::Integer, stage::Integer = 2)
    sample!(scenarios(structure, stage), sampler, n)
    return nothing
end
# ========================== #


text:
module Architectures

export
    @hascuda,
    AbstractArchitecture, AbstractCPUArchitecture, AbstractGPUArchitecture, CPU, GPU,
    device, architecture, array_type, arch_array

using CUDA

using KernelAbstractions

"""
    AbstractArchitecture

Abstract supertype for architectures supported by Oceananigans.
"""
abstract type AbstractArchitecture end

"""
    AbstractCPUArchitecture

Abstract supertype for CPU architectures supported by Oceananigans.
"""
abstract type AbstractCPUArchitecture <: AbstractArchitecture end

"""
    AbstractGPUArchitecture

Abstract supertype for GPU architectures supported by Oceananigans.
"""
abstract type AbstractGPUArchitecture <: AbstractArchitecture end

"""
    CPU <: AbstractArchitecture

Run Oceananigans on one CPU node. Uses multiple threads if the environment
variable `JULIA_NUM_THREADS` is set.
"""
struct CPU <: AbstractCPUArchitecture end

"""
    GPU <: AbstractArchitecture

Run Oceananigans on a single NVIDIA CUDA GPU.
"""
struct GPU <: AbstractGPUArchitecture end

"""
    @hascuda expr

A macro to compile and execute `expr` only if CUDA is installed and available. Generally used to
wrap expressions that can only be compiled if `CuArrays` and `CUDAnative` can be loaded.
"""
macro hascuda(expr)
    return has_cuda() ? :($(esc(expr))) : :(nothing)
end

device(::AbstractCPUArchitecture) = KernelAbstractions.CPU()
device(::AbstractGPUArchitecture) = KernelAbstractions.CUDADevice()

architecture(::Number)  = nothing
architecture(::Array)   = CPU()
architecture(::CuArray) = GPU()

array_type(::CPU) = Array
array_type(::GPU) = CuArray

arch_array(::AbstractCPUArchitecture, A::Array) = A
arch_array(::AbstractCPUArchitecture, A::CuArray) = Array(A)
arch_array(::AbstractGPUArchitecture, A::Array) = CuArray(A)
arch_array(::AbstractGPUArchitecture, A::CuArray) = A

end


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

using Test, Random
using DelimitedFiles

isequaldlm(m1, m2, t) = isequal(m1, m2) && (eltype(m1) == eltype(m2) == t)

@testset "readdlm" begin
    @test isequaldlm(readdlm(IOBuffer("1\t2\n3\t4\n5\t6\n")), [1. 2; 3 4; 5 6], Float64)
    @test isequaldlm(readdlm(IOBuffer("1\t2\n3\t4\n5\t6\n"), Int), [1 2; 3 4; 5 6], Int)
    @test isequaldlm(readdlm(IOBuffer("1,22222222222222222222222222222222222222,0x3,10e6\n2000.1,true,false,-10.34"), ',', Any),
        reshape(Any[1,2000.1,Float64(22222222222222222222222222222222222222),true,0x3,false,10e6,-10.34], 2, 4), Any)

    @test isequaldlm(readdlm(IOBuffer("-9223355253176920979,9223355253176920979"), ',', Int64), Int64[-9223355253176920979  9223355253176920979], Int64)

    @test size(readdlm(IOBuffer("1,2,3,4"), ',')) == (1,4)
    @test size(readdlm(IOBuffer("1,2,3,"), ',')) == (1,4)
    @test size(readdlm(IOBuffer("1,2,3,4\n"), ',')) == (1,4)
    @test size(readdlm(IOBuffer("1,2,3,\n"), ',')) == (1,4)
    @test size(readdlm(IOBuffer("1,2,3,4\n1,2,3,4"), ',')) == (2,4)
    @test size(readdlm(IOBuffer("1,2,3,4\n1,2,3,"), ',')) == (2,4)
    @test size(readdlm(IOBuffer("1,2,3,4\n1,2,3"), ',')) == (2,4)

    @test size(readdlm(IOBuffer("1,2,3,4\r\n"), ',')) == (1,4)
    @test size(readdlm(IOBuffer("1,2,3,4\r\n1,2,3\r\n"), ',')) == (2,4)
    @test size(readdlm(IOBuffer("1,2,3,4\r\n1,2,3,4\r\n"), ',')) == (2,4)
    @test size(readdlm(IOBuffer("1,2,3,\"4\"\r\n1,2,3,4\r\n"), ',')) == (2,4)

    @test size(readdlm(IOBuffer("1 2 3 4\n1 2 3"))) == (2,4)
    @test size(readdlm(IOBuffer("1\t2 3 4\n1 2 3"))) == (2,4)
    @test size(readdlm(IOBuffer("1\t 2 3 4\n1 2 3"))) == (2,4)
    @test size(readdlm(IOBuffer("1\t 2 3 4\n1 2 3\n"))) == (2,4)
    @test size(readdlm(IOBuffer("1,,2,3,4\n1,2,3\n"), ',')) == (2,5)

    let result1 = reshape(Any["", "", "", "", "", "", 1.0, 1.0, "", "", "", "", "", 1.0, 2.0, "", 3.0, "", "", "", "", "", 4.0, "", "", ""], 2, 13),
        result2 = reshape(Any[1.0, 1.0, 2.0, 1.0, 3.0, "", 4.0, ""], 2, 4)

        @test isequaldlm(readdlm(IOBuffer(",,,1,,,,2,3,,,4,\n,,,1,,,1\n"), ','), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("   1    2 3   4 \n   1   1\n")), result2, Any)
        @test isequaldlm(readdlm(IOBuffer("   1    2 3   4 \n   1   1\n"), ' '), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1 2\n3 4 \n")), [[1.0, 3.0] [2.0, 4.0]], Float64)
    end

    let result1 = reshape(Any["", "", "", "", "", "", "भारत", 1.0, "", "", "", "", "", 1.0, 2.0, "", 3.0, "", "", "", "", "", 4.0, "", "", ""], 2, 13)
        @test isequaldlm(readdlm(IOBuffer(",,,भारत,,,,2,3,,,4,\n,,,1,,,1\n"), ',') , result1, Any)
    end

    let result1 = reshape(Any[1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, ""], 2, 4)
        @test isequaldlm(readdlm(IOBuffer("1\t 2 3 4\n1 2 3")), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1\t 2 3 4\n1 2 3 ")), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1\t 2 3 4\n1 2 3\n")), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1,2,3,4\n1,2,3\n"), ','), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1,2,3,4\n1,2,3"), ','), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1,2,3,4\r\n1,2,3\r\n"), ','), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("1,2,3,\"4\"\r\n1,2,3\r\n"), ','), result1, Any)
    end

    let result1 = reshape(Any["abc", "hello", "def,ghi", " \"quote\" ", "new\nline", "world"], 2, 3),
        result2 = reshape(Any["abc", "line\"", "\"hello\"", "\"def", "", "\" \"\"quote\"\" \"", "ghi\"", "", "world", "\"new", "", ""], 3, 4)

        @test isequaldlm(readdlm(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ','), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ',', quotes=false), result2, Any)
    end

    let result1 = reshape(Any["t", "c", "", "c"], 2, 2),
        result2 = reshape(Any["t", "\"c", "t", "c"], 2, 2)
        @test isequaldlm(readdlm(IOBuffer("t  \n\"c\" c")), result1, Any)
        @test isequaldlm(readdlm(IOBuffer("t t \n\"\"\"c\" c")), result2, Any)
    end

    @test isequaldlm(readdlm(IOBuffer("\n1,2,3\n4,5,6\n\n\n"), ',', skipblanks=false),
                     reshape(Any["",1.0,4.0,"","","",2.0,5.0,"","","",3.0,6.0,"",""], 5, 3), Any)
    @test isequaldlm(readdlm(IOBuffer("\n1,2,3\n4,5,6\n\n\n"), ',', skipblanks=true), reshape([1.0,4.0,2.0,5.0,3.0,6.0], 2, 3), Float64)
    @test isequaldlm(readdlm(IOBuffer("1,2\n\n4,5"), ',', skipblanks=false), reshape(Any[1.0,"",4.0,2.0,"",5.0], 3, 2), Any)
    @test isequaldlm(readdlm(IOBuffer("1,2\n\n4,5"), ',', skipblanks=true), reshape([1.0,4.0,2.0,5.0], 2, 2), Float64)

    let x = bitrand(5, 10), io = IOBuffer()
        writedlm(io, x)
        seek(io, 0)
        @test readdlm(io, Bool) == x
    end

    let x = [1,2,3], y = [4,5,6], io = IOBuffer()
        writedlm(io, zip(x,y), ",  ")
        seek(io, 0)
        @test readdlm(io, ',') == [x y]
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm(io, x, ", ")
        seek(io, 0)
        @test read(io, String) == "0.1, 0.3, 0.5\n"
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm(io, x, ", ")
        seek(io, 0)
        @test readdlm(io, ',') == [0.1 0.3 0.5]
    end

    let x = ["abc", "def\"ghi", "jk\nl"], y = [1, ",", "\"quoted\""], io = IOBuffer()
        writedlm(io, zip(x,y), ',')
        seek(io, 0)
        @test readdlm(io, ',') == [x y]
    end

    let x = ["a" "b"; "d" ""], io = IOBuffer()
        writedlm(io, x)
        seek(io, 0)
        @test readdlm(io) == x
    end

    let x = ["\"hello\"", "world\""], io = IOBuffer()
        writedlm(io, x, quotes=false)
        @test String(take!(io)) == "\"hello\"\nworld\"\n"

        writedlm(io, x)
        @test String(take!(io)) == "\"\"\"hello\"\"\"\n\"world\"\"\"\n"
    end
end

@testset "comments" begin
    @test isequaldlm(readdlm(IOBuffer("#this is comment\n1,2,3\n#one more comment\n4,5,6"), ',', comments=true), [1. 2. 3.;4. 5. 6.], Float64)
    @test isequaldlm(readdlm(IOBuffer("#this is \n#comment\n1,2,3\n#one more \n#comment\n4,5,6"), ',', comments=true), [1. 2. 3.;4. 5. 6.], Float64)
    @test isequaldlm(readdlm(IOBuffer("1,2,#3\n4,5,6"), ',', comments=true), [1. 2. "";4. 5. 6.], Any)
    @test isequaldlm(readdlm(IOBuffer("1#,2,3\n4,5,6"), ',', comments=true), [1. "" "";4. 5. 6.], Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,\"#3\"\n4,5,6"), ',', comments=true), [1. 2. "#3";4. 5. 6.], Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,3\n #with leading whitespace\n4,5,6"), ',', comments=true), [1. 2. 3.;" " "" "";4. 5. 6.], Any)
end

@testset "without comments" begin
    @test isequaldlm(readdlm(IOBuffer("1,2,#3\n4,5,6"), ','), [1. 2. "#3";4. 5. 6.], Any)
    @test isequaldlm(readdlm(IOBuffer("1#,2,3\n4,5,6"), ','), ["1#" 2. 3.;4. 5. 6.], Any)
    @test isequaldlm(readdlm(IOBuffer("1,2,\"#3\"\n4,5,6"), ','), [1. 2. "#3";4. 5. 6.], Any)
end

@testset "skipstart" begin
    x = ["a" "b" "c"; "d" "e" "f"; "g" "h" "i"; "A" "B" "C"; 1 2 3; 4 5 6; 7 8 9]
    io = IOBuffer()

    writedlm(io, x, quotes=false)
    seek(io, 0)
    (data, hdr) = readdlm(io, header=true, skipstart=3)
    @test data == [1 2 3; 4 5 6; 7 8 9]
    @test hdr == ["A" "B" "C"]

    x = ["a" "b" "\nc"; "d" "\ne" "f"; "g" "h" "i\n"; "A" "B" "C"; 1 2 3; 4 5 6; 7 8 9]
    io = IOBuffer()

    writedlm(io, x, quotes=true)
    seek(io, 0)
    (data, hdr) = readdlm(io, header=true, skipstart=6)
    @test data == [1 2 3; 4 5 6; 7 8 9]
    @test hdr == ["A" "B" "C"]

    io = IOBuffer()
    writedlm(io, x, quotes=false)
    seek(io, 0)
    (data, hdr) = readdlm(io, header=true, skipstart=6)
    @test data == [1 2 3; 4 5 6; 7 8 9]
    @test hdr == ["A" "B" "C"]
end

@testset "i18n" begin
    # source: http://www.i18nguy.com/unicode/unicode-example-utf8.zip
    let i18n_data = ["Origin (English)", "Name (English)", "Origin (Native)", "Name (Native)",
            "Australia", "Nicole Kidman", "Australia", "Nicole Kidman",
            "Austria", "Johann Strauss", "Österreich", "Johann Strauß",
            "Belgium (Flemish)", "Rene Magritte", "België", "René Magritte",
            "Belgium (French)", "Rene Magritte", "Belgique", "René Magritte",
            "Belgium (German)", "Rene Magritte", "Belgien", "René Magritte",
            "Bhutan", "Gonpo Dorji", "འབྲུག་ཡུལ།", "མགོན་པོ་རྡོ་རྗེ།",
            "Canada", "Celine Dion", "Canada", "Céline Dion",
            "Canada - Nunavut (Inuktitut)", "Susan Aglukark", "ᓄᓇᕗᒻᒥᐅᑦ", "ᓱᓴᓐ ᐊᒡᓗᒃᑲᖅ",
            "Democratic People's Rep. of Korea", "LEE Sol-Hee", "조선 민주주의 인민 공화국", "이설희",
            "Denmark", "Soren Hauch-Fausboll", "Danmark", "Søren Hauch-Fausbøll",
            "Denmark", "Soren Kierkegaard", "Danmark", "Søren Kierkegård",
            "Egypt", "Abdel Halim Hafez", "ﻣﺼﺮ", "ﻋﺑﺪﺍﻠﺣﻟﻳﻢ ﺤﺎﻓﻅ",
            "Egypt", "Om Kolthoum", "ﻣﺼﺮ", "ﺃﻡ ﻛﻟﺛﻭﻡ",
            "Eritrea", "Berhane Zeray", "ብርሃነ ዘርኣይ", "ኤርትራ",
            "Ethiopia", "Haile Gebreselassie", "ኃይሌ ገብረሥላሴ", "ኢትዮጵያ",
            "France", "Gerard Depardieu", "France", "Gérard Depardieu",
            "France", "Jean Reno", "France", "Jean Réno",
            "France", "Camille Saint-Saens", "France", "Camille Saint-Saëns",
            "France", "Mylene Demongeot", "France", "Mylène Demongeot",
            "France", "Francois Truffaut", "France", "François Truffaut",
            "France (Braille)", "Louis Braille", "⠋⠗⠁⠝⠉⠑", "⠇⠕⠥⠊⠎⠀<BR>⠃⠗⠁⠊⠇⠇⠑",
            "Georgia", "Eduard Shevardnadze", "საქართველო", "ედუარდ შევარდნაძე",
            "Germany", "Rudi Voeller", "Deutschland", "Rudi Völler",
            "Germany", "Walter Schultheiss", "Deutschland", "Walter Schultheiß",
            "Greece", "Giorgos Dalaras", "Ελλάς", "Γιώργος Νταλάρας",
            "Iceland", "Bjork Gudmundsdottir", "Ísland", "Björk Guðmundsdóttir",
            "India (Hindi)", "Madhuri Dixit", "भारत", "माधुरी दिछित",
            "Ireland", "Sinead O'Connor", "Éire", "Sinéad O'Connor",
            "Israel", "Yehoram Gaon", "ישראל", "יהורם גאון",
            "Italy", "Fabrizio DeAndre", "Italia", "Fabrizio De André",
            "Japan", "KUBOTA Toshinobu", "日本", "久保田    利伸",
            "Japan", "HAYASHIBARA Megumi", "日本", "林原 めぐみ",
            "Japan", "Mori Ogai", "日本", "森鷗外",
            "Japan", "Tex Texin", "日本", "テクス テクサン",
            "Norway", "Tor Age Bringsvaerd", "Noreg", "Tor Åge Bringsværd",
            "Pakistan (Urdu)", "Nusrat Fatah Ali Khan", "پاکستان", "نصرت فتح علی خان",
            "People's Rep. of China", "ZHANG Ziyi", "中国", "章子怡",
            "People's Rep. of China", "WONG Faye", "中国", "王菲",
            "Poland", "Lech Walesa", "Polska", "Lech Wałęsa",
            "Puerto Rico", "Olga Tanon", "Puerto Rico", "Olga Tañón",
            "Rep. of China", "Hsu Chi", "臺灣", "舒淇",
            "Rep. of China", "Ang Lee", "臺灣", "李安",
            "Rep. of Korea", "AHN Sung-Gi", "대한민국", "안성기",
            "Rep. of Korea", "SHIM Eun-Ha", "대한민국", "심은하",
            "Russia", "Mikhail Gorbachev", "Россия", "Михаил Горбачёв",
            "Russia", "Boris Grebenshchikov", "Россия", "Борис Гребенщиков",
            "Slovenia", "\"Frane \"\"Jezek\"\" Milcinski", "Slovenija", "Frane Milčinski - Ježek",
            "Syracuse (Sicily)", "Archimedes", "Συρακούσα", "Ἀρχιμήδης",
            "Thailand", "Thongchai McIntai", "ประเทศไทย", "ธงไชย แม็คอินไตย์",
            "U.S.A.", "Brad Pitt", "U.S.A.", "Brad Pitt",
            "Yugoslavia (Cyrillic)", "Djordje Balasevic", "Југославија", "Ђорђе Балашевић",
            "Yugoslavia (Latin)", "Djordje Balasevic", "Jugoslavija", "Đorđe Balašević"]

        i18n_arr = permutedims(reshape(i18n_data, 4, Int(floor(length(i18n_data)/4))), [2, 1])
        i18n_buff = PipeBuffer()
        writedlm(i18n_buff, i18n_arr, ',')
        @test i18n_arr == readdlm(i18n_buff, ',')

        hdr = i18n_arr[1:1, :]
        data = i18n_arr[2:end, :]
        writedlm(i18n_buff, i18n_arr, ',')
        @test (data, hdr) == readdlm(i18n_buff, ',', header=true)

        writedlm(i18n_buff, i18n_arr, '\t')
        @test (data, hdr) == readdlm(i18n_buff, '\t', header=true)
    end
end

@testset "issue #13028" begin
    for data in ["A B C", "A B C\n"]
        data,hdr = readdlm(IOBuffer(data), header=true)
        @test hdr == AbstractString["A" "B" "C"]
        @test data == Matrix{Float64}(undef, 0, 3)
    end
end

# fix #13179 parsing unicode lines with default delmiters
@test isequaldlm(readdlm(IOBuffer("# Should ignore this π\n1\tα\n2\tβ\n"), comments=true), Any[1 "α"; 2 "β"], Any)

# BigInt parser
let data = "1 2 3"
    readdlm(IOBuffer(data), ' ', BigInt) == BigInt[1 2 3]
end

@testset "show with MIME types" begin
    @test sprint(show, "text/csv", [1 2; 3 4]) == "1,2\n3,4\n"

    for writefunc in ((io,x) -> show(io, "text/csv", x),
                      (io,x) -> invoke(writedlm, Tuple{IO,Any,Any}, io, x, ","))
        # iterable collections of iterable rows:
        let x = [(1,2), (3,4)], io = IOBuffer()
            writefunc(io, x)
            seek(io, 0)
            @test readdlm(io, ',') == [1 2; 3 4]
        end
        # vectors of strings:
        let x = ["foo", "bar"], io = IOBuffer()
            writefunc(io, x)
            seek(io, 0)
            @test vec(readdlm(io, ',')) == x
        end
    end
end

# Test that we can read a write protected file
let fn = tempname()
    open(fn, "w") do f
        write(f, "Julia")
    end
    chmod(fn, 0o444)
    readdlm(fn)[] == "Julia"
    rm(fn)
end

# issue #21180
let data = "\"721\",\"1438\",\"1439\",\"…\",\"1\""
    @test readdlm(IOBuffer(data), ',') == Any[721  1438  1439  "…"  1]
end

# issue #21207
let data = "\"1\",\"灣\"\"灣灣灣灣\",\"3\""
    @test readdlm(IOBuffer(data), ',') == Any[1 "灣\"灣灣灣灣" 3]
end

# reading from a byte array (#16731)
let data = Vector{UInt8}("1,2,3\n4,5,6"), origdata = copy(data)
    @test readdlm(data, ',') == [1 2 3; 4 5 6]
    @test data == origdata
end

# issue #11484: useful error message for invalid readdlm filepath arguments
@test_throws ArgumentError readdlm(tempdir())

# displaying as text/csv
let d = TextDisplay(IOBuffer())
    display(d, "text/csv", [3 1 4])
    @test String(take!(d.io)) == "3,1,4\n"
end

@testset "complex" begin
    @test readdlm(IOBuffer("3+4im, 4+5im"), ',', Complex{Int}) == [3+4im 4+5im]
end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
module MyExample

    include("extra_file.jl")

    export my_f
end


text:
# Autogenerated wrapper script for Minuit2_jll for i686-linux-musl-cxx03
export libMinuit2, libMinuit2Math

JLLWrappers.@generate_wrapper_header("Minuit2")
JLLWrappers.@declare_library_product(libMinuit2, "libMinuit2.so")
JLLWrappers.@declare_library_product(libMinuit2Math, "libMinuit2Math.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libMinuit2,
        "lib/libMinuit2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libMinuit2Math,
        "lib/libMinuit2Math.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Parts of this code were taken / derived from NetworkX. See LICENSE for
# licensing details.

"""
    pagerank(g, α=0.85, n=100, ϵ=1.0e-6)

Calculate the [PageRank](https://en.wikipedia.org/wiki/PageRank) of the
directed graph `g` parameterized by damping factor `α`, number of
iterations `n`, and convergence threshold `ϵ`. Return a vector representing
the centrality calculated for each node in `g`, or an error if convergence
is not reached within `n` iterations.
"""
function pagerank end

@traitfn function pagerank(g::::IsDirected, α=0.85, n=100, ϵ=1.0e-6)
    A = adjacency_matrix(g,:in,Float64)
    S = vec(sum(A,1))
    S = 1./S
    S[find(S .== Inf)]=0.0
    M = A' # need a separate line due to bug #17456 in julia
    M = (Diagonal(S) * M)'
    N = Int(nv(g))
    x = fill(1.0/N, N)
    p = fill(1.0/N, N)
    dangling_weights = p
    is_dangling = find(S .== 0)

    for _ in 1:n
        xlast = x
        x = α * (M * x + sum(x[is_dangling]) * dangling_weights) + (1 - α) * p
        err = sum(abs, (x - xlast))
        if (err < N * ϵ)
            return x
        end
    end
    error("Pagerank did not converge after $n iterations.")
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# Combined (lower and upper) bounds type for the RockSample problem
include("rockSampleParticleLB.jl")
include("rockSampleFringeUB.jl")
include("../../upperBound/upperBoundNonStochastic.jl")

import DESPOT: bounds, init_bounds, default_action

mutable struct RockSampleBounds
    lb::RockSampleParticleLB
    ub::UpperBoundNonStochastic{RockSampleState,RockSampleAction,RockSampleObs}

    function RockSampleBounds(pomdp::RockSample)
        this = new()
        this.lb = RockSampleParticleLB(pomdp)
        this.ub = UpperBoundNonStochastic{RockSampleState,RockSampleAction,RockSampleObs}(pomdp)

        return this
    end
end

function bounds(b::RockSampleBounds,
       pomdp::RockSample,
       particles::Vector{DESPOTParticle{RockSampleState}},
       config::DESPOTConfig)

    ubound::Float64 = upper_bound(b.ub, pomdp, particles, config)
    lbound::Float64 = lower_bound(b.lb, pomdp, particles, b.ub.upper_bound_act, config)

    return lbound, ubound
end

function init_bounds(bounds::RockSampleBounds, pomdp::RockSample, config::DESPOTConfig)
    init_bound(bounds.ub, pomdp, config)
    # lower bound init not currently needed
end

function default_action(b::RockSampleBounds, pomdp::RockSample, particles::Vector, config::DESPOTConfig)
    bounds(b, pomdp, particles, config)
    return b.lb.best_action
end


text:
@testset "Memento.TestUtils" begin
    @testset "@test_log" begin
        logger = getlogger("test_log")
        msg = "Hello!"
        setlevel!(logger, "info")
        @test_log(logger, "info", msg, Memento.info(logger, msg))

        # Partial matches
        msg = "Hello World!"
        @test_log(logger, "info", "Hello", Memento.info(logger, msg))
        @test_log(logger, "info", r"Hello*", Memento.info(logger, msg))
        @test_log(logger, "info", ("Hello", "World"), Memento.info(logger, msg))
        @test_log(logger, "info", x -> x == msg, Memento.info(logger, msg))
    end

    @testset "@test_nolog" begin
        logger = getlogger("test_nolog")
        setlevel!(logger, "info")

        msg = "Hello!"
        different_msg = "Goodbye"

        @test_nolog(logger, "info", msg, nothing)
        @test_nolog(logger, "info", msg, Memento.debug(logger, msg))
        @test_nolog(logger, "info", msg, Memento.info(logger, different_msg))
        @test_nolog(logger, "info", msg, Memento.warn(logger, different_msg))
        @test_nolog(logger, "info", r".*foo.*", Memento.info(logger, "baz"))
    end

    @testset "@test_warn" begin
        logger = getlogger("test_warn")
        msg = "Hello!"
        @test_warn(logger, msg, Memento.warn(logger, "Hello!"))
    end

    @testset "@test_throws" begin
        logger = getlogger("test_log")
        @test_throws(logger, ErrorException, error(logger, "Error!"))
    end
end


text:
// Testing the return checker

int f () {
   if (true)
     return 0;
   else
     {}
}

int g () {
  if (false)
      {}
  else
      return 0;
}

void p () {}

int main() {
  p();
  return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
module FirebaseyTestPlace end
@testset "Firebasey" begin
    Core.eval(FirebaseyTestPlace, quote
        using Test
        include("../src/webserver/Firebasey.jl")
    end)
end

module DataUrlTestPlace end
@testset "DataUrl" begin
    Core.eval(DataUrlTestPlace, quote
        using Test
        include("../src/webserver/data_url.jl")
    end)
end


text:
using BlockArrays, LinearAlgebra, Base64

struct PartiallyImplementedBlockVector <: AbstractBlockArray{Float64,1} end

@testset "partially implemented block array" begin
    # the error thrown before was incorrect
    A = PartiallyImplementedBlockVector()
    @test_throws(ErrorException("getblock for PartiallyImplementedBlockVector is not implemented"),
                getblock(A, 1))
    @test_throws(ErrorException("getblock! for PartiallyImplementedBlockVector is not implemented"),
        getblock!(zeros(5), A, Block(1)))
    @test_throws(ErrorException("setblock! for PartiallyImplementedBlockVector is not implemented"),
        BlockArrays.setblock!(A, zeros(5), Block(1)))
    @test_throws(ErrorException( "axes for PartiallyImplementedBlockVector is not implemented"),
        BlockArrays.axes(A))
end

@testset "Array block interface" begin
    @test 1[Block()] == 1

    A = randn(5)
    @test blocksize(A) == (1,)
    @test A[Block(1)] == A
    view(A,Block(1))[1] = 2
    @test A[1] == 2
    @test_throws BlockBoundsError A[Block(2)]

    A = randn(5,5)
    @test A[Block(1,1)] == A
end

@testset "Triangular/Symmetric/Hermitian block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (1:4))
    A .= reshape(1:length(A), size(A))

    @test blocksize(UpperTriangular(A)) == blocksize(Symmetric(A)) == blocksize(A)
    @test UpperTriangular(A)[Block(2,2)] == UpperTriangular(A[2:3,2:3])
    @test UpperTriangular(A)[Block(2,3)] == A[2:3,4:6]
    @test UpperTriangular(A)[Block(3,2)] == zeros(3,2)
    @test Symmetric(A)[Block(2,2)] == Symmetric(A[2:3,2:3])
    @test Symmetric(A)[Block(2,3)] == A[2:3,4:6]
    @test Symmetric(A)[Block(3,2)] == transpose(A[2:3,4:6])
    @test Hermitian(A)[Block(2,2)] == Hermitian(A[2:3,2:3])
    @test Hermitian(A)[Block(2,3)] == A[2:3,4:6]
    @test Hermitian(A)[Block(3,2)] == A[2:3,4:6]'
    if VERSION ≥ v"1.2"
        @test stringmime("text/plain", UpperTriangular(A)) == "10×10 UpperTriangular{Complex{Float64},PseudoBlockArray{Complex{Float64},2,Array{Complex{Float64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:10×1:1:10:\n 1.0+0.0im  │  11.0+0.0im  21.0+0.0im  │  31.0+0.0im  41.0+0.0im  51.0+0.0im  │  61.0+0.0im  71.0+0.0im  81.0+0.0im   91.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │  12.0+0.0im  22.0+0.0im  │  32.0+0.0im  42.0+0.0im  52.0+0.0im  │  62.0+0.0im  72.0+0.0im  82.0+0.0im   92.0+0.0im\n     ⋅      │       ⋅      23.0+0.0im  │  33.0+0.0im  43.0+0.0im  53.0+0.0im  │  63.0+0.0im  73.0+0.0im  83.0+0.0im   93.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │  34.0+0.0im  44.0+0.0im  54.0+0.0im  │  64.0+0.0im  74.0+0.0im  84.0+0.0im   94.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅      45.0+0.0im  55.0+0.0im  │  65.0+0.0im  75.0+0.0im  85.0+0.0im   95.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅      56.0+0.0im  │  66.0+0.0im  76.0+0.0im  86.0+0.0im   96.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │  67.0+0.0im  77.0+0.0im  87.0+0.0im   97.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅      78.0+0.0im  88.0+0.0im   98.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅      89.0+0.0im   99.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅           ⋅      100.0+0.0im"
    end
end

@testset "Adjoint/Transpose block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (2:5))
    A .= randn.() .+ randn.().*im

    @test blocksize(A') == (4,4)
    @test blocksize(Transpose(A)) == (4,4)

    @test A'[Block(2,2)] == A[Block(2,2)]' == A[2:3,3:5]'
    @test transpose(A)[Block(2,2)] == transpose(A[2:3,3:5])
    @test A'[Block(2,3)] == A[Block(3,2)]'
    @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

    @test BlockArray(A') == A'
    @test BlockArray(transpose(A)) == transpose(A)
end

@testset "Diagonal BlockArray" begin
    A = mortar(Diagonal(fill([1 2],2)))
    @test A isa BlockMatrix{Int,Diagonal{Matrix{Int}, Vector{Matrix{Int}}}}
    @test A[Block(1,2)] == [0 0]
    @test_throws BlockBoundsError A[Block(1,3)]
    @test A == [1 2 0 0; 0 0 1 2]
end

@testset "non-standard block axes" begin
    A = BlockArray([1 2; 3 4], Fill(1,2),Fill(1,2))
    @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},NTuple{2,BlockedUnitRange{StepRange{Int64,Int64}}}}
    A = BlockArray([1 2; 3 4], Fill(1,2),[1,1])
    @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},Tuple{BlockedUnitRange{StepRange{Int64,Int64}},BlockedUnitRange{Vector{Int}}}}
end


text:
using Documenter, GPUifyLoops

makedocs(
    modules = [GPUifyLoops],
    format = :html,
    sitename = "GPUifyLoops.jl",
    pages = [
        "Home"    => "index.md",
    ],
    doctest = true
)


text:
using LoggingFacilities
using Test

using Logging
using JSON
using Dates

@testset "LoggingFacilities.jl" begin

    log1 = (level = Logging.Info, message = "hello", kwargs = ())
    log2 = (level = Logging.Info, message = "", kwargs = (:a => 1,))
    log3 = (level = Logging.Info, message = "hello", kwargs = (:a => 1, :b => 2))

    @testset "Inject" begin

        @testset "Message Location" begin
            # test all locations
            @test inject(log1, BeginningMessageLocation(), 1).message == "1 hello"
            @test inject(log1, EndingMessageLocation(), 1).message    == "hello 1"
            @test inject(log1, WholeMessageLocation(), 1).message     == "1"

            # message separator for beginning/ending message location
            @test inject(log1, BeginningMessageLocation(), 1; sep="").message == "1hello"
            @test inject(log1, EndingMessageLocation(), 1; sep="").message    == "hello1"

            # callable
            @test inject(log1, WholeMessageLocation(), () -> 2).message == "2"
        end

        @testset "Kwargs Location" begin
            # inject into empty kwargs (non-callable and callable)
            @test inject(log1, KwargsLocation(), [:a => 1]).kwargs == (:a => 1,)
            @test inject(log1, KwargsLocation(), () -> [:a => 1]).kwargs == (:a => 1,)

            # inject into non-empty kwargs
            @test inject(log2, KwargsLocation(), [:b => 2]).kwargs |> length == 2

            # merged with same key (expectation: does not dedup)
            @test inject(log2, KwargsLocation(), [:a => 2]).kwargs |> length == 2
        end

        @testset "Level Location" begin
            @test inject(log1, LevelLocation(), Logging.Error).level == Logging.Error
            @test inject(log1, LevelLocation(), () -> Logging.Error).level == Logging.Error
        end
    end

    @testset "Remove" begin
        @test remove(log1, MessageProperty()).message == ""
        @test remove(log1, KwargsProperty()).kwargs |> length == 0
        @test_throws ErrorException remove(log1, LevelProperty())
    end

    @testset "Migrate" begin
        # message => kwargs
        @test migrate(log1, MessageProperty(), KwargsProperty()).kwargs[1] |> first == :message
        @test migrate(log1, MessageProperty(), KwargsProperty()).kwargs[1] |> last  == "hello"

        # kwargs => message
        @test migrate(log2, KwargsProperty(), MessageProperty()).message == "a=1"
        @test migrate(log2, KwargsProperty(), MessageProperty()).message == "a=1"
        @test migrate(log3, KwargsProperty(), MessageProperty()).message == "hello a=1 b=2"
        @test migrate(log3, KwargsProperty(), MessageProperty(); sep=":").message == "hello a:1 b:2"
        @test migrate(log3, KwargsProperty(), MessageProperty(); divider="_").message == "hello a=1_b=2"
        @test migrate(log3, KwargsProperty(), MessageProperty(); prepend="|").message == "hello|a=1 b=2"
        @test migrate(log3, KwargsProperty(), MessageProperty(); transform = string).message == "hello (:a => 1, :b => 2)"

        # level => message
        @test migrate(log1, LevelProperty(), MessageProperty()).message == "Info hello"
        @test migrate(log1, LevelProperty(), MessageProperty(); location = BeginningMessageLocation()).message == "Info hello"
        @test migrate(log1, LevelProperty(), MessageProperty(); location = EndingMessageLocation()).message == "hello Info"
        @test migrate(log1, LevelProperty(), MessageProperty(); transform = (level) -> uppercase(string(level))).message == "INFO hello"

        # level => kwargs
        @test migrate(log1, LevelProperty(), KwargsProperty()).kwargs == (:level => Logging.Info,)
        @test migrate(log1, LevelProperty(), KwargsProperty(); label = :_level).kwargs == (:_level => Logging.Info,)
        @test migrate(log1, LevelProperty(), KwargsProperty(); transform = string).kwargs == (:level => "Info",)
    end

    @testset "Mutate" begin
        escalate_info_to_warn(log) = log.level == Logging.Info ? Logging.Warn : log.level
        @test mutate(log1, LevelProperty(); transform = escalate_info_to_warn).level == Logging.Warn

        quote_me(log) = "\"$(log.message)\""
        @test mutate(log1, MessageProperty(); transform = quote_me).message == "\"hello\""

        replace_b_with_c(log) = tuple(collect(((k == :b ? :c : k) => v) for (k,v) in log.kwargs)...)
        @test mutate(log3, KwargsProperty(); transform = replace_b_with_c).kwargs == (:a => 1, :c => 2)
    end

    # The following tests must use `current_logger` such that the results can be collected
    # by `Test.collect_test_logs`.

    @testset "One Line Logger" begin
        logs, value = Test.collect_test_logs() do
            with_logger(OneLineTransformerLogger(current_logger())) do
                x = 1
                y = 2
                @info "hello" x y
            end
        end
        @test logs[1].message == "hello x=1 y=2"
        @test logs[1].kwargs |> length == 0
        @test logs[1].level == Logging.Info
    end

    @testset "Timestamp Logger" begin
        logs, value = Test.collect_test_logs() do
            with_logger(TimestampTransformerLogger(current_logger(), BeginningMessageLocation();
                                                   format = "yyyy-mm-dd")) do
                @info "hello"
            end
            with_logger(TimestampTransformerLogger(current_logger(), EndingMessageLocation();
                                                   format = "yyyy-mm-dd")) do
                @info "hello"
            end
            with_logger(TimestampTransformerLogger(current_logger(), KwargsLocation();
                                                   format = "yyyy-mm-dd")) do
                @info "hello"
            end
        end
        is_date_format(s) = match(r"\d\d\d\d-\d\d-\d\d", s) !== nothing
        @test logs[1].message[1:10]      |> is_date_format
        @test logs[2].message[end-9:end] |> is_date_format
        @test first(logs[3].kwargs) |> first == :timestamp
        @test first(logs[3].kwargs) |> last |> is_date_format
    end

    @testset "JSON Logger" begin
        logs, value = Test.collect_test_logs() do
            with_logger(() -> @info("hello"), JSONTransformerLogger(current_logger()))
            with_logger(() -> @info("hello"), JSONTransformerLogger(current_logger(); level_label = :LEVEL))
            with_logger(() -> @info("hello"), JSONTransformerLogger(current_logger(); message_label = :MSG))
            with_logger(() -> @info("hello"), JSONTransformerLogger(current_logger(); indent = 2))
        end
        # parsable
        foreach(log -> @test_nowarn(JSON.Parser.parse(log.message)), logs)

        let d = JSON.Parser.parse(logs[1].message)
            @test d["message"] == "hello"
            @test d["level"] == "Info"
        end

        let d = JSON.Parser.parse(logs[2].message)
            @test d["message"] == "hello"
            @test d["LEVEL"] == "Info"
        end

        let d = JSON.Parser.parse(logs[3].message)
            @test d["MSG"] == "hello"
            @test d["level"] == "Info"
        end

        @test split(logs[4].message, "\n") |> length == 4
    end

    @testset "Color message" begin
        colors = Dict(Logging.Info => ColorSpec(:red, false))
        logs, value = Test.collect_test_logs() do
            clogger = ColorMessageTransformerLogger(colors)(current_logger())
            with_logger(() -> @info("hello"), clogger)
        end
        @test length(logs[1].message) > 5    # since it has additional chars for switching colors
    end

    @testset "Fixed width message" begin
        logs, value = Test.collect_test_logs() do
            with_logger(() -> @info("hello"), FixedMessageWidthTransformerLogger(current_logger(), 30))
        end
        @test length(logs[1].message) == 30
    end

    # Composition
    @testset "Compose transformers" begin
        logs, value = Test.collect_test_logs() do
            logger = compose(
                current_logger(),
                TimestampTransformerLogger(BeginningMessageLocation()),
                OneLineTransformerLogger
            )
            with_logger(logger) do
                x = 1
                @info "hello" x
            end
        end
        @test match(r"\d\d\d\d.* hello x=1", logs[1].message) !== nothing
    end

    # Sinks
    @testset "MessageOnlyLogger" begin
        io = IOBuffer()
        with_logger(MessageOnlyLogger(io)) do
            x = 1
            @info "hello"
        end
        @test String(take!(io)) == "hello\n"  # no level, no kwargs, contains newline
    end

end