text:
using SQLite
using Test, Dates, Random, WeakRefStrings, Tables, DataFrames

import Base: +, ==

dbfile = joinpath(dirname(pathof(SQLite)),"../test/Chinook_Sqlite.sqlite")
dbfile2 = joinpath(dirname(pathof(SQLite)),"../test/test.sqlite")
cp(dbfile, dbfile2; force=true)
db = SQLite.DB(dbfile2)

# regular SQLite tests
ds = SQLite.Query(db, "SELECT name FROM sqlite_master WHERE type='table';") |> columntable
@test length(ds) == 1
@test keys(ds) == (:name,)
@test length(ds.name) == 11

results1 = SQLite.tables(db, columntable)
@test isequal(ds, results1)

results = SQLite.Query(db,"SELECT * FROM Employee;") |> DataFrame
@test size(results) == (8, 15)

SQLite.Query(db,"SELECT * FROM Album;")
SQLite.Query(db,"SELECT a.*, b.AlbumId
	FROM Artist a
	LEFT OUTER JOIN Album b ON b.ArtistId = a.ArtistId
	ORDER BY name;")

r = SQLite.execute!(db,"create table temp as select * from album")
r = SQLite.Query(db,"select * from temp limit 10") |> DataFrame
@test size(r) == (10,3)
SQLite.execute!(db,"alter table temp add column colyear int")
SQLite.execute!(db,"update temp set colyear = 2014")
r = SQLite.Query(db,"select * from temp limit 10") |> DataFrame
@test size(r) == (10,4)
@test all(Bool[x == 2014 for x in r[4]])
SQLite.execute!(db,"alter table temp add column dates blob")
stmt = SQLite.Stmt(db,"update temp set dates = ?")
SQLite.bind!(stmt,1,Dates.Date(2014,1,1))
SQLite.execute!(stmt)
finalize(stmt); stmt = nothing; GC.gc()
r = SQLite.Query(db,"select * from temp limit 10") |> DataFrame
@test size(r) == (10,5)
@test typeof(r[5][1]) == Date
@test all(Bool[x == Date(2014,1,1) for x in r[5]])
SQLite.execute!(db,"drop table temp")

dt = DataFrame([1.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 1.0])
tablename = dt |> SQLite.load!(db, "temp")
r = SQLite.Query(db,"select * from $tablename") |> DataFrame
@test size(r) == (5, 5)
@test names(r) == [:x1, :x2, :x3, :x4, :x5]
SQLite.drop!(db,"$tablename")

dt = DataFrame(zeros(5, 5))
tablename = dt |> SQLite.load!(db, "temp")
r = SQLite.Query(db, "select * from $tablename") |> DataFrame
@test size(r) == (5,5)
@test all([i for i in r[1]] .== 0.0)
@test all([typeof(i) for i in r[1]] .== Float64)
SQLite.drop!(db, "$tablename")

dt = DataFrame(zeros(5, 5))
tablename = dt |> SQLite.load!(db, "temp")
r = SQLite.Query(db, "select * from $tablename") |> DataFrame
@test size(r) == (5,5)
@test all([i for i in r[1]] .== 0)
@test all([typeof(i) for i in r[1]] .== Float64)

dt = DataFrame(ones(Int, 5, 5))
tablename = dt |> SQLite.load!(db, "temp")
r = SQLite.Query(db, "select * from $tablename") |> DataFrame
@test size(r) == (10,5)
@test r[1] == [0,0,0,0,0,1,1,1,1,1]
@test all([typeof(i) for i in r[1]] .== Float64)
SQLite.drop!(db, "$tablename")

rng = Dates.Date(2013):Dates.Day(1):Dates.Date(2013,1,5)
dt = DataFrame(i=collect(rng), j=collect(rng))
tablename = dt |> SQLite.load!(db, "temp")
r = SQLite.Query(db, "select * from $tablename") |> DataFrame
@test size(r) == (5,2)
@test all([i for i in r[1]] .== collect(rng))
@test all([typeof(i) for i in r[1]] .== Dates.Date)
SQLite.drop!(db, "$tablename")

SQLite.Query(db, "CREATE TABLE temp AS SELECT * FROM Album")
r = SQLite.Query(db, "SELECT * FROM temp LIMIT ?"; values=[3]) |> DataFrame
@test size(r) == (3,3)
r = SQLite.Query(db, "SELECT * FROM temp WHERE Title LIKE ?"; values=["%time%"]) |> DataFrame
@test r[1] == [76, 111, 187]
SQLite.Query(db, "INSERT INTO temp VALUES (?1, ?3, ?2)"; values=[0,0,"Test Album"])
r = SQLite.Query(db, "SELECT * FROM temp WHERE AlbumId = 0") |> DataFrame
@test r[1][1] === 0
@test r[2][1] == "Test Album"
@test r[3][1] === 0
SQLite.drop!(db, "temp")

SQLite.Query(db, "CREATE TABLE temp AS SELECT * FROM Album")
r = SQLite.Query(db, "SELECT * FROM temp LIMIT :a"; values=Dict(:a => 3)) |> DataFrame
@test size(r) == (3,3)
r = SQLite.Query(db, "SELECT * FROM temp WHERE Title LIKE @word"; values=Dict(:word => "%time%")) |> DataFrame
@test r[1] == [76, 111, 187]
SQLite.Query(db, "INSERT INTO temp VALUES (@lid, :title, \$rid)"; values=Dict(:rid => 0, :lid => 0, :title => "Test Album"))
r = SQLite.Query(db, "SELECT * FROM temp WHERE AlbumId = 0") |> DataFrame
@test r[1][1] === 0
@test r[2][1] == "Test Album"
@test r[3][1] === 0
SQLite.drop!(db, "temp")

register(db, SQLite.regexp, nargs=2, name="regexp")
r = SQLite.Query(db, SQLite.@sr_str("SELECT LastName FROM Employee WHERE BirthDate REGEXP '^\\d{4}-08'")) |> DataFrame
@test r[1][1] == "Peacock"

triple(x) = 3x
@test_throws AssertionError SQLite.register(db, triple, nargs=186)
SQLite.register(db, triple, nargs=1)
r = SQLite.Query(db, "SELECT triple(Total) FROM Invoice ORDER BY InvoiceId LIMIT 5") |> DataFrame
s = SQLite.Query(db, "SELECT Total FROM Invoice ORDER BY InvoiceId LIMIT 5") |> DataFrame
for (i, j) in zip(r[1], s[1])
    @test abs(i - 3*j) < 0.02
end

SQLite.@register db function add4(q)
    q+4
end
r = SQLite.Query(db, "SELECT add4(AlbumId) FROM Album") |> DataFrame
s = SQLite.Query(db, "SELECT AlbumId FROM Album") |> DataFrame
@test r[1][1] == s[1][1] + 4

SQLite.@register db mult(args...) = *(args...)
r = SQLite.Query(db, "SELECT Milliseconds, Bytes FROM Track") |> DataFrame
s = SQLite.Query(db, "SELECT mult(Milliseconds, Bytes) FROM Track") |> DataFrame
@test (r[1][1] * r[2][1]) == s[1][1]
t = SQLite.Query(db, "SELECT mult(Milliseconds, Bytes, 3, 4) FROM Track") |> DataFrame
@test (r[1][1] * r[2][1] * 3 * 4) == t[1][1]

SQLite.@register db sin
u = SQLite.Query(db, "select sin(milliseconds) from track limit 5") |> DataFrame
@test all(-1 .< convert(Vector{Float64},u[1]) .< 1)

SQLite.register(db, hypot; nargs=2, name="hypotenuse")
v = SQLite.Query(db, "select hypotenuse(Milliseconds,bytes) from track limit 5") |> DataFrame
@test [round(Int,i) for i in v[1]] == [11175621,5521062,3997652,4339106,6301714]

SQLite.@register db str2arr(s) = Vector{UInt8}(s)
r = SQLite.Query(db, "SELECT str2arr(LastName) FROM Employee LIMIT 2") |> DataFrame
@test r[1][2] == UInt8[0x45,0x64,0x77,0x61,0x72,0x64,0x73]

SQLite.@register db big
r = SQLite.Query(db, "SELECT big(5)") |> DataFrame
@test r[1][1] == big(5)

doublesum_step(persist, current) = persist + current
doublesum_final(persist) = 2 * persist
SQLite.register(db, 0, doublesum_step, doublesum_final, name="doublesum")
r = SQLite.Query(db, "SELECT doublesum(UnitPrice) FROM Track") |> DataFrame
s = SQLite.Query(db, "SELECT UnitPrice FROM Track") |> DataFrame
@test abs(r[1][1] - 2*sum(convert(Vector{Float64},s[1]))) < 0.02

mycount(p, c) = p + 1
SQLite.register(db, 0, mycount)
r = SQLite.Query(db, "SELECT mycount(TrackId) FROM PlaylistTrack") |> DataFrame
s = SQLite.Query(db, "SELECT count(TrackId) FROM PlaylistTrack") |> DataFrame
@test r[1][1] == s[1][1]

bigsum(p, c) = p + big(c)
SQLite.register(db, big(0), bigsum)
r = SQLite.Query(db, "SELECT bigsum(TrackId) FROM PlaylistTrack") |> DataFrame
s = SQLite.Query(db, "SELECT TrackId FROM PlaylistTrack") |> DataFrame
# @test r[1][1] == big(sum(convert(Vector{Int},s[1])))

SQLite.Query(db, "CREATE TABLE points (x INT, y INT, z INT)")
SQLite.Query(db, "INSERT INTO points VALUES (?, ?, ?)"; values=[1, 2, 3])
SQLite.Query(db, "INSERT INTO points VALUES (?, ?, ?)"; values=[4, 5, 6])
SQLite.Query(db, "INSERT INTO points VALUES (?, ?, ?)"; values=[7, 8, 9])
mutable struct Point3D{T<:Number}
    x::T
    y::T
    z::T
end
==(a::Point3D, b::Point3D) = a.x == b.x && a.y == b.y && a.z == b.z
+(a::Point3D, b::Point3D) = Point3D(a.x + b.x, a.y + b.y, a.z + b.z)
sumpoint(p::Point3D, x, y, z) = p + Point3D(x, y, z)
SQLite.register(db, Point3D(0, 0, 0), sumpoint)
r = SQLite.Query(db, "SELECT sumpoint(x, y, z) FROM points") |> DataFrame
@test r[1][1] == Point3D(12, 15, 18)
SQLite.drop!(db, "points")

db2 = SQLite.DB()
SQLite.Query(db2, "CREATE TABLE tab1 (r REAL, s INT)")

@test_throws SQLite.SQLiteException SQLite.drop!(db2, "nonexistant")
# should not throw anything
SQLite.drop!(db2, "nonexistant", ifexists=true)
# should drop "tab2"
SQLite.drop!(db2, "tab2", ifexists=true)
@test !in("tab2", SQLite.tables(db2)[1])

SQLite.drop!(db, "sqlite_stat1")
@test size(SQLite.tables(db)) == (11,1)

finalize(db); db = nothing; GC.gc(); GC.gc();

#Test removeduplicates!
db = SQLite.DB() #In case the order of tests is changed
ints = Int64[1,1,2,2,3]
strs = String["A", "A", "B", "C", "C"]
dt = DataFrame(ints=ints, strs=strs)
tablename = dt |> SQLite.load!(db, "temp")
SQLite.removeduplicates!(db, "temp", ["ints", "strs"]) #New format
dt3 = SQLite.Query(db, "Select * from temp") |> DataFrame
@test dt3[1][1] == 1
@test dt3[2][1] == "A"
@test dt3[1][2] == 2
@test dt3[2][2] == "B"
@test dt3[1][3] == 2
@test dt3[2][3] == "C"

# issue #104
db = SQLite.DB() #In case the order of tests is changed
SQLite.execute!(db, "CREATE TABLE IF NOT EXISTS tbl(a  INTEGER);")
stmt = SQLite.Stmt(db, "INSERT INTO tbl (a) VALUES (@a);")
SQLite.bind!(stmt, "@a", 1)

binddb = SQLite.DB()
SQLite.Query(binddb, "CREATE TABLE temp (n NULL, i6 INT, f REAL, s TEXT, a BLOB)")
SQLite.Query(binddb, "INSERT INTO temp VALUES (?1, ?2, ?3, ?4, ?5)"; values=Any[missing, convert(Int64,6), 6.4, "some text", b"bytearray"])
r = SQLite.Query(binddb, "SELECT * FROM temp") |> DataFrame
@test isa(r[1][1], Missing)
@test isa(r[2][1], Int)
@test isa(r[3][1], Float64)
@test isa(r[4][1], AbstractString)
@test isa(r[5][1], Base.CodeUnits)
SQLite.Query(binddb, "CREATE TABLE blobtest (a BLOB, b BLOB)")
SQLite.Query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", b"b"])
SQLite.Query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", BigInt(2)])
mutable struct Point{T}
    x::T
    y::T
end
==(a::Point, b::Point) = a.x == b.x && a.y == b.y
p1 = Point(1, 2)
p2 = Point(1.3, 2.4)
SQLite.Query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", p1])
SQLite.Query(binddb, "INSERT INTO blobtest VALUES (?1, ?2)"; values=Any[b"a", p2])
r = SQLite.Query(binddb, "SELECT * FROM blobtest"; stricttypes=false) |> DataFrame
for value in r[1]
    @test value == b"a"
end
@test r[2][3] == p1
@test r[2][4] == p2
############################################

#test for #158
@test_throws SQLite.SQLiteException SQLite.DB("nonexistentdir/not_there.db")

#test for #180 (Query)
param = "Hello!"
query = SQLite.Query(SQLite.DB(), "SELECT ?1 UNION ALL SELECT ?1", values = Any[param])
param = "x"
for row in query
    @test row[1] == "Hello!"
    GC.gc() # this must NOT garbage collect the "Hello!" bound value
end

#test for #180 (bind! and clear!)
params = tuple("string", UInt8[1, 2, 3]) # parameter types that can be finalized
wkdict = WeakKeyDict{Any, Any}(param => 1 for param in params)
stmt = SQLite.Stmt(SQLite.DB(), "SELECT ?, ?")
SQLite.bind!(stmt, params)
params = "x"
GC.gc() # this MUST NOT garbage collect any of the bound values
@test length(wkdict) == 2
SQLite.clear!(stmt)
GC.gc() # this will garbage collect the no longer bound values
@test isempty(wkdict)


text:
using PyPlot
using LinearAlgebra

export plot_state_space, plot_quadrotor, plot_target, plot_states,
  plot_trajectory, plot_robot, plot_observation, visualize_filter,
  visualize_filters

const object_scale = 13^2
const default_linewidth = 2.0
const default_cmap = "viridis"

function plot_state_space(g::Grid; color=:k)
  states = get_states(g)
  xs = [x[1] for x in states]
  ys = [x[2] for x in states]

  scatter(xs, ys, color=color)
end

# Plots a visualization of a quadrotor
rotate(theta) = [cos(theta) -sin(theta); sin(theta) cos(theta)]
function circle(p; radius=1, scale=1, dth=0.1)
  scale*hcat(map(x->p+[radius*cos(x);radius*sin(x)], 0:dth:(2*pi)+dth)...)
end
plot_quadrotor(p::State; kwargs...) = plot_quadrotor([p...]; kwargs...)
function plot_quadrotor(p::Vector; scale=0.20,
                        color="blue",
                        linewidth=default_linewidth,
                        theta=pi/6,
                        alpha=1.0)
  ret = []

  rotor_radius = 0.5

  rot = rotate(theta)

  circle_ps = rot*scale*[[1;0] [0;1] [-1;0] [0;-1]]

  p3 = 0.0
  if length(p) == 3
    p3 = p[3]
  end

  p2 = p[1:2]

  # Plot X shape
  l1 = [p2+scale*rot[:,1] p2-scale*rot[:,1]]
  l2 = [p2+scale*rot[:,2] p2-scale*rot[:,2]]
  append!(ret, plot(l1[1,:][:], l1[2,:][:], p3*ones(2), color=color, linewidth=linewidth,
                  alpha=alpha))
  append!(ret, plot(l2[1,:][:], l2[2,:][:], p3*ones(2), color=color, linewidth=linewidth,
                  alpha = alpha))

  # Plot and fill circles
  for ii = 1:size(circle_ps,2)
    c = circle(p2+circle_ps[:,ii], radius=scale*rotor_radius)
    append!(ret, plot(c[1,:][:], c[2,:][:], p3*ones(size(c,2)), color=color,
                    linestyle="--", linewidth=0.5linewidth, alpha=alpha))
    append!(ret, fill(c[1,:][:], c[2,:][:], p3*ones(size(c,2)), color=:black,
                    alpha=0.3alpha))
  end

  ret
end

function plot_target(p; color="red", kwargs...)
  [scatter(p[1], p[2]; color=color, marker="^", edgecolors=:k,
           s=object_scale, kwargs...)]
end

# When a robot or target crosses a boundary of the grid, that creates a
# discontinuity.
#
# Cut the trajectory into continuous parts and add segments going off the grid
# to visualize the wrap around.
function wrap_discontinuities(states)
  state_ranges = continuous_ranges(states)
  ranges = convert(Array{Array{Tuple{Float64,Float64}}}, state_ranges)

  # Add segments going off the grid to each discontinuity
  for ii = 1:length(ranges) - 1
    disc_start = last(ranges[ii])
    disc_end = first(ranges[ii+1])

    dir = disc_end .- disc_start
    off_grid = -0.5 .* dir ./ norm(dir)

    push!(ranges[ii], disc_start .+ off_grid)
    insert!(ranges[ii+1], 1, disc_end .- off_grid)
  end

  ranges
end

function plot_states(states::Array{State}; linewidth=default_linewidth,
                     kwargs...)
  vcat(map(wrap_discontinuities(states)) do range
    xs = [x[1] for x in range]
    ys = [x[2] for x in range]

    plot(xs, ys; linewidth=linewidth, kwargs...)
  end...)
end

function plot_trajectory(states; color=:red, linewidth=default_linewidth,
                         kwargs...)

  first = states[1]
  last = states[end]

  # start end
  ret = []
  push!(ret, scatter(first[1], first[2], color=color, marker="X", edgecolors=:k,
                     s=object_scale))
  push!(ret, scatter(last[1], last[2], color=color, marker="^", edgecolors=:k,
                     s=object_scale))

  range_plots = map(wrap_discontinuities(states)) do range
    xs = [x[1] for x in range]
    ys = [x[2] for x in range]
    plot(xs, ys; color=color, linewidth=default_linewidth, kwargs...)
  end

  # return values are arrays of plot objects, concatenate
  vcat(ret, range_plots...)
end

function plot_robot(state; color=:blue, kwargs...)
  [scatter(state[1], state[2], color=color, marker="v", edgecolors=:k,
           s=object_scale)]
end

function plot_observation(state, range; linestyle="-", kwargs...)
  plot_circle([state...]; radius=range, linestyle=linestyle, kwargs...)
end

# Visualize a discrete histogram filter
function visualize_filter(filter::AbstractFilter; kwargs...)
  visualize_filters([filter]; kwargs...)
end
function visualize_filters(filters::Vector{<:AbstractFilter};
                           show_colorbar = false,
                           alpha=0.75)
  if length(filters) == 0
    error("No filters to visualize")
  end

  range = get_range(filters[1])

  # Adjust bounds because indices are in centers
  limits = [range[1][[1,end]]'; range[2][[1,end]]'] + repeat([-0.5 0.5], 2)

  data = sum(get_data, filters)

  # Note: visualize filter does not normalize by default
  visualize_pdf(Matrix(data), show_colorbar = show_colorbar, limits=limits,
               alpha=alpha, cmap=default_cmap)
end

ext: jl

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
### Format registry infrastructure
@compat abstract type OS end
@compat abstract type Unix <: OS end
immutable Windows <: OS end
immutable OSX <: Unix end
immutable Linux <: Unix end

immutable LOAD end
immutable SAVE end

split_predicates(list) = filter(x-> x <: OS, list), filter(x-> !(x <: OS), list)
applies_to_os(os::Vector) = isempty(os) || any(applies_to_os, os)
applies_to_os{O <: OS}(os::Type{O})              = false

applies_to_os{U <: Unix}(os::Type{U}) = is_unix()
applies_to_os(os::Type{Windows}) = is_windows()
applies_to_os(os::Type{OSX}) = is_apple()
applies_to_os(os::Type{Linux}) = is_linux()

function add_loadsave(format, predicates)
    library = shift!(predicates)
    os, loadsave = split_predicates(predicates)
    if applies_to_os(os)
        if isempty(loadsave) || (LOAD in loadsave)
            add_loader(format, library)
        end
        if isempty(loadsave) || (SAVE in loadsave)
            add_saver(format, library)
        end
    end
end

"""
`DataFormat{sym}()` indicates a known binary or text format of kind `sym`,
where `sym` is always a symbol. For example, a .csv file might have
`DataFormat{:CSV}()`.

An easy way to write `DataFormat{:CSV}` is `format"CSV"`.
"""
immutable DataFormat{sym} end

macro format_str(s)
    :(DataFormat{$(Expr(:quote, Symbol(s)))})
end

const unknown_df = DataFormat{:UNKNOWN}

"""
`unknown(f)` returns true if the format of `f` is unknown.
"""
unknown(::Type{format"UNKNOWN"})    = true
unknown{sym}(::Type{DataFormat{sym}}) = false

const ext2sym    = Dict{String, Union{Symbol,Vector{Symbol}}}()
const magic_list = Vector{Pair}(0)    # sorted, see magic_cmp below
const sym2info   = Dict{Symbol,Any}()  # Symbol=>(magic, extension)
const magic_func = Vector{Pair}(0)    # for formats with complex magic #s

function add_format(fmt, magic, extension, load_save_libraries...)
    add_format(fmt, magic, extension)
    for library in load_save_libraries
        add_loadsave(fmt, library)
    end
    fmt
end

"""
`add_format(fmt, magic, extention)` registers a new `DataFormat`.
For example:

    add_format(format"PNG", (UInt8[0x4d,0x4d,0x00,0x2b], UInt8[0x49,0x49,0x2a,0x00]), [".tiff", ".tif"])
    add_format(format"PNG", [0x89,0x50,0x4e,0x47,0x0d,0x0a,0x1a,0x0a], ".png")
    add_format(format"NRRD", "NRRD", [".nrrd",".nhdr"])

Note that extensions, magic numbers, and format-identifiers are case-sensitive.
"""
function add_format{sym}(fmt::Type{DataFormat{sym}}, magic::Union{Tuple,AbstractVector,String}, extension)
    haskey(sym2info, sym) && error("format ", fmt, " is already registered")
    m = canonicalize_magic(magic)
    rng = searchsorted(magic_list, m, lt=magic_cmp)
    if !isempty(m) && !isempty(rng)
        error("magic bytes ", m, " are already registered")
    end
    insert!(magic_list, first(rng), Pair(m, sym))  # m=>sym in 0.4
    sym2info[sym] = (m, extension)
    add_extension(extension, sym)
    fmt
end

# for multiple magic bytes
function add_format{sym, T <: Vector{UInt8}, N}(fmt::Type{DataFormat{sym}}, magics::NTuple{N, T}, extension)
    haskey(sym2info, sym) && error("format ", fmt, " is already registered")
    magics = map(canonicalize_magic, magics)
    for magic in magics
        rng = searchsorted(magic_list, magic, lt=magic_cmp)
        if !isempty(magic) && !isempty(rng)
            error("magic bytes ", magic, " are already registered")
        end
        insert!(magic_list, first(rng), Pair(magic, sym))  # m=>sym in 0.4
    end
    sym2info[sym] = (magics, extension)
    add_extension(extension, sym)
    fmt
end

# For when "magic" is supplied as a function (see the HDF5 example in
# registry.jl)
function add_format{sym}(fmt::Type{DataFormat{sym}}, magic, extension)
    haskey(sym2info, sym) && error("format ", fmt, " is already registered")
    push!(magic_func, Pair(magic,sym))  # magic=>sym in 0.4
    sym2info[sym] = (magic, extension)
    add_extension(extension, sym)
    fmt
end

"""
`del_format(fmt::DataFormat)` deletes `fmt` from the format registry.
"""
function del_format{sym}(fmt::Type{DataFormat{sym}})
    magic, extension = sym2info[sym]
    del_magic(magic, sym)
    delete!(sym2info, sym)
    del_extension(extension)
    nothing
end

# Deletes mutliple magic bytes
del_magic(magic::Tuple, sym) = for m in magic
    del_magic(m, sym)
end
# Deletes single magic bytes
function del_magic{N}(magic::NTuple{N, UInt8}, sym)
    rng = searchsorted(magic_list, magic, lt=magic_cmp)
    if length(magic) == 0
        fullrng = rng
        found = false
        for idx in fullrng
            if last(magic_list[idx]) == sym
                rng = idx:idx
                found = true
                break
            end
        end
        found || error("format ", sym, " not found")
    end
    @assert length(rng) == 1
    deleteat!(magic_list, first(rng))
    nothing
end

function del_magic(magic::Function, sym)
    deleteat!(magic_func, findfirst(magic_func, Pair(magic,sym)))
    nothing
end

"""
`info(fmt)` returns the magic bytes/extension information for
`DataFormat` `fmt`.
"""
Base.info{sym}(::Type{DataFormat{sym}}) = sym2info[sym]

canonicalize_magic{N}(m::NTuple{N,UInt8}) = m
canonicalize_magic(m::AbstractVector{UInt8}) = tuple(m...)
canonicalize_magic(m::String) = canonicalize_magic(Vector{UInt8}(m))

function add_extension(ext::String, sym)
    if haskey(ext2sym, ext)
        v = ext2sym[ext]
        if isa(v, Symbol)
            ext2sym[ext] = Symbol[v, sym]
        else
            push!(ext2sym[ext], sym)
        end
        return
    end
    ext2sym[ext] = sym
end
function add_extension(ext::Union{Array,Tuple}, sym)
    for e in ext
        add_extension(e, sym)
    end
end

del_extension(ext::String) = delete!(ext2sym, ext)
function del_extension(ext::Union{Array,Tuple})
    for e in ext
        del_extension(e)
    end
end

# magic_cmp results in magic_list being sorted in order of increasing
# length(magic), then (among tuples with the same length) in
# dictionary order. This ordering has the advantage that you can
# incrementally read bytes from the stream without worrying that
# you'll encounter an EOF yet still have potential matches later in
# the list.
function magic_cmp(p::Pair, t::Tuple)
    pt = first(p)
    lp, lt = length(pt), length(t)
    lp < lt && return true
    lp > lt && return false
    pt < t
end
function magic_cmp(t::Tuple, p::Pair)
    pt = first(p)
    lp, lt = length(pt), length(t)
    lt < lp && return true
    lt > lp && return false
    t < pt
end

@compat abstract type Formatted{F<:DataFormat} end  # A specific file or stream

"""
`File(fmt, filename)` indicates that `filename` is a file of known
DataFormat `fmt`.  For example, `File{fmtpng}(filename)` would indicate a PNG
file.
"""
immutable File{F<:DataFormat} <: Formatted{F}
    filename::String
end
File{sym}(fmt::Type{DataFormat{sym}}, filename) = File{fmt}(filename)

"""
`filename(file)` returns the filename associated with `File` `file`.
"""
filename(f::File) = f.filename

"""
`file_extension(file)` returns the file extension associated with `File` `file`.
"""
file_extension(f::File) = splitext(filename(f))[2]

"""
`Stream(fmt, io, [filename])` indicates that the stream `io` is
written in known `Format`.  For example, `Stream{PNG}(io)` would
indicate PNG format.  If known, the optional `filename` argument can
be used to improve error messages, etc.
"""
immutable Stream{F<:DataFormat,IOtype<:IO} <: Formatted{F}
    io::IOtype
    filename::Nullable{String}
end

Stream{F<:DataFormat}(::Type{F}, io::IO) = Stream{F,typeof(io)}(io, Nullable{String}())
Stream{F<:DataFormat}(::Type{F}, io::IO, filename::AbstractString) = Stream{F,typeof(io)}(io,String(filename))
Stream{F<:DataFormat}(::Type{F}, io::IO, filename) = Stream{F,typeof(io)}(io,filename)
Stream{F}(file::File{F}, io::IO) = Stream{F,typeof(io)}(io,filename(file))

"`stream(s)` returns the stream associated with `Stream` `s`"
stream(s::Stream) = s.io

"""
`filename(stream)` returns a nullable-string of the filename
associated with `Stream` `stream`.
"""
filename(s::Stream) = s.filename

"""
`file_extension(file)` returns a nullable-string for the file extension associated with `Stream` `stream`.
"""
function file_extension(f::Stream)
    isnull(filename(f)) && return filename(f)
    splitext(get(filename(f)))[2]
end

# Note this closes the stream. It's useful when you've opened
# the file to check the magic bytes, but don't want to leave
# a dangling stream.
function file!{F}(strm::Stream{F})
    f = filename(strm)
    if isnull(f)
        error("filename unknown")
    end
    close(strm.io)
    File{F}(get(f))
end

# Implement standard I/O operations for File and Stream
@inline function Base.open{F<:DataFormat}(file::File{F}, args...)
    fn = filename(file)
    Stream(F, open(fn, args...), abspath(fn))
end
Base.close(s::Stream) = close(stream(s))

Base.position(s::Stream) = position(stream(s))
Base.seek(s::Stream, offset::Integer) = (seek(stream(s), offset); s)
Base.seekstart(s::Stream) = (seekstart(stream(s)); s)
Base.seekend(s::Stream) = (seekend(stream(s)); s)
Base.skip(s::Stream, offset::Integer) = (skip(stream(s), offset); s)
Base.eof(s::Stream) = eof(stream(s))

@inline Base.read(s::Stream, args...)  = read(stream(s), args...)
Base.read!(s::Stream, array::Array) = read!(stream(s), array)
@inline Base.write(s::Stream, args...) = write(stream(s), args...)
# Note: we can't sensibly support the all keyword. If you need that,
# call read(stream(s), ...; all=value) manually
Base.readbytes!(s::Stream, b) = readbytes!(stream(s), b)
Base.readbytes!(s::Stream, b, nb) = readbytes!(stream(s), b, nb)
Base.read(s::Stream) = read(stream(s))
Base.read(s::Stream, nb) = read(stream(s), nb)
Base.flush(s::Stream) = flush(stream(s))

Base.isreadonly(s::Stream) = isreadonly(stream(s))
Base.isopen(s::Stream) = isopen(stream(s))

"`magic(fmt)` returns the magic bytes of format `fmt`"
magic{F<:DataFormat}(fmt::Type{F}) = UInt8[info(fmt)[1]...]

"""
`skipmagic(s)` sets the position of `Stream` `s` to be just after the magic bytes.
For a plain IO object, you can use `skipmagic(io, fmt)`.
"""
skipmagic{F}(s::Stream{F}) = (skipmagic(stream(s), F); s)
function skipmagic{sym}(io, fmt::Type{DataFormat{sym}})
    magic, _ = sym2info[sym]
    skipmagic(io, magic)
    nothing
end
skipmagic(io, magic::Function) = nothing
skipmagic{N}(io, magic::NTuple{N,UInt8}) = seek(io, length(magic))
function skipmagic(io, magic::Tuple)
    lengths = map(length, magic)
    all(x->lengths[1] == x, lengths) && return seek(io, lengths[1]) # it doesn't matter what magic bytes get skipped as they all have the same length
    magic = [magic...]
    sort!(magic, lt=(a,b)-> length(a)>= length(b)) # start with longest first, to avoid overlapping magic bytes
    seekend(io)
    len = position(io)
    seekstart(io)
    filter!(x-> length(x) <= len, magic) # throw out magic bytes that are longer than IO
    tmp = read(io, length(first(magic))) # now, first is both the longest and guaranteed to fit into io, so we can just read the bytes
    for m in magic
        if magic_equal(m, tmp)
            seek(io, length(m))
            return nothing
        end
    end
    error("tried to skip magic bytes of an IO that does not contain the magic bytes of the format. IO: $io")
end
function magic_equal(magic, buffer)
    for (i,elem) in enumerate(magic)
        buffer[i] != elem && return false
    end
    true
end

unknown{F}(::File{F}) = unknown(F)
unknown{F}(::Stream{F}) = unknown(F)

"""
`query(filename)` returns a `File` object with information about the
format inferred from the file's extension and/or magic bytes."""
function query(filename::AbstractString)
    _, ext = splitext(filename)
    if haskey(ext2sym, ext)
        sym = ext2sym[ext]
        no_magic = !hasmagic(sym)
        if lensym(sym) == 1 && (no_magic || !isfile(filename)) # we only found one candidate and there is no magic bytes, or no file, trust the extension
            return File{DataFormat{sym}}(filename)
        elseif !isfile(filename) && lensym(sym) > 1
            return File{DataFormat{sym[1]}}(filename)
        end
        if no_magic && !hasfunction(sym)
            error("Some formats with extension ", ext, " have no magic bytes; use `File{format\"FMT\"}(filename)` to resolve the ambiguity.")
        end
    end
    !isfile(filename) && return File{unknown_df}(filename) # (no extension || no magic byte) && no file
    # Otherwise, check the magic bytes
    file!(query(open(filename), abspath(filename)))
end

lensym(s::Symbol) = 1
lensym(v::Vector) = length(v)

hasmagic(s::Symbol) = hasmagic(sym2info[s][1])
hasmagic(v::Vector) = any(hasmagic, v)

hasmagic(t::Tuple) = !isempty(t)
hasmagic(::Any) = false   # for when magic is a function

hasfunction(s::Symbol) = hasfunction(sym2info[s][1])
hasfunction(v::Vector) = any(hasfunction, v)
hasfunction(s::Any) = true #has function
hasfunction(s::Tuple) = false #has magic

"""
`query(io, [filename])` returns a `Stream` object with information about the
format inferred from the magic bytes."""
query(io::IO, filename) = query(io, Nullable(String(filename)))

function query(io::IO, filename::Nullable{String}=Nullable{String}())
    magic = Vector{UInt8}(0)
    pos = position(io)
    for p in magic_list
        m = first(p)
        length(m) == 0 && continue
        while length(m) > length(magic)
            if eof(io)
                seek(io, pos)
                return Stream{unknown_df,typeof(io)}(io, filename)
            end
            push!(magic, read(io, UInt8))
        end
        if iter_eq(magic, m)
            seek(io, pos)
            return Stream{DataFormat{last(p)},typeof(io)}(io, filename)
        end
    end
    if seekable(io)
        for p in magic_func
            seek(io, pos)
            f = first(p)
            if f(io)
                return Stream{DataFormat{last(p)},typeof(io)}(seek(io, pos), filename)
            end
        end
        seek(io, pos)
    end
    Stream{unknown_df,typeof(io)}(io, filename)
end

seekable(io::IOBuffer) = io.seekable
seekable(::IOStream) = true
seekable(::Any) = false

function iter_eq(A, B)
    length(A) == length(B) || return false
    i,j = 1,1
    for _=1:length(A)
        a=A[i]; b=B[j]
        a == b && (i+=1; j+=1; continue)
        a == UInt32('\r') && (i+=1; continue) # this seems like the shadiest solution to deal with windows \r\n
        b == UInt32('\r') && (j+=1; continue)
        return false #now both must be unequal, and no \r windows excemption any more
    end
    true
end


text:
#### Beta diversity calculation functions

## Function to calculate beta diversity statistics
"""
    betadiv(Y::Matrix)

Computes the beta diversity statistics from the community matrix `Y` based on
Legendre & De Cáceres (2013). This function was based on the `beta.div` R
function from the supplementary material of that paper, and results were also
tested against the `beta.div` function from `adespatial`.
"""
function betadiv(Y::Matrix)
    @assert !any(isnothing, Y) "Y must only contain observed sites and no nothing values"
    # S -> squared deviations from column mean
    S = (Y .- mean(Y; dims=1)) .^ 2.0
    # SStotal -> total sum of squares
    SStotal = sum(S)
    # BDtotal -> index of beta diversity, unbiased & comparable estimator of Var(Y)
    BDtotal = SStotal / (size(Y, 1) - 1)
    # SSj -> sum of squares for species j
    SSj = sum(S; dims=1)
    # SCBDj -> species contribution to beta diversity (species j, relative)
    SCBDj = SSj ./ SStotal
    # SSi -> sum of squares for site i
    SSi = sum(S; dims=2)
    # LCBD -> local contribution to beta diversity (site i, relative)
    LCBDi = SSi ./ SStotal
    # Combine results in tuple
    res = (
        S=S, SStotal=SStotal, BDtotal=BDtotal, SSj=SSj, SCBDj=SCBDj, SSi=SSi, LCBDi=LCBDi
    )
    return res
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
"""
    λmin,xmin = minimize_geig_bracketing(A,B,C,λ₀,x₀)

Minimize the (maximum) generalized eigenvalue for a pair of matrices `A`, `B` under the additional constraint that `B`≻0 and `C`≻0.

All the three matrices are actually affine functions of a vector `x`, that is, `A(x) = A₀+A₁x₁+…+Aₙxₙ` (and the same holds for `B(x)` and `C(x)`). They are entered as (1D) arrays of matrices.

An initial guesses for the upper bound on generalized eigenvalues `λ₀` and the solution vector `x₀` must be provided.

# Example
```julia
A₀ = [1.0 2.0; 2.0 3.0]
A₁ = [1.0 0.0; 0.0 1.0]
A₂ = [4.0 5.0; 5.0 6.0]
A = [A₀,A₁,A₂]

B₀ = [1.0 5.0; 5.0 5.0]
B₁ = [2.0 3.0; 3.0 4.0]
B₂ = [5.0 1.0; 1.0 6.0]
B = [B₀,B₁,B₂]

C₀ = [1.0 3.0; 3.0 1.0]
C₁ = [2.0 0.0; 0.0 3.0]
C₂ = [4.0 5.0; 5.0 1.0]
C = [C₀,C₁,C₂]
```
"""
function minimize_geig_bracketing(A,B,C,λ₀,x₀)
   n = 2
   x = Variable(n)
   A = a[1] + x[1]*a[2] + x[2]*a[3]
   B = b[1] + x[1]*b[2] + x[2]*b[3]
   C = c[1] + x[1]*c[2] + x[2]*c[3]
   for λ = 5.0:-0.1:0.7
       constr1 = (λ*B-A ⪰ 0)
       constr2 = (B ⪰ 0)
       constr3 = (C ⪰ 0)
       problem = satisfy([constr1,constr2,constr3])
       solve!(problem,SCS.Optimizer(verbose=false))
       @show λ
       @show evaluate(x)
   end
end


text:
module Ecto

__precompile__(true)

include("exports.jl")
include("reflection.jl")
include("mirrortypes.jl")
include("schema.jl")
include("changeset.jl")
include("query.jl")
include("repo.jl")
include("english.jl")

end # module Ecto


text:
"""
    SimpleANN <: AbstractNeuralNetwork{Int,Int,1}

A simple, 1-layer, feed-forward artificial neural network with perceptron nodes.

    SimpleANN(ws::AbstractArray{Float64,2}, ts::AbstractVector{Float64})

Create a SimpleANN with weights `ws` and thresholds `ts`.

    SimpleANN(n::Int,m::Int) = SimpleANN(zeros(m,n), zeros(m))

Create a SimpleANN with `n` inputs and `m` outputs. All weights and thresholds are
initialized to zero.
"""
struct SimpleANN <: AbstractNeuralNetwork{Int,Int,1}
    insize :: Int
    outsize :: Int
    weights :: Array{Float64,2}
    thresholds :: Vector{Float64}
    function SimpleANN(ws::AbstractArray{Float64,2}, ts::AbstractVector{Float64})
        if any(size(ws).==0)
            throw(ArgumentError("weights matrix; invalid size"))
        elseif length(ts) != size(ws,1)
            throw(DimensionMismatch("weights and thresholds"))
        end
        m, n = size(ws)
        new(n, m, ws, ts)
    end
end

SimpleANN(n::Int,m::Int) = SimpleANN(zeros(m,n), zeros(m))

insize(net::SimpleANN) = (net.insize,)
outsize(net::SimpleANN) = (net.outsize,)

function fire!(net::SimpleANN, input::AbstractVector{Int},
               output::AbstractVector{Int})
    if net.insize != length(input)
        throw(DimensionMismatch("invalid input"))
    elseif net.outsize != length(output)
        throw(DimensionMismatch("invalid output"))
    end
    temp = net.weights * input - net.thresholds
    for i in 1:net.outsize
        output[i] = (temp[i] > 0.0)
    end
    output
end


text:
module Syntax

using ..OmegaCore.Util: mapf
using ..OmegaCore.Var: pw, liftapply
import ..OmegaCore: AbstractVariable
export @joint, @~, @uid

"Reduces to `@uid ~ expr`"
macro ~(ex)
  esc(:(~(@uid, ($ex))))
end

macro ~(args...)
  ids = args[1:end-1]
  ex = args[end]
  esc(:(~((@uid, $(ids...)), ($ex))))
end

"autotomatically generated id"
macro uid()
  rand(Int)
end

"""
Random variable over named tuple using varnames as keys

`@joint randvar1 randvar2 ...`

```julia
using OmegaCore
a = 1 ~ StdNormal()
b = 2 ~ StdNormal()
c = a >=ₚ b
randsample(@joint a b c)
=> (Ashoots = false, Bshoots = false)
```
"""
macro joint(args::Symbol...)
  esc(:(ω -> NamedTuple{$args}(OmegaCore.Util.mapf(ω, tuple($(args...))))))
end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
    aggstat(X::Union{AbstractMatrix, AbstractVector}; group, fun = mean)
Compute the mean (or other statistic) of each column of `X`, by group.
* `X` : A matrix (n, p) or vector (n).
* `group` : A vector (n) defining the groups.
* `fun` : Function to compute.
"""
function aggstat(X::Union{AbstractMatrix, AbstractVector}; group,
        fun = mean)
    X = ensure_mat(X)
    group = vec(group)
    q = size(X, 2)
    ztab = tab(group)
    lev = ztab.keys
    nlev = length(lev)
    ni = collect(values(ztab))
    zX = similar(X, nlev, q)
    for i in 1:nlev
        s = group .== lev[i]
        zX[i, :] .= vec(fun(vrow(X, s), dims = 1))
    end
    (X = zX, lev = lev, ni = ni)
end

"""
    aggstat(X::DataFrame; var_nam, group_nam, fun = mean)
Compute the mean (or other statistic) of selected columns of `X`, by group.
* `X` : A dataframe.
* `var_nam` : Names (vector) of the variables to summarize.
* `group_nam` : Names (vector) of the group variables to consider.
* `fun` : Function given the statistic to compute.

Variables defined in `var_nam` and `group_nam` must be columns of `X`.
"""
function aggstat(X::DataFrame; var_nam, group_nam, fun = mean)
    gdf = groupby(X, group_nam)
    combine(gdf, var_nam .=> fun, renamecols = false)
end

"""
    center(X, v)
Center each column of `X`.
* `X` : Data.
* `v` : Centering factors.
"""
function center(X, v)
    M = copy(X)
    center!(M, v)
    M
end

function center!(X, v)
    p = size(X, 2)
    @inbounds for j = 1:p
        X[:, j] .= vcol(X, j) .- v[j]
    end
end

"""
    checkdupl(X; digits = 3)
Finding replicated rows in a dataset.
* `X` : A dataset.
* `digits` : Nb. digits used to round `X`.

## Examples
```julia
X = rand(5, 3)
Z = vcat(X, X[1:3, :])
checkdupl(X)
checkdupl(Z)

M = hcat(X, fill(missing, 5))
Z = vcat(M, M[1:3, :])
checkdupl(M)
checkdupl(Z)
```
"""
function checkdupl(X; digits = 3)
    X = round.(ensure_mat(X), digits = digits)
    n = size(X, 1)
    rownum1 = []
    rownum2 = []
    @inbounds for i = 1:n
        @inbounds for j = 1:n
            res = isequal(vrow(X, i), vrow(X, j))
            if res
                push!(rownum1, i)
                push!(rownum2, j)
            end
        end
    end
    u = findall(rownum1 .!= rownum2)
    res = DataFrame((rownum1 = rownum1[u], rownum2 = rownum2[u]))
    k = length(u)
    if k > 0
        res = res[1:Int64(k / 2), :]
    end
    res
end

"""
    checkmiss(X)
Finding rows with missing data in a dataset.
* `X` : A dataset.

## Examples
```julia
X = rand(5, 4)
zX = hcat(rand(2, 3), fill(missing, 2))
Z = vcat(X, zX)
checkmiss(X)
checkmiss(Z)
```
"""
function checkmiss(X)
    X = ensure_mat(X)
    z = vec(sum(ismissing.(X); dims = 2))
    u = findall(z .> 0)
    DataFrame((rownum = u, nmissing = z[u]))
end

"""
    colmeans(X)
    colmeans(X, w)
Compute the mean of each column of `X`.
* `X` : Data.
* `w` : Weights of the observations.

Return a vector.

For a true weighted mean, `w` must preliminary be normalized to sum to 1.
"""
colmeans(X) = vec(Statistics.mean(X; dims = 1))
colmeans(X, w) = vec(w' * ensure_mat(X))

"""
    colvars(X)
    colvars(X, w)
Compute the (uncorrected) variance of each column of `X`.
* `X` : Data.
* `w` : Weights of the observations.

Return a vector.

**Note:** For a true weighted variance, `w` must preliminary be normalized to sum to 1.
"""
colvars(X) = vec(Statistics.var(X; corrected = false, dims = 1))
function colvars(X, w)
    p = size(X, 2)
    z = colmeans(X, w)
    @inbounds for j = 1:p
        z[j] = dot(view(w, :), (vcol(X, j) .- z[j]).^2)
    end
    z
end
colstds(X) = sqrt.(colvars(X))
colstds(X, w) = sqrt.(colvars(X, w))

"""
    dummy(y)
Examples
≡≡≡≡≡≡≡≡≡≡
y = ["d", "a", "b", "c", "b", "c"]
#y =  rand(1:3, 7)
dummy(y)
"""
function dummy(y)
    ztab = tab(y)
    lev = ztab.keys
    nlev = length(lev)
    ni = collect(values(ztab))
    Y = BitArray(undef, length(y), nlev)
    for i = 1:nlev
        Y[:, i] = y .== lev[i]
    end
    Y = Float64.(Y)
    (Y = Y, lev = lev, ni = ni)
end

function dummy2(y)
    z = tab(y)
    lev = z.keys
    nlev = length(lev)
    ni = collect(values(tab(y)))
    res = list(nlev, BitVector)
    for i = 1:nlev
        res[i] = y .== lev[i]
    end
    Y = reduce(hcat, res)
    (Y = Y, lev = lev, ni = ni)
end

"""
    ensure_df(X)
Reshape `X` to a dataframe if necessary.
"""
ensure_df(X::DataFrame) = X
ensure_df(X::AbstractVector) = DataFrame([X], :auto)
ensure_df(X::AbstractMatrix) = DataFrame(X, :auto)

"""
    ensure_mat(X)
Reshape `X` to a matrix if necessary.
"""
ensure_mat(X::AbstractMatrix) = X
ensure_mat(X::AbstractVector) = reshape(X, :, 1)
ensure_mat(X::Number) = reshape([X], 1, 1)
ensure_mat(X::LinearAlgebra.Adjoint) = Matrix(X)
ensure_mat(X::DataFrame) = Matrix(X)

"""
    findmax_cla(x, weights = nothing)
Find the most occurent level in `x`.
"""
function findmax_cla(x, weights = nothing)
    isnothing(weights) ? weights = ones(length(x)) : nothing
    res = aggstat(weights; group = x, fun = sum)
    res.lev[argmax(res.res)]   # if equal, argmax takes the first
end

"""
    iqr(x)
Compute the interquartile interval (IQR).
"""
iqr(x) = quantile(x, .75) - quantile(x, .25)

"""
    list(n::Integer)
Create a Vector{Any}(nothing, n).

`isnothing(object, i)` can be used to check if cell i is empty.

## Examples
```julia
list(5)
```
"""
list(n::Integer) = Vector{Any}(nothing, n)

"""
    list(n::Integer, type)
Create a Vector{type}(undef, n).

`isassigned(object, i)` can be used to check if cell i is empty.

## Examples
```julia
list(5, Float64)
list(5, Array{Float64})
list(5, Matrix{Float64})
```
"""
list(n::Integer, type) = Vector{type}(undef, n)

"""
    mad(x)
Compute the median absolute deviation (MAD),
adjusted by a factor (1.4826) for asymptotically normal consistency.
"""
mad(x) = 1.4826 * median(abs.(x .- median(x)))

"""
    matcov(X)
    matcov(X, w)
Compute the (uncorrected) covariance matrix of `X`.
* `X` : Data.
* `w` : Weights of the observations.

Uncorrected covariance matrix of the columns of `X`.

**Note:** For true weighted covariances, `w` must preliminary be normalized to sum to 1.
"""
matcov(X) = Statistics.cov(ensure_mat(X); corrected = false)

function matcov(X, w)
    X = ensure_mat(X)
    xmeans = colmeans(X, w)
    X = center(X, xmeans)
    z = Diagonal(sqrt.(w)) * X
    z' * z
end

"""
    mweights(w)
Return a vector of weights that sums to 1.
"""
mweights(w) = w / sum(w)

"""
    pnames(x)
Return the names of the elements of `x`.
"""
pnames(x) = propertynames(x)

"""
    psize(x)
Return the type and size of `x`.
"""
function  psize(x)
    println(typeof(x))
    println(size(x))
end

"""
    recodcat2num(x; start = 1)
Recode a categorical variable to a numeric variable
* `x` : Variable to recode.
* `start` : Numeric value that will be set to the first category.

The codes correspond to the sorted categories.

## Examples
```julia
x = ["b", "a", "b"]
zx = recodcat2num(x)
[x zx]
recodcat2num(x; start = 0)
recodcat2num([25, 1, 25])
```
"""
function recodcat2num(x; start = 1)
    z = dummy(x).Y
    ncla = size(z, 2)
    u = z .* collect(start:(start + ncla - 1))'
    u = sum(u; dims = 2)  ;
    u = vec(u)
end

"""
    recodnum2cla(x, q)
Recode a continuous variable to classes
* `x` : Variable to recode.
* `q` : Values separating the classes.
## Examples
```julia
x = [collect(1:10); 8.1 ; 3.1]
q = [3; 8] ;
zx = recodnum2cla(x, q)
[x zx]
probs = [.33; .66]
q = Statistics.quantile(x, probs)
zx = recodnum2cla(x, q) ;
[x zx]
```
"""
function recodnum2cla(x, q)
    zx = similar(x)
    q = sort(q)
    for i = 1:length(x)
        k = 1
        for j = 1:length(q)
            x[i] > q[j] ? k = k + 1 : nothing
        end
        zx[i] = k
    end
    zx
end

"""
    replacebylev(x, lev)
Replaces the elements of x
by levels of corresponding order.

* `x` : Vector of values to replace.
* `lev` : Vector containing the levels.

The number of levels in `x` and `lev` must be the same.
Before replacement, the levels in `x` and `lev` are sorted.

## Examples
```julia
x = [10, 4, 3, 3, 4, 4]
lev = ["B"; "C"; "AA"]
[x replacebylev(x, lev)]
[string.(x) replacebylev(string.(x), lev)]

x = [10, 4, 3, 3, 4, 4]
lev = [3; 0; -1]
[x replacebylev(x, lev)]
```
"""
function replacebylev(x, lev)
    m = length(x)
    lev = sort(lev)
    nlev = length(lev)
    x_lev = tab(x).keys
    v = similar(lev, m)
    @inbounds for i = 1:nlev
        u = findall(x .== x_lev[i])
        v[u] .= lev[i]
    end
    v
end

"""
    rmcol(X, s)
Remove the columns of `X` having indexes `s`.
## Examples
```julia
X = rand(5, 3) ;
rmcol(X, 1:2)
rmcol(X, [1, 3])
```
"""
rmcol(X::AbstractMatrix, s::Union{Vector, Number}) = X[:, setdiff(1:end, Int64.(s))]
rmcol(x::Vector, s::Union{Vector, Number}) = x[setdiff(1:end, Int64.(s))]

"""
    rmrow(X, s)
Remove the rows of `X` having indexes `s`.
## Examples
```julia
X = rand(5, 2) ;
rmrow(X, 2:3)
rmrow(X, [1, 4])
```
"""
rmrow(X::AbstractMatrix, s::Union{Vector, Number}) = X[setdiff(1:end, Int64.(s)), :]
rmrow(x::Vector, s::Union{Vector, Number}) = x[setdiff(1:end, Int64.(s))]

"""
    scale(X, v)
Scale each column of `X`.
* `X` : Data.
* `v` : Scaling factors.
"""
function scale(X, v)
    M = copy(X)
    scale!(M, v)
    M
end

function scale!(X, v)
    p = size(X, 2)
    @inbounds for j = 1:p
        X[:, j] .= vcol(X, j) ./ v[j]
    end
end

# Below: Much slower and requires more memories
scale2(X, v) = mapslices(function f(x) ; x ./ v ; end, X, dims = 2)

"""
    sourcedir(path)
Include all the files contained in a directory.
"""
function sourcedir(path)
    z = readdir(path)  ## List of files in path
    n = length(z)
    for i in 1:n
        include(string(path, "/", z[i]))
    end
end

"""
    summ(X; digits = 3)
Summarize a dataset (or a variable).
"""
function summ(X; digits = 3)
    X = ensure_df(X)
    res = DataFrames.describe(X, :mean, :min, :max, :nmissing) ;
    insertcols!(res, 5, :n => size(X, 1) .- res.nmissing)
    for j = 2:4
        z = vcol(res, j)
        s = findall(isa.(z, Float64))
        res[s, j] .= round.(res[s, j], digits = digits)
        end
    (res = res, ntot = size(X, 1))
end

"""
    summ(X, group; digits = 1)
Summarize a dataset (or a variable), by group.
* `X` : Dataset (n, p) or (n,).
* `group` : A vector (n,) defing the groups.
"""
function summ(X, group; digits = 1)
    zgroup = sort(unique(group))
    for i = 1:length(zgroup)
        u = findall(group .== zgroup[i])
        z = X[u, :]
        res = summ(z).res
        println("Group: ", zgroup[i])
        println(res)
        println("") ; println("")
        #println(repeat("-", 70))
    end
end

"""
    tab(x)
Univariate tabulation.
* `x` : Univariate class membership.

In the output, the levels in `x` are sorted.
Levels and values can be get by `tab(x).keys` and `tab(x).vals`.
"""
function tab(x)
    x = vec(x)
    sort(StatsBase.countmap(x))
end

"""
    tabnum(x)
Univariate tabulation (only integer classes).
* `x` : Univariate class membership.

In the output, the levels in `x` are sorted.
"""
function tabnum(x)
    x = vec(x)
    lev = sort(unique(x))
    cnt = StatsBase.counts(x)
    cnt = cnt[cnt .> 0]
    (cnt = cnt, lev = lev)
end

"""
    vrow(X::Matrix, j)
    vrow(x::Vector, i)
    vcol(X::Matrix, j)
View of the i-th row(s) or j-th column(s) of a matrix `X`,
or of the i-th element(s) of vector `x`.
"""
vrow(X, i) = view(X, i, :)
vrow(x::Vector, i) = view(x, i)
vcol(X, j) = view(X, :, j)
vcol(x::Vector, i) = view(x, i)


text:
module ExtensibleMCMCPlots

    using ExtensibleMCMC
    const eMCMC = ExtensibleMCMC
    using RecipesBase, Plots
    include("online_chains.jl")
    include("callbacks.jl")
    export PlottingCallback

end # module


text:
struct Ellipse
    x::Number
    y::Number

    rh::Number
    rv::Number
end

function circumference(ellipse::Ellipse)
    return pi * (3 * (ellipse.rh + ellipse.rv) - sqrt(10 * ellipse.rh * ellipse.rv + 3 * (ellipse.rh^2 + ellipse.rv^2)))
end

Base.:(==)(lhs::Ellipse, rhs::Ellipse) = lhs.x == rhs.x && lhs.y == rhs.y && lhs.rv == rhs.rv && lhs.rh == rhs.rh

ellipse_t(theta, e::Ellipse) = atan(e.rh / e.rv * tan(theta))

ellipse_f(t, e::Ellipse) = sqrt(e.rh^2 * sin(t)^2 + e.rv^2 * cos(t)^2)

# Calculates step size based on estimated arc length at current section compared to the 'optimum'.
# https://math.stackexchange.com/questions/433094/how-to-determine-the-arc-length-of-ellipse
function pointsOnEllipse(ellipse::Ellipse; filled::Bool=false, checks::Int=max(180, ceil(Int, circumference(ellipse) / 2)))
    res = Set{Tuple{Number, Number}}()
    step = 2 * pi / checks

    if checks > 0
        theta = -pi / 2
        while theta <= pi / 2
            ox = cos(theta) * ellipse.rh
            oy = sin(theta) * ellipse.rv

            push!(res, (round(Int, ellipse.x + ox), round(Int, ellipse.y + oy)))
            push!(res, (round(Int, ellipse.x - ox), round(Int, ellipse.y + oy)))

            if filled
                for i in floor(Int, -ox + 1):ceil(Int, ox - 1)
                    push!(res, (round(Int, ellipse.x + i), round(Int, ellipse.y + oy)))
                end
            end

            t = ellipse_t(theta, ellipse)
            f = ellipse_f(t, ellipse)
            mr = step * f
            theta += min(step / mr, step)
        end
    end

    return collect(res)
end


text:
__precompile__()

module Napier

export loge, logn, ln

function loge(n::Real, li::Real, ls::Real)
  if abs(li-ls) <= 0.000001
    (li+ls)/2.0
  elseif (exp(li)-n)*(exp((li+ls)/2.0)-n) < 0
    loge(n,li,(li+ls)/2.0)
  else
    loge(n,(li+ls)/2.0,ls)
  end
end

function logn(n::Real)
  if n <= 0
    throw(DomainError())
  elseif n < 1
    loge(n, 0, -n-80)
  elseif n == 1
    0
  else
    loge(n, 0, n)
  end
end

const ln = logn

end # module


text:
#
# rubberindex.jl -
#
# Rubber-index for Julia arrays.
#

"""
    colons(n)

yields a `n`-tuple of colons `:` (a.k.a. `Colon()`).

When `n` is known at compile time, it is faster to call:

    colons(Val(n))

This method is suitable to extract sub-arrays of build views when some kind of
rubber index is needed.  For instance:

    slice(A::AbstractArray{T,N}, i::Integer) where {T,N} =
        A[colons(Val(N-1))..., i]

defines a function that returns the `i`-th slice of `A` assuming index `i`
refers the last index of `A`.  Using the rubber-index `..`, a shorter
definition is:

    slice(A::AbstractArray, i) = A[.., i]

which is also able to deal with multiple trailing indices if `i` is a
`CartesianIndex`.

See also: `..`, [`RubberIndex`](@ref).

"""
colons(n::Integer) =
    (n ==  0 ? () :
     n ==  1 ? (:,) :
     n ==  2 ? (:,:,) :
     n ==  3 ? (:,:,:,) :
     n ==  4 ? (:,:,:,:,) :
     n ==  5 ? (:,:,:,:,:,) :
     n ==  6 ? (:,:,:,:,:,:,) :
     n ==  7 ? (:,:,:,:,:,:,:,) :
     n ==  8 ? (:,:,:,:,:,:,:,:,) :
     n ==  9 ? (:,:,:,:,:,:,:,:,:,) :
     n == 10 ? (:,:,:,:,:,:,:,:,:,:,) :
     _colons(n))

function _colons(n::Integer)
    n ≥ 0 || bad_ndims(n)
    return ([Colon() for i in 1:n]...,)
end

@noinline bad_ndims(n::Integer) =
    throw(ArgumentError(string("number of dimensions should be ≥ 0, got ", n)))

colons(v::Val{N}) where {N} = ntuple(colon, v)

# Just yields a colon whatever the argument.
colon(x) = Colon()

"""

`RubberIndex` is the singleron type that represents any number of indices.  The
constant `..` is defined as `RubberIndex()` and can be used in array indexation
to left and/or right justify the other indices.  For instance, assuming `A` is
a `3×4×5×6` array, then all the following equalities hold:

    A[..]           == A[:,:,:,:]
    A[..,3]         == A[:,:,:,3]
    A[2,..]         == A[2,:,:,:]
    A[..,2:4,5]     == A[:,:,2:4,5]
    A[2:3,..,1,2:4] == A[2:3,:,1,2:4]

As you can see, the advantage of the rubber index `..` is that it automatically
expands as the number of colons needed to have the correct number of indices.
The expressions are also more readable.

The rubber index may also be used for setting values.  For instance:

    A[..] .= 1         # to fill A with ones
    A[..,3] = A[..,2]  # to copy A[:,:,:,2] in A[:,:,:,3]
    A[..,3] .= A[..,2] # idem but faster
    A[2,..] = A[3,..]  # to copy A[3,:,:,:] in A[2,:,:,:]
    A[..,2:4,5] .= 7   # to set all elements in A[:,:,2:4,5] to 7

Leading/trailing indices may be specified as Cartesian indices (of type
`CartesianIndex`).

!!! warning
    There are two known limitations:
    1. The `end` reserved word can only be used in intervals specified *before*
       the rubber index but not *after*.  This limitation is due to the Julia
       parser cannot be avoided.
    2. At most 9 indices can be specified before the rubber index.  This
       can be extended by editing the source code.

See also: [`colons`](@ref).

""" RubberIndex

struct RubberIndex end

const .. = RubberIndex()
const … = .. # FIXME: should be deprecated

# Quickly get the tuple inside a Cartesian index.
to_tuple(I::CartesianIndex) = I.I

# Grow tuple of colons until tuple of axes is empty.
@inline growcolons(colons, inds::Tuple{}) = colons
@inline growcolons(colons, inds::Tuple) = growcolons((colons..., :), tail(inds))

# Drop as many colons as there are specified indices.
@inline dropcolons(colons, I::Tuple{}) = colons
@inline dropcolons(colons::Tuple{}, ::Tuple{}) = ()
@noinline dropcolons(colons::Tuple{}, I::Tuple) =
    throw(ArgumentError("too many indices specifed"))
@inline dropcolons(colons, I::Tuple) =
    dropcolons(tail(colons), tail(I))
@inline dropcolons(colons, I::Tuple{CartesianIndex, Vararg}) =
    dropcolons(dropcolons(colons, to_tuple(I[1])), tail(I))
@noinline dropcolons(colons, I::Tuple{RubberIndex, Vararg}) =
    throw(ArgumentError("more than one rubber index specified"))

@inline function to_indices(A, inds, I::Tuple{RubberIndex, Vararg})
    # Align the remaining indices to the tail of the `inds`.  First
    # `growcolons` is called to build a tuple of as many as colons as the
    # number of remaining axes. Second, `dropcolons` is used to drop as
    # many colons as the number of specified indices (taking into account
    # Cartesian indices).
    colons = dropcolons(growcolons((), inds), tail(I))
    to_indices(A, inds, (colons..., tail(I)...))
end

# avoid copying if indexing with .. alone, see
# https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/214
@inline Base.getindex(A::AbstractArray, ::RubberIndex) = A

# The following is needed to allow for statements like `A[..] .+= expr` to
# work properly.
dotview(A::AbstractArray{T,N}, ::RubberIndex) where {T,N} =
    dotview(A, ntuple(colon, Val{N}())...)


text:
using RLCore
using Test
using Random

mutable struct TestAgent <: AbstractAgent
    cur_s::Array{Float64, 1}
end
TestAgent(size) = TestAgent(zeros(size))

function RLCore.start!(agent::TestAgent, s)
    agent.cur_s .= s
    a = rand(1:4)
end

function RLCore.start!(agent::TestAgent, s, rng::Random.AbstractRNG)
    agent.cur_s .= s
    a = rand(rng, 1:4)
end

function RLCore.step!(agent::TestAgent, s, r, t)
    agent.cur_s .= s
    a = rand(1:4)
end

function RLCore.step!(agent::TestAgent, s, r, t, rng::Random.AbstractRNG)
    agent.cur_s .= s
    a = rand(rng, 1:4)
end

function test_agent()

    @testset "Test Agent GlobalRNG" begin
        agent = TestAgent(2)
        Random.seed!(1)
        s = rand(2)
        @test start!(agent, s) == 4
        @test all(agent.cur_s .== s)
        s = rand(2)
        @test step!(agent, s, 0, false) == 3
        @test all(agent.cur_s .== s)
    end

    @testset "Test Agent local RNG" begin
        agent = TestAgent(2)
        rng = MersenneTwister(1)
        s = rand(rng, 2)
        @test start!(agent, s, rng) == 4
        @test all(agent.cur_s .== s)
        s = rand(rng, 2)
        @test step!(agent, s, 0, false, rng) == 3
        @test all(agent.cur_s .== s)
    end

end


text:
@testset "Induction Motor" begin
    #valid (non-default) A and B
    im = SingleCageInductionMachine(
        name = "init",
        R_s = 0.0,
        R_r = 0.018,
        X_ls = 0.1,
        X_lr = 0.18,
        X_m = 3.2,
        H = 0.5,
        A = 0.2,
        B = 0.0,
        base_power = 100.0,
    )
    @test im isa PowerSystems.Component

    #invalid A, B or C
    @test_throws ErrorException SingleCageInductionMachine(
        name = "init",
        R_s = 0.0,
        R_r = 0.018,
        X_ls = 0.1,
        X_lr = 0.18,
        X_m = 3.2,
        H = 0.5,
        A = 0.2,
        B = 0.9,
        base_power = 100.0,
    )

    sys = System(100)
    bus = Bus(nothing)
    add_component!(sys, bus)
    static_load = PowerLoad(nothing)
    add_component!(sys, static_load)
    add_component!(sys, im, static_load)
    IMs = collect(get_components(SingleCageInductionMachine, sys))
    @test length(IMs) == 1

    im = SimplifiedSingleCageInductionMachine(
        name = "init",
        R_s = 0.0,
        R_r = 0.018,
        X_ls = 0.1,
        X_lr = 0.18,
        X_m = 3.2,
        H = 0.5,
        A = 0.2,
        B = 0.0,
        base_power = 100.0,
    )
    @test im isa PowerSystems.Component

    #invalid A, B or C
    @test_throws ErrorException SimplifiedSingleCageInductionMachine(
        name = "init",
        R_s = 0.0,
        R_r = 0.018,
        X_ls = 0.1,
        X_lr = 0.18,
        X_m = 3.2,
        H = 0.5,
        A = 0.2,
        B = 0.9,
        base_power = 100.0,
    )

    sys = System(100)
    bus = Bus(nothing)
    add_component!(sys, bus)
    static_load = PowerLoad(nothing)
    add_component!(sys, static_load)
    add_component!(sys, im, static_load)
    IMs = collect(get_components(SimplifiedSingleCageInductionMachine, sys))
    @test length(IMs) == 1
end


text:
###Importing dependencies
using Pkg
Pkg.activate("../")
using DataFrames, MClassification
include("../src/io.jl")

###Setting path of the datasets and getting the their names
path = "../datasets/sinthetic/"
datasets = split(read(`ls $path`, String), "\n")[1:end-1]

dataset_results = DataFrame()
for dataset in datasets
    n_instances = 0
    hits = 0
    misses = 0
    dataset_file = open(path * dataset , "r")

    classifier = MClassification.fit(dataset_file, 150, 0.1)

    while true
        sample = load_sample(dataset_file)
        if sample == ""
            break
        end
        n_instances += 1

        new_instance = sample[1:end-1]
        real_label = sample[end]

        predicted_label = MClassification.predict(classifier, new_instance)

        if real_label == predicted_label
            hits += 1
        else
            misses += 1
        end
    end

    x = hits / n_instances

    dataset_result = DataFrame(dataset =  dataset, n_initial_instances = 150 ,n_samples = n_instances, accuracy = round(x * 100, digits=2, RoundDown))
    append!(dataset_results, dataset_result)
end
println(dataset_results)


text:
#= from Nemo.jl
function show(io::IO, x::arb)
  d = ceil(parent(x).prec * 0.30102999566398119521)
  cstr = ccall((:arb_get_str, :libarb), Ptr{UInt8}, (Ref{arb}, Int, UInt),
                                                  x, Int(d), UInt(0))
  print(io, unsafe_string(cstr))
  ccall((:flint_free, :libflint), Void, (Ptr{UInt8},), cstr)
end
=#

function ArbFloat{P}(s::AbstractString) where {P}
    return ArbFloat{P}(BigFloat(s))
end

function ArbReal{P}(s::AbstractString) where {P}
    return ArbReal{P}(BigFloat(s))
end

function ArbFloat(s::AbstractString)
    P = DEFAULT_PRECISION[]
    return ArbFloat{P}(s)
end

function ArbReal(s::AbstractString)
    P = DEFAULT_PRECISION[]
    return ArbReal{P}(s)
end

function ArbComplex(s::AbstractString)
    P = DEFAULT_PRECISION[]
    x = parse(Complex{BigFloat}, s)
    return ArbComplex{P}(x)
end

function ArbComplex{P}(s::AbstractString) where {P}
    x = parse(Complex{BigFloat}, s)
    return ArbComplex{P}(x)
end

function Base.tryparse(::Type{ArbFloat}, s::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(BigFloat, s, base=base)
    P = DEFAULT_PRECISION[]
    return ifelse(isnothing(x), nothing, ArbFloat{P}(x))
end

function Base.tryparse(::Type{ArbFloat{P}}, s::S; base::Int=10) where {P, S<:AbstractString}
    x = tryparse(BigFloat, s, base=base)
    return ifelse(isnothing(x), nothing, ArbFloat{P}(x))
end

function Base.tryparse(::Type{ArbReal}, s::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(BigFloat, s, base=base)
    P = DEFAULT_PRECISION[]
    return ifelse(isnothing(x), nothing, ArbReal{P}(x))
end

function Base.tryparse(::Type{ArbReal{P}}, s::S; base::Int=10) where {P, S<:AbstractString}
    x = tryparse(BigFloat, s, base=base)
    return ifelse(isnothing(x), nothing, ArbReal{P}(x))
end

function Base.tryparse(::Type{ArbComplex}, s::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(Complex{BigFloat}, s, base=base)
    P = DEFAULT_PRECISION[]
    return ifelse(isnothing(x), nothing, ArbComplex{P}(x))
end

function Base.tryparse(::Type{ArbComplex{P}}, s::S; base::Int=10) where {P, S<:AbstractString}
    x = tryparse(Complex{BigFloat}, s, base=base)
    return ifelse(isnothing(x), nothing, ArbComplex{P}(x))
end


text:
# Advent of code 2021, day 5

function multicrossings(fname, countdiagonals)
    grid = zeros(Int, 1000, 1000)
    for txtsegment in [split(s, r"[\D]+") for s in strip.(readlines(fname)) if !isempty(s)]
        x1, y1, x2, y2 = map(s -> parse(Int, s), txtsegment)
        xsign, ysign = sign(x2 - x1), sign(y2 - y1)
        if xsign * ysign == 0 || countdiagonals
            for p in [
                [x1 + xsign * i, y1 + ysign * i] for
                i in 0:abs(xsign == 0 ? y1 - y2 : x1 - x2)
            ]
                grid[first(p), last(p)] += 1
            end
        end
    end
    return count(grid .> 1)
end

println("Part 1: ", multicrossings("AoCdata/AoC_2021_day5.txt", false))
println("Part 2: ", multicrossings("AoCdata/AoC_2021_day5.txt", true))


text:
using Dates
using TimeSeries
using WaterSystems

DayAhead  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime("1/1/2024  23:00:00", "d/m/y  H:M:S"))

junctions = [   Junction(name="J2",elevation=20, head = 20.0, minimum_pressure=1, coordinates = (lat = 0.0,lon = 0.0)),
                Junction(name="J3",elevation=18, head = 20.0, minimum_pressure=1, coordinates = (lat = 0.0,lon = 0.0)),
                Junction(name="J4",elevation=10, head = 20.0, minimum_pressure=1, coordinates = (lat = 0.0,lon = 0.0)),
                Junction(name="T1",elevation=320,head = 20.0, minimum_pressure=1, coordinates = (lat = 0.0,lon = 0.0)),
                Junction(name="R1",elevation=0,  head = 50.0, minimum_pressure=1, coordinates = (lat = 0.0,lon = 0.0))
            ];

tanks = [ RoundTank("T1", junctions[4], (min = 0.0, max = 2500.0), 10.0, 500.0, pi*25, 200/(10*pi), (min=0.0, max=1000/(10*pi)))];

reservoirs = [SourceReservoir("R1",junctions[1],0)];
nodes = vcat(junctions,tanks,reservoirs)

pipes = [   ReversibleFlowPipe("P1", (from=junctions[4],to=junctions[3]), 0.6, 100.0, 200.0, 2.0, 10.0, 1),
            StandardPositiveFlowPipe("P2",(from=junctions[2],to=junctions[2]),0.6,150.0,200.0,5.0,30.0,1)
        ];

valves = [  PressureReducingValve("V1", (from=junctions[1],to=junctions[2]), "open", 1.0, 1.0)
        ];

elec_price = rand(24)

pumps = [ConstSpeedPump("PMP1",(from=junctions[5],to=junctions[4]), 1,
                        [(0.0,350.0),(400.0,328.0),(800.0,295.0),(1400.0,260.0),(1600.0,245.0),(1800.0,225.0),(2200.0,153.0)],
                        0.85,TimeSeries.TimeArray(DayAhead, elec_price))]

demands = [ WaterDemand("D1", junctions[3], true, 10.0, TimeSeries.TimeArray(DayAhead, rand(24)), TimeSeries.TimeArray(DayAhead, rand(24)))];
#demand = TimeSeries.TimeArray(DayAhead, rand(length(DayAhead))), demandforcast = TimeSeries.TimeArray(DayAhead, rand(length(DayAhead)))

simulations = Simulation(24, 1.0, 24, DateTime(2013,7,10), DateTime(2013,7,11))
links = vcat(pipes,valves, pumps)
link_classes = LinkClasses(links, pipes, pumps, valves)
node_classes = NodeClasses(junctions, junctions[1:3], tanks, reservoirs)
net = Network(junctions,links)
system = WaterSystem(node_classes, link_classes, demands)


text:
using NeuroAnalysis,YAML,Plots,Interact

cd(@__DIR__)
## Generate hartley subspace
dk=0.2
kbegin=0.2
kend=6.6
# add 4 uniform gray as blanks to match 4 phases of each grating with the same mean luminance
nblank=4

hs = hartleysubspace(;kbegin,kend,dk,phase=0.0,shape=:circle,addhalfcycle=true,blank=(kx=0.0,ky=0.0,phase=0.375),nblank)
ss = map(i->cas2sin(i...),hs)

# check gratings
@manipulate for i in 1:length(ss)
    Gray.(grating(θ=ss[i].θ,sf=ss[i].f,phase=ss[i].phase,size=(5,5),ppd=30))
end

# save condition
hartleyconditions = Dict(:Ori=>map(i->rad2deg(mod(i.θ,2π)),ss),:SpatialFreq=>map(i->i.f,ss),:SpatialPhase=>map(i->i.phase,ss))
title = "Hartley_k[$kbegin,$kend]_dk[$dk]_nblank[$nblank]"
YAML.write_file("$title.yaml",hartleyconditions)

max_stars_repo_name: VLAB-ExperimentSystem/VLab

max_issues_repo_name: VLAB-ExperimentSystem/VLab

max_forks_repo_name: VLAB-ExperimentSystem/VLab


text:
#=
Fast Type-I Discrete Cosine Transform and Chebyshev polynomial interpolation
=#
module Fct

#fast type I DCT:
#(fct(x_0,...,x_N))_n == sum_{0<=m<=N}x_m*cos(pi*m*n/N)
#with 1/2 weight given to first and last summands.
#WARNING: assumes N is power of 2.
function fct(x::AbstractVector,T::Type=Complex{Float64})::AbstractVector{T}
    fct!(x,fct_alloc(x,T)...)
end
function fct_alloc(x::AbstractVector,T::Type=Complex{Float64})
    N = length(x)-1
    twoN = 2N
    Vector{T}(undef,twoN), Vector{T}(undef,twoN), Vector{T}(undef,twoN)
end
function fct!(dct_in::AbstractVector, #input
              ifft_in::AbstractVector, #workspace: ifft input
              ifft_out::AbstractVector, #output is view(ifft_out,1:N+1)
              ifft_work::AbstractVector, #workspace: ifft workspace
              N=length(dct_in)-1)
    ifft_in[1] = dct_in[1]
    i,j = 2,2N
    while i<=N
        z = dct_in[i]
        ifft_in[i] = z
        ifft_in[j] = z
        i+=1
        j-=1
    end
    ifft_in[N+1] = dct_in[N+1]
    ifft!(ifft_in,ifft_out,ifft_work,2N)
    dct_out = view(ifft_out,1:N+1)
    dct_out .*= 0.5
end

#inverse of fct and dct
function idct(x::AbstractVector,T::Type=Complex{Float64})::AbstractVector{T}
    idct!(x,fct_alloc(x,T)...)
end
function idct!(idct_in::AbstractVector, #input
               ifft_in::AbstractVector, #workspace: ifft input
               ifft_out::AbstractVector, #output is view(ifft_out,1:N+1)
               ifft_work::AbstractVector, #workspace: ifft workspace
               n=length(idct_in)-1)
    idct_out = fct!(idct_in,ifft_in,ifft_out,ifft_work,n)
    T = eltype(ifft_work)
    r = real(T)(2)/n
    idct_out .*= r
end

ifft_cache=Dict{Pair{Int64,Type},Vector}()

#fast inverse DFT without 1/N factor
#WARNING: assumes N is power of 2.
#(ifft(x_0,...,x_{N-1}))_n = sum_{0<=m<N}x_m*exp(2*pi*m*n/N)
function ifft(x::AbstractVector,T::Type=Complex{Float64})
    N=length(x)
    y=Vector{T}(undef,N) #output
    z=Vector{T}(undef,N) #work space
    ifft!(x,y,z)
end
function ifft!(input::AbstractVector,output::AbstractVector,workspace::AbstractVector,
               N=length(input))
    T=eltype(workspace)
    key = N=>T
    if key in keys(ifft_cache)
        exps = ifft_cache[key]
    else
        exps = ifft_cache[key] = exp_populate(N,T)
    end
    ifft_helper!(input,output,workspace,exps,0,1,N)
end

#x is input; y is output; z is work space
#assumes length(x)==length(y)==length(exps)<=length(z)
#input/output indices: 1+offset,1+offset+stride,...,1+offset+(len-1)*stride
function ifft_helper!(x::AbstractVector,
                      y::AbstractVector,
                      z::AbstractVector,
                      exps::AbstractVector,
                      offset::Integer,
                      stride::Integer,
                      len::Integer)
    if len==1
        y[1+offset]=x[1+offset]
        return y
    end
    doublestride=stride<<1
    halflen=len>>1
    ifft_helper!(x,y,z,exps,offset,doublestride,halflen)
    ifft_helper!(x,y,z,exps,offset+stride,doublestride,halflen)
    i=1
    j=1+offset
    k=1
    while k<=halflen
        z[k]=y[j]+y[j+stride]*exps[i]
        i+=stride
        j+=doublestride
        k+=1
    end
    j=1+offset
    while k<=len
        z[k]=y[j]+y[j+stride]*exps[i]
        i+=stride
        j+=doublestride
        k+=1
    end
    #println()
    i=1+offset
    k=1
    while k<=len
        y[i]=z[k]
        i+=stride
        k+=1
    end
    y
end

#output: exps[1+m]==exp(2im*pi*m/N) for 0<=m<N
#N can be any positive integer.
function exp_populate(N::Integer,T::Type=Complex{Float64})
    exps=Vector{T}(undef,N)
    w = 2im*T(pi)/N
    for m in 0:N-1
        exps[m+1] = exp(w*m)
    end
    exps
end

#=
Chebyshev polynomial interpolation of f on [-1,1] intended for Clenshaw-Curtis quadrature.
On input y_0,...,y_N, outputs c_0,...,c_N such that y_n == sum_m c_m cos(m*n*pi/N).
y_m should be f(cos(m*pi/N)).
WARNING: N must be a power of 2.
=#
function chebcoeffs(samples::AbstractVector,T::Type=promote_type(Complex{Bool},eltype(samples)))
    chebcoeffs!(samples,fct_alloc(samples,T)...)
end
function chebcoeffs!(samples::AbstractVector,
                     ifft_in::AbstractVector,
                     ifft_out::AbstractVector,#output is view(ifft_out,1:length(samples))
                     ifft_work::AbstractVector,
                     n=length(samples)-1)
    coeffs = idct!(samples,ifft_in,ifft_out,ifft_work,n)
    coeffs[1] *= 0.5
    coeffs[end] *= 0.5
    coeffs
end

function chebsample(f,n::Integer,base=0,stride=1;T::Type=Complex{Float64})
    samples = Vector{T}(undef,length(base:stride:n))
    chebsample!(f,samples,n,base,stride)
end
function chebsample!(f,samples,n::Integer=length(samples)-1,base=0,stride=1)
    p = real(eltype(samples))(pi)
    i=1
    for j in base:stride:n
        samples[i] = f(cos(p*j/n))
        i += 1
    end
    samples
end

function shiftedchebsample(center,radius,f1,f2,n::Integer,base=0,stride=1;T::Type=Complex{Float64})
    samples = Vector{T}(undef,length(base:stride:n))
    shiftedchebsample!(center,radius,f1,f2,samples,n,base,stride)
end
function shiftedchebsample!(center,radius,f1,f2,samples,n::Integer=length(samples)-1,base=0,stride=1)
    p = real(eltype(samples))(pi)
    i=1
    for j in base:stride:n
        x = center + radius*cos(p*j/n)
        samples[i] = f1(x)*f2(x)
        i += 1
    end
    samples
end

function tonechebsample(cfreq,center,radius,f1,f2,n::Integer,base=0,stride=1;T::Type=Complex{Float64})
    samples = Vector{T}(undef,length(base:stride:n))
    tonechebsample!(cfreq,center,radius,f1,f2,samples,n,base,stride)
end
function tonechebsample!(cfreq,center,radius,f1,f2,samples,n::Integer=length(samples)-1,base=0,stride=1)
    p = real(eltype(samples))(pi)
    i=1
    for j in base:stride:n
        x = cos(p*j/n)
        y = center + radius*x
        z = x * -cfreq
        s,c = sincos(z)
        samples[i] = f1(y)*f2(y)*complex(c,s)
        i += 1
    end
    samples
end

function chirpchebsample(cfreq,chirp,center,radius,f1,f2,n::Integer,base=0,stride=1;T::Type=Complex{Float64})
    samples = Vector{T}(undef,length(base:stride:n))
    chirpchebsample!(cfreq,chirp,center,radius,f1,f2,samples,n,base,stride)
end
function chirpchebsample!(cfreq,chirp,center,radius,f1,f2,samples,n::Integer=length(samples)-1,base=0,stride=1)
    p = real(eltype(samples))(pi)
    i=1
    for j in base:stride:n
        x = cos(p*j/n)
        y = center + radius*x
        z = -x*(cfreq + x*chirp)
        s,c = sincos(z)
        samples[i] = f1(y)*f2(y)*complex(c,s)
        i += 1
    end
    samples
end

#double Chebyshev degree
function doublesample(f,oldsamples)
    N = length(oldsamples)-1
    T = eltype(oldsamples)
    newsamples = chebsample(f,2N,1,2;T=T)
    samples = Vector{T}(undef,2N+1)
    samples[1] = oldsamples[1]
    for n in 1:N
        samples[2n] = newsamples[n]
        samples[2n+1] = oldsamples[n+1]
    end
    samples
end

function shifteddoublesample(center,radius,f1,f2,oldsamples)
    N = length(oldsamples)-1
    T = eltype(oldsamples)
    newsamples = shiftedchebsample(center,radius,f1,f2,2N,1,2;T=T)
    samples = Vector{T}(undef,2N+1)
    samples[1] = oldsamples[1]
    for n in 1:N
        samples[2n] = newsamples[n]
        samples[2n+1] = oldsamples[n+1]
    end
    samples
end

function chirpdoublesample(cfreq,chirp,center,radius,f1,f2,oldsamples)
    N = length(oldsamples)-1
    T = eltype(oldsamples)
    newsamples = chirpchebsample(cfreq,chirp,center,radius,f1,f2,2N,1,2;T=T)
    samples = Vector{T}(undef,2N+1)
    samples[1] = oldsamples[1]
    for n in 1:N
        samples[2n] = newsamples[n]
        samples[2n+1] = oldsamples[n+1]
    end
    samples
end

function doublesample!(f,samples,N)
    for n in N:-1:1
        samples[2n+1] = samples[n+1]
    end
    chebsample!(f,view(samples,2:2:2N),2N,1,2)
    samples
end

function shifteddoublesample!(center,radius,f1,f2,samples,N)
    for n in N:-1:1
        samples[2n+1] = samples[n+1]
    end
    shiftedchebsample!(center,radius,f1,f2,view(samples,2:2:2N),2N,1,2)
    samples
end

function tonedoublesample!(cfreq,center,radius,f1,f2,samples,N)
    for n in N:-1:1
        samples[2n+1] = samples[n+1]
    end
    tonechebsample!(cfreq,center,radius,f1,f2,view(samples,2:2:2N),2N,1,2)
    samples
end

end #module


text:
"""
    support(cvxpoly, dir)

Compute the point of contact between a convex object and its
supporting hyperplane defined by the given normal direction as
a StaticArray.
"""
function support(vertices::SMatrix, dir::AbstractVector)
    @inbounds vertices[:, argmax(vertices'*dir)]
end

# TODO support for ordered vertices / simplices in a convex hull
#  struct CvxHull{A<:SMatrix}
#      pts::A
#  end
#  function support(vertices::CvxHull{SMatrix{2, M}}, dir::SVector{2, T}) where {M, T}
#      li = 1; ri = M
#      left = vertices[:, 1]'*dir
#      right = vertices[:, M]'*dir
#      cnt = 10
#      while ri - li > 2
#          mi = div(li + ri, 2)
#          mid = vertices[:, mi]'*dir
#          if left < right
#              left = mid
#              li = mi
#          else
#              right = mid
#              ri = mi
#          end
#          cnt -= 1
#      end
#      if left < right
#          return vertices[:, ri]
#      else
#          return vertices[:, li]
#      end
#  end


text:
using DSP, Random
export AutoRegressive

"""
    AutoRegressive(n, method = LPCLevinson())

Autoregressive surrogates of order-`n`. The autoregressive coefficients `φ` are estimated
using `DSP.lpc(x, n, method)`, and thus see the documentation of DSP.jl for possible
`method`s.

While these surrogates are obviously suited to test the null hypothesis whether the data
are coming from a autoregressive process, the Fourier Transform-based surrogates are
probably a better option. The current method is more like an explicit way to
produce surrogates for the same hypothesis by fitting a model.
It can be used as convient way to estimate
autoregressive coefficients and automatically generate surrogates based on them.

The coefficients φ of the autoregressive fit can be found by doing
```julia
sg = surrogenerator(x, AutoRegressive(n))
φ = sg.init.φ
```
"""
struct AutoRegressive{M} <: Surrogate
    n::Int
    m::M
end
AutoRegressive(n::Int) = AutoRegressive(n, LPCLevinson())

function surrogenerator(x, method::AutoRegressive)
    φ, e = lpc(x, method.n, method.m)
    init = (d = Normal(0, std(x)), φ = φ)
    return SurrogateGenerator(method, x, init)
end

function (sg::SurrogateGenerator{<:AutoRegressive})()
    N, φ, d = length(sg.x), sg.init.φ, sg.init.d
    return autoregressive(d, N, φ)
end

function autoregressive(d::Normal{T}, N, φ) where {T}
    f = length(φ)
    z = zeros(T, N+f)
    @inbounds for i in 1:f; z[i] = rand(d); end
    @inbounds for i in f+1:length(z)
        s = zero(T)
        for j in 1:f
            s += φ[j]*z[i-j]
        end
        s += rand(d)
        z[i] = s
    end
    return view(z, f+1:N+f)
end


text:
using Test, Scanf

@testset "Scanf" begin

    # pointer
    @testset "%p" begin
        r, p = @scanf("0x42      ", "%4p", Ptr{Nothing})
        @test p == Ptr{Nothing}(UInt(0x42))
    end

    # whitespace is squeezed
    @testset "whitespace" begin
        f = Scanf.format"  "
        @test length(f.formats) == 1
    end

    # %% is maintained before WS
    @testset "%% before WS" begin
        f = Scanf.format"%% "
        @test length(f.formats) == 2
        @test f.formats[1] isa Scanf.LiteralSpec
        @test f.formats[2] isa Scanf.Spec{Val{Char(Scanf.SKIP)}}
    end

    # floating point numbers
    @testset "float to $T" for T in (Float64, Float32, Float16, BigFloat)
        f1 = Scanf.Format("%g")
        f2 = Scanf.Format("%4e")
        @testset "$res" for (inp, res) in (
            ("-13", big"-13.0"),
            ("99.76 ", big"99.76"),
            ("42.11f3", big"42.11e3"),
            ("-12.7491234797912347971e-4", big"-12.7491234797912347971e-4"),
            ("0x1.4711p2", 0x1.4711p2)
        )
            @test scanf(inp, f1, T) == (1, T(res))
            @test scanf(inp, f2, T) == (1, parse(T, inp[1:min(4, length(inp))]))
        end
    end

    # floating point nan and infinity
    @testset "float to $T" for T in (Float64, Float32, Float16, BigFloat)
        f = Scanf.format"%f%n"
        @testset "$inp" for (inp, res, rr, nn) in (
            ("InfX", T(Inf), 2, 3),
            ("-InFX", T(-Inf), 2, 4),
            ("infiNITYX", T(Inf), 2, 8),
            ("NANX", T(NaN), 2, 3),
            ("-nanX", T(-NaN), 2, 4),
            ("infiX", T(0.0), 0, 0),
            ("naX", T(0.0), 0, 0),
        )
            io = IOBuffer(inp)
            r, x, n = scanf(io, f, T(0.0), 0)
            @test r == rr
            @test n == nn
            @test x === res || T <: BigFloat && (x == res || isnan(x) && isnan(res))
            @test peek(io, Char) == 'X'
        end
    end

    @testset "marginal floats" begin
    end

    @testset "incomplete floats $inp" for inp in ("Z", "0xZ", ".Z", "e", "+e", "-e", "1eZ", "0E-Z")
        io = IOBuffer(inp)
        @test ( @scanf io "%f" 0.0 ) == (0, 0.0)
        @test peek(io, Char) == inp[end]
    end

    # integers
    @testset "integer %i to $T" for T in (UInt64, UInt32, UInt16)
        f1 = Scanf.Format("%i")
        @testset "$res" for res in ("17", "4711", "0x7ffe")
            @test scanf(res, f1, T) == (1, parse(T, res))
        end
    end
    sigreals = (Int64, Int32, Int16, BigInt, Float64, Float32, Float16, BigFloat)
    @testset "integer %i to $T" for T in sigreals
        f1 = Scanf.Format("%i")
        @testset "$res" for res in ("+17", "-4711", "0x7ffe", "-0x7ffe")
            @test scanf(res, f1, T) == (1, parse(T, res))
        end
    end
    @testset "integer %i octal to $T " for T in sigreals
        @testset "$res" for (inp, res) in (("077", 63), ("-077", -63))
            @test scanf(inp, Scanf.format"%i", T) == (1, res)
        end
    end

    @testset "incomplete integers $inp" for inp in ("Z", "0xZ", )
        io = IOBuffer(inp)
        @test ( @scanf io "%i" 0 ) == (0, 0)
        @test peek(io) == UInt8('Z')
    end
    @testset "incomplete pointers $inp $arg" for inp in ("Z", "0Z", "0XZ"), arg in (Ptr{String}, Ptr{String}(0))
        io = IOBuffer(inp)
        @test ( @scanf io "%p" arg) == (0, Ptr{String}(0))
        @test peek(io) == UInt8('Z')
    end

    @testset "integer %i octal input" for T in (Int64,)
        @test scanf("0377", Scanf.Format("%i"), T) == (1, 255)
    end
    @testset "integer %i negated input" for T in (UInt64, UInt16)
        @test scanf("-99", Scanf.Format("%i"), T) == (1, -T(99))
    end

    @testset "integer %o to $T" for T in (Int64, UInt64, Int32, UInt32, Int16, UInt16)
        f1 = Scanf.Format("%o")
        @testset "$res" for res in ("17", "4711", "0377")
            @test scanf(res, f1, T) == (1, parse(T, res, base = 8))
        end
    end

    @testset "integer %x to $T" for T in (Int64, UInt64, Int32, UInt32, Int16, UInt16)
        f1 = Scanf.Format("%x")
        @testset "$res" for res in ("0x4711",)
            @test scanf(res, f1, T) == (1, parse(T, res))
        end
    end

    @testset "convert integer to other type" begin
        f = Scanf.format"%i"
        @test scanf("  14", f, String) == (1, "14")
    end

    @testset "%i follow up" for (inp, rr, a, b) in [
        ("Z", 0, 0, ""),
        ("0", 1, 0, ""),
        ("0x", 0, 0, ""),
        ("0xZ", 0, 0, ""),
    ]
        f = Scanf.format"%i%s"
        r, x, y = scanf(inp, f, 42, String)
        @test r == rr
        @test r < 1 || a == x
    end

    # character sets
    @testset "character set [" begin
        f1 = Scanf.Format("%10[a-c ]")
        f2 = Scanf.Format("%10[abc ]")
        f3 = Scanf.Format("%10[abcx-y ]")
        res = " abbcbbcbadabbdbbabbbann"
        @test scanf(res, f1, String) == (1, res[1:10])
        @test scanf(res, f2, "") == (1, res[1:10])
        @test scanf(res, f3, "") == (1, res[1:10])
    end

    @testset "character set ^" begin
        f1 = Scanf.Format("%10[^A-A]")
        f2 = Scanf.Format("%10[^A-B]")
        f3 = Scanf.Format("%10[^A-BX]")
        res = " abbcb Abadabbdbbabbbann"
        @test scanf(res, f1, "") == (1, res[1:7])
        @test scanf(res, f2, String) == (1, res[1:7])
        @test scanf(res, f3, "") == (1, res[1:7])
    end

    @testset "string to Char" begin
        @test @scanf("äbc", "%s", Char[]) == (1, ['ä', 'b', 'c'])
        @test @scanf("äbc", "%[a-zäöü]", Char) == (1, 'ä')
        @test @scanf("Äbc", "%[a-zäöü]", Char[]) == (0, Char[])
        @test @scanf("Äbc", "%[a-zäöü]", Char) == (0, '\0')
    end

    @testset "many arguments" begin
        f = Scanf.Format("%d%d%d%d%d%d%d%d%d%d ")
        @test scanf("1 2 3 4 5 6 7 8 9 10", f, fill(Int, 10)...)[1] == 10
    end

    @testset "single characters" begin
        @test @scanf("a%bX", "a%%b%c", Char) == (1, 'X')
        @test @scanf("a%bX", "a%%b%c", String) == (1, "X")
    end

    @testset "multiple characters" begin
        v = Char[]
        r, cc = @scanf("abXYZ", "ab%3c", v)
        @test ['X', 'Y', 'Z'] == cc == v
        r, cc = @scanf("abXYZ", "ab%3c", Char[])
        @test ['X', 'Y', 'Z'] == cc
        r, cc = @scanf("abXYZ", "ab%3c", Char)
        @test 'X' == cc
        r, cc = @scanf("abXYZ", "ab%3c", String)
        @test "XYZ" == cc
    end

    # string
    @testset "strings" begin
        r, a, b = @scanf("Hällo\u1680heimør", "%s%s", String, String)
        @test r == 2
        @test "Hällo" == a
        @test "heimør" == b
    end

    # position
    @testset "%n" begin
        @test @scanf("aäc 15 16  \n", "aäc %*hhd %*Ld%n", 0) == (1, 10)
        @test @scanf("abc", "abd%n", Int) == (0, 0)
    end

    # basics
    @testset "basics" begin
        @test_throws ArgumentError try
            @eval @scanf(1, 2, 3)
        catch ex
            rethrow(ex.error)
        end
        @test @scanf("%15", "%%%d", Int) == (1, 15)
        @test_throws ArgumentError Scanf.Format("")
        @test_throws ArgumentError Scanf.Format("%")
        @test_throws ArgumentError Scanf.Format("%l")
        @test_throws ArgumentError Scanf.Format("%hh")
        @test_throws ArgumentError Scanf.Format("%hhU")
        @test_throws ArgumentError Scanf.Format("%+")
        @test_throws ArgumentError Scanf.Format("%.")
        @test_throws ArgumentError Scanf.Format("%.0")
        @test Scanf.Format("%%").formats[1].string == "%"
        @test_throws ArgumentError Scanf.Format("%y%d")
        @test_throws ArgumentError Scanf.Format("%\u00d0%d")
        @test_throws ArgumentError Scanf.Format("%\u0f00%d")
        @test_throws ArgumentError Scanf.Format("%\U0001ffff%d")
        @test_throws ArgumentError Scanf.Format("%[")
        @test_throws ArgumentError Scanf.Format("%[^")
        @test_throws ArgumentError Scanf.Format("%[]")
        @test_throws ArgumentError Scanf.Format("%[^]")
        @test Scanf.Format("%[A-A]").formats[1].set == "A"
        f = Scanf.Format("%d ")
        @test Scanf.Format(f) === f
        @test Scanf.format"%d " == f
        @test_throws ArgumentError scanf("", f)

        toobig_argnum = (Int64(typemax(Scanf.ARGNUM_TYPE)) + 1)
        @test Scanf.Format("%d"^(toobig_argnum - 1)) isa Scanf.Format
        @test_throws ArgumentError Scanf.Format("%d"^toobig_argnum)

        toobig_width = (Int64(typemax(Scanf.WIDTH_TYPE)) + 1)
        @test Scanf.Format("%$(toobig_width - 1)d") isa Scanf.Format
        @test_throws ArgumentError Scanf.Format("%$(toobig_width)d")
    end

    @testset "default arguments" begin
        def = (1, 2.0, "x", 'y')
        r, i, f, s, c = @scanf("42X", "%i%f%s%c", def...)
        @test r == 1
        @test (i, f, s, c) == (42, def[2:end]...)
    end

    # failing literal match consumes all matching characters
    @testset "literal string" begin
        io = IOBuffer("AÖÜdef")
        @test scanf(io, Scanf.format"AÖÜDEF")[1] == 0
        @test read(io, String) == "def"
    end

    @testset "%n" begin
        @test @scanf("1234", "%3d4%n", Int, 0) == (2, 123, 4)
        @test @scanf("😉", "%s%n", String, 0) == (2, "😉", 1)
        @test @scanf("1234 ", "%s%n", String, 0) == (2, "1234", 4)
    end

    @testset "show specifiers" begin
        f = Scanf.format"%dABC%[a-cA-C]%[^]a-cx] %[]B-Aa-zC-]%[B-A]"
        @test sprint(show, f.formats) ==
            "(%d, \"ABC\", %[a-cA-C], %[^]xa-c], %*_, %[]Ca-z-], %[1-0])"
        @test sprint(show, f) ==
            "Scanf.format\"%dABC%[a-cA-C]%[^]a-cx] %[]B-Aa-zC-]%[B-A]\""
    end

    @testset "scanf from stdin" begin
        f1 = () -> scanf(Scanf.format"abc%i", 0)
        f2 = () -> @scanf "abc%i" 0

        file = tempname(cleanup = true)
        write(file, "abc42")

        io = open(file)
        @test redirect_stdin(f1, io) == (1, 42)
        close(io)

        io = open(file)
        @test redirect_stdin(f2, io) == (1, 42)
        close(io)
    end

    @testset "overflow $T" for T in (Float64, Float32, Float16, BigFloat)
        @testset "$T $res" for (input, res) in [
            ("+1.55e9999999999999999999", Inf),
            ("-1.55e+9999999999999999999", -Inf),
            ("1.5e-9999999999999999999", 0.0),
            ("-1.5e-9999999999999999999", -0.0),
        ]
            @test ((r, a) = @scanf(input, "%f", oneunit(T)); r == 1) && isequal(a, T(res))
        end
    end
    @testset "overflow $T" for T in (Int128, Int64, Int32, Int16, Int8)
        str = "270141183460469231731687303715884105727"
        @testset "$T $input" for input in [str, "-" * str]
            @test @scanf(input, "%d", T) == (0, 0)
        end
    end
    @testset "overflow $T" for T in (UInt128, UInt64, UInt32, UInt16, UInt8)
        str = "470141183460469231731687303715884105727"
        @testset "$T $input" for input in [str, "-" * str]
            @test @scanf(input, "%d", T) == (0, 0)
        end
    end
    @testset "overflow Ptr" begin
        T = Ptr{String}
        input = "0x12341234abcdabcd0"
        @test @scanf(input, "%p", T) == (0, T(0))
    end

    @testset "convert to different type" begin
        inp = "-1.5e+3"
        f = Scanf.format"%f"
        @test scanf(inp, f, "") == (1, inp)
        @test scanf(inp, f, Int) == (1, -1500)
    end
    @testset "convert string to $T" for T in (Int, Float32)
        f = Scanf.format"%s"
        @test scanf("-1.5f+3", f, Float64) == (1, -1500.0)
        @test scanf("-1.5f+3", f, Int) == (0, 0)
        @test scanf("-077", f, Int) == (1, -77) # like %d
        @test scanf("-0x77", f, Int) == (1, -119) # like %x
    end

end # @testset "Scanf"


text:
# f_n = \left( 1 + \exp\left[ \epsilon_{n} - \mu \right]/kY \right)^{-1}
# kT = smearing width
function smear_FD( ev::Array{Float64,1}, efermi::Float64, kT::Float64; is_spinpol=false )
    x = (ev[:] .- efermi)/kT
    #
    Nstates = length(ev)
    f = zeros(Nstates)
    if is_spinpol
        for ist = 1:Nstates
            f[ist] = 1.0/( 1.0 + exp(x[ist]))
        end
    else
        for ist = 1:Nstates
            f[ist] = 2.0/( 1.0 + exp(x[ist]))
        end
    end
    return f
end

function smear_FD( ev::Float64, efermi::Float64, kT::Float64; is_spinpol=false )
    x = (ev - efermi)/kT
    if is_spinpol
        f = 1.0/( 1.0 + exp(x))
    else
        f = 2.0/( 1.0 + exp(x))
    end
    return f
end