text:
@testset "Indexing" begin
    @test maxtransformlevels(randn(4,2), 1) == 2
    @test maxtransformlevels(randn(4,2), 2) == 1
    @test_throws AssertionError maxtransformlevels(randn(4,2), 3)

    Xw = collect(1:12) |> x -> reshape(x,(4,3))
    @test getbasiscoef(Xw, maketree(4,2,:dwt)) == [9,10,7,8]
    @test getbasiscoef(Xw, maketree(4,2,:full)) == [9,10,11,12]
    @test_throws AssertionError getbasiscoef(Xw, BitVector([0,1,0]))
    @test_throws AssertionError getbasiscoef(Xw, BitVector([1,1,1,1]))
    @test_throws AssertionError getbasiscoef(Xw, maketree(8,2,:dwt))
    @test_throws AssertionError getbasiscoef(zeros(4,4), maketree(4,2,:dwt))
    @test_throws ArgumentError getbasiscoef(zeros(4,2), maketree(4,2,:dwt))

    Xw = collect(1:24) |> x -> reshape(x,(4,3,2))
    tree₁ = maketree(4,2,:dwt)
    tree₂ = maketree(4,2,:full)
    @test getbasiscoefall(Xw, tree₁) == [9 21;10 22;7 19;8 20]
    @test getbasiscoefall(Xw, tree₂) == [9 21;10 22;11 23;12 24]
    @test getbasiscoefall(Xw, [tree₁ tree₂]) == [9 21;10 22;7 23;8 24]
    @test_throws AssertionError getbasiscoefall(Xw, [tree₁ tree₂ tree₁])
    @test_throws AssertionError getbasiscoefall(Xw, BitVector([0,1,0]))
    @test_throws AssertionError getbasiscoefall(Xw, BitVector([1,1,1,1]))
    @test_throws AssertionError getbasiscoefall(Xw, maketree(8,2,:dwt))
    @test_throws AssertionError getbasiscoefall(zeros(4,4,2), maketree(4,2,:dwt))
    @test_throws ArgumentError getbasiscoefall(zeros(4,2,2), maketree(4,2,:dwt))
    @test_throws AssertionError getbasiscoefall(Xw, [tree₁ tree₂ tree₁])
    @test_throws AssertionError getbasiscoefall(Xw, BitMatrix([0 0;1 0;0 1]))
    @test_throws AssertionError getbasiscoefall(Xw, BitMatrix([1 1;1 1;1 1;1 1]))
    @test_throws AssertionError getbasiscoefall(Xw, [maketree(8,2,:dwt) maketree(8,2,:dwt)])
    @test_throws AssertionError getbasiscoefall(zeros(4,4,2), BitMatrix([1 1;1 1;1 1]))
    @test_throws ArgumentError getbasiscoefall(zeros(4,2,2), BitMatrix([1 1;1 1;1 1]))

    @test nodelength(8,2) == 2

    @test Utils.main2depthshift(10,4) == [0,0,2,2,10]
    @test Utils.main2depthshift(5,5) == [0,1,1,5,5,5]
    @test_throws AssertionError Utils.main2depthshift(8,3)
    @test_throws AssertionError Utils.main2depthshift(8,2)

    # TODO: @test packet()

    tree = maketree(4,1,:dwt)
    @test coarsestscalingrange(zeros(4), tree) == 1:2
    @test coarsestscalingrange(4, tree) == 1:2
    @test coarsestscalingrange(zeros(4,3), tree, true) == (1:4, 2)
    @test coarsestscalingrange(4, tree, true) == (1:4, 2)
    @test_throws AssertionError coarsestscalingrange(5, tree, true)
    @test_throws AssertionError coarsestscalingrange(5, tree, false)

    @test finestdetailrange(zeros(4), tree) == 3:4
    @test finestdetailrange(4, tree) == 3:4
    @test finestdetailrange(zeros(4,3), tree, true) == (1:4, 3)
    @test finestdetailrange(4, tree, true) == (1:4, 3)
    @test_throws AssertionError finestdetailrange(5, tree, true)
    @test_throws AssertionError finestdetailrange(5, tree, false)

    @test getrowrange(8,2) == 1:4
    @test getrowrange(8,3) == 1:4
    @test getrowrange(8,4) == 5:8
    @test getrowrange(8,5) == 5:8
    @test_throws AssertionError getrowrange(8,86)

    @test getcolrange(8,2) == 1:4
    @test getcolrange(8,3) == 5:8
    @test getcolrange(8,4) == 1:4
    @test getcolrange(8,5) == 5:8
    @test_throws AssertionError getcolrange(8,86)
end

@testset "Traverse Trees" begin
    @test isvalidtree(zeros(4,4), maketree(4,4,2,:full))
    @test isvalidtree(zeros(4,4), maketree(4,4,2,:dwt))
    @test !isvalidtree(zeros(4,4), BitVector([0,1,0,0,0]))
    @test !isvalidtree(zeros(4,4), trues(4))

    @test getchildindex(3,:left) == 6
    @test getchildindex(3,:right) == 7
    @test getchildindex(3,:topleft) == 10
    @test getchildindex(3,:topright) == 11
    @test getchildindex(3,:bottomleft) == 12
    @test getchildindex(3,:bottomright) == 13
    @test_throws AssertionError getchildindex(3,:fail)

    @test getparentindex(4,:binary) == 2
    @test getparentindex(5,:binary) == 2
    @test getparentindex(10,:quad) == 3
    @test getparentindex(11,:quad) == 3
    @test getparentindex(12,:quad) == 3
    @test getparentindex(13,:quad) == 3
    @test_throws AssertionError getparentindex(15, :fail)

    qleaf = falses(21); qleaf[[3,4,5,6,7,8,9]] .= 1
    @test getleaf(maketree(4,2,:dwt),:binary) == BitVector([0,0,1,1,1,0,0])
    @test getleaf(maketree(4,4,2,:dwt),:quad) == qleaf
    @test_throws AssertionError getleaf(maketree(4,2,:dwt),:fail)
    @test_throws AssertionError getleaf(maketree(4,4,2,:dwt),:binary)
    @test_throws AssertionError getleaf(BitVector([0,1,0]),:binary)
    @test_throws AssertionError getleaf(BitVector([0,1]),:binary)
    @test_throws AssertionError getleaf(maketree(4,2,:dwt),:quad)
    @test_throws AssertionError getleaf(BitVector([0,1,1,1,1]),:quad)
    @test_throws AssertionError getleaf(BitVector([0,1]),:quad)

    @test maketree(zeros(4,4)) == trues(5)
    @test maketree(zeros(4,4),:dwt) == BitVector([1,1,0,0,0])
    @test maketree(4,4,2) == trues(5)
    @test maketree(4,4,2,:dwt) == BitVector([1,1,0,0,0])
    @test_throws AssertionError maketree(4,4,3,:dwt)
    @test_throws AssertionError maketree(4,4,2,:fail)

    @test getdepth(5,:binary) == 2
    @test getdepth(5,:quad) == 1
    @test_throws AssertionError getdepth(5,:fail)
    @test_throws AssertionError getdepth(0,:binary)

    @test gettreelength(8) == 7
    @test gettreelength(8,8) == 21
    @test gettreelength(8,16) == 21
end

@testset "Metrics" begin
    x₀ = ones(5)
    x = 2*ones(5)
    @test relativenorm(x, x₀) == 1
    @test relativenorm(x, x₀, 1.0) == 1
    @test psnr(x, x₀) == 0
    @test snr(x, x₀) == 0
    @test ssim(x, x₀) == assess_ssim(x, x₀)
end

@testset "Dataset" begin
    x = [1, 0, 0, 0]
    @test duplicatesignals(x, 2, 1) == [1 0; 0 1; 0 0; 0 0]
    @test duplicatesignals(x, 2, 1, true) != duplicatesignals(x, 2, 1)
    @test duplicatesignals(x, 2, 1, true, 0.5) != duplicatesignals(x, 2, 1)

    @test length(generatesignals(:blocks, 5)) == 32
    @test length(generatesignals(:bumps, 5)) == 32
    @test length(generatesignals(:doppler, 5)) == 32
    @test length(generatesignals(:heavisine, 5)) == 32
    @test length(generatesignals(:quadchirp, 5)) == 32
    @test length(generatesignals(:mishmash, 5)) == 32
    @test_throws ArgumentError generatesignals(:fail, 5)

    @test typeof(ClassData(:tri, 5, 5, 5)) == ClassData
    @test typeof(ClassData(:cbf, 5, 5, 5)) == ClassData
    @test_throws ArgumentError ClassData(:fail, 5, 5, 5)

    @test size(generateclassdata(ClassData(:tri, 5, 5, 6))[1]) == (32,16)
    @test size(generateclassdata(ClassData(:cbf, 5, 5, 5))[1]) == (128,15)
    @test_nowarn generateclassdata(ClassData(:tri, 5, 5, 5), true)
end


text:
using AlgamesDriving
using Algames
using StaticArrays
using LinearAlgebra
using MeshCat

# Generate Autonomous Driving Scenario

# Create Roadway
roadway_opts = FourIntersectionRoadwayOptions()
roadway = build_roadway(roadway_opts)

# Create Dynamics Model
p = 3
model = UnicycleGame(p=p)

players = Vector{Player}(undef, p)
# For each player specify the start and goal state
# west_straight
players[1] = Player(model, roadway.lane[1],
    Q=Diagonal(SVector{model.ni[1]}([1,1,0.1,1])),
    x0=VehicleState(-2.2, -0.1, 0.0, 0.9),
    xf=VehicleState( 0.5, -0.1, 0.0, 0.9))

# west_right
players[2] = Player(model, roadway.lane[2],
    Q=Diagonal(SVector{model.ni[2]}([1,1,0.1,1])),
    x0=VehicleState(-1.8, -0.1, 0.0, 0.9),
    xf=VehicleState(-0.1, -1.0, -π/2, 0.9))

# west_left
players[3] = Player(model, roadway.lane[3],
    Q=Diagonal(SVector{model.ni[3]}([1,1,0.1,1])),
    x0=VehicleState(-1.4, -0.1, 0.0, 0.9),
    xf=VehicleState( 0.1,  1.0, π/2, 0.9))

# # east_straight
# players[1] = Player(model, roadway.lane[4],
#     Q=Diagonal(SVector{model.ni[1]}([1,1,0.1,1])),
#     x0=VehicleState(2.2, 0.1, 3.14, 0.9),
#     xf=VehicleState( -0.5, 0.1, 3.14, 0.9))
#
# # east_right
# players[2] = Player(model, roadway.lane[5],
#     Q=Diagonal(SVector{model.ni[2]}([1,1,0.1,1])),
#     x0=VehicleState(1.8, 0.1, 3.14, 0.9),
#     xf=VehicleState(0.1, 1.0, π/2, 0.9))
#
# # east_left
# players[3] = Player(model, roadway.lane[6],
#     Q=Diagonal(SVector{model.ni[3]}([1,1,0.1,1])),
#     x0=VehicleState(1.4, 0.1, 3.14, 0.9),
#     xf=VehicleState( -0.1, -1.0, 3*π/2, 0.9))

# # north_straight
# players[1] = Player(model, roadway.lane[7],
#     Q=Diagonal(SVector{model.ni[1]}([1,1,0.1,1])),
#     x0=VehicleState(-0.1, 2.2, -1.57, 0.9),
#     xf=VehicleState( -0.1, -0.5, -π/2, 0.9))

# # north_right
# players[2] = Player(model, roadway.lane[8],
#     Q=Diagonal(SVector{model.ni[2]}([1,1,0.1,1])),
#     x0=VehicleState(-0.1, 1.8, -1.57, 0.9),
#     xf=VehicleState(-1.0, 0.1, -π, 0.9))

# # north_left
# players[3] = Player(model, roadway.lane[9],
#     Q=Diagonal(SVector{model.ni[3]}([1,1,0.1,1])),
#     x0=VehicleState(-0.1, 1.4, -1.57, 0.9),
#     xf=VehicleState( 1.0, -0.1, 0., 0.9))

# # south_straight
# players[1] = Player(model, roadway.lane[10],
#     Q=Diagonal(SVector{model.ni[1]}([1,1,0.1,1])),
#     x0=VehicleState(0.1, -2.2, 1.57, 0.9),
#     xf=VehicleState( 0.1, 0.5, π/2, 0.9))

# # south_right
# players[2] = Player(model, roadway.lane[11],
#     Q=Diagonal(SVector{model.ni[2]}([1,1,0.1,1])),
#     x0=VehicleState(0.1, -1.8, 1.57, 0.9),
#     xf=VehicleState(1.0, -0.1, 0., 0.9))

# # south_left
# players[3] = Player(model, roadway.lane[12],
#     Q=Diagonal(SVector{model.ni[3]}([1,1,0.1,1])),
#     x0=VehicleState(0.1, -1.4, 1.57, 0.9),
#     xf=VehicleState( -1.0, 0.1, 3.14, 0.9))

# Create Scenario
sce = Scenario(model, roadway, players)

# # Create Dynamic Game Problem
N = 40 # number of time steps
dt = 0.1 # time step
solver_opts = Options() # solver options
prob = GameProblem(N, dt, sce, solver_opts) # autonomous driving problem
@time newton_solve!(prob) # solve the problem

# Plotting vehicles' trajectories and the solver's progress
# plot!(model, prob.pdtraj.pr)
# plot!(prob.stats)

# Visualize Scenario
vis = Visualizer()
open(vis)
# Build environment
set_scenario!(vis, sce)
# Set camera
# set_camera_birdseye!(vis, height=7)
set_env!(vis, VehicleState(0.,0.,0.,0.))
# Animate the vehicles
set_traj!(vis, model, sce, prob.pdtraj.pr)

ext: jl


text:
using CircGeometry

function create_circle()
    vf = 0.4
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    radius = 1.5
    center = Point(-0.5,1.0)
    outline = OutlineCircle(radius,center)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(
        outline,material,between_buffer;log=true)

    write_circ("circle.circ",ps)
    save_image("circle.svg",ps,outline)
end

function create_rectangle()
    vf = 0.4
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    p1 = Point(-1.0,0.0)
    p2 = Point(1.0,2.0)
    outline = OutlineRectangle(p1,p2)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(
        outline,material,between_buffer;log=true)

    write_circ("rectangle.circ",ps)
    save_image("rectangle.svg",ps,outline)
end

function create_car()
    vf = 0.4
    n_bodies = 400
    material = MaterialParameters(vf,n_bodies)

    polygon = csv_to_polygon("car.txt")
    outline = OutlinePolygon(polygon)

    between_buffer = compute_between_buffer(outline,material)
    ps = generate_porous_structure(
        outline,material,between_buffer;log=true)

    write_circ("car.circ",ps)
    save_image("car.svg",ps,outline)
end

function run_all_examples()
    println("running circle example...")
    create_circle()
    println("running rectangle example...")
    create_rectangle()
    println("running car (polygon) example...")
    create_car()
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
"""
computeElemsByMEBI
This function computes the vector of the element indexes for a given vector of MEBI parameter
"""
function computeElemsByMEBI( MEBIVals, MEBIVec )

    # get the relevant BCs definitions
    MEBIUnique = sort( unique( MEBIVals ) )     # get and sort them
    MEBIUnique[1]==0 && popfirst!( MEBIUnique ) # and remove zeros

    elements = []
    # loop
    print(MEBIUnique)
    for MEBIval in MEBIUnique
        print("\nMEBI:", MEBIval, "\n" )
        # find the elements with the current MEBIval
        aux = findall( x->x==MEBIval, MEBIVec )

        # if there are any, add them to the list
        length(aux) > 0 && push!( elements, aux )
    end
    return elements
end

function computeNeumDofs( nNodes )

    neumDofs = nodes2dofs( (1:nNodes), 6 )

    neumDofs[2:2:end] .= 0

    return neumDofs

end

"""
nodes2dofs computes the vector of dofs for an input vector of nodes
"""
function nodes2dofs( nodes , degreespernode )
    n    = length(nodes);
    dofs = zeros( Int, n*degreespernode ) ;
    for i in (1:n);
        print("\ni:",i,"\n")
        print( (i-1)*degreespernode .+ Vector(1:degreespernode) )
        print((degreespernode*(nodes[i]-1)).+Vector(1:degreespernode))
        dofs[ (i-1)*degreespernode .+ Vector(1:degreespernode) ] = (degreespernode*(nodes[i]-1)).+Vector(1:degreespernode)  ;
    end
    return dofs;
end


text:
module AdapterBase # Octo

using ..Octo
using .Octo.Queryable: Structured, FromItem, SubQuery
using .Octo: SQLElement, SQLAlias, SQLExtract, SQLFunctionName, SQLFunction, Field, Predicate, Raw, Enclosed, PlaceHolder, SQLKeyword, KeywordAllKeyword, VectorOfTuples
using .Octo.Schema
using .Octo.DBMS
using .Octo.DBMS: AbstractDatabase
using .Octo: Year, Month, Day, Hour, Minute, Second, CompoundPeriod, DatePeriod, TimePeriod, DateTime, format
using .Octo: @sql_keywords, @sql_functions, db_keywords, db_functionnames

const current = Dict{Symbol,AbstractDatabase}(
    :database => DBMS.SQL()
)

@sql_keywords  ADD ALL ALTER AND AS ASC BEGIN BETWEEN BIGINT BY COMMIT COLUMN CONSTRAINT CREATE DATABASE DEFAULT DELETE DESC DISTINCT DROP EXCEPT EXECUTE EXISTS FOREIGN FROM FULL GROUP
@sql_keywords  HAVING IF IN INDEX INNER INSERT INTERSECT INTO IS JOIN KEY LEFT LIKE LIMIT NULL OFF OFFSET ON OR ORDER OUTER OVER
@sql_keywords  PARTITION PREPARE PRIMARY RECURSIVE REFERENCES RELEASE RIGHT ROLLBACK SAVEPOINT SELECT SET TABLE TO TRANSACTION TRIGGER UNION UPDATE USING VALUES WHERE WITH

# dates
@sql_keywords EXTRACT YEAR MONTH DAY HOUR MINUTE SECOND TIMESTAMP INTERVAL

# @sql_keywords ANY (Julia TypeVar)

# aggregates
@sql_functions AVG COUNT EVERY MAX MIN NOT SOME SUM

# rankings
@sql_functions DENSE_RANK RANK ROW_NUMBER

# Colored SQL statements
module Beautiful # Octo.AdapterBase

struct ElementStyle
    color::Symbol
    bold::Bool
    ElementStyle(color::Symbol, bold::Bool=false) = new(color, bold)
end

struct Element
    style::ElementStyle
    body
end

struct Container
    elements::Vector{Union{Element,Container}}
    sep::String
end

end # Octo.AdapterBase.Beautiful

const SqlPartElement = Beautiful.Element
const SqlPart        = Beautiful.Container
using .Beautiful: ElementStyle

const style_normal                  = ElementStyle(:normal)
const style_subquery_select         = ElementStyle(:light_green, true)
const style_subquery_non_select     = ElementStyle(:light_blue, true)
const style_placeholder             = ElementStyle(:green, true)
const style_keyword                 = ElementStyle(:cyan)
const style_functionname            = ElementStyle(:cyan)
const style_function                = ElementStyle(:yellow)
const style_string                  = ElementStyle(:light_magenta)
const style_dates                   = ElementStyle(:light_green)
const style_predicate_enclosed      = ElementStyle(:normal)
const style_table_name              = ElementStyle(:normal, false)
const style_table_alias             = ElementStyle(:normal, true)
const style_field_fromclause_alias  = ElementStyle(:normal)
const style_field_fromclause_dot    = ElementStyle(:normal)
const style_field_fromclause_name   = ElementStyle(:normal)
const style_field_subquery_alias    = ElementStyle(:light_green)
const style_field_subquery_dot      = ElementStyle(:normal)
const style_field_subquery_name     = ElementStyle(:normal)
const style_vector_of_tuples        = ElementStyle(:yellow)

# sqlrepr -> SqlPartElement

sqlrepr(::DB where DB<:AbstractDatabase, sym::Symbol)::SqlPartElement        = SqlPartElement(style_normal, sym)
sqlrepr(::DB where DB<:AbstractDatabase, num::Number)::SqlPartElement        = SqlPartElement(style_normal, num)
sqlrepr(::DB where DB<:AbstractDatabase, f::Function)::SqlPartElement        = SqlPartElement(style_normal, f)
sqlrepr(::DB where DB<:AbstractDatabase, h::PlaceHolder)::SqlPartElement     = SqlPartElement(style_placeholder, h.body)
sqlrepr(::DB where DB<:AbstractDatabase, el::SQLKeyword)::SqlPartElement     = SqlPartElement(style_keyword, el.name)
sqlrepr(::DB where DB<:AbstractDatabase, var::TypeVar)::SqlPartElement       = SqlPartElement(style_keyword, string(var)) # ANY
sqlrepr(::DB where DB<:AbstractDatabase, f::SQLFunctionName)::SqlPartElement = SqlPartElement(style_functionname, f.name)
sqlrepr(::DB where DB<:AbstractDatabase, T::Type)::SqlPartElement            = SqlPartElement(style_table_name, _table_name_of(T))

function sqlrepr(::DB where DB<:AbstractDatabase, M::Type{PlaceHolder})::SqlPartElement
    SqlPartElement(style_placeholder, '?')
end

# sqlrepr -> SqlPart

function enclosed_part(style, body::SqlPart)::SqlPart
    SqlPart([
        SqlPartElement(style, '('),
        body,
        SqlPartElement(style, ')')], "")
end

function sqlrepr(::DB where DB<:AbstractDatabase, str::String)::SqlPart
    quot = SqlPartElement(style_string, "'")
    SqlPart([quot, SqlPartElement(style_string, str), quot], "")
end

function sqlrepr(::DB where DB<:AbstractDatabase, field::Field)::SqlPart
    if field.clause isa Nothing ||
       field.clause.__octo_alias isa Nothing
        SqlPart([SqlPartElement(style_normal, field.name)], "")
    else
        if field.clause isa SubQuery
            (style_field_alias, style_field_dot, style_field_name) = (style_field_subquery_alias, style_field_subquery_dot, style_field_subquery_name)
        else # FromItem
            (style_field_alias, style_field_dot, style_field_name) = (style_field_fromclause_alias, style_field_fromclause_dot, style_field_fromclause_name)
        end
        SqlPart([
            SqlPartElement(style_field_alias, field.clause.__octo_alias),
            SqlPartElement(style_field_dot, '.'),
            SqlPartElement(style_field_name, field.name)], "")
    end
end

function sqlrepr(::DB where DB<:AbstractDatabase, raw::Raw)::SqlPart
    lines = split(raw.string, '\n'; keepempty=true)
    parts = []
    for (line_idx, line) in enumerate(lines)
        words = split(line, ' '; keepempty=true)
        for (word_idx, word) in enumerate(words)
            chars = collect(word)
            if all(isuppercase, chars)
                if word in db_keywords
                    push!(parts, SqlPartElement(style_keyword, word))
                elseif word in db_functionnames
                    push!(parts, SqlPartElement(style_functionname, word))
                else
                    push!(parts, SqlPartElement(style_normal, word))
                end
            elseif all(x -> isuppercase(x) || x in ('(', ','), chars)
                if word[1:end-1] in db_keywords
                    push!(parts, SqlPartElement(style_keyword, word[1:end-1]))
                    push!(parts, SqlPartElement(style_normal, ","))
                elseif word[1:end-1] in db_functionnames
                    push!(parts, SqlPartElement(style_functionname, word[1:end-1]))
                    push!(parts, SqlPartElement(style_normal, ","))
                else
                    push!(parts, SqlPartElement(style_normal, word))
                end
            else
                push!(parts, SqlPartElement(style_normal, word))
            end
            word_idx != length(words) && push!(parts, SqlPartElement(style_normal, ' '))
        end
        line_idx != length(lines) && push!(parts, SqlPartElement(style_normal, '\n'))
    end
    SqlPart(parts, "")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, el::KeywordAllKeyword)::SqlPart
    els = [el.left, *, el.right]
    SqlPart(sqlrepr.(Ref(db), els), " ")
end

function both_isa((a,b), T::Type)::Bool
    a isa T && b isa T
end

function subquery_startswith_select(subquery::SubQuery)::Bool
    first(subquery.__octo_query) === SELECT
end

function _subquery_predicate_side(db::DB where DB<:AbstractDatabase, side)
    if side isa SubQuery && side.__octo_alias isa Symbol
        if subquery_startswith_select(side)
            sqlrepr(db, side)
        else
            SqlPartElement(style_subquery_non_select, side.__octo_alias)
        end
    else
        sqlrepr(db, side)
    end
end

function sqlrepr(db::DB where DB<:AbstractDatabase, pred::Predicate)::SqlPart
    if ==(pred.func, ==)
        op = :(=)
    else
        op = pred.func
    end
    (left, right) = _subquery_predicate_side.(Ref(db), (pred.left, pred.right))
    parts = [left, sqlrepr(db, op), right]
    body = SqlPart(parts, " ")
    if op in (+, -) && both_isa((pred.left, pred.right), Union{Field, SQLFunction})
        enclosed_part(style_predicate_enclosed, body)
    else
        body
    end
end

function _table_name_of(M::Type)::String # throw Schema.TableNameError
    Tname = Base.typename(M)
    if haskey(Schema.tables, Tname)
        table = Schema.tables[Tname]
        table[:table_name]
    else
        name = nameof(M)
        throw(Schema.TableNameError("""Provide schema table_name by `Schema.model($name, table_name="tbl")`"""))
    end
end

function _sqlrepr(db::DB where DB<:AbstractDatabase, clause::FromItem; with_as::Bool)::SqlPart
    table_name = _table_name_of(clause.__octo_model)
    part = SqlPartElement(style_table_name, table_name)
    if clause.__octo_alias isa Nothing
        parts = [part]
    else
        if table_name == String(clause.__octo_alias)
            parts = [part]
        else
            parts = [part, (with_as ? [sqlrepr(db, AS)] : [])..., SqlPartElement(style_table_alias, clause.__octo_alias)]
        end
    end
    SqlPart(parts, " ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, clause::FromItem)::SqlPart
    _sqlrepr(db, clause; with_as=true)
end

function sqlrepr(db::DB where DB<:AbstractDatabase, subquery::SubQuery)::SqlPart
    query = subquery.__octo_query
    body = SqlPart(vcat(sqlrepr.(Ref(db), query)...), " ")
    style = subquery_startswith_select(subquery) ? style_subquery_select : style_subquery_non_select
    part = enclosed_part(style, body)
    if subquery.__octo_alias isa Nothing
        part
    else
        SqlPart([
            part,
            sqlrepr(db, AS),
            SqlPartElement(style, subquery.__octo_alias),
        ], " ")
    end
end

function sqlrepr(db::DB where DB<:AbstractDatabase, tup::NamedTuple)::SqlPart
    els = map(kv -> Predicate(==, kv.first, kv.second), collect(pairs(tup)))
    SqlPart(sqlrepr.(Ref(db), els), ", ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, enclosed::Enclosed)::SqlPart
    body = SqlPart(sqlrepr.(Ref(db), enclosed.values), ", ")
    enclosed_part(style_normal, body)
end

function sqlrepr(db::DB where DB<:AbstractDatabase, vec_of_tuples::VectorOfTuples)::SqlPart
    parts = map(x -> enclosed_part(style_vector_of_tuples, SqlPart(sqlrepr.(Ref(db), collect(x)), ", ")), vec_of_tuples.values)
    SqlPart(parts, ", ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, a::SQLAlias)::SqlPart
    els = [a.field, AS, a.alias]
    SqlPart(sqlrepr.(Ref(db), els), " ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, extract::SQLExtract)::SqlPart
    body = SqlPart(sqlrepr.(Ref(db), [extract.field, FROM, extract.from]), " ")
    SqlPart([
        sqlrepr(db, EXTRACT),
        enclosed_part(style_normal, body)
    ], "")
end

sqlrepr(db::DB where DB<:AbstractDatabase, ::Type{Year})::SqlPartElement   = sqlrepr(db, YEAR)
sqlrepr(db::DB where DB<:AbstractDatabase, ::Type{Month})::SqlPartElement  = sqlrepr(db, MONTH)
sqlrepr(db::DB where DB<:AbstractDatabase, ::Type{Day})::SqlPartElement    = sqlrepr(db, DAY)
sqlrepr(db::DB where DB<:AbstractDatabase, ::Type{Hour})::SqlPartElement   = sqlrepr(db, HOUR)
sqlrepr(db::DB where DB<:AbstractDatabase, ::Type{Minute})::SqlPartElement = sqlrepr(db, MINUTE)
sqlrepr(db::DB where DB<:AbstractDatabase, ::Type{Second})::SqlPartElement = sqlrepr(db, SECOND)

function sqlrepr(db::DB where DB<:AbstractDatabase, dt::DateTime)::SqlPart
    str = format(dt, "yyyy-mm-dd HH:MM:SS")
    quot = SqlPartElement(style_dates, "'")
    SqlPart([
        sqlrepr(db, TIMESTAMP),
        SqlPart([quot, SqlPartElement(style_dates, str), quot], "")
    ], " ")
end

function compound_period_string(x::CompoundPeriod)
    if isempty(x.periods)
        return "empty period"
    else
        s = ""
        for p in x.periods
            s *= string(' ', p)
        end
        return s[2:end]
    end
end

function sqlrepr(db::DB where DB<:AbstractDatabase, period::Union{DatePeriod, TimePeriod})::SqlPart
    str = string(period)
    quot = SqlPartElement(style_dates, "'")
    SqlPart([
        sqlrepr(db, INTERVAL),
        SqlPart([quot, SqlPartElement(style_dates, str), quot], "")
    ], " ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, period::CompoundPeriod)::SqlPart
    str = compound_period_string(period)
    quot = SqlPartElement(style_dates, "'")
    SqlPart([
        sqlrepr(db, INTERVAL),
        SqlPart([quot, SqlPartElement(style_dates, str), quot], "")
    ], " ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, f::SQLFunction)::SqlPart
    SqlPart([
        SqlPartElement(style_function, f.name),
        sqlrepr(db, Enclosed(collect(f.fields))),
        ], "")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, tup::Tuple)::SqlPart
    parts = []
    for el in tup
        push!(parts, sqlrepr(db, el))
    end
    SqlPart(parts, ", ")
end

function sqlrepr(db::DB where DB<:AbstractDatabase, query::Structured)::SqlPart
    nth = 1
    els = []
    prev = nothing
    for el in query
        if el isa Type{PlaceHolder}
            push!(els, _placeholder(db, nth))
            nth += 1
        elseif el isa Predicate && el.right isa Type{PlaceHolder}
            push!(els, Predicate(el.func, el.left, _placeholder(db, nth)))
            nth += 1
        elseif el isa SubQuery
            if prev === AS
                push!(els, typeof(el)(el.__octo_query, nothing))
            elseif prev === JOIN
                push!(els, el.__octo_alias)
            else
                push!(els, el)
            end
        elseif el isa Tuple
            if isa(prev, FromItem)
                if el isa NTuple{N, Field} where N
                    push!(els, Enclosed(map(field -> Field(nothing, field.name), collect(el))))
                else
                    push!(els, Enclosed(collect(el)))
                end
            elseif prev === IN ||
                   prev === OVER ||
                   prev === USING ||
                   prev === VALUES
                push!(els, Enclosed(collect(el)))
            else
                push!(els, el)
            end
        elseif prev === USING
            push!(els, Enclosed([el]))
        else
            push!(els, el)
        end
        prev = el
    end
    SqlPart(sqlrepr.(Ref(db), els), " ")
end

function joinpart(part::SqlPart)::String
    join(map(part.elements) do el
        if el isa SqlPart
            joinpart(el)
        elseif el isa SqlPartElement
            el.body
        end
    end, part.sep)
end

function printpart(io::IO, el::SqlPartElement)
    printstyled(io, el.body; color=el.style.color, bold=el.style.bold)
end

function printpart(io::IO, part::SqlPart)
    for (idx, el) in enumerate(part.elements)
        printpart(io, el)
        length(part.elements) == idx || print(io, part.sep)
    end
end

# _to_sql

function _to_sql(db::DB where DB<:AbstractDatabase, element::Union{E,Structured} where E<:SQLElement)::String
    joinpart(sqlrepr(db, element))
end

# _placeholder, _placeholders

function _placeholder(db::DB where DB<:AbstractDatabase, nth::Int)
    PlaceHolder("?")
end

function _placeholders(db::DB where DB<:AbstractDatabase, dims::Int)
    Enclosed(fill(_placeholder(db, 1), dims))
end

# _show

function Base.show(io::IO, mime::MIME"text/plain", element::Union{E,Structured} where E<:SQLElement)
    db = current[:database]    # to be changed by Repo.connect
    _show(io, mime, db, element)
end

function _show(io::IO, ::MIME"text/plain", db::DB where DB<:AbstractDatabase, element::E where E<:SQLElement)
    printstyled(io, nameof(typeof(element)), color=:underline)
    print(io, ' ')
    printpart(io, sqlrepr(db, element))
end

function _show(io::IO, ::MIME"text/plain", db::DB where DB<:AbstractDatabase, query::Structured)
    if any(x -> x isa SQLElement, query)
        printpart(io, sqlrepr(db, query))
    else
        Base.show(io, query)
    end
end

end # module Octo.AdapterBase


text:
using PtFEM
using Test

# write your own tests here

code_tests = [
  "p41.1.jl",
  "p41.2.jl",
  "p42.2.jl",
  "p43.1b.jl",
  "p44.2.jl",
  "p45.2.jl",
  "p46.2.jl",
  "p47.1.jl",

  "p51.1.jl",
  "p51.2.jl",
  "p51.3.jl",
  "p51.4.jl",
  "p51.5.jl",
  "p52.1.jl",
  "p53.1.jl",
  "p54.1.jl",
  "p54.2.jl",
  "p55.1.jl",
  "p56.1.jl",

  "p61.1.jl",
  "p62.1.jl",
  "p63.1.jl"
]

println("\n\nRunning PtFEM/PtFEM.jl tests:\n\n")

@testset "PtFEM.jl" begin
  for test in code_tests
      println("\n\n  * $(test) *\n")
      include("test_"*test)
  end
end


text:
module Intersection

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

export intersection

Genie.Renderer.Html.register_normal_element("q__intersection", context = @__MODULE__)

function intersection(args...;
                      wrap::Function = StippleUI.DEFAULT_WRAPPER,
                      kwargs...)

  wrap() do
    q__intersection(args...; attributes([kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
export StepperBS

# @code_warntype ✓
function stepsize_tup(H::Real, ns::NTuple{KMAX,Int}) where KMAX
    if @generated
        a = Vector{Expr}(undef, KMAX)
        @inbounds for k in 1:KMAX
            a[k] = :(H / ns[$k])
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:tuple, a...))
        end
    else
        return ntuple(k -> H / ns[k], KMAX)
    end
end

# @code_warntype ✓
function stepinv2_tup(ns::NTuple{KMAX,Int}) where KMAX
    if @generated
        a = Vector{Expr}(undef, KMAX)
        @inbounds for k in 1:KMAX
            a[k] = :(4 / abs2(ns[$k]))
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:tuple, a...))
        end
    else
        return ntuple(k -> 4 / abs2(ns[k]), KMAX)
    end
end

# @code_warntype ✓
function multistepBS!(ηt::VecI, t0::Real, hk::Real, y0::VecI, diff!::Function, nk::Int, ηm::VecB, ηn::VecB; fargs...)
    h2 = 2.0 * hk
    copy!(ηm, y0)
    diff!(ηn, t0, ηm; fargs...)
    xpby!(ηm, hk, ηn)
    one2n = eachindex(ηt)

    tj = t0 + hk; j = 1
    while j < nk
        diff!(ηt, tj, ηn; fargs...)
        xpby!(ηm, h2, ηt)
        @simd for i in one2n
            @inbounds ηm[i], ηn[i] = ηn[i], ηt[i]
        end
        tj += hk; j += 1
    end

    diff!(ηt, tj, ηn; fargs...)
    @simd for i in one2n
        @inbounds ηt[i] = 0.5 * (ηn[i] + ηm[i] + hk * ηt[i])
    end
end

# @code_warntype ✓
function extrap2zero!(Ts::MatO{T}, xs::NTuple, ηt::VecI) where T
    zeroT = zero(T)
    @inbounds for i in eachindex(ηt)
        Ts[i,1], Ts[i,2] = ηt[i], zeroT
    end
end
# @code_warntype ✓
function extrap2zero!(Ts::MatO{T}, xs::NTuple, ηt::VecI, k::Int) where T
    isone(k) && return extrap2zero!(Ts, xs, ηt)
    zeroT = zero(T)
    one2n = eachindex(ηt)
    for j in 1:k-1
        @simd for i in one2n
            @inbounds Ts[i,j] -= ηt[i]
        end
    end
    @simd for i in one2n
        @inbounds Ts[i,k] = zeroT
    end
    for j in reverse(1:k-1)
        @inbounds xratio = xs[k] / xs[j] - 1.0
        @inbounds for i in one2n
            Ts[i,j] += (Ts[i,j] - Ts[i,j+1]) / xratio
        end
    end
    for j in eachindex(1:k)
        @simd for i in one2n
            @inbounds Ts[i,j] += ηt[i]
        end
    end
end

# @code_warntype ✓
function estimate_err(Tsol::VecI, Tref::VecI, y0::VecI, atol::Real, rtol::Real)
    e = 0.0 # error
    @inbounds for i in eachindex(y0)
        s  = atol + rtol * max(abs(y0[i]), abs(Tref[i])) # scale
        e += abs2((Tsol[i] - Tref[i]) / s)
    end
    return sqrt(e / length(y0))
end

struct StepperBS{KMAX}
    ns::NTuple{KMAX,Int}
    hs::NTuple{KMAX,Float64}
    xs::NTuple{KMAX,Float64}
    ηt::Vector{Float64}
    ηm::Vector{Float64}
    ηn::Vector{Float64}
    Ts::Matrix{Float64}

    StepperBS{13}(H::Real, ydim::Int) = StepperBS(H, ydim, (2,  6, 10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50))
    StepperBS{ 9}(H::Real, ydim::Int) = StepperBS(H, ydim, (2,  8, 14, 20, 26, 32, 38, 44, 50))
    StepperBS{ 7}(H::Real, ydim::Int) = StepperBS(H, ydim, (2, 10, 18, 26, 34, 42, 50))

    function StepperBS(H::Real, ydim::Int, ns::NTuple{KMAX,Int}) where KMAX
        hs = stepsize_tup(H, ns)
        xs = stepinv2_tup(ns)
        ηt = Vector{Float64}(undef, ydim)
        ηm = Vector{Float64}(undef, ydim)
        ηn = Vector{Float64}(undef, ydim)
        Ts = Matrix{Float64}(undef, ydim, KMAX)
        return new{KMAX}(ns, hs, xs, ηt, ηm, ηn, Ts)
    end
end

# @code_warntype ✓
function (stbs::StepperBS{KMAX})(ys::VecO{T}, t0::Real, y0::VecI, f!::Function; disp::Bool=false, fargs...) where {KMAX,T<:Real}
    ns = stbs.ns
    hs = stbs.hs
    xs = stbs.xs
    ηt = stbs.ηt
    ηm = stbs.ηm
    ηn = stbs.ηn
    Ts = stbs.Ts

    atol = 1.0e-12
    rtol = 1.0e-12
    preE = Inf # previous error

    @nview Ts Tsol Tref # Ts[solution], Ts[reference]
    @inbounds for k in 1:KMAX
        multistepBS!(ηt, t0, hs[k], y0, f!, ns[k], ηm, ηn; fargs...)
        extrap2zero!(Ts, xs, ηt, k)

        nowE = estimate_err(Tsol, Tref, y0, atol, rtol)
        if nowE > preE || iszero(nowE)
            disp && println("converged error @ k = $k: $nowE")
            break
        end
        preE = nowE
    end
    @simd for i in eachindex(ys)
        @inbounds ys[i] = Tsol[i]
    end
end


text:
"""
A structure that stores information created by make_cue_matrix:
C is the cue matrix;
f2i is a dictionary returning the indices for features;
i2f is a dictionary returning the features for indices;
gold_ind is a list of indices of gold paths;
A is the adjacency matrix;
grams is the number of grams for cues;
target_col is the column name for target strings;
tokenized is whether the dataset target is tokenized;
sep_token is the separator;
keep_sep is whether to keep separators in cues;
start_end_token is the start and end token in boundary cues.
"""
struct Cue_Matrix_Struct
    C::Union{Matrix,SparseMatrixCSC}
    f2i::Dict
    i2f::Dict
    gold_ind::Vector{Vector{Int64}}
    A::SparseMatrixCSC
    grams::Int64
    target_col::Union{Symbol, String}
    tokenized::Bool
    sep_token::Union{String, Nothing}
    keep_sep::Bool
    start_end_token::String
end

"""
Construct cue matrix.
"""
function make_cue_matrix end

"""
Construct cue matrix where combined features and adjacencies for both training datasets and validation datasets.
"""
function make_combined_cue_matrix end

"""
Given a list of string tokens, extract their n-grams.
"""
function make_ngrams end

"""
    make_cue_matrix(data::DataFrame)

Make the cue matrix for training datasets and corresponding indices as well as the adjacency matrix
and gold paths given a dataset in a form of dataframe.

# Obligatory Arguments
- `data::DataFrame`: the dataset

# Optional Arguments
- `grams::Int64=3`: the number of grams for cues
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings
- `tokenized::Bool=false`:if true, the dataset target is assumed to be tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `keep_sep::Bool=false`: if true, keep separators in cues
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# make cue matrix without tokenization
cue_obj_train = JudiLing.make_cue_matrix(
     latin_train,
    grams=3,
    target_col=:Word,
    tokenized=false,
    sep_token="-",
    start_end_token="#",
    keep_sep=false,
    verbose=false
    )

# make cue matrix with tokenization
cue_obj_train = JudiLing.make_cue_matrix(
    french_train,
    grams=3,
    target_col=:Syllables,
    tokenized=true,
    sep_token="-",
    start_end_token="#",
    keep_sep=true,
    verbose=false
    )
```
"""
function make_cue_matrix(
    data::DataFrame;
    grams = 3,
    target_col = :Words,
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    start_end_token = "#",
    verbose = false,
)

    # split tokens from words or other columns
    if tokenized && !isnothing(sep_token)
        tokens = split.(data[:, target_col], sep_token)
    else
        tokens = split.(data[:, target_col], "")
    end

    # making ngrams from tokens
    # make_ngrams function are below
    ngrams = make_ngrams.(tokens, grams, keep_sep, sep_token, start_end_token)

    # find all unique ngrams features
    ngrams_features = unique(vcat(ngrams...))

    f2i = Dict(v => i for (i, v) in enumerate(ngrams_features))
    i2f = Dict(i => v for (i, v) in enumerate(ngrams_features))

    n_f = sum([length(v) for v in ngrams])

    m = size(data, 1)
    n = length(ngrams_features)
    I = zeros(Int64, n_f)
    J = zeros(Int64, n_f)
    V = ones(Int64, n_f)

    A = [Int64[] for i = 1:length(ngrams_features)]

    cnt = 0
    for (i, v) in enumerate(ngrams)
        last = 0
        for (j, f) in enumerate(v)
            cnt += 1
            I[cnt] = i
            fi = f2i[f]
            J[cnt] = fi
            if j == 1
                last = fi
            else
                push!(A[last], fi)
                last = fi
            end
        end
    end

    cue = sparse(I, J, V, m, n, *)

    ngrams_ind = [[f2i[x] for x in y] for y in ngrams]

    verbose && println("making adjacency matrix...")
    A = [sort(unique(i)) for i in A]
    n_adj = sum(length.(A))
    I = zeros(Int64, n_adj)
    J = zeros(Int64, n_adj)
    V = ones(Int64, n_adj)

    cnt = 0
    iter = enumerate(A)
    if verbose
        pb = Progress(length(A))
    end
    for (i, v) in iter
        for j in v
            cnt += 1
            I[cnt] = i
            J[cnt] = j
        end
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    A = sparse(I, J, V, length(f2i), length(f2i))

    Cue_Matrix_Struct(cue, f2i, i2f, ngrams_ind, A, grams, target_col,
        tokenized, sep_token, keep_sep, start_end_token)
end

"""
    make_cue_matrix(data::DataFrame, cue_obj::Cue_Matrix_Struct)

Make the cue matrix for validation datasets and corresponding indices as well as the adjacency matrix
and gold paths given a dataset in a form of dataframe.

# Obligatory Arguments
- `data::DataFrame`: the dataset
- `cue_obj::Cue_Matrix_Struct`: training cue object

# Optional Arguments
- `grams::Int64=3`: the number of grams for cues
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings
- `tokenized::Bool=false`:if true, the dataset target is assumed to be tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `keep_sep::Bool=false`: if true, keep separators in cues
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# make cue matrix without tokenization
cue_obj_val = JudiLing.make_cue_matrix(
  latin_val,
  cue_obj_train,
  grams=3,
  target_col=:Word,
  tokenized=false,
  sep_token="-",
  keep_sep=false,
  start_end_token="#",
  verbose=false
  )

# make cue matrix with tokenization
cue_obj_val = JudiLing.make_cue_matrix(
    french_val,
    cue_obj_train,
    grams=3,
    target_col=:Syllables,
    tokenized=true,
    sep_token="-",
    keep_sep=true,
    start_end_token="#",
    verbose=false
    )
```
"""
function make_cue_matrix(
    data::DataFrame,
    cue_obj::Cue_Matrix_Struct;
    grams = 3,
    target_col = "Words",
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    start_end_token = "#",
    verbose = false,
)

    # split tokens from words or other columns
    if tokenized && !isnothing(sep_token)
        tokens = split.(data[:, target_col], sep_token)
    else
        tokens = split.(data[:, target_col], "")
    end

    # making ngrams from tokens
    # make_ngrams function are below
    ngrams = make_ngrams.(tokens, grams, keep_sep, sep_token, start_end_token)

    f2i = cue_obj.f2i
    i2f = cue_obj.i2f

    n_f = sum([length(v) for v in ngrams])

    m = size(data, 1)
    n = length(f2i)
    I = zeros(Int64, n_f)
    J = zeros(Int64, n_f)
    V = ones(Int64, n_f)

    cnt = 0
    for (i, v) in enumerate(ngrams)
        for (j, f) in enumerate(v)
            cnt += 1
            I[cnt] = i
            J[cnt] = f2i[f]
        end
    end

    cue = sparse(I, J, V, m, n, *)
    ngrams_ind = [[f2i[x] for x in y] for y in ngrams]

    Cue_Matrix_Struct(cue, f2i, i2f, ngrams_ind, cue_obj.A, grams, target_col,
        tokenized, sep_token, keep_sep, start_end_token)
end

"""
    make_cue_matrix(data_train::DataFrame, data_val::DataFrame)

Make the cue matrix for traiing and validation datasets at the same time.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset

# Optional Arguments
- `grams::Int64=3`: the number of grams for cues
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings
- `tokenized::Bool=false`:if true, the dataset target is assumed to be tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `keep_sep::Bool=false`: if true, keep separators in cues
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# make cue matrix without tokenization
cue_obj_train, cue_obj_val = JudiLing.make_cue_matrix(
    latin_train,
    latin_val,
    grams=3,
    target_col=:Word,
    tokenized=false,
    keep_sep=false
    )

# make cue matrix with tokenization
cue_obj_train, cue_obj_val = JudiLing.make_cue_matrix(
    french_train,
    french_val,
    grams=3,
    target_col=:Syllables,
    tokenized=true,
    sep_token="-",
    keep_sep=true,
    start_end_token="#",
    verbose=false
    )
```
"""
function make_cue_matrix(
    data_train::DataFrame,
    data_val::DataFrame;
    grams = 3,
    target_col = "Words",
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    start_end_token = "#",
    verbose = false,
)

    cue_obj_train = make_cue_matrix(
        data_train,
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose,
    )

    cue_obj_val = make_cue_matrix(
        data_val,
        cue_obj_train,
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose,
    )

    cue_obj_train, cue_obj_val
end

"""
  make_cue_matrix(data::DataFrame, pyndl_weights::Pyndl_Weight_Struct)

Make the cue matrix for pyndl mode.
"""
function make_cue_matrix(
    data::DataFrame,
    pyndl_weights::Pyndl_Weight_Struct;
    grams = 3,
    target_col = "Words",
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    start_end_token = "#",
    verbose = false,
)

    # split tokens from words or other columns
    if tokenized && !isnothing(sep_token)
        tokens = split.(data[:, target_col], sep_token)
    else
        tokens = split.(data[:, target_col], "")
    end

    # making ngrams from tokens
    # make_ngrams function are below
    ngrams = make_ngrams.(tokens, grams, keep_sep, sep_token, start_end_token)

    # find all unique ngrams features
    ngrams_features = unique(vcat(ngrams...))

    f2i = Dict(v => i for (i, v) in enumerate(pyndl_weights.cues))
    i2f = Dict(i => v for (i, v) in enumerate(pyndl_weights.cues))

    n_f = sum([length(v) for v in ngrams])

    m = size(data, 1)
    n = length(ngrams_features)
    I = zeros(Int64, n_f)
    J = zeros(Int64, n_f)
    V = ones(Int64, n_f)

    A = [Int64[] for i = 1:length(ngrams_features)]

    cnt = 0
    for (i, v) in enumerate(ngrams)
        last = 0
        for (j, f) in enumerate(v)
            cnt += 1
            I[cnt] = i
            fi = f2i[f]
            J[cnt] = fi
            if j == 1
                last = fi
            else
                push!(A[last], fi)
                last = fi
            end
        end
    end

    cue = sparse(I, J, V, m, n, *)

    ngrams_ind = [[f2i[x] for x in y] for y in ngrams]

    verbose && println("making adjacency matrix...")
    A = [sort(unique(i)) for i in A]
    n_adj = sum(length.(A))
    I = zeros(Int64, n_adj)
    J = zeros(Int64, n_adj)
    V = ones(Int64, n_adj)

    cnt = 0
    iter = enumerate(A)
    if verbose
        pb = Progress(length(A))
    end
    for (i, v) in iter
        for j in v
            cnt += 1
            I[cnt] = i
            J[cnt] = j
        end
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    A = sparse(I, J, V, length(f2i), length(f2i))

    Cue_Matrix_Struct(cue, f2i, i2f, ngrams_ind, A, grams, target_col,
        tokenized, sep_token, keep_sep, start_end_token)
end

"""
    make_combined_cue_matrix(data_train, data_val)

Make the cue matrix for training and validation datasets at the same time, where the features and adjacencies are combined.

# Obligatory Arguments
- `data_train::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset

# Optional Arguments
- `grams::Int64=3`: the number of grams for cues
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings
- `tokenized::Bool=false`:if true, the dataset target is assumed to be tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `keep_sep::Bool=false`: if true, keep separators in cues
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# make cue matrix without tokenization
cue_obj_train, cue_obj_val = JudiLing.make_combined_cue_matrix(
    latin_train,
    latin_val,
    grams=3,
    target_col=:Word,
    tokenized=false,
    keep_sep=false
    )

# make cue matrix with tokenization
cue_obj_train, cue_obj_val = JudiLing.make_combined_cue_matrix(
    french_train,
    french_val,
    grams=3,
    target_col=:Syllables,
    tokenized=true,
    sep_token="-",
    keep_sep=true,
    start_end_token="#",
    verbose=false
    )
```
"""
function make_combined_cue_matrix(
    data_train,
    data_val;
    grams = 3,
    target_col = "Words",
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    start_end_token = "#",
    verbose = false,
)

    data_combined = copy(data_train)
    append!(data_combined, data_val)

    cue_obj_combined = make_cue_matrix(
        data_combined,
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose,
    )

    cue_obj_train = make_cue_matrix(
        data_train,
        cue_obj_combined,
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose,
    )

    cue_obj_val = make_cue_matrix(
        data_val,
        cue_obj_combined,
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose,
    )

    cue_obj_train, cue_obj_val
end

"""
    make_ngrams(tokens, grams, keep_sep, sep_token, start_end_token)

Given a list of string tokens return a list of all n-grams for these tokens.
"""
function make_ngrams(
    tokens,
    grams,
    keep_sep,
    sep_token,
    start_end_token,
)

    push!(pushfirst!(tokens, start_end_token), start_end_token)
    if keep_sep
        # collect ngrams
        ngrams =
            join.(
                collect(zip((Iterators.drop(tokens, k) for k = 0:grams-1)...)),
                sep_token,
            )
    else
        ngrams =
            join.(
                collect(zip((Iterators.drop(tokens, k) for k = 0:grams-1)...)),
                "",
            )
    end

    ngrams
end


text:
# Flags
const RSB_FLAG_NOFLAGS                                 = 0x000000
const RSB_FLAG_IDENTICAL_FLAGS                         = RSB_FLAG_NOFLAGS
const RSB_FLAG_FORTRAN_INDICES_INTERFACE               = 0x000001
const RSB_FLAG_C_INDICES_INTERFACE                     = 0x000000
const RSB_FLAG_USE_HALFWORD_INDICES                    = 0x000002
const RSB_FLAG_WANT_ROW_MAJOR_ORDER                    = 0x000000
const RSB_FLAG_WANT_COLUMN_MAJOR_ORDER                 = 0x4000000
const RSB_FLAG_SORTED_INPUT                            = 0x000004
const RSB_FLAG_TRIANGULAR                              = 0x000008
const RSB_FLAG_LOWER                                   = 0x000010
const RSB_FLAG_UPPER                                   = 0x000020
const RSB_FLAG_UNIT_DIAG_IMPLICIT                      = 0x000040
const RSB_FLAG_WANT_COO_STORAGE                        = 0x000100
const RSB_FLAG_DUPLICATES_KEEP_LAST                    = 0x000000
const RSB_FLAG_DUPLICATES_DEFAULT_HANDLE               = 0x000000
const RSB_FLAG_DUPLICATES_SUM                          = 0x000200
const RSB_FLAG_DISCARD_ZEROS                           = 0x000400
const RSB_FLAG_QUAD_PARTITIONING                       = 0x002000
const RSB_FLAG_WANT_BCSS_STORAGE                       = 0x004000
const RSB_FLAG_ASSEMBLED_IN_COO_ARRAYS                 = 0x040000
const RSB_FLAG_EXPERIMENTAL_IN_PLACE_PERMUTATION_SORT  = 0x080000
const RSB_FLAG_SYMMETRIC                               = 0x400000
const RSB_FLAG_HERMITIAN                               = 0x800000
const RSB_FLAG_RECURSIVE_MORE_LEAVES_THAN_THREADS      = 0x1000000
const RSB_FLAG_LOWER_HERMITIAN          = (RSB_FLAG_HERMITIAN | RSB_FLAG_LOWER)
const RSB_FLAG_UPPER_HERMITIAN          = (RSB_FLAG_HERMITIAN | RSB_FLAG_UPPER)
const RSB_FLAG_LOWER_TRIANGULAR         = (RSB_FLAG_TRIANGULAR | RSB_FLAG_LOWER)
const RSB_FLAG_UPPER_TRIANGULAR         = (RSB_FLAG_TRIANGULAR | RSB_FLAG_UPPER)
const RSB_FLAG_LOWER_SYMMETRIC          = (RSB_FLAG_SYMMETRIC | RSB_FLAG_LOWER)
const RSB_FLAG_DIAGONAL                 = (RSB_FLAG_UPPER | RSB_FLAG_LOWER)
const RSB_FLAG_UPPER_SYMMETRIC          = (RSB_FLAG_SYMMETRIC | RSB_FLAG_UPPER)
const RSB_FLAG_RECURSIVE_SUBDIVIDE_MORE_ON_DIAG        = 0x8000000
const RSB_FLAG_EXTERNALLY_ALLOCATED_ARRAYS             = 0x40000000
const RSB_FLAG_USE_CSR_RESERVED                        = 0x200000
const RSB_FLAG_USE_HALFWORD_INDICES_CSR      = (RSB_FLAG_USE_HALFWORD_INDICES|RSB_FLAG_USE_CSR_RESERVED)
const RSB_FLAG_USE_HALFWORD_INDICES_COO      = (RSB_FLAG_USE_HALFWORD_INDICES|RSB_FLAG_WANT_COO_STORAGE)
const RSB_FLAG_MUTUALLY_EXCLUSIVE_SWITCHES   = (RSB_FLAG_USE_HALFWORD_INDICES_COO|RSB_FLAG_USE_HALFWORD_INDICES_CSR)
const RSB_FLAG_DEFAULT_STORAGE_FLAGS         = (RSB_FLAG_WANT_BCSS_STORAGE|RSB_FLAG_WANT_COO_STORAGE)
const RSB_FLAG_DEFAULT_COO_MATRIX_FLAGS      = RSB_FLAG_WANT_COO_STORAGE
const RSB_FLAG_DEFAULT_CSR_MATRIX_FLAGS      = RSB_FLAG_WANT_BCSS_STORAGE
const RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS      = (RSB_FLAG_QUAD_PARTITIONING|RSB_FLAG_USE_HALFWORD_INDICES|RSB_FLAG_WANT_COO_STORAGE|RSB_FLAG_WANT_BCSS_STORAGE)
const RSB_FLAG_DEFAULT_MATRIX_FLAGS          = RSB_FLAG_DEFAULT_RSB_MATRIX_FLAGS

# Transpose flags
const RSB_TRANSPOSITION_N = 0x4E
const RSB_TRANSPOSITION_T = 0x54
const RSB_TRANSPOSITION_C = 0x43

# Info flags
const RSB_MIF_INDEX_STORAGE_IN_BYTES__TO__SIZE_T             = 0x00000001
const RSB_MIF_INDEX_STORAGE_IN_BYTES_PER_NNZ__TO__RSB_REAL_T = 0x00000002
const RSB_MIF_MATRIX_ROWS__TO__RSB_COO_INDEX_T               = 0x00000004
const RSB_MIF_MATRIX_COLS__TO__RSB_COO_INDEX_T               = 0x00000008
const RSB_MIF_MATRIX_NNZ__TO__RSB_NNZ_INDEX_T                = 0x00000010
const RSB_MIF_TOTAL_SIZE__TO__SIZE_T                         = 0x00000020
const RSB_MIF_MATRIX_FLAGS__TO__RSB_FLAGS_T                  = 0x00000040
const RSB_MIF_MATRIX_TYPECODE__TO__RSB_TYPE_T                = 0x00000080
const RSB_MIF_MATRIX_INFO__TO__CHAR_P                        = 0x00000100
const RSB_MIF_LEAVES_COUNT__TO__RSB_BLK_INDEX_T              = 0x00000200

# Vector extraction flags
const RSB_EXTF_NORM_ONE  = 0x00001001
const RSB_EXTF_NORM_TWO  = 0x00001002
const RSB_EXTF_NORM_INF  = 0x00001003
const RSB_EXTF_SUMS_ROW  = 0x00001004
const RSB_EXTF_SUMS_COL  = 0x00001005
const RSB_EXTF_ASUMS_ROW = 0x00001006
const RSB_EXTF_ASUMS_COL = 0x00001007
const RSB_EXTF_DIAG      = 0x00000004

# Errors
const RSB_ERR_NO_ERROR = 0x000

const RSB_DEFAULT_ROW_BLOCKING = 1
const RSB_DEFAULT_COL_BLOCKING = 1
const RSB_DEFAULT_BLOCKING     = 1


text:
abstract type AbstractFont end

include("fonts/Terminus_32_16.jl")

const FONTS = [
               TERMINUS_32_16,
              ]

struct Character{I, C <: AbstractChar, F <: AbstractFont} <: AbstractShape
    position::Point{I}
    char::C
    font::F
end

has_char(font::Terminus_32_16, char) = isascii(char) && isprint(char)

function get_bitmap(font::Terminus_32_16, char::Char)
    bitmap = font.bitmap

    codepoint_begin = codepoint(' ')
    k = codepoint(char) - codepoint_begin + one(codepoint_begin)

    char_bitmap = @view bitmap[:, :, k]

    return char_bitmap
end

is_valid(shape::Character) = has_char(shape.font, shape.char)

get_i_min(shape::Character) = shape.position.i
get_i_max(shape::Character) = shape.position.i + size(shape.font.bitmap, 1) - one(shape.position.i)

get_j_min(shape::Character) = shape.position.j
get_j_max(shape::Character) = shape.position.j + size(shape.font.bitmap, 2) - one(shape.position.j)

function draw!(image::AbstractMatrix, shape::Character{I, C, Terminus_32_16} where {I, C}, color)
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    char = shape.char
    font = shape.font

    if is_outbounds(image, shape)
        return nothing
    end

    if char == ' '
        return nothing
    else
        bitmap_shape = Bitmap(position, get_bitmap(font, char))
        draw!(put_pixel_inbounds!, image, clip(image, bitmap_shape), color)
    end

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::Character{I, C, Terminus_32_16} where {I, C}, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    char = shape.char
    font = shape.font

    if char == ' '
        return nothing
    else
        bitmap_shape = Bitmap(position, get_bitmap(font, char))
        draw!(f, image, bitmap_shape, color)
    end

    return nothing
end


text:
#=
maximize expected returns subject to risk constraints

TODO
- add entropic ball constraint using entropy cone
- describe formulation and options
=#

using SparseArrays

struct PortfolioNative{T <: Real} <: ExampleInstanceNative{T}
    num_stocks::Int
    epinormeucl_constr::Bool # add L2 ball constraints
    epinorminf_constrs::Bool # add Linfty ball constraints
    use_epinorminf::Bool # use epinorminf cone, else nonnegative cones
end

function build(inst::PortfolioNative{T}) where {T <: Real}
    @assert xor(inst.epinormeucl_constr, inst.epinorminf_constrs)
    num_stocks = inst.num_stocks

    returns = rand(T, num_stocks)
    sigma_half = randn(T, num_stocks, num_stocks)
    x = randn(T, num_stocks)
    x ./= norm(x)
    gamma = sum(abs, sigma_half * x) / sqrt(T(num_stocks))

    c = -returns
    # investments add to one, nonnegativity
    A = ones(T, 1, num_stocks)
    G = sparse(-one(T) * I, num_stocks, num_stocks)
    b = T[1]
    h = zeros(T, num_stocks)
    cones = Cones.Cone{T}[Cones.Nonnegative{T}(num_stocks)]
    cone_offset = num_stocks

    function add_ball_constr(cone, gamma_new)
        G = vcat(G, spzeros(T, 1, num_stocks), -sigma_half)
        h_risk = vcat(gamma_new, zeros(T, num_stocks))
        h = vcat(h, h_risk)
        push!(cones, cone)
        cone_offset += num_stocks + 1
    end

    if inst.epinormeucl_constr
        add_ball_constr(Cones.EpiNormEucl{T}(num_stocks + 1), gamma)
    end

    if inst.epinorminf_constrs
        if inst.use_epinorminf
            add_ball_constr(Cones.EpiNormInf{T, T}(num_stocks + 1,
                use_dual = true), gamma * sqrt(T(num_stocks)))
            add_ball_constr(Cones.EpiNormInf{T, T}(num_stocks + 1), gamma)
        else
            c = vcat(c, zeros(T, 2 * num_stocks))
            A = [
                A    spzeros(T, 1, 2 * num_stocks);
                sigma_half    -I    I;
                ]
            padding = spzeros(T, num_stocks, 2 * num_stocks)
            G = [
                G    spzeros(T, size(G, 1), 2 * num_stocks);
                spzeros(T, 2 * num_stocks, num_stocks)    -I;
                spzeros(T, 1, num_stocks)    ones(T, 1, 2 * num_stocks);
                sigma_half    padding;
                -sigma_half    padding;
                ]
            b = vcat(b, zeros(T, num_stocks))
            h = vcat(h, zeros(T, 2 * num_stocks), gamma * sqrt(T(num_stocks)),
                gamma * ones(T, 2 * num_stocks))
            push!(cones, Cones.Nonnegative{T}(4 * num_stocks + 1))
            cone_offset += 4 * num_stocks + 1
        end
    end

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end


text:
# Autogenerated wrapper script for xdm_jll for i686-linux-musl
export SpiritCommon, SpiritExprCommon, XdmRapidXmlReader

using boostpython_jll
JLLWrappers.@generate_wrapper_header("xdm")
JLLWrappers.@declare_library_product(SpiritCommon, "SpiritCommon.so")
JLLWrappers.@declare_library_product(SpiritExprCommon, "SpiritExprCommon.so")
JLLWrappers.@declare_library_product(XdmRapidXmlReader, "XdmRapidXmlReader.so")
function __init__()
    JLLWrappers.@generate_init_header(boostpython_jll)
    JLLWrappers.@init_library_product(
        SpiritCommon,
        "xdm_bundle/SpiritCommon.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        SpiritExprCommon,
        "xdm_bundle/SpiritExprCommon.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        XdmRapidXmlReader,
        "xdm_bundle/XdmRapidXmlReader.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function solve()
    x = 0
    cs = []
    n = parse(Int, readline())
    for i=1:n
        a, b = [parse(Int, x) for x in split(readline())]
        x += a
        push!(cs, 2a+b)
    end
    sort!(cs)
    ret = 0
    while x ≥ 0
        ret += 1
        c = pop!(cs)
        x -= c
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
@testset "StudentTEffectSize" begin
    @testset "One group" begin
        # This example is taken from 5.2 in the Gitbook "Bayesian Statistics":
        # https://statswithr.github.io/book/hypothesis-testing-with-normal-populations.html
        # R code result:
        # ## Single numerical variable
        # n = 10, y-bar = 0.0804, s = 0.0523
        # (Using Zellner-Siow Cauchy prior:  mu ~ C(0, 1*sigma)
        # (Using Jeffreys prior: p(sigma^2) = 1/sigma^2
        #
        # Hypotheses:
        # H1: mu = 0 versus H2: mu != 0
        # Priors:
        # P(H1) = 0.5 , P(H2) = 0.5
        # Results:
        # BF[H2:H1] = 50.7757
        # P(H1|data) = 0.0193  P(H2|data) = 0.9807
        n = 10
        xbar = 0.0804
        s = 0.0523

        normalstat = NormalStatistics(n=n, meanx=xbar, sdx=0.0523)
        stats = StudentTStatistics(normalstat)
        model = StudentTEffectSize(r=1.0)

        bf = bayesfactor(model, stats)
        @test isapprox(bf, 50.6, rtol=0.01)
    end

    @testset "Two Samples" begin
        # Hypotheses:
        # H1: mu_mature mom  = mu_younger mom
        # H2: mu_mature mom != mu_younger mom
        #
        # Priors: P(H1) = 0.5  P(H2) = 0.5
        #
        # Results:
        # BF[H1:H2] = 5.7162
        # P(H1|data) = 0.8511
        # P(H2|data) = 0.1489
        #
        # Posterior summaries for under H2:
        # 95% Cred. Int.: (-4.3386 , 0.8633)

        # statistics from the nc dataset in R's `statsr` package
        # Group 1: mature mom
        # Group 2: young mom
        normalstats = TwoNormalStatistics(
            NormalStatistics(meanx=28.8, sdx=13.5, n=133),
            NormalStatistics(meanx=30.6, sdx=14.3, n=867)
        )
        model = StudentTEffectSize(r=1.0)

        # pooled case
        stats_pooled = StudentTStatistics(normalstats, pooled=true)
        bf12_pooled = 1/bayesfactor(model, stats_pooled)
        @test isapprox(bf12_pooled, 5.45, rtol=0.001)

        # welch case
        stats_welch = StudentTStatistics(normalstats, pooled=false)
        bf12_welch = 1/bayesfactor(model, stats_welch)
        @test isapprox(bf12_welch, 5.03,  rtol=0.001)
    end

    @testset "One sample: Sleep data" begin
        # This example uses the sleep data and "BayesFactor" package in R:
        # https://richarddmorey.github.io/BayesFactor/
        # Environment:
        # R: 4.0.4
        # BayesFactor:0.9.12-4.2
        # "sleep" data in R have records of 10 paired observations
        sleepdata = [-1.2, -2.4, -1.3, -1.3, 0.0, -1.0, -1.8, -0.8, -4.6, -1.4]

        normalstats = NormalStatistics(sleepdata)
        model = StudentTEffectSize()
        stats = StudentTStatistics(normalstats)

        # test the calculation for one-sample t-statistics
        @test isapprox(stats.t,  -4.0621, rtol=0.001)
        @test stats.n == length(sleepdata)
        @test stats.dof == length(sleepdata) - 1

        # test bayes factor calculation
        model = StudentTEffectSize(r=sqrt(2)/2)
        bf = bayesfactor(model, stats)

        @test isapprox(bf, 17.259, rtol=0.001)
    end

    @testset "Two samples: Chicken weights" begin
        # Chicken Weights by Feed Type in R:
        # Environment:
        # R: 4.0.4
        # BayesFactor:0.9.12-4.2
        #
        #> t.test(weight ~ feed, data=chickwts, var.eq=TRUE)
        # 	Two Sample t-test
        #
        # data:  weight by feed
        # t = -2.934, df = 20, p-value = 0.008205
        # alternative hypothesis: true difference in means is not equal to 0
        # 95 percent confidence interval:
        #  -100.17618  -16.92382
        # sample estimates:
        # mean in group horsebean   mean in group linseed
        #                  160.20                  218.75
        #
        # >bf = ttestBF(
        #     x=(chickwts%>%filter(feed=="horsebean"))$weight,
        #     y=(chickwts%>%filter(feed=="linseed"))$weight
        #     )
        #
        # Bayes factor analysis
        # --------------
        # [1] Alt., r=0.707 : 5.975741 ±0%
        #
        # Against denominator:
        # Null, mu1-mu2 = 0
        # ---
        # Bayes factor type: BFindepSample, JZS
        horsebean=[179, 160, 136, 227, 217, 168, 108, 124, 143, 140]
        linseed=[309, 229, 181, 141, 260, 203, 148, 169, 213, 257, 244, 271]

        stats1 = NormalStatistics(horsebean)
        stats2 = NormalStatistics(linseed)
        # Here we assume equal variance
        tstat = StudentTStatistics(TwoNormalStatistics(stats1, stats2))

        @test isapprox(tstat.t, -2.934, rtol=0.001)
        @test tstat.dof == 20

        model = StudentTEffectSize(r=sqrt(2)/2)
        bf = bayesfactor(model, tstat)
        @test isapprox(bf, 5.9757,  rtol=0.001)
    end

end


text:
##########################################################################
# Network definition
type InterferingBroadcastChannel{System_t <: System} <: CanonicalNetwork
    MSs::Vector{CanonicalMS}
    BSs::Vector{CanonicalBS}

    system::System_t
    num_MSs_per_cell::Int
    alpha::Float64
    aux_network_params::AuxNetworkParams

    assignment::Assignment
end

# Convenience constructor without network params and assignments
InterferingBroadcastChannel(MSs, BSs, system, num_MSs_per_cell, alpha) =
    InterferingBroadcastChannel(MSs, BSs, system, num_MSs_per_cell, alpha, AuxNetworkParams(), Assignment())

Base.show(io::IO, x::InterferingBroadcastChannel) =
    print(io, "IBC(I = $(length(x.BSs)), Kc = $(x.num_MSs_per_cell), α = $(x.alpha))")
Base.showcompact(io::IO, x::InterferingBroadcastChannel) =
    print(io, "IBC($(length(x.BSs)), $(x.num_MSs_per_cell), $(x.alpha))")

function setup_interfering_broadcast_channel(
    num_BSs, num_MSs_per_cell, num_MS_antennas, num_BS_antennas;
    system = SinglecarrierSystem(),
    alpha = 1.,
    transmit_power = 1., transmit_powers = transmit_power*ones(Float64, num_BSs),
    user_priority = 1., user_priorities = user_priority*ones(Float64, num_BSs*num_MSs_per_cell),
    num_streams = 1, num_streamss = num_streams*ones(Int, num_BSs*num_MSs_per_cell),
    receiver_noise_power = 1., receiver_noise_powers = receiver_noise_power*ones(Float64, num_BSs*num_MSs_per_cell))

    isa(num_MS_antennas, Vector) || (num_MS_antennas = num_MS_antennas*ones(Int, num_BSs*num_MSs_per_cell))
    isa(num_BS_antennas, Vector) || (num_BS_antennas = num_BS_antennas*ones(Int, num_BSs))

    BSs = [ CanonicalBS(num_BS_antennas[i], transmit_powers[i]) for i = 1:num_BSs ]
    MSs = [ CanonicalMS(num_MS_antennas[k], user_priorities[k], num_streamss[k], receiver_noise_powers[k]) for k = 1:num_BSs*num_MSs_per_cell ]

    InterferingBroadcastChannel(MSs, BSs, system, num_MSs_per_cell, alpha)
end

##########################################################################
# Standard cell assignment functions
function IDCellAssignment!(channel, network::InterferingBroadcastChannel)
    Kc = network.num_MSs_per_cell; I = get_num_BSs(network)
    cell_assignment = Array(Int, I*Kc)

    for i = 1:I
        cell_assignment[(i-1)*Kc+1:i*Kc] = i
    end

    network.assignment = Assignment(cell_assignment, I)

    return AssignmentResults()
end

# We don't have large scale fading in this network.
LargeScaleFadingCellAssignment!(channel, network::InterferingBroadcastChannel) =
    IDCellAssignment!(channel, network)

##########################################################################
# Simulation functions
draw_user_drop!(network::InterferingBroadcastChannel) = nothing

function draw_channel{System_t <: SinglecarrierSystem}(network::InterferingBroadcastChannel{System_t})
    I = get_num_BSs(network); Kc = network.num_MSs_per_cell
    Ns = get_num_MS_antennas(network); Ms = get_num_BS_antennas(network)

    coefs = Array(Matrix{Complex128}, I*Kc, I)
    large_scale_fading_factor = Array(Float64, I*Kc, I)

    for k = 1:I*Kc
        serving_id = div(k - 1, Kc) + 1

        for i = 1:I
            if i == serving_id
                large_scale_fading_factor[k,i] = 1
            else
                large_scale_fading_factor[k,i] = sqrt(network.alpha)
            end

            # Apply scale factors
            coefs[k,i] = large_scale_fading_factor[k,i]*(1/sqrt(2))*(randn(Ns[k], Ms[i]) + im*randn(Ns[k], Ms[i]))
        end
    end

    return SinglecarrierChannel(coefs, Ns, Ms, I*Kc, I, large_scale_fading_factor)
end

function draw_channel{System_t <: MulticarrierSystem}(network::InterferingBroadcastChannel{System_t})
    I = get_num_BSs(network); Kc = network.num_MSs_per_cell; Lc = system.num_subcarriers
    Ns = get_num_MS_antennas(network); Ms = get_num_BS_antennas(network)

    coefs = Array(Array{Complex128, 3}, I*Kc, I)
    for k = 1:I*Kc
        serving_id = div(k - 1, Kc) + 1

        for i = 1:I
            for l = 1:Lc
                if i == serving_id
                    scale_factor = 1
                else
                    scale_factor = sqrt(network.alpha)
                end

                coefs[k,i,l] = scale_factor*(1/sqrt(2))*(randn(Ns[k], Ms[i]) + im*randn(Ns[k], Ms[i]))
            end
        end
    end

    return MulticarrierChannel(coefs, Ns, Ms, I*Kc, I, Lc)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "ZstdBuilder"
version = v"1.3.5"

# Collection of sources required to build ZstdBuilder
sources = [
    "https://github.com/facebook/zstd/archive/v$(version).tar.gz" =>
    "d6e1559e4cdb7c4226767d4ddc990bff5f9aab77085ff0d0490c828b025e2eea",
    "https://github.com/facebook/zstd/releases/download/v$(version)/zstd-v$(version)-win64.zip" =>
    "5dabee0ab14a92d73226cf8cb90d09a17b96b77542dffd173bdeeb3c3a3bab4a",
]

# Bash recipe for building across all platforms
script = raw"""
if [ $target = "x86_64-w64-mingw32" ]; then
    mkdir -p ${WORKSPACE}/destdir/bin/
    cp dll/libzstd.dll ${WORKSPACE}/destdir/bin/
else
    cd $WORKSPACE/srcdir/zstd-*
    make prefix=${prefix} install
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl, :eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libzstd", :libzstd)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using PDSampler, Test, Random, LinearAlgebra, Distributions

@testset "gaussian"    begin include("gaussian_test.jl")   end
@testset "t-gaussian"  begin include("tgaussian_test.jl")  end
@testset "logreg"      begin include("logreg_test.jl")     end
@testset "pmf"         begin include("pmf_test.jl")        end

@testset "geometry"    begin include("geometry_test.jl")   end
@testset "ippsampler"  begin include("ippsampler_test.jl") end
@testset "kernels"     begin include("kernels_test.jl")    end
@testset "path"        begin include("path_test.jl")       end
@testset "simulate"    begin include("simulate_test.jl")   end

@testset "l_event"     begin include("local_event_test.jl")        end
@testset "l_fgraph"    begin include("local_factorgraph_test.jl")  end
@testset "l_sim"       begin include("local_simulate_test.jl")     end

@testset "ex_gbps1"    begin include("ex_gbps1.jl") end
@testset "ex_genbps"   begin include("ex_genbps.jl") end

@testset "ex_lbps1"    begin include("ex_lbps1.jl") end


text:
using MCMCChains
using FFTW

using Random
using Statistics
using Test

@testset "ESS per second" begin
    c1 = Chains(randn(100,5, 1), info = (start_time=time(), stop_time = time()+1))
    c2 = Chains(randn(100,5, 1), info = (start_time=time()+1, stop_time = time()+2))
    c = chainscat(c1, c2)

    wall = MCMCChains.wall_duration(c)
    compute = MCMCChains.compute_duration(c)

    @test round(wall, digits=1) ≈ round(c2.info.stop_time - c1.info.start_time, digits=1)
    @test compute ≈ (MCMCChains.compute_duration(c1) + MCMCChains.compute_duration(c2))

    s = ess(c)
    @test length(s[:,:ess_per_sec]) == 5
    @test all(map(!ismissing, s[:,:ess_per_sec]))
end

@testset "copy and split" begin
    # check a matrix with even number of rows
    x = rand(50, 20)

    # check incompatible sizes
    @test_throws DimensionMismatch MCMCChains.copyto_split!(similar(x, 25, 20), x)
    @test_throws DimensionMismatch MCMCChains.copyto_split!(similar(x, 50, 40), x)

    y = similar(x, 25, 40)
    MCMCChains.copyto_split!(y, x)
    @test reshape(y, size(x)) == x

    # check a matrix with odd number of rows
    x = rand(51, 20)

    # check incompatible sizes
    @test_throws DimensionMismatch MCMCChains.copyto_split!(similar(x, 25, 20), x)
    @test_throws DimensionMismatch MCMCChains.copyto_split!(similar(x, 51, 40), x)

    MCMCChains.copyto_split!(y, x)
    @test reshape(y, 50, 20) == x[vcat(1:25, 27:51), :]
end

@testset "ESS and R̂ (IID samples)" begin
    Random.seed!(20)

    rawx = randn(10_000, 40, 10)

    # Repeat tests with different scales
    for scale in (1, 50, 100)
        x = scale * rawx

        ess_standard, rhat_standard = MCMCChains.ess_rhat(x)
        ess_standard2, rhat_standard2 = MCMCChains.ess_rhat(x; method = ESSMethod())
        ess_fft, rhat_fft = MCMCChains.ess_rhat(x; method = FFTESSMethod())
        ess_bda, rhat_bda = MCMCChains.ess_rhat(x; method = BDAESSMethod())

        # check that we get (roughly) the same results
        @test ess_standard == ess_standard2
        @test ess_standard ≈ ess_fft
        @test rhat_standard == rhat_standard2 == rhat_fft == rhat_bda

        # check that the estimates are reasonable
        @test all(x -> isapprox(x, 100_000; atol = 5_000), ess_standard)
        @test all(x -> isapprox(x, 100_000; atol = 5_000), ess_bda)
        @test all(x -> isapprox(x, 1; atol = 0.1), rhat_standard)

        # BDA method fluctuates more
        @test var(ess_standard) < var(ess_bda)
    end
end

@testset "ESS and R̂ (identical samples)" begin
    x = ones(10_000, 40, 10)

    ess_standard, rhat_standard = MCMCChains.ess_rhat(x)
    ess_standard2, rhat_standard2 = MCMCChains.ess_rhat(x; method = ESSMethod())
    ess_fft, rhat_fft = MCMCChains.ess_rhat(x; method = FFTESSMethod())
    ess_bda, rhat_bda = MCMCChains.ess_rhat(x; method = BDAESSMethod())

    # check that the estimates are all NaN
    for ess in (ess_standard, ess_standard2, ess_fft, ess_bda)
        @test all(isnan, ess)
    end
    for rhat in (rhat_standard, rhat_standard2, rhat_fft, rhat_bda)
        @test all(isnan, rhat)
    end
end

@testset "ESS and R̂ (chains)" begin
    x = rand(10_000, 40, 10)
    chain = Chains(x)

    for method in (ESSMethod(), FFTESSMethod(), BDAESSMethod())
        # analyze chain
        ess_df = ess(chain; method = method)

        # analyze array
        ess_array, rhat_array = MCMCChains.ess_rhat(x; method = method)

        @test ess_df[:,2] == ess_array
        @test ess_df[:,3] == rhat_array
    end
end

@testset "ESS and R̂ (single sample)" begin # check that issue #137 is fixed
    val = rand(1, 5, 3)
    chain = Chains(val, ["a", "b", "c", "d", "e"])

    for method in (ESSMethod(), FFTESSMethod(), BDAESSMethod())
        # analyze chain
        ess_df = ess(chain; method = method)

        # analyze array
        ess_array, rhat_array = MCMCChains.ess_rhat(val; method = method)

        @test ismissing(ess_df[:,2][1]) # since min(maxlag, niter - 1) = 0
        @test ismissing(ess_df[:,3][1])
    end
end


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

using Dierckx #for interpolation

"""
    𝚿 = numComp(Ψ, t, fs)
    𝚿 = numComp(Ψ, t)
    𝚿 = numComp(Ψ, fs)
    𝚿 = numComp(Ψ)

Create a 'numComp' consisting of a complex-valed signal `Ψ`, and time index `t`, and sampling frequency `fs`.

# Examples
```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), t, fs )
```
Another way to call 'numComp' by providing complex-valed signal `Ψ`, and sampling frequency `fs`.

```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
```
Another way to call 'numComp' by providing complex-valed signal `Ψ`, and time index `t`.

```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), t )
```
Another convenient way to call 'numComp' by providing complex-valed signal `Ψ`.

```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t) )
```
"""
struct numComp
  Ψ::Vector{ComplexF64}
  t::Vector{Float64}
  fs::Real
  s::Spline1D
  σ::Spline1D
end

# CONSTRUCTIONS
numComp(Ψ::Vector{ComplexF64}, t::Vector{Float64}, fs::Real) = numComp(Ψ, t, fs, Spline1D(t , real(Ψ)), Spline1D(t , imag(Ψ)))
numComp(Ψ::Vector{ComplexF64}, fs::Real) = numComp(Ψ, collect(0:length(Ψ)-1)/fs, fs)
numComp(Ψ::Vector{ComplexF64}, t::Vector{Float64}) = numComp(Ψ, t, 1/(t[2]-t[1]) )
numComp(Ψ::Vector{ComplexF64}, t::StepRangeLen) = numComp(Ψ, collect(t) )
numComp(Ψ::Vector{ComplexF64}) = numComp(Ψ, collect(0.0:1.0:length(Ψ)-1), 1.0)
numComp(Ψ::Vector{ComplexF64}, t::StepRangeLen, fs::Real) = numComp(Ψ, collect(t), fs)
numComp(Ψ::Vector{ComplexF64}, t::UnitRange, fs::Real) = numComp(Ψ, collect(t), fs)

# METHODS
function (𝚿::numComp)(t::Real)
  return ifelse(𝚿.t[1]<=t<=𝚿.t[end], 𝚿.s(t)+1im*𝚿.σ(t), NaN)
end
function (𝚿::numComp)(t::Vector{<:Real})
  return [ifelse(𝚿.t[1]<=tₖ<=𝚿.t[end], 𝚿.s(tₖ)+1im*𝚿.σ(tₖ), NaN) for tₖ∈t]
end
function (𝚿::numComp)(t::StepRangeLen)
  return [ifelse(𝚿.t[1]<=tₖ<=𝚿.t[end], 𝚿.s(tₖ)+1im*𝚿.σ(tₖ), NaN) for tₖ∈t]
end
function (𝚿::numComp)(t::UnitRange)
  return [ifelse(𝚿.t[1]<=tₖ<=𝚿.t[end], 𝚿.s(tₖ)+1im*𝚿.σ(tₖ), NaN) for tₖ∈t]
end

# DISPLAY
function display𝚿(x)
  text = ""
  T = typeof(x)
  for (name, typ) in zip(fieldnames(T), T.types)
      text *= "\n"*" "^2*"$name is $typ"
  end
  return text
end
Base.show(io::IO, x::numComp) = print(io, "numerical AM--FM component$(display𝚿(x))")

"""
    numTriplet

The default way to create a *numerical triplet* 'numTriplet' is by calling 'AMFMdemod()' on a *numerical component* 'numComp'.

# Examples
```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
𝐂₀ = AMFMdemod(𝚿₀)
```
"""
struct numTriplet
  Ψ::Vector{ComplexF64}
  t::Vector{Float64}
  fs::Float64
  a::Vector{Float64}
  ω::Vector{Float64}
  s::Vector{Float64}
  σ::Vector{Float64}
  θ::Vector{Float64}
end

# DISPLAY
function display𝐂(x)
  text = ""
  T = typeof(x)
  for (name, typ) in zip(fieldnames(T), T.types)
      text *= "\n"*" "^2*"$name is $typ"
  end
  return text
end
Base.show(io::IO, x::numTriplet) = print(io, "numerical canonical triplet$(display𝐂(x))")


text:
__precompile__()
module Utils

import GZip: gzopen, gzdopen
using ArgParse
using Bio.Seq
using SeqHax.ProgressLoggers

export  add_common_args,
        add_paired_args,
        open_zipped,
        foreach_read,
        foreach_readpair

function add_common_args(argparser)
    @add_arg_table argparser begin
        "--quiet", "-q"
            help = "Don't print any logging to STDERR"
            action = :store_true
    end
end

function add_paired_args(argparser)
    @add_arg_table argparser begin
        "--paired", "-p"
            help = "Process reads in paired mode"
            action = :store_true
    end
end

"Opens a file or stream, using gzip or not"
function open_zipped(file_or_io, mode="r", ziplevel::Int=0)
    stream = ifelse(startswith(mode, "w"), STDOUT, STDIN)
    zmode = ifelse(startswith(mode, "r"), "r", "w$ziplevel")
    if typeof(file_or_io) <: IO
        stream = file_or_io
    end
    if file_or_io == "-" || typeof(file_or_io) <: IO
        if ziplevel > 0
            return gzdopen(stream, zmode)
        else
            return stream
        end
    else
        if ziplevel > 0 || startswith(mode, "r")
            return gzopen(file_or_io, zmode)
        else
            return open(file_or_io, mode)
        end
    end
end

function foreach_read(fn::Function, readfile::AbstractString;
                      quiet::Bool=false)
    log = ProgressLogger(ifelse(quiet, 0, 100000), "sequences")
    nread = 0
    fp = open_zipped(readfile, "r")
    stream = open(fp, FASTQ, Seq.SANGER_QUAL_ENCODING)
    try
        read = eltype(stream)()
        while !eof(stream)
            update!(log, nread)
            read!(stream, read)
            fn(nread += 1, read)
        end
        return nread
    catch e
        throw(e)
    finally
        close(stream)
        flush!(log, nread)
    end
end

function foreach_readpair(fn::Function, readfile::AbstractString;
                          quiet::Bool=false)
    log = ProgressLogger(ifelse(quiet, 0, 100000), "reads")
    nread = 0
    fp = open_zipped(readfile, "r")
    stream = open(fp, FASTQ, Seq.SANGER_QUAL_ENCODING)
    try
        r1 = eltype(stream)()
        r2 = eltype(stream)()
        while !eof(stream)
            update!(log, nread)
            read!(stream, r1)
            if eof(stream)
                println(STDERR, "Non-paired file (unexpected eof after ",
                        "$nread reads")
                break
            end
            read!(stream, r2)
            fn(nread += 2, r1, r2)
        end
        return nread
    catch e
        throw(e)
    finally
        close(stream)
        flush!(log, nread)
    end
end

function foreach_readpair(fn::Function, readfile1::AbstractString,
                          readfile2::AbstractString; quiet::Bool=false)
    log = ProgressLogger(ifelse(quiet, 0, 100000), "reads")
    nread = 0
    r1fp = open_zipped(readfile1, "r")
    r2fp = open_zipped(readfile2, "r")
    r1s = open(r1fp, FASTQ, Seq.SANGER_QUAL_ENCODING)
    r2s = open(r2fp, FASTQ, Seq.SANGER_QUAL_ENCODING)
    try
        r1 = eltype(r1s)()
        r2 = eltype(r1s)() # deliberately using r1s here to keep them the same
        while !eof(r1s)
            update!(log, nread)
            read!(r1s, r1)
            if eof(r2s)
                println(STDERR, "Non-paired file (unexpected eof in R2 after ",
                        "$nread reads")
                break
            end
            read!(r2s, r2)
            fn(nread += 2, r1, r2)
        end
        if !eof(r2s)
            println(STDERR, "Non-paired file (unexpected eof in R1 after ",
                    "$nread reads")
        end
        return nread
    catch e
        throw(e)
    finally
        close(r1s)
        close(r2s)
        flush!(log, nread)
    end
end

end # module Utils