text:
################# GtkReactive time widgets ###############
import GtkReactive: timewidget, InputWidget

immutable TimeWidget2{T <: Dates.TimeType} <: InputWidget{T}
    signal::Signal{T}
    widget::Gtk.GtkFrame
end

"""
    timewidget(time)

Return a time widget that includes the `Time` and a `GtkFrame` with the hour, minute, and
second widgets in it. You can specify the specific `GtkFrame` widget (useful when using the `Gtk.Builder` and `glade`). Time is guaranteed to be positive.
"""
function timewidget(t1::Dates.Time; widget=nothing, signal=nothing)
    const zerotime = Dates.Time(0,0,0)
    b = Gtk.GtkBuilder(filename=joinpath(@__DIR__, "time.glade"))
    if signal == nothing
        signal = Signal(t1)
    end
    S = map(signal) do x
        (Dates.Second(x), x)
    end
    M = map(S) do x
        x = last(x)
        (Dates.Minute(x), x)
    end
    H = map(M) do x
        x = last(x)
        (Dates.Hour(x), x)
    end
    t2 = map(last, H)
    bind!(signal, t2)
    Sint = Signal(Dates.value(first(value(S))))
    Ssb = spinbutton(-1:60, widget=b["second"], signal=Sint)
    foreach(Sint) do x
        Δ = Dates.Second(x) - first(value(S))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Second(new_t)
        push!(S, (new_x, new_t))
    end
    Sint2 = map(src -> Dates.value(Dates.Second(src)), t2)
    Sint3 = droprepeats(Sint2)
    bind!(Sint, Sint3, false)
    Mint = Signal(Dates.value(first(value(M))))
    Msb = spinbutton(-1:60, widget=b["minute"], signal=Mint)
    foreach(Mint) do x
        Δ = Dates.Minute(x) - first(value(M))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Minute(new_t)
        push!(M, (new_x, new_t))
    end
    Mint2 = map(src -> Dates.value(Dates.Minute(src)), t2)
    Mint3 = droprepeats(Mint2)
    bind!(Mint, Mint3, false)
    Hint = Signal(Dates.value(first(value(H))))
    Hsb = spinbutton(0:23, widget=b["hour"], signal=Hint)
    foreach(Hint) do x
        Δ = Dates.Hour(x) - first(value(H))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Hour(new_t)
        push!(H, (new_x, new_t))
    end
    Hint2 = map(src -> Dates.value(Dates.Hour(src)), t2)
    Hint3 = droprepeats(Hint2)
    bind!(Hint, Hint3, false)

    if widget == nothing
        return TimeWidget2(signal, b["frame"])
    else
        push!(widget, b["frame"])
        return TimeWidget2(signal, widget)
    end
end

"""
    datetimewidget(datetime)

Return a datetime widget that includes the `DateTime` and a `GtkBox` with the
year, month, day, hour, minute, and second widgets in it. You can specify the
specific `SpinButton` widgets for the hour, minute, and second (useful when using
`Gtk.Builder` and `glade`). Date and time are guaranteed to be positive.
"""
function datetimewidget(t1::DateTime; widget=nothing, signal=nothing)
    const zerotime = DateTime(0,1,1,0,0,0)
    b = Gtk.GtkBuilder(filename=joinpath(@__DIR__, "datetime.glade"))
    # t1 = eps(t0) < Dates.Second(1) ? round(t0, Dates.Second(1)) : t0
    if signal == nothing
        signal = Signal(t1)
    end
    S = map(signal) do x
        (Dates.Second(x), x)
    end
    M = map(S) do x
        x = last(x)
        (Dates.Minute(x), x)
    end
    H = map(M) do x
        x = last(x)
        (Dates.Hour(x), x)
    end
    d = map(H) do x
        x = last(x)
        (Dates.Day(x), x)
    end
    m = map(d) do x
        x = last(x)
        (Dates.Month(x), x)
    end
    y = map(m) do x
        x = last(x)
        (Dates.Year(x), x)
    end
    t2 = map(last, y)
    bind!(signal, t2)
    Sint = Signal(Dates.value(first(value(S))))
    Ssb = spinbutton(-1:60, widget=b["second"], signal=Sint)
    foreach(Sint) do x
        Δ = Dates.Second(x) - first(value(S))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Second(new_t)
        push!(S, (new_x, new_t))
    end
    Sint2 = map(src -> Dates.value(Dates.Second(src)), t2)
    Sint3 = droprepeats(Sint2)
    bind!(Sint, Sint3, false)
    Mint = Signal(Dates.value(first(value(M))))
    Msb = spinbutton(-1:60, widget=b["minute"], signal=Mint)
    foreach(Mint) do x
        Δ = Dates.Minute(x) - first(value(M))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Minute(new_t)
        push!(M, (new_x, new_t))
    end
    Mint2 = map(src -> Dates.value(Dates.Minute(src)), t2)
    Mint3 = droprepeats(Mint2)
    bind!(Mint, Mint3, false)
    Hint = Signal(Dates.value(first(value(H))))
    Hsb = spinbutton(-1:24, widget=b["hour"], signal=Hint)
    foreach(Hint) do x
        Δ = Dates.Hour(x) - first(value(H))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Hour(new_t)
        push!(H, (new_x, new_t))
    end
    Hint2 = map(src -> Dates.value(Dates.Hour(src)), t2)
    Hint3 = droprepeats(Hint2)
    bind!(Hint, Hint3, false)
    dint = Signal(Dates.value(first(value(d))))
    dsb = spinbutton(-1:32, widget=b["day"], signal=dint)
    foreach(dint) do x
        Δ = Dates.Day(x) - first(value(d))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Day(new_t)
        push!(d, (new_x, new_t))
    end
    dint2 = map(src -> Dates.value(Dates.Day(src)), t2)
    dint3 = droprepeats(dint2)
    bind!(dint, dint3, false)
    mint = Signal(Dates.value(first(value(m))))
    msb = spinbutton(-1:13, widget=b["month"], signal=mint)
    foreach(mint) do x
        Δ = Dates.Month(x) - first(value(m))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Month(new_t)
        push!(m, (new_x, new_t))
    end
    mint2 = map(src -> Dates.value(Dates.Month(src)), t2)
    mint3 = droprepeats(mint2)
    bind!(mint, mint3, false)
    yint = Signal(Dates.value(first(value(y))))
    ysb = spinbutton(-1:10000, widget=b["year"], signal=yint)
    foreach(yint) do x
        Δ = Dates.Year(x) - first(value(y))
        new_t = value(signal) + Δ
        new_t = new_t < zerotime ? zerotime : new_t
        new_x = Dates.Year(new_t)
        push!(y, (new_x, new_t))
    end
    yint2 = map(src -> Dates.value(Dates.Year(src)), t2)
    yint3 = droprepeats(yint2)
    bind!(yint, yint3, false)

    if widget == nothing
        return TimeWidget2(signal, b["frame"])
    else
        push!(widget, b["frame"])
        return TimeWidget2(signal, widget)
    end
end

######################## ProgressBar #############################

immutable ProgressBar <: GtkReactive.Widget
    signal::Signal{Int}
    widget::Gtk.GtkProgressBar
    preserved::Vector{Any}

    function (::Type{ProgressBar})(signal::Signal{Int}, widget, preserved)
        obj = new(signal, widget, preserved)
        GtkReactive.gc_preserve(widget, obj)
        obj
    end
end

"""
    progressbar(n; widget=nothing, signal=nothing)

Create a progressbar displaying the progress of a process; push to the widget new
iteration integers to update the progressbar. Optionally specify:
  - the GtkProgressBar `widget` (by default, creates a new one)
  - the (Reactive.jl) `signal` coupled to this progressbar (by default, creates a new signal)
"""
function progressbar(n::Int;
               widget=nothing,
               signal=nothing,
               syncsig=true,
               own=nothing)
    signalin = signal
    signal, value = GtkReactive.init_wsigval(Int, signal, GtkReactive.value(signal))
    if own == nothing
        own = signal != signalin
    end
    if widget == nothing
        widget = Gtk.GtkProgressBar()
    else
        setproperty!(widget, :fraction, value)
    end
    preserved = []
    if syncsig
        push!(preserved, map(signal) do val
            setproperty!(widget, :fraction, val/n)
        end)
    end
    if own
        Gtk.ondestroy(widget, preserved)
    end
    ProgressBar(signal, widget, preserved)
end

########################## ProgressBar ############################

immutable ProgressBar{T<:Number} <: InputWidget{T}
    signal::Signal{T}
    widget::GtkProgressBarLeaf
    id::Culong
    preserved::Vector

    function (::Type{ProgressBar{T}}){T}(signal::Signal{T}, widget, id, preserved)
        obj = new{T}(signal, widget, id, preserved)
        gc_preserve(widget, obj)
        obj
    end
end
ProgressBar{T}(signal::Signal{T}, widget::GtkProgressBarLeaf, id, preserved) =
    ProgressBar{T}(signal, widget, id, preserved)

progressbar(signal::Signal, widget::GtkProgressBarLeaf, id, preserved = []) =
    ProgressBar(signal, widget, id, preserved)

range2fraction(range::Range{T}, i::T) where T<:Number = (i - first(range) + step(range))/step(range)/length(range)

"""
    progressbar(range; widget=nothing, value=nothing, signal=nothing, orientation="horizontal")

Create a progressbar widget with the specified `range`. Optionally provide:
  - the GtkProgressBar `widget` (by default, creates a new one)
  - the starting `value` (defaults to the median of `range`)
  - the (Reactive.jl) `signal` coupled to this progressbar (by default, creates a new signal)
  - the `orientation` of the progressbar.
"""
function progressbar{T}(range::Range{T};
                   widget=nothing,
                   value=nothing,
                   signal=nothing,
                   orientation="horizontal",
                   syncsig=true,
                   own=nothing)
    signalin = signal
    signal, value = init_wsigval(T, signal, value; default=first(range))
    if own == nothing
        own = signal != signalin
    end
    if widget == nothing
        widget = GtkProgressBar()
        Gtk.G_.size_request(widget, 200, -1)
    else
        adj = Gtk.Adjustment(widget)
        Gtk.G_.fraction(adj, range2fraction(range, value))
    end
    Gtk.G_.fraction(widget, range2fraction(range, value))

    ## widget -> signal
    # id = signal_connect(widget, :value_changed) do w
        # push!(signal, defaultgetter(w))
    # end

    ## signal -> widget
    preserved = []
    if syncsig
        push!(preserved, init_signal2widget(widget, id, map(x -> range2fraction(range, x), signal)))
    end
    if own
        ondestroy(widget, preserved)
    end

    ProgressBar(signal, widget, id, preserved)
end

using Gtk.ShortNames, GtkReactive
w = Window("BeetleWay")
p = progressbar(2:2:12)
push!(w, p)
showall(w)

push!(p, 7)


text:
contact_control_path = "/home/taylor/Research/ContactControl.jl/src"

using Parameters

# Utilities
include(joinpath(contact_control_path, "utils.jl"))

# Solver
include(joinpath(contact_control_path, "solver/cones.jl"))
include(joinpath(contact_control_path, "solver/interior_point.jl"))
include(joinpath(contact_control_path, "solver/lu.jl"))

n = 3
nz = 3 * 7
nθ = 3 * 3

function residual(z, θ, κ)
    u = z[1:n]
    s1 = z[n .+ (1:n)]
    s2 = z[n + n .+ (1:n)]
    y1 = z[n + n + n .+ (1:n)]
    y2 = z[n + n + n + n .+ (1:n)]
    z1 = z[n + n + n + n + n .+ (1:n)]
    z2 = z[n + n + n + n + n + n .+ (1:n)]

    ū = θ[1:n]
    ul = θ[n .+ (1:n)]
    uu = θ[n + n .+ (1:n)]

    [
     (u - ū) - y1 + y2;
     -y1 - z1;
     -y2 - z2;
     (uu - u) - s1;
     (u - ul) - s2;
     s1 .* z1 .- κ;
     s2 .* z2 .- κ;
    ]
end

@variables z[1:nz], θ[1:nθ], κ[1:1]

r = residual(z, θ, κ)
r .= simplify(r)
r_func = eval(Symbolics.build_function(r, z, θ, κ)[2])

rz = Symbolics.jacobian(r, z)
rz = simplify.(rz)
rz_func = eval(Symbolics.build_function(rz, z, θ)[2])

rθ = Symbolics.jacobian(r, θ)
rθ = simplify.(rθ)
rθ_func = eval(Symbolics.build_function(rθ, z, θ)[2])

idx_ineq = collect([(n .+ (1:(n + n)))..., (n + n + n + n + n .+ (1:(n + n)))...])

opts = InteriorPointOptions(
  κ_init = 0.1,
  κ_tol = 1.0e-5,
  r_tol = 1.0e-8,
  diff_sol = true)

ū = [0.0, 0.0, 0.0]
ul = [-1.0, -1.0, -1.0]
uu = [1.0, 1.0, 1.0]

z0 = 0.1 * ones(nz)
z0[1:n] = copy(ū)
θ0 = [ū; ul; uu]

# solver
ip = interior_point(z0, θ0,
  r! = r_func, rz! = rz_func,
  rz = similar(rz, Float64),
  rθ! = rθ_func,
  rθ = similar(rθ, Float64),
  idx_ineq = idx_ineq,
  opts = opts)

interior_point_solve!(ip)

ip.z[1:n]

ext: jl

max_issues_repo_name: thowell/DirectMotionPlanning


text:
abstract type SeedBlk end

# [201] Murdock Event Detection Blockette (60 bytes)
mutable struct Blk201 <: SeedBlk
  sig::Array{Float32,1}
  flags::UInt8
  t::Array{Int32,1}
  det::String
  snr::Array{UInt8,1} # includes spots for lookback and pick algorithm
  Blk201() = new(zeros(Float32, 3),
                 0x00,
                 Array{Int32,1}(undef,7),
                 "None",
                 zeros(UInt8,8))
end

#  [500] Timing Blockette (200 bytes)
mutable struct Blk500 <: SeedBlk
  vco_correction::Float32
  t::Array{Int32,1}
  μsec::Int8
  reception_quality::Int8
  exception_count::UInt32
  exception_type::String
  clock_model::String
  clock_status::String
  Blk500() = new(0.0f0,
                 Array{Int32,1}(undef,7),
                 Int8(0),
                 Int8(0),
                 0x0000,
                 "",
                 "",
                 "")
end

# [2000] Variable Length Opaque Data Blockette
mutable struct Blk2000 <: SeedBlk
  n::UInt32
  NB::UInt16
  os::UInt16
  flag::UInt8
  hdr::Array{UInt8,1}
  data::Array{UInt8,1}
  Blk2000() = new(
                  0x00000000,
                  0x0000,
                  0x0000,
                  0x00,
                  Array{UInt8,1}(undef,0),
                  Array{UInt8,1}(undef,0)
                  )
end

# Calibration blockettes: [300], [310], [320], [390]
mutable struct BlkCalib <: SeedBlk
  t::Array{Int32,1}
  n::UInt8
  flags::UInt8
  dur1::UInt32
  dur2::UInt32
  amplitude::Float32
  period::Float32
  channel::Array{UInt8,1}
  ref::UInt32
  coupling::Array{UInt8,1}
  rolloff::Array{UInt8,1}
  noise::Array{UInt8,1}
  BlkCalib() = new( Array{Int32,1}(undef,7),
                    0x00,
                    0x00,
                    0x00000000,
                    0x00000000,
                    zero(Float32),
                    zero(Float32),
                    Array{UInt8,1}(undef, 3),
                    0x00000000,
                    Array{UInt8,1}(undef, 12),
                    Array{UInt8,1}(undef, 12),
                    Array{UInt8,1}(undef, 8)
                  )
end

mutable struct SeisIOBuf
  fmt::UInt8    # 0x0a
  nx::UInt16    # 0x1000
  wo::UInt8     # 0x01
  tc::Int32     # Time correction
  swap::Bool
  calibs::Array{UInt16,1}

  # hdr:
  seq::Vector{UInt8}
  hdr::Vector{UInt8}
  hdr_old::Vector{UInt8}
  id::Array{UInt8,1}
  id_str::String
  #=  order of u16:
  1 year
  2 jdy
  3 record start time
  4 beginning of data
  5 first blockette
  6 position of next blockette (relative to record begin) =#

  # Values read from file/stream
  dt::Float64
  r1::Int16
  r2::Int16
  r1_old::Int16
  r2_old::Int16
  Δ::Int64
  xs::Bool
  k::Int
  n::UInt16

  # Data-related
  x::Array{Float32,1}         # Data buffer
  buf::Array{UInt8,1}         # Buffer for reading UInt8 data
  uint16_buf::Array{UInt16,1} # Used in SEG Y and SEED
  uint32_buf::Array{UInt32,1} # Unsigned 32-bit steim-encoded data
  int16_buf::Array{Int16,1}   # Buffer for int16s in headers
  int32_buf::Array{Int32,1}   # Buffer for int32s in headers
  int64_buf::Array{Int64,1}   # Buffer for int64s; used in rseis/wseis
  sac_fv::Array{Float32,1}
  sac_iv::Array{Int32,1}
  sac_cv::Array{UInt8,1}
  sac_dv::Array{Float64,1}
  flags::Array{UInt8,1}       # SEED flags (Stored as four UInt8s)

  # For parsing dates
  date_buf::Array{Int32,1}

  # For mini-SEED data header
  dh_arr::Array{UInt8, 1}
  dh_buf::Base.GenericIOBuffer{Array{UInt8,1}}

  # Blockette containers
  B201::Blk201
  B500::Blk500
  B2000::Blk2000
  Calib::BlkCalib

  function SeisIOBuf()
    new(0xff,                               # fmt
        0x1000,                             # nx
        0x01,                               # wo
        zero(Int32),                        # tc
        false,                              # swap

        # SEED calibiration blockettes (all use the same parser)
        UInt16[0x012c, 0x0136, 0x0140, 0x0186],

        # SEED header
        Array{UInt8,1}(undef,8),            # seq::Vector{UInt8}
        Array{UInt8,1}(undef,12),           # hdr::Vector{UInt8}
        Array{UInt8,1}(undef,12),           # hdr_old::Vector{UInt8}
        Array{UInt8,1}(undef,15),           # id::Array{UInt8,1}
        "",                                 # id_str::String

        # computed
        0.0,                                # dt::Float64
        zero(Int16),                        # r1::Int16
        zero(Int16),                        # r2::Int16
        zero(Int16),                        # r1_old::Int16
        zero(Int16),                        # r2_old::Int16
        0,                                  # Δ::Int64
        false,                              # xs::Bool
        0,                                  # k::Int
        0x0000,                             # n::UInt16

        # data-related arrays
        Array{Float32,1}(undef, 65535),     # x::Array{Float32,1}
        Array{UInt8,1}(undef, 65535),       # buf::Array{UInt8,1}
        Array{UInt16,1}(undef, 6),          # uint16_buf::Array{UInt16,1}
        Array{UInt32,1}(undef, 16384),      # uint32_buf::Array{UInt32,1}
        Array{Int16,1}(undef, 62),          # int16_buf::Array{Int16,1}
        Array{Int32,1}(undef, 100),         # int32_buf::Array{Int32,1}
        Array{Int64,1}(undef, 8),           # int64_buf::Array{Int64,1}
        Array{Float32,1}(undef, 70),        # sac_fv::Array{Float32,1}
        Array{Int32,1}(undef, 40),          # sac_iv::Array{Int32,1}
        Array{UInt8,1}(undef, 192),         # sac_cv::Array{UInt8,1}
        Array{Float64,1}(undef,22),         # sac_dv::Array{Float64,1}
        Array{UInt8,1}(undef, 4),           # flags::Array{UInt8,1}

        # dedicated array for parsing dates
        Array{Int32,1}(undef, 7),           # date_buf::Array{Int32,1}

        # dedicated array for mini-SEED data header
        Array{UInt8,1}(undef, 48),          # dh_arr::Array{UInt8,1}
        IOBuffer(),                         # dh_buf::Base.GenericIOBuffer{Array{UInt8,1}}

        # blockette fields
        Blk201(),                           # Blk201:: Blk201
        Blk500(),                           # Blk500:: Blk500
        Blk2000(),                          # Blk2000:: Blk2000
        BlkCalib()                          # Calib:: BlkCalib
        )
  end
end

const BUF = SeisIOBuf()
BUF.dh_buf = IOBuffer(BUF.dh_arr)


text:
#! /usr/bin/env julia

using Pkg

using soapysdr_jll

using Clang.Generators

include_dir = joinpath(soapysdr_jll.artifact_dir, "include") |> normpath
clang_dir = joinpath(include_dir, "clang-c")

options = load_options(joinpath(@__DIR__, "generator.toml"))

@show options

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$include_dir")

headers = [joinpath(include_dir, "SoapySDR", header) for header in readdir(joinpath(include_dir, "SoapySDR")) if endswith(header, ".h")]
# there is also an experimental `detect_headers` function for auto-detecting top-level headers in the directory
# headers = detect_headers(clang_dir, args)

filter!(s -> basename(s) != "Logger.h", headers) # ignore Logger as it hangs the build

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx)


text:
using Klimakoffer

Klimakoffer.upscale_world(joinpath(@__DIR__,"..","input","world"), "The_World128x65.dat",150)

upscaled = Klimakoffer.read_geography(joinpath(@__DIR__,"..","input","world","The_World150x76.dat"), 150, 76)

ref = Klimakoffer.read_geography(joinpath(@__DIR__,"..","examples","reference_data","world","ref_The_World150x76.dat"), 150, 76)

diff = ref-upscaled

res = sum(abs.(diff))


text:
# TESTS - the classical test statistics for linear restrictions
#
# written by M. Creel, March 10, 1998.  michael.creel@uab.es
# converted to Octave 13 Jan. 2003
# converted to Julia 25 Aug. 2017
#
# purpose: calculates qF, Wald, Score and Likelihood Ratio tests for
# linear model						 y=XB+e
# 									 e~N(0,sig^2*I_n)
# subject to linear restrictions  	 RB=r
#
# format: [qF, W, LR, S] = TestStatistics(y,x,R,r)
#
# inputs:
# 		 y: nx1 dependent variable
# 		 x: nxk regressor matrix
# 		 R: R above, a qxk matrix
# 		 r: r above, a qx1 vector
#
# returns:
#        qF: the qF statistic
# 		 W: the Wald statistic
# 		 S: the score statistic
# 		 LR: the likelihood ratio statistic
using Distributions
function TestStatistics(y, x, R, r)

    n,k = size(x)
    q = size(R,1)
	b = x\y
	xx_inv = inv(x'*x)
	P_inv = inv(R*xx_inv*R')
	b_r = b - xx_inv*R'*P_inv*(R*b-r)
    e = y - x*b
    ess = (e'*e)[1,1]
    e_r = y - x*b_r
    ess_r = (e_r' * e_r)[1,1]
    sigsqhat_ols = ess/(n-k)
    sigsqhat_mle = ess/(n)
    sigsqhat_mle_r = ess_r/(n)
    # F-test
    F = (ess_r-ess)/(q*sigsqhat_ols)
    # Wald test (uses unrestricted model's est. of sig^2
    W = (R*b-r)'*P_inv*(R*b-r)/sigsqhat_mle
    # Score test (uses restricted model's est. of sig^2
    P_x = x * xx_inv * x'
    S = e_r' * P_x * e_r/(sigsqhat_mle_r)
    # LR test
    lnl = -n/2*log(2*pi) - n/2*log(sigsqhat_mle) - e' * e/(2*sigsqhat_mle)
    lnl_r = -n/2*log(2*pi) - n/2*log(sigsqhat_mle_r) - e_r' * e_r/(2*sigsqhat_mle_r)
    LR = 2.0*(lnl-lnl_r)
	tests = [q*F[1], W[1], LR[1], S[1]]
	pvalues = ccdf(Chisq(q), tests)
    tests = [tests pvalues]
	TESTS = ["qF","Wald","LR","Score"]
	labels = ["Value","p-value"]
    prettyprint(tests, labels, TESTS)
    return F, W, LR, S
end


text:
# this is a modified version of SavingCallback from
# https://github.com/SciML/DiffEqCallbacks.jl/blob/master/src/saving.jl

#=
"""
    SavedValues{tType<:Real, output_valsType}
A struct used to save values of the time in `t::Vector{tType}` and
additional values in `output_vals::Vector{output_valsType}`.
"""
struct SavedValues{tType, output_valsType}
    t::Vector{tType}
    output_vals::Vector{output_valsType}
end

"""
    SavedValues(tType::DataType, output_valsType::DataType)
Return `SavedValues{tType, output_valsType}` with empty storage vectors.
"""
function SavedValues(::Type{tType}, ::Type{output_valsType}) where {tType,output_valsType}
    SavedValues{tType, output_valsType}(Vector{tType}(), Vector{output_valsType}())
end

function Base.show(io::IO, saved_values::SavedValues)
    tType = eltype(saved_values.t)
    output_valsType = eltype(saved_values.output_vals)
    print(io, "SavedValues{tType=", tType, ", output_valsType=", output_valsType, "}",
                "\nt:\n", saved_values.t, "\noutput_vals:\n", saved_values.output_vals)
end

@recipe function plot(saved_values::SavedValues)
    DiffEqArray(saved_values.t, saved_values.output_vals)
end
=#

mutable struct SavingEvent{SaveFunc, SavedValues, saveatType, saveatCacheType}
  save_func::SaveFunc
  saved_values::SavedValues
  saveat::saveatType
  saveat_cache::saveatCacheType
  save_everystep::Bool
  save_start::Bool
  save_end::Bool
  save_scope::Bool
  saveiter::Int
end

function (affect!::SavingEvent)(integrator,force_save = false; scope = :ode_)

  just_saved = false
  # see OrdinaryDiffEq.jl -> integrator_utils.jl, function output_valsues!
  while !isempty(affect!.saveat) && integrator.tdir * first(affect!.saveat) <= integrator.tdir * integrator.t # Perform saveat
      affect!.saveiter += 1
      curt = pop!(affect!.saveat) # current time

      if curt != integrator.t # If <t, interpolate
          if typeof(integrator) <: OrdinaryDiffEq.ODEIntegrator
              # Expand lazy dense for interpolation
              DiffEqBase.addsteps!(integrator)
          end
          if !DiffEqBase.isinplace(integrator.sol.prob)
              curu = integrator(curt)
          else
              curu = first(get_tmp_cache(integrator))
              integrator(curu,curt) # inplace since save_func allocates
          end
          copyat_or_push!(affect!.saved_values.t, affect!.saveiter, curt)
          affect!.save_scope && copyat_or_push!(affect!.saved_values.scope, affect!.saveiter, scope, Val{false})
          copyat_or_push!(affect!.saved_values.u, affect!.saveiter,
                          affect!.save_func(curu, curt, integrator),Val{false})
      else # ==t, just save
          just_saved = true
          copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
          affect!.save_scope && copyat_or_push!(affect!.saved_values.scope, affect!.saveiter, scope, Val{false})
          copyat_or_push!(affect!.saved_values.u, affect!.saveiter, affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
      end
  end
  if !just_saved &&
      affect!.save_everystep || force_save ||
      (affect!.save_end && integrator.t == integrator.sol.prob.tspan[end])

      affect!.saveiter += 1
      copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
      affect!.save_scope && copyat_or_push!(affect!.saved_values.scope, affect!.saveiter, scope, Val{false})
      copyat_or_push!(affect!.saved_values.u, affect!.saveiter, affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
  end
  u_modified!(integrator, false)
end

function saving_initialize(cb, u, t, integrator)
  if cb.affect!.saveiter != 0
      if integrator.tdir > 0
          cb.affect!.saveat = BinaryMinHeap(cb.affect!.saveat_cache)
      else
          cb.affect!.saveat = BinaryMaxHeap(cb.affect!.saveat_cache)
      end
      cb.affect!.saveiter = 0
  end
  clear_savings(cb.affect!.saved_values)
  cb.affect!.save_start && save_timepoint!(integrator, :start_)
end

function saving_wrapper(save_func, saved_values::SavedValues;
                      saveat=Vector{eltype(saved_values.t)}(),
                      save_everystep=isempty(saveat),
                      save_start = save_everystep || isempty(saveat) || saveat isa Number,
                      save_end = save_everystep || isempty(saveat) || saveat isa Number,
                      save_scope = true,
                      tdir=1)
  # saveat conversions, see OrdinaryDiffEq.jl -> integrators/type.jl
  saveat_vec = collect(saveat)
  if tdir > 0
      saveat_internal = BinaryMinHeap(saveat_vec)
  else
      saveat_internal = BinaryMaxHeap(saveat_vec)
  end
  affect! = SavingEvent(save_func, saved_values, saveat_internal, saveat_vec, save_everystep, save_start, save_end, save_scope, 0)
  condtion = (u, t, integrator) -> true
  DiscreteCallback(condtion, affect!;
                   initialize = saving_initialize,
                   save_positions=(false,false))
end

# force saving at current timepoint
function save_timepoint!(integrator::DiffEqBase.AbstractODEIntegrator, scope=:ode_)

  affect! = integrator.opts.callback.discrete_callbacks[1].affect!
  affect!.saveiter += 1
  copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
  affect!.save_scope && copyat_or_push!(affect!.saved_values.scope, affect!.saveiter, scope, Val{false})
  copyat_or_push!(affect!.saved_values.u, affect!.saveiter, affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
end

# check saveat points before applying a callback
function save_after_step!(integrator, scope = :ode_)
  affect! = integrator.opts.callback.discrete_callbacks[1].affect!
  while !isempty(affect!.saveat) && integrator.tdir * first(affect!.saveat) <= integrator.tdir * integrator.t # Perform saveat
    affect!.saveiter += 1
    curt = pop!(affect!.saveat) # current time
    curu = integrator(curt)

    if curt != integrator.t # If <t, interpolate
      if typeof(integrator) <: OrdinaryDiffEq.ODEIntegrator
          # Expand lazy dense for interpolation
          DiffEqBase.addsteps!(integrator)
      end
      if !DiffEqBase.isinplace(integrator.sol.prob)
          curu = integrator(curt)
      else
          curu = first(get_tmp_cache(integrator))
          integrator(curu,curt) # inplace since save_func allocates
      end
      copyat_or_push!(affect!.saved_values.t, affect!.saveiter, curt)
      affect!.save_scope && copyat_or_push!(affect!.saved_values.scope, affect!.saveiter, scope, Val{false})
      copyat_or_push!(affect!.saved_values.u, affect!.saveiter,
                      affect!.save_func(curu, curt, integrator),Val{false})
    else # ==t, just save
      copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
      affect!.save_scope && copyat_or_push!(affect!.saved_values.scope, affect!.saveiter, scope, Val{false})
      copyat_or_push!(affect!.saved_values.u, affect!.saveiter, affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
    end
  end
end


text:
using Luxor
const JULIA_RED = (0.796, 0.235, 0.2)
const JULIA_PURPLE = (0.584, 0.345, 0.698)
const JULIA_BLUE = (0.251, 0.388, 0.847)
const JULIA_GREEN = (0.22, 0.598, 0.149)

const JULIA_LIGHT_RED = (0.828125, 0.40234375, 0.390625)
const JULIA_LIGHT_PURPLE = (0.66796875, 0.47265625, 0.7421875)
const JULIA_LIGHT_GREEN = (0.375, 0.671875, 0.31640625)

function halfring(w, h, c, bk = "white")
  @layer begin
    box(Point(-256, -256), Point(256, 1), :clip)
    sethue(JULIA_LIGHT_RED)
    ellipse(O, w, h, :fill)
    sethue(c)
    setline(5)
    ellipse(O, w, h, :stroke)
  end
  @layer begin
    box(Point(-256, -256), Point(256, 2), :clip)
    sethue(bk)
    ellipse(O, 0.8w, 0.8h, :fill)
    sethue(c)
    setline(3)
    ellipse(O, 0.8w, 0.8h, :stroke)
  end
end

function planet(w, h, deg, rc, pc, bk = "white")
  C = - Point(0.25h,0.35h)
  L = (0.9, 0.6, 1, 1)
  @layer (rotate(deg2rad(deg)); halfring(w, h, rc, bk))
  @layer begin
    sethue(JULIA_LIGHT_PURPLE)
    circle(O, h, :fill)
    sethue(pc)
    setline(5)
    circle(O, h, :stroke)
  end
  @layer (rotate(-π + deg2rad(deg)); halfring(w, h, rc, bk))
  @layer begin
    rotate(-π + deg2rad(deg))
    ellipse(O, 0.8w, 0.8h, :clip)
    @layer begin
      sethue(JULIA_LIGHT_PURPLE)
      circle(O, h, :fill)
      sethue(pc)
      setline(5)
      circle(O, h, :stroke)
    end
  end
  @layer totem(h, JULIA_PURPLE)
end

function logotext(h)
  @layer begin
    setline(0.5)
    fontsize(57)
    fontface("Tamil MN")
    textoutlines("Transformers.jl", Point(-2.2h, 1.8h), :path, halign=:left)
    fillpreserve()
  end
end

function totem(h, pc)
  @layer begin
    sethue(pc)
    setline(3)
    #left up
    C = O - 0.3h
    r = 0.2h
    circle(C, r, :stroke)
    line(C - Point(r, 0), C - Point(2r, 0), :stroke)
    sector(C, 2r, 2r, π/6, 1.2π, :stroke)
    line(C + 4r/3 + 1, C + 2r +r/3, :stroke)
    C = O + 0.4h + Point(0.2h, -0.2h)
  end
end

function moon(h, gc)
  @layer begin
    P = O + 0.3h + Point(0, 0.35h)
    l = 0.15h
    sethue(JULIA_LIGHT_GREEN)
    setline(3)
    circle(P, l, :fill)
    sethue(gc)
    circle(P, l, :stroke)
    setline(2)
    circle(P - Point(-0.04h, 0.03h), 0.03h, :stroke)
  end
end

function logo(h, r, rc, pc, gc, text)
  @layer planet(5h, h, r, rc, pc, "white")
  @layer moon(h, gc)
  text && @layer logotext(h)
end

const h = 100
const r = 30

@svg begin
  scale(1.15)
  translate(Point(0, -25))
  logo(h, r, JULIA_RED, JULIA_PURPLE, JULIA_GREEN, true)
end 512 400 "transformerslogo.svg"

@png begin
  scale(1.15)
  translate(Point(0, -25))
  logo(h, r, JULIA_RED, JULIA_PURPLE, JULIA_GREEN, true)
end 512 400 "transformerslogo.png"

@svg begin
  scale(1.15)
  logo(h, r, JULIA_RED, JULIA_PURPLE, JULIA_GREEN, false)
end 512 320 "logo.svg"

@png begin
  scale(1.15)
  logo(h, r, JULIA_RED, JULIA_PURPLE, JULIA_GREEN, false)
end 512 320 "logo.png"

@png begin
  scale(2.2)
  logo(h, r, JULIA_RED, JULIA_PURPLE, JULIA_GREEN, false)
end 1024 640 "logolarge.png"


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ e2d165ea-8efe-11eb-0505-65311ef915fa
using Pkg

# ╔═╡ c41ad652-8efe-11eb-2bec-45044027da71
using CSV, DataFrames, Plots, StatsPlots, EmpiricalCDFs, StatsBase, ShiftedArrays, Statistics, Distributions, Dates

# ╔═╡ 70a19b4e-8f05-11eb-224b-5df1b5308b30
Pkg.add("ShiftedArrays")

# ╔═╡ 84c098cc-8f06-11eb-29b3-e11faff4f67e
pyplot()

# ╔═╡ 1b2f8ff0-8eff-11eb-181e-3542d87b88fd
btcdata = DataFrame(
	CSV.File("BTC-USD.csv",
		# types=Dict("Open" => Float16,
		# 	 "High" => Float16,
		# 	 "Low" => Float16,
		# 	 "Close" => Float16,
		# 	 "Adj Close" => Float16,
		# 	 "Volume" => Int)
	)
)

# ╔═╡ cd0240d6-8f10-11eb-30c0-f72bc355d709
describe(btcdata)

# ╔═╡ e6ea09ea-8eff-11eb-2e0c-7d3e532fef5b
@df btcdata plot(:Date, :Close)

# ╔═╡ 67ad1e5e-8f01-11eb-239a-877da6996b3a
@df btcdata plot(:Date, :Close, yaxis=:log)

# ╔═╡ 841df9f2-8f01-11eb-10af-fde225b0ed04
@df btcdata histogram(:Close)

# ╔═╡ 0818e5bc-8f07-11eb-3b0f-4712461dfc62
returns = [btcdata[i, :Close] - btcdata[i-1,:Close] for i=2:length(btcdata[:Close]) ]

# ╔═╡ 21c6e3fc-8f06-11eb-3e81-652aad55e40e
logreturns = [ log(btcdata[i, :Close] / btcdata[i-1,:Close]) for i=2:length(btcdata[:Close]) ]

# ╔═╡ 1293b38c-8f07-11eb-3fac-814f22d8ecd8
plot(
	histogram(returns),
	histogram(rand(Normal(mean(returns), std(returns)), length(returns)), xlims=(minimum(returns),maximum(returns))),
	layout=(2,1)
)

# ╔═╡ fbad7c0c-8f10-11eb-2eac-f3588fdce9fb
histogram([returns, rand(Normal(mean(returns), std(returns)), length(returns))], bins=-10^4:200:10^4)

# ╔═╡ 37abb3d4-8f06-11eb-1f55-098daa698f59
histogram(logreturns)

# ╔═╡ 75b8244c-8f08-11eb-2413-f9e77b9c0cec
function empiricalcdf(v)
	empcdf = EmpiricalCDF()
	append!(empcdf, v)
	sort!(empcdf)
	return empcdf
end

# ╔═╡ b5023930-8f03-11eb-3e0d-c598e1526937
begin
	lr_close_empcdf = empiricalcdf(logreturns)
	r_close_empcdf = empiricalcdf(returns)
	close_empcdf = empiricalcdf(btcdata[:Close])
end

# ╔═╡ 2c76370e-8f07-11eb-1b20-41861bebfe81
plot(
	plot(returns, r_close_empcdf(returns), seriestype=:scatter, label="returns"),
	plot(logreturns, lr_close_empcdf(logreturns), seriestype=:scatter, label="logreturns"),
	plot(btcdata[:Close], close_empcdf(btcdata[:Close]), seriestype=:scatter, label="close"),
	layout=(1,3)
)

# ╔═╡ 39c023b6-8f07-11eb-3b5c-7341e02b7e90
plot(
	returns,
	r_close_empcdf(returns),
	seriestype=:scatter,
	yaxis=:log,
	xaxis=:log,
	xlims=(1,10^5)
)

# ╔═╡ 11458512-8f04-11eb-08c4-733a2ae79b86
plot(
	logreturns,
	lr_close_empcdf(logreturns),
	seriestype=:scatter,
	yaxis=:log,
	xaxis=:log,
	xlims=(0.01,0.5)
)

# ╔═╡ Cell order:
# ╠═e2d165ea-8efe-11eb-0505-65311ef915fa
# ╠═70a19b4e-8f05-11eb-224b-5df1b5308b30
# ╠═c41ad652-8efe-11eb-2bec-45044027da71
# ╠═84c098cc-8f06-11eb-29b3-e11faff4f67e
# ╠═1b2f8ff0-8eff-11eb-181e-3542d87b88fd
# ╠═cd0240d6-8f10-11eb-30c0-f72bc355d709
# ╠═e6ea09ea-8eff-11eb-2e0c-7d3e532fef5b
# ╠═67ad1e5e-8f01-11eb-239a-877da6996b3a
# ╠═841df9f2-8f01-11eb-10af-fde225b0ed04
# ╠═0818e5bc-8f07-11eb-3b0f-4712461dfc62
# ╠═21c6e3fc-8f06-11eb-3e81-652aad55e40e
# ╠═1293b38c-8f07-11eb-3fac-814f22d8ecd8
# ╠═fbad7c0c-8f10-11eb-2eac-f3588fdce9fb
# ╠═37abb3d4-8f06-11eb-1f55-098daa698f59
# ╠═75b8244c-8f08-11eb-2413-f9e77b9c0cec
# ╠═b5023930-8f03-11eb-3e0d-c598e1526937
# ╠═2c76370e-8f07-11eb-1b20-41861bebfe81
# ╠═39c023b6-8f07-11eb-3b5c-7341e02b7e90
# ╠═11458512-8f04-11eb-08c4-733a2ae79b86


text:
# Auxiliary functions for working with vector of vectors
function vector_alloc(T, M, N)
    v = Vector{Vector{T}}(undef, N)
    for i in eachindex(v)
        v[i] = zeros(T, M)
    end
    v
end

function flatten_vector!(dest::T, src::Vector{T}) where T<:AbstractArray
    N = length(src)
    M = length(src[1])
    for i in eachindex(src)
        dest[((i-1)*M)+1:i*M] = src[i]
    end
end

function nest_vector!(dest::Vector{T}, src::T) where T<:AbstractArray
    M = length(dest[1])
    for i in eachindex(dest)
        copyto!(dest[i], src[(M*(i-1))+1:(M*i)])
    end
end


text:
module MultivariateBases

import MutableArithmetics
const MA = MutableArithmetics

using MultivariatePolynomials
const MP = MultivariatePolynomials

export AbstractPolynomialBasis
export maxdegree_basis, basis_covering_monomials, empty_basis
include("interface.jl")

export AbstractMonomialBasis, MonomialBasis, ScaledMonomialBasis
include("monomial.jl")
include("scaled.jl")

export FixedPolynomialBasis, AbstractMultipleOrthogonalBasis, ProbabilistsHermiteBasis, PhysicistsHermiteBasis, LaguerreBasis
export AbstractGegenbauerBasis, LegendreBasis, ChebyshevBasis, ChebyshevBasisFirstKind, ChebyshevBasisSecondKind
export univariate_orthogonal_basis, reccurence_first_coef, reccurence_second_coef, reccurence_third_coef, reccurence_deno_coef
include("fixed.jl")
include("orthogonal.jl")
include("hermite.jl")
include("laguerre.jl")
include("legendre.jl")
include("chebyshev.jl")

end # module

max_stars_repo_name: JuliaAlgebra/PolynomialBases

max_issues_repo_name: JuliaAlgebra/PolynomialBases

max_forks_repo_name: JuliaAlgebra/PolynomialBases


text:
using ArrayLayouts, Random, FillArrays, Test, SparseArrays, Base64
import ArrayLayouts: MemoryLayout, @_layoutlmul, triangulardata

Random.seed!(0)

include("test_layouts.jl")
include("test_muladd.jl")
include("test_ldiv.jl")
include("test_layoutarray.jl")
include("test_cumsum.jl")


text:
function _extract_lp_solution_sub(graph, x)
    # Same as for SpanningTree, but for a natively directed graph (no check
    # for reverse edges).
    solution = Edge{Int}[]
    for e in edges(graph)
        if JuMP.value(x[e]) >= 0.5
            push!(solution, e)
        end
    end
    return solution
end

function _extract_lp_solution(i::ElementaryPathInstance{Int, Maximise}, x)
    return _extract_lp_solution_sub(i.graph, x)
end

function _extract_lp_solution(
    i::MinimumBudget{ElementaryPathInstance{Int, Maximise}, Int},
    x,
)
    return _extract_lp_solution_sub(i.instance.graph, x)
end


text:
struct Configserror <:Exception
    msg::String
end

function deleteenvkey!(key::String)
    if haskey(ENV, key)
        delete!(ENV, key)
    end
end

function parseenvkey(key::String, value)
    haskey(ENV, key) ? ENV[key] : value
end

function parsecustomenv!(tree::Dict, key::String, value)
    haskey(ENV, value) ? (tree[key] = ENV[value]) : delete!(tree, key)
end

function parsecustomenv!(tree::Dict)
    for (key, value) in tree
        value isa Dict ? parsecustomenv!(tree[key]) : parsecustomenv!(tree, key, value)
    end
end

tosymbol(value) = value
tosymbol(dict::Dict) = Dict{Symbol, Any}(Symbol(key) => value for (key, value) in dict)

makeimmutable(value) = value
function makeimmutable(array::Array)
    shadow = Array{Any, 1}()
    for value in array
        push!(shadow, makeimmutable(value))
    end
    tuple(shadow...)
end

function makeimmutable(dict::Dict)
    shadow = Dict()
    dict = tosymbol(dict)
    for (key, value) in dict
        shadow[key] = makeimmutable(dict[key])
    end
    (; shadow...)
end

function override!(baseconf::Dict, newconf::Dict)
    for (key, value) in newconf
        if value isa Dict
            if haskey(baseconf, key) && baseconf[key] isa Dict
                override!(baseconf[key], newconf[key])
            else
                baseconf[key] = Dict{String, Any}();
                override!(baseconf[key], newconf[key])
            end
        else
            baseconf[key] = value
        end
    end
end

function getfiles(path::String, retry::Bool = false)
    try
        readdir(path)
    catch err
        if retry
            throw(Configserror("no such config directory: " * path))
        else
            getfiles(joinpath(pwd(), path), true)
        end
    end
end

function pathtodict(path::String, value)::Dict{String, Any}
    path === "" && throw(Configserror("a path is required to set a config"))
    subpaths = split(path, ".")
    base = Dict{String, Any}()
    ref = base
    for i in eachindex(subpaths)
        subpath = subpaths[i]
        if length(subpaths) === i
            ref[subpath] = value
            return base
        else
            ref[subpath] = Dict{String, Any}()
        end
        ref = ref[subpath]
    end
end

max_stars_repo_name: tisztamo/Configs

max_issues_repo_name: tisztamo/Configs

max_forks_repo_name: tisztamo/Configs


text:
@testset "Corpus" begin

    sample_text1 = "This is a string"
    sample_text2 = "This is also a string"
    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    sd = StringDocument(sample_text1)
    fd = FileDocument(sample_file)
    td = TokenDocument(sample_text1)
    ngd = NGramDocument(sample_text1)

    crps = Corpus([sd, fd, td, ngd])
    @test typeof(crps) <: Corpus{<:String, <:AbstractDocument}

    crps2 = Corpus([ngd, ngd])
    update_inverse_index!(crps2)
    @test typeof(crps2) <: Corpus{<:String, <:NGramDocument}

    documents(crps)

    for doc in crps
    	@test isa(doc, AbstractDocument)
    end

    @test isempty(lexicon(crps))
    @test isempty(inverse_index(crps))
    for ngram_complexity in [1,2]
        update_lexicon!(crps, ngram_complexity)
        @test !isempty(lexicon(crps))
        @test lexicon(crps) == create_lexicon(crps, ngram_complexity)

        update_inverse_index!(crps, ngram_complexity)
        @test !isempty(inverse_index(crps))
        @test inverse_index(crps) == create_inverse_index(crps, ngram_complexity)
    end

    @test hash_function(hash_function(crps)) === hash
    hash_function!(crps, TextHashFunction(hash, 10))
    @test hash_function(hash_function(crps)) === hash
    @test cardinality(hash_function(crps)) == 10

    # Indexing
    @test crps[1] isa StringDocument
    @test text(crps[1]) == "This is a string"
    @test crps["string"] == [1, 3, 4]
    @test crps2["string"] == [1, 2]

    # Directory corpus
    function generate_dir_corpus()
        # Make directories
        tmp_path = tempdir()
        data_path = abspath(joinpath(tmp_path, "StringAnalysis"))
        mkpath(data_path)
        directories = ["XXX","XXX/YYY","ZZZ"]
        for (i, dir) in enumerate(directories)
            file_path = joinpath(data_path, dir)
            mkpath(file_path)
            filename = joinpath(file_path, "file_$i.txt")
            open(filename, "w") do fid
                text = """
                       This is what
                       And this is who
                       One makes cat
                       The other boo.
                       ---
                       Filename: $(filename)
                       """
                write(fid, text)
            end
        end
        return data_path, directories
    end

    path, dirs = generate_dir_corpus()
    crps = DirectoryCorpus(path)
    standardize!(crps, StringDocument{String})
    for doc in crps
        # The file path should be in the metadata
        # and the text data.
        if !Sys.isapple()  # apparently fails on Apple boxes
            @test occursin(StringAnalysis.name(doc), text(doc)) == true
        end
    end

    rm(path, recursive=true, force=true)

    # Cotainer treatment
    doc = StringDocument("~pushed~")
    n = length(crps)
    push!(crps, doc);
    @test length(crps) == n+1
    pdoc = pop!(crps)
    @test length(crps) == n
    @test text(pdoc) == text(doc)
    pushfirst!(crps, doc)
    @test length(crps) == n+1
    @test text(crps[1]) == text(doc)
    pdoc = popfirst!(crps)
    @test length(crps) == n
    @test text(pdoc) == text(doc)
    pos = 3
    insert!(crps, pos, pdoc)
    @test length(crps) == n+1
    @test text(crps[3]) == text(pdoc)
    deleteat!(crps, 3)
    @test length(crps) == n
end


text:
module MLBase

    using Reexport
    using Iterators
    using Compat
    using Compat: view
    @reexport using StatsBase

    import Base: length, show, keys, precision, length, getindex
    import Base: start, next, done
    import Base.Order: lt, Ordering, ForwardOrdering, ReverseOrdering, Forward, Reverse
    import StatsBase: RealVector, IntegerVector, RealMatrix, IntegerMatrix, RealArray

    export

    # reexport from Base.Order
    Forward, Reverse,

    # utils
    repeach,        # repeat each element in a vector
    repeachcol,     # repeat each column in a matrix
    repeachrow,     # repeat each row in a matrix

    # classification
    LabelMap,       # a type to represent a label map

    classify,       # predict class label(s) based on score values
    classify!,      # inplace version of classify
    classify_withscore,     # classify with additional output of best score
    classify_withscores,    # classify with additional output of best scores
    classify_withscores!,   # inplace version of classify_withscores
    labelmap,       # construct a label map from a list of labels
    labelencode,    # encode a sequence of discrete values using a label map
    labeldecode,    # decode the label to the associated discrete value
    groupindices,   # grouped indices based on labels

    # crossval
    CrossValGenerator,  # abstract base class for all cross-validation plans
    Kfold,              # K-fold cross validation plan
    StratifiedKfold,        # stratified K-fold cross validation plan
    LOOCV,              # leave-one-out cross validation plan
    RandomSub,              # repeated random subsampling cross validation
    StratifiedRandomSub,    # stratified repeated random subsampling

    cross_validate,     # perform cross-validation

    # perfeval
    ROCNums,        # A class to capture ROC numbers

    correctrate,    # compute correct rate of predictions
    errorrate,      # compute error rate of predictions
    confusmat,      # compute confusion matrix
    counthits,      # count the number of hits
    hitrate,        # compute hit-rate of ranked lists at a specific rank
    hitrates,       # compute hit-rate of ranked lists at multiple ranks
    roc,            # compute roc numbers from predictions
    true_positive,      # number of true positives
    true_negative,      # number of true negatives
    false_positive,     # number of false positives
    false_negative,     # number of false negatives
    true_positive_rate,     # rate of true positives
    true_negative_rate,     # rate of true negatives
    false_positive_rate,    # rate of false positives
    false_negative_rate,    # rate of false negatives
    recall,             # recall computed from ROCNums
    sensitivity,        # sensitivity computed from ROCNums
    specificity,        # specificity computed from ROCNums
    precision,          # precision computed from ROCNums
    positive_predictive_value,  # positive predictive value computed from ROCNums
    negative_predictive_value,  # negative predictive value computed from ROCNums
    false_discovery_rate,       # false discovery rate computed from ROCNums
    accuracy,           # accuracy computed from ROCNums
    f1score,            # F1-score computed from ROCNums

    # modeltune
    gridtune           # grid-based model tuning (search best config)

    # source files

    include("utils.jl")
    include("classification.jl")
    include("crossval.jl")
    include("perfeval.jl")
    include("modeltune.jl")

    include("deprecates.jl")
end


text:
# Autogenerated wrapper script for SuperLU_jll for x86_64-linux-musl
export libsuperlu

using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("SuperLU")
JLLWrappers.@declare_library_product(libsuperlu, "libsuperlu.so.5")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libsuperlu,
        "lib/libsuperlu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DifferentialEquations
using ModelingToolkit
using OrdinaryDiffEq
using DataFrames
using DataFrames
using StatsPlots
using BenchmarkTools

@parameters t β c γ
@variables S(t) I(t) R(t)
@derivatives D'~t
N=S+I+R # This is recognized as a derived variable
eqs = [D(S) ~ -β*c*I/N*S,
       D(I) ~ β*c*I/N*S-γ*I,
       D(R) ~ γ*I];

sys = ODESystem(eqs)
sys = ode_order_lowering(sys);

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;

u0 = [S => 990.0,
      I => 10.0,
      R => 0.0];

p = [β=>0.05,
     c=>10.0,
     γ=>0.25];

prob_ode = ODEProblem(sys,u0,tspan,p;jac=true)

sol_ode = solve(prob_ode);

df_ode = DataFrame(sol_ode(t)')
df_ode[!,:t] = t;

@df df_ode plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_ode)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: ChrisRackauckas/sir-julia

max_issues_repo_name: ChrisRackauckas/sir-julia

max_forks_repo_name: ChrisRackauckas/sir-julia


text:
"""
---
Apply OQAM demapping to incoming OQAM matrix qamMat of size nbSubcarriers x 2nbSymb
# --- Syntax
qamMat = oqamDemapping(qamMat)
# --- Input parameters
- oqamMat : OQAM matrix (pure real and pure imag. alterns) [Array{Complex{Float64,nbSubcarriers,nbSymb*2}}]
# --- Output parameters
- qamMat : Output complex qam Matrix [Array{Complex{Float64,nbSubcarriers,nbSymb}}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function oqamDemapping(oQamMat)
	# --- Size of incoming matrix
	(nbSubcarriers,nbOQAM)	  = size(oQamMat);
	nbSymb	  = nbOQAM ÷ 2;
	# Real and imag split:
	qamMat	  = zeros(Complex{Float64},nbSubcarriers,nbSymb);
	# --- Init matrix
	for iS = 1 : 1 : nbSubcarriers
		# Mask to keep real or imag part
		pattern		= (-1im)^( mod( mod(iS,2) + 1,2));
		pattern2	= (-1im)^( mod( mod(iS,2) + 2,2));
		# 1 // 1i shift
		scaleReIm	= (1im)^( mod( mod(iS,2) + 1,2));
		scaleReIm2	= (1im)^( mod( mod(iS,2) + 2,2));
		for iN = 1 : 1 : nbSymb
			# OQAM processing
			qamMat[iS,iN]	= scaleReIm*real(pattern*oQamMat[iS,2*(iN-1)+1])+scaleReIm2*real(pattern2*oQamMat[iS,2*(iN-1)+2]);
		end
	end
	return qamMat;
end

"""
---
Demodulate FBMC waveform based on the dual operation of fbmcSigGen

FBMC is parametrized by its FFT size, its cyclic prefix length (in samples) and the allocSubcarriers vector
# Syntax
sigId	= fbmcSigDecode(sigRx,nFFT,K,allocSubcarriers)
# ---  Input parameters
- sigRx	  : Time domain FBMC signal  [Array{Complex{Float64},nbEch}]
- nFFT	  : FFT size [Int]
- K		  : Overlapping factor [Int]
allocSubcarriers : Vector of index of allocated subcarriers [Array{Int,nbSubcarriers}]
# ---  Output parameters
- qamMat  : Time frequency matrix : [Array{Complex{Float64},nbSubcarriers,nbSymb}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function fbmcSigDecode(sigRx,nFFT,K,allocatedSubcarriers)
	# Set 4 core for FFT computation
	FFTW.set_num_threads(4)
	# --- Getting symbol size
	sizeSymb	   = nFFT * K;
	# --- Init FBMC filter
	p			   = getFBMCFilter(K,nFFT);
	# --- Getting number of symbol to decode
	# Here, symbol are QAM symbols (need to demodulate twice of that)
	nbSymbDecode   = Int(floor(length(sigRx)/ nFFT - K + 1/2));
	sBeforeFFT	   = zeros(Complex{Float64},nFFT,2*nbSymbDecode);
	# --- Iterative decoding with overlap and sum undoing
	for n = 1 : 1 : 2*nbSymbDecode
		# --- Received signal and Rx apodisation
		# Output is a KN vector multiplied by the filter
		yTmp	= p  .* sigRx[(n-1)*nFFT÷2 .+ collect(1:sizeSymb)];
		# From this convert to N sequence. All samples delayed by nFFT are sumed
		sBeforeFFT[:,n] = sum(reshape(yTmp,nFFT,K),dims=2);
	end
	# FFT
	oQamMat	= fft(sBeforeFFT,1);
	# --- OQAM to QAM conversion
	qamMat	= oqamDemapping(oQamMat[allocatedSubcarriers,:]);
end

# MD for input structure
function fbmcSigDecode(sigChan,fbmc::StrucFBMC)
	return fbmcSigDecode(sigChan,fbmc.nFFT,fbmc.K,fbmc.allocatedSubcarriers);
end


text:
# implementation/strides.jl
#
# Implements the stride calculations of the various problems

@generated function add_strides{N}(dims::NTuple{N,Int}, stridesA::NTuple{N,Int}, stridesC::NTuple{N,Int})
    minstridesex = Expr(:tuple,[:(min(stridesA[$d],stridesC[$d])) for d = 1:N]...)
    quote
        minstrides = $minstridesex
        p = sortperm(collect(minstrides))
        dims = _permute(dims, p)
        stridesA = _permute(stridesA, p)
        stridesC = _permute(stridesC, p)
        minstrides = _permute(minstrides, p)

        return dims, stridesA, stridesC, minstrides
    end
end

@generated function trace_strides{NA,NC}(dims::NTuple{NA,Int}, stridesA::NTuple{NA,Int}, stridesC::NTuple{NC,Int})
    M = div(NA-NC,2)
    dimsex = Expr(:tuple,[:(dims[$d]) for d=1:(NC+M)]...)
    stridesAex = Expr(:tuple,[:(stridesA[$d]) for d = 1:NC]...,[:(stridesA[$(NC+d)]+stridesA[$(NC+M+d)]) for d = 1:M]...)
    stridesCex = Expr(:tuple,[:(stridesC[$d]) for d = 1:NC]...,[0 for d = 1:M]...)
    minstridesex = Expr(:tuple,[:(min(stridesA[$d],stridesC[$d])) for d = 1:NC]...,[:(stridesA[$(NC+d)]+stridesA[$(NC+M+d)]) for d = 1:M]...)
    quote
        minstrides = $minstridesex
        p = sortperm(collect(minstrides))
        newdims = _permute($dimsex, p)
        newstridesA = _permute($stridesAex, p)
        newstridesC = _permute($stridesCex, p)
        minstrides = _permute(minstrides, p)

        return newdims, newstridesA, newstridesC, minstrides
    end
end

@generated function contract_strides{NA, NB, NC}(dimsA::NTuple{NA, Int}, dimsB::NTuple{NB, Int},
    stridesA::NTuple{NA, Int}, stridesB::NTuple{NB, Int}, stridesC::NTuple{NC, Int})
    meta = Expr(:meta, :inline)
    cN = div(NA+NB-NC, 2)
    oNA = NA - cN
    oNB = NB - cN

    dimsex = Expr(:tuple, [:(dimsA[$d]) for d = 1:oNA]..., [:(dimsB[$d]) for d = 1:oNB]..., [:(dimsA[$(oNA+d)]) for d = 1:cN]...)

    stridesAex = Expr(:tuple, [:(stridesA[$d]) for d = 1:oNA]..., [0 for d = 1:oNB]..., [:(stridesA[$(oNA+d)]) for d = 1:cN]...)
    stridesBex = Expr(:tuple, [0 for d = 1:oNA]..., [:(stridesB[$d]) for d = 1:oNB]..., [:(stridesB[$(oNB+d)]) for d = 1:cN]...)
    stridesCex = Expr(:tuple, [:(stridesC[$d]) for d = 1:(oNA+oNB)]..., [0 for d = 1:cN]...)

    minstridesex = Expr(:tuple, [:(min(stridesA[$d], stridesC[$d])) for d = 1:oNA]...,
    [:(min(stridesB[$d], stridesC[$(oNA+d)])) for d = 1:oNB]...,
    [:(min(stridesA[$(oNA+d)], stridesB[$(oNB+d)])) for d = 1:cN]...)
    quote
        $meta
        minstrides = $minstridesex
        p = sortperm(collect(minstrides))
        dims = _permute($dimsex, p)
        stridesA = _permute($stridesAex, p)
        stridesB = _permute($stridesBex, p)
        stridesC = _permute($stridesCex, p)
        minstrides = _permute(minstrides, p)

        return dims, stridesA, stridesB, stridesC, minstrides
    end
end


text:
# This file is a part of HistTools.jl, licensed under the MIT License (MIT).

# TODO: make it more general
"""
    reshape(h, range)

Change histogram format (min, max, bin size).
"""
function reshape(h::Histogram{<:Real,1}, format::AbstractRange; isdensity::Bool=false)
    hr = fit(Histogram,
             StatsBase.midpoints(h.edges[1]),
             Weights(h.weights),
             format)
    isdensity && (hr.weights .*= step(h.edges[1])/step(format))
    hr
end


text:
# Makes VRT of all the geotiffs in a folder and in a region.  One vrt per UTM zone.

using VAWTools
import ArchGDAL
const AG = ArchGDAL

regions = [11]
dir = "../output/rasters-fit_sigma/production-with-dhdt/"

for region in [3,4,7,11,14][end:end]
    @show region
    fls = Dict() # UTM zones => list of files
    for f in readdir(dir)
        ff, ext = splitext(f)
        ext!=".tif" && continue
        reg, nr = parse.(Int, split(split(ff, '-')[end], '.') )
        reg!=region && continue
        # get all projections
        #p = readstring(`gdalsrsinfo -o proj4 $(joinpath(dir,f))`)
        # Faster with ArchGDAL?
        p = AG.registerdrivers() do
            AG.read(joinpath(dir,f)) do dataset
                proj4 = strip(AG.toPROJ4(AG.importWKT(AG.getproj(dataset))))
            end
        end
        zone = [[nr,parse(Int, split(s,'=')[2])] for s in split(p) if split(s, '=')[1]=="+zone"][1]
        push!(get!(fls, zone[2], []), joinpath(dir,f))
    end

    for (z,fs) in fls
        @show z
        @show vrt = joinpath(dir, "composite-$(VAWTools.int2str2(region))-zone_$z.vrt")
        run(`gdalbuildvrt $vrt $(fs...)`)
    end
end


text:
global DISABLESTBPRTLINES = false

function togglePrtStbLines()
  global DISABLESTBPRTLINES
  DISABLESTBPRTLINES = !DISABLESTBPRTLINES
end

function plotLsrScanFeats(br::Array{Float64,2})
  Cart = zeros(size(br))
  Cart[:,1] = br[:,2].*cos(br[:,1])
  Cart[:,2] = br[:,2].*sin(br[:,1])
  plot(x=Cart[:,1],y=Cart[:,2],Geom.point,
  Guide.xticks(ticks=collect(-60:10:60)),
  Guide.yticks(ticks=collect(0:10:80)))
end

function drawFeatTrackers(trkrs::Dict{Int,Feature}, bfts::Array{Float64,2})
  musX = Float64[]
  varX = Float64[]
  musY = Float64[]
  varY = Float64[]
  allPtsX = Float64[]
  allPtsY = Float64[]

  for ftr in trkrs
    pts = getPoints(ftr[2].bel)
    allPtsX = [allPtsX; vec(pts[1,:])]
    allPtsY = [allPtsY; vec(pts[2,:])]

    push!(musX, Base.mean(vec(pts[1,:])))
    push!(varX, Base.std(vec(pts[1,:])))
    push!(musY, Base.mean(vec(pts[2,:])))
    push!(varY, Base.std(vec(pts[2,:])))
  end

  X = Float64[]
  Y = Float64[]

  if size(bfts,2) > 0
    if bfts[1,1] != 0.0 && bfts[2,1] != 0.0 && bfts[3,1] != 0.0
      for i in 1:size(bfts,2)
          u, R = p2c(vec(bfts[:,i]))
          push!(X, u[1])
          push!(Y, u[2])
      end
    end
  end

  # Guide.yticks(ticks=collect(-60:10:60)),
  # Guide.xticks(ticks=collect(0:10:80))
  p = plot(layer(x=musX, y=musY, Geom.point, Theme(default_color=colorant"red")),
  layer(x=allPtsX, y=allPtsY, Geom.histogram2d),
  Guide.yticks(ticks=collect(-70:10:70)),
  Guide.xticks(ticks=collect(-40:10:80)))
  for i in 1:length(X)
    push!(p.layers, Gadfly.layer(x=[0.0;X[i]], y=[0.0;Y[i]], Geom.line, Gadfly.Theme(default_color=colorant"magenta"))[1])
  end
  p
end

# moved to RobotUtils.jl
# lsrBR(a) = [a[2,:];a[1,:]]';

function saveImgSeq(d::Dict{Int,Array{Float64,2}}; from::Int=1,to::Int=10,step::Int=1)
  for i in from:step:to
    p = plotLsrScanFeats(lsrBR(d[i]));
    Gadfly.draw(PNG(string("imgs/img",i,".png"),25cm,25cm),p)
  end
  nothing
end

# --------------------------------------------------------------
# transfered in from IncrementalInference

## TODO -- you were here with port starboard lines
function stbPrtLineLayers!(pl, Xpp, Ypp, Thpp; l::Float64=5.0)
    if DISABLESTBPRTLINES
      return nothing
    end

    lnstpr = [0.0;l;0.0]
    lnstpg = [0.0;-l;0.0]

    Rd  =SE2(lnstpr)
    Gr = SE2(lnstpg)

    for i in 1:length(Xpp)
      lnstt = [Xpp[i];Ypp[i];Thpp[i]]
      Ps = SE2(lnstt)
      lnr = se2vee(Ps*Rd)
      lng = se2vee(Ps*Gr)
      xsr = [Xpp[i];lnr[1]]
      ysr = [Ypp[i];lnr[2]]
      xsg = [Xpp[i];lng[1]]
      ysg = [Ypp[i];lng[2]]

      push!(pl.layers, layer(x=xsr, y=ysr, Geom.path(), Gadfly.Theme(default_color=colorant"red", line_width=1.5pt))[1] )
      push!(pl.layers, layer(x=xsg, y=ysg, Geom.path(), Gadfly.Theme(default_color=colorant"green", line_width=1.5pt))[1] )
    end
    nothing
end

# function lblsFromTo(from,to)
#   lbls=String[]
#   [push!(lbls, "$(i)") for i in from:to]
#   return lbls
# end

function drawPoses(fg::FactorGraph; from::Int=0,to::Int=99999999,
                    meanmax=:max, lbls=true, drawhist=true,
                    spscale::Float64=5.0,
                    api::DataLayerAPI=IncrementalInference.localapi  )
    #Gadfly.set_default_plot_size(20cm, 30cm)
    Xp,Yp = get2DPoseSamples(fg, from=from, to=to)
    Xpp = Float64[]; Ypp=Float64[]; Thpp=Float64[]; LBLS=String[];
    if meanmax == :mean
      Xpp,Ypp, Thpp, LBLS = get2DPoseMeans(fg, from=from, to=to, api=api)
    elseif meanmax == :max
      Xpp,Ypp, Thpp, LBLS = get2DPoseMax(fg, from=from, to=to, api=api)
    end

    # lbls = lblsFromTo(1,length(Xpp))
    psplt = Union{}
    if lbls
      psplt = Gadfly.plot(
      Gadfly.layer(x=Xpp,y=Ypp,label=LBLS,Geom.path(), Theme(line_width=1pt), Geom.label)
      )
    else
      psplt = Gadfly.plot(
      Gadfly.layer(x=Xpp,y=Ypp,Geom.path(), Theme(line_width=1pt))
      )
    end
    stbPrtLineLayers!(psplt, Xpp, Ypp, Thpp, l=spscale)
    if drawhist
      push!(psplt.layers,  Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)[1] )#(xbincount=100, ybincount=100))
    end
    psplt
end

function drawLandms(fg::FactorGraph;
              from::Int=0, to::Int=99999999, minnei::Int=0,
              meanmax=:max,
              lbls=true,showmm=false,drawhist=true,
              c="red",
              MM=Union{},
              api::DataLayerAPI=IncrementalInference.localapi  )
    #Gadfly.set_default_plot_size(20cm, 30cm)
    Xp,Yp = get2DLandmSamples(fg, from=from, to=to)
    Xpp = Float64[]; Ypp=Float64[]; Thpp=Float64[]; lblstags=String[];
    if meanmax==:mean
      Xpp,Ypp, t, lbltags = get2DLandmMeans(fg, from=from, to=to, api=api)
    elseif meanmax==:max
      Xpp,Ypp, t, lbltags = get2DLandmMax(fg, from=from, to=to,showmm=showmm,MM=MM, api=api)
    end

    if lbls
      psplt = Gadfly.plot(
      Gadfly.layer(x=Xpp,y=Ypp, label=lbltags, Geom.point, Theme(line_width=1pt, default_color=parse(Colorant,c), point_size=1pt), Geom.label)
      # ,Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)#(xbincount=100, ybincount=100)
      )
    else
      psplt = Gadfly.plot(
      Gadfly.layer(x=Xpp,y=Ypp, Geom.point, Theme(line_width=1pt, default_color=parse(Colorant,c), point_size=1pt))
      )
    end

    if drawhist
      push!(psplt.layers, Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)[1])#(xbincount=100, ybincount=100)
    end

    psplt
end

function drawPosesLandms(fg::FactorGraph;
                    from::Int=0, to::Int=99999999, minnei::Int=0,
                    meanmax=:max,lbls=true,drawhist=true, MM=Union{}, showmm=true,
                    spscale::Float64=5.0,
                    api::DataLayerAPI=IncrementalInference.localapi  )
  p = drawPoses(fg, from=from,to=to,meanmax=meanmax,lbls=lbls,drawhist=drawhist, spscale=spscale, api=api)
  pl = drawLandms(fg, from=from, to=to, minnei=minnei,lbls=lbls,drawhist=drawhist, MM=MM, showmm=showmm, api=api)
  for l in pl.layers
    push!(p.layers, l)
  end
  return p
end

function drawSubmaps(fgl::FactorGraph, fromto::Array{Int,2};
                    m1hist=false,m2hist=false,m3hist=false, showmm=false, MM=Union{},
                    api::DataLayerAPI=IncrementalInference.localapi )
  p = drawLandms(fgl, from=fromto[1,1], to=fromto[1,2], drawhist=m1hist, showmm=showmm, MM=MM, api=api)
  if size(fromto,1) >1
    p2 = drawLandms(fgl, from=fromto[2,1], to=fromto[2,2], drawhist=m2hist,c="blue", showmm=showmm, MM=MM, api=api)
    for l in p2.layers
      push!(p.layers, l)
    end
  end
  if size(fromto,1) >2
    p3 = drawLandms(fgl, from=fromto[3,1], to=fromto[3,2], drawhist=m3hist,c="magenta", showmm=showmm, MM=MM, api=api)
    for l in p3.layers
      push!(p.layers, l)
    end
  end
  return p
end

function drawSubmaps(fgl::FactorGraph, fromto::Array{Int,1}; spread::Int=25,
                    m1hist=false,m2hist=false,m3hist=false, showmm=false, MM=Union{})
  ft = zeros(Int,length(fromto),2)
  for i in 1:length(fromto)
    ft[i,1] = fromto[i]-spread; ft[i,2] = fromto[i]+spread;
  end
  drawSubmaps(fgl, ft, m1hist=m1hist, m2hist=m2hist, m3hist=m3hist, showmm=showmm, MM=MM)
end

# function getKDEMax(p::BallTreeDensity;N=200)
#   m = zeros(p.bt.dims)
#   for i in 1:p.bt.dims
#     mm = marginal(p,[i])
#     rangeV = getKDERange(mm)
#     X = linspace(rangeV[1],rangeV[2],N)
#     yV = evaluateDualTree(mm,X)
#     m[i] = X[findfirst(yV,maximum(yV))]
#   end
#   return m
# end

function investigatePoseKDE(p::BallTreeDensity, p0::BallTreeDensity)
    # co = ["black"; "blue"]
    # h = Union{}
    # x = plotKDE([marginal(p,[1]); marginal(p0,[1])], c=co )
    # y = plotKDE([marginal(p,[2]); marginal(p0,[2])], c=co )
    # if p.bt.dims >= 3
    #   th = plotKDE([marginal(p,[3]); marginal(p0,[3])], c=co )
    #   h = hstack(x,y,th)
    # else
    #   h = hstack(x,y)
    # end
    #
    # return h
    return investigateMultidimKDE(p, p0)
end

function investigatePoseKDE(p::Array{BallTreeDensity,1})
    # co = ["black"; "blue"; "green"; "red"; "magenta"; "cyan"; "cyan1"; "cyan2";
    # "magenta"; "cyan"; "cyan1"; "cyan2"; "magenta"; "cyan"; "cyan1"; "cyan2"; "magenta";
    # "cyan"; "cyan1"; "cyan2"; "magenta"; "cyan"; "cyan1"; "cyan2"]
    # # compute all the marginals
    # Pm = Array{Array{BallTreeDensity,1},1}()
    # push!(Pm,stackMarginals(p,1)) #[marginal(p[1],[1]); marginal(p[2],[1])]
    # push!(Pm,stackMarginals(p,2)) #[marginal(p[1],[2]); marginal(p[2],[2])]
    #
    # h = Union{}
    # x = plotKDE(Pm[1], c=co )
    # y = plotKDE(Pm[2], c=co )
    # if p[1].bt.dims >= 3
    #   #Pm3 = [marginal(p[1],[3]); marginal(p[2],[3])]
    #   push!(Pm,stackMarginals(p,3)) # [marginal(p[1],[3]); marginal(p[2],[3])]
    #   th = plotKDE(Pm[3], c=co )
    #   h = hstack(x,y,th)
    # else
    #   h = hstack(x,y)
    # end
    # return h
    return investigateMultidimKDE(p)
end

function investigatePoseKDE(p::BallTreeDensity)
    # x = plotKDE(marginal(p,[1]) )
    # y = plotKDE(marginal(p,[2]) )
    # if p.bt.dims >= 3
    #   th = plotKDE(marginal(p,[3]) )
    #   return hstack(x,y,th)
    # end
    # return hstack(x,y)
    return investigateMultidimKDE(p)
end

function drawMarginalContour(fgl::FactorGraph, lbl::String;
    xmin=-150,xmax=150,ymin=-150,ymax=150,n=200,
    api::DataLayerAPI=IncrementalInference.localapi )
  #
  p = getVertKDE(fgl,lbl, api=api)  # p = getKDE(getVert(fgl,lbl))
  Gadfly.plot(z=(x,y)->evaluateDualTree(p,vectoarr2([x,y]))[1],
    x=linspace(xmin,xmax,n),
    y=linspace(ymin,ymax,n),
    Geom.contour,
    Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
    Guide.title(lbl)
  )
end

function accumulateMarginalContours(fgl, order;
    xmin=-150,xmax=150,ymin=-150,ymax=150,n=200,
    api::DataLayerAPI=IncrementalInference.localapi )
  #
  pl = drawMarginalContour(fgl, order[1],xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,n=n, api=api)
  pl2 = nothing
  PL = []
  for or in order[1:end]
    pl2 = drawMarginalContour(fgl, or, xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,n=n, api=api)
    push!(PL, pl2)
    push!(pl.layers, pl2.layers[1])
  end
  return pl, PL
end

function plotPose3Pairs(fgl::FactorGraph, sym::Symbol; fill::Bool=true)
  p1= plotKDE(fgl, :x1, dims=[1;2], fill=fill)
  p2 = plotKDE(fgl, :x1, dims=[6;3], fill=fill)
  p3 = plotKDE(fgl, :x1, dims=[4;5], fill=fill)
  Gadfly.draw(PDF("/tmp/RoMEvstackPose3.pdf",15cm, 20cm), vstack(p1,p2,p3) )
  @async run(`evince /tmp/RoMEvstackPose3.pdf`)
  nothing
end

#


text:
import Base.push!

mutable struct HistoryProvider
    history::Array{String, 1}
    filename::String
    cur_idx::Int
    search_results::Array{Int, 1}
    prefix::String
    idx_search::Int
    worker_idx::Int

    HistoryProvider(w_idx::Int) = new(String[""], nothing, 0, 0, [], "", 1, widx)
    HistoryProvider(h::Array{String, 1}, hf, cidx::Int, w_idx::Int) = new(h, hf, cidx, [], "", 1, w_idx)
end

function setup_history(w_idx::Int)
    #load history, etc
    h = HistoryProvider(String["x = pi"],
        joinpath(HOMEDIR, "..", "config", "history_" * string(w_idx)), 1, w_idx)

    if isfile(h.filename)
        h.history = parse_history(h)
        h.cur_idx = length(h.history)+1
    else
        f = open(h.filename, "w")
        close(f)
    end
    return h
end

function setup_history(name)
    #load history, etc
    h = HistoryProvider(String[""], joinpath(HOMEDIR, "..", "config", "history_" * name), 1, 1)

    if isfile(h.filename)
        h.history = parse_history(h)
        h.cur_idx = length(h.history)+1
    else
        f = open(h.filename, "w")
        close(f)
    end
    return h
end

function push!(h::HistoryProvider, str)
    isempty(strip(str)) && return
    push!(h.history, str)

    if isfile(h.filename)
        f = open(h.filename, "a")
        str = "
# _history_entry_
$str"
        write(f, str)
        close(f)
    else
        write(console, "unable to open history file " * h.filename)
    end
end
## use JSON ?
function parse_history(h::HistoryProvider)

    f = open(h.filename, "r")
    lines = readlines(f)
    close(f)

    out = Array{String}(undef, 0)
    current_command = ""
    for line in lines
        if match(r"^# _history_entry_", line) != nothing
            current_command != "" && push!(out, current_command)
            current_command = ""
        else
            current_command = string(current_command, line)
        end
    end
    current_command != "" && push!(out, current_command)
    return out
end

function search(h::HistoryProvider, prefix)

    idx = Array{Int}(undef, 0)
    for i = length(h.history):-1:1
        startswith(h.history[i], prefix) && push!(idx, i)
    end
    h.search_results = idx
    h.prefix = prefix

    return idx
end

function history_up(h::HistoryProvider, prefix, cmd)

    if length(prefix) > 0

        if prefix != h.prefix #new search

            results = search(h, prefix)
            isempty(results) && return false

            h.idx_search = 1
            h.cur_idx = results[h.idx_search]

        else  #we already searched but want to see next result

            h.idx_search = min(h.idx_search+1, length(h.search_results))
            h.cur_idx =  length(h.search_results) > 0 ? h.search_results[h.idx_search] : h.cur_idx
        end
    else
        cmd == "" && seek_end(h) #we go back to the end of the list
        history_move(h, -1)
    end
    return true
end

function history_down(h::HistoryProvider, prefix, cmd)

        if length(prefix) > 0
            h.idx_search = h.idx_search-1
            if h.idx_search <= 0
                seek_end(h)
                return
            end
            h.cur_idx =  length(h.search_results) > 0 ? h.search_results[h.idx_search] : h.cur_idx
        else
            history_move(h, +1)
        end
end

##
function history_move(h::HistoryProvider, m::Int)
    h.cur_idx = clamp(h.cur_idx+m, 1, length(h.history)+1) #+1 is the empty state when we are at the end of history and press down
end

history_get_current(h::HistoryProvider) = h.cur_idx == length(h.history)+1 ? "" : h.history[h.cur_idx]

function seek_end(h::HistoryProvider)
    h.cur_idx = length(h.history)+1
end


text:
using WordCloud
using Random

background = loadmask(pkgdir(WordCloud)*"/res/butterfly.png")
istrans = c->maximum(c[1:3])*(c[4]/255)<128
mask = WordCloud.imagemask(background, istrans)
showmask(background, mask, highlight=(1,0,0,0.7))
#md# `showmask` might be helpful to find a proper `istrans` function
words = [randstring(1) for i in 1:600]
weights = randexp(length(words)) .+ 1

wc = wordcloud(
    words, weights,
    mask = background,
    colors = "LimeGreen",
    angles = -30,
    density = 0.45,
    transparent = istrans,
    spacing = 1,
) |> generate!;
#md# ## average style
recolor!(wc, style=:average)
avgimg = paint(wc, background=loadmask(background, color=0.99))
#md# ## clipping style
recolor!(wc, style=:clipping)
clipimg = paint(wc, background=loadmask(background, color=0.99))
#md# ## blending style
recolor!(wc, style=:reset)
recolor!(wc, style=:blending, alpha=0.5) #blending with origin color - LimeGreen
blendimg = paint(wc, background=loadmask(background, color=0.99))
#md# ## mix style
#md# styles can also be mixed
# setcolors!(wc, :, "LimeGreen")
recolor!(wc, style=:reset)
recolor!(wc, 1:3:length(words), style=:average) #vector index is ok
recolor!(wc, 2:3:length(words), style=:clipping)
recolor!(wc, 3:3:length(words), style=:blending)
setcolors!(wc, 200:250, "black")
recolor!(wc, 200:250, style=:reset)
setcolors!(wc, 1, "black")
recolor!(wc, 1, style=:reset) #single index is ok
mixstyleimg = paint(wc, background=loadmask(background, color=0.99))
#md#
h, w = size(avgimg)
lw = 21
lc = eltype(avgimg)(parsecolor(0.1))
vbar = zeros(eltype(avgimg), (h, lw))
hbar = zeros(eltype(avgimg), (lw, 2w+lw))
vbar[:, lw÷2+1] .= lc
hbar[lw÷2+1, :] .= lc
image = [avgimg vbar clipimg; hbar; blendimg vbar mixstyleimg]
println("results are saved to recolor.png")
WordCloud.save("recolor.png", image)
wc
#eval# runexample(:recolor)
#md# ![recolor](recolor.png)


text:
d = Uniform([1])

@test rand(d) == 1
@test rand(MersenneTwister(4), d) == 1
@test collect(support(d)) == [1]
@test Random.gentype(d) == typeof(1)
@test Random.gentype(typeof(d)) == typeof(1)
@test pdf(d, 0) == 0.0
@test pdf(d, 1) == 1.0
@test mode(d) == 1
@test mean(d) == 1
@test typeof(mean(d)) == typeof(mean([1]))
@test collect(weighted_iterator(d)) == [1=>1.0]

@test sprint((io,d)->show(io,MIME("text/plain"),d), d) == sprint((io,d)->showdistribution(io,d,title="Uniform distribution"), d)

d2 = Uniform((:symbol,))
@test rand(d2) == :symbol
@test rand(MersenneTwister(4), d2) == :symbol
@test collect(support(d2)) == [:symbol]
@test Random.gentype(d2) == typeof(:symbol)
@test Random.gentype(typeof(d2)) == typeof(:symbol)
@test pdf(d2, :another) == 0.0
@test pdf(d2, :symbol) == 1.0
@test mode(d2) == :symbol
@test collect(weighted_iterator(d2)) == [:symbol=>1.0]

# uniqueness test
@test_throws ErrorException Uniform((:symbol, :symbol))

d3 = UnsafeUniform([1])

@test rand(d3) == 1
@test rand(MersenneTwister(4), d3) == 1
@test collect(support(d3)) == [1]
@test Random.gentype(d3) == typeof(1)
@test Random.gentype(typeof(d3)) == typeof(1)
@test pdf(d3, 1) == 1.0
@test mean(d3) == 1
@test mode(d3) == 1
@test typeof(mean(d3)) == typeof(mean([1]))
@test collect(weighted_iterator(d3)) == [1=>1.0]

@test sprint((io,d3)->show(io,MIME("text/plain"),d3), d3) == sprint((io,d3)->showdistribution(io,d3,title="UnsafeUniform distribution"), d3)

d4 = UnsafeUniform((:symbol,))
@test rand(d4) == :symbol
@test rand(MersenneTwister(4), d4) == :symbol
@test collect(support(d4)) == [:symbol]
@test Random.gentype(d4) == typeof(:symbol)
@test Random.gentype(typeof(d4)) == typeof(:symbol)
# @test pdf(d4, :another) == 0.0 # this will not work
@test pdf(d4, :symbol) == 1.0
@test mode(d4) == :symbol
@test collect(weighted_iterator(d4)) == [:symbol=>1.0]


text:
# BSE kernel for the u channel
function compute_u_BSE!(
    Λ    :: Float64,
    buff :: Matrix{Float64},
    v    :: Float64,
    dv   :: Float64,
    u    :: Float64,
    vu   :: Float64,
    vup  :: Float64,
    r    :: Reduced_lattice,
    m    :: Mesh,
    a    :: Action_su2,
    temp :: Array{Float64, 3}
    )    :: Nothing

    # get propagator
    p = -get_propagator(Λ, v - 0.5 * u, v + 0.5 * u, m, a)

    # get buffers for left vertex
    bs1 = get_buffer_s(v + vu, 0.5 * (u - v + vu), 0.5 * (-u - v + vu), m)
    bt1 = get_buffer_t(v - vu, 0.5 * (u + v + vu), 0.5 * (-u + v + vu), m)
    bu1 = get_buffer_empty()

    # get buffers for right vertex
    bs2 = get_buffer_s( v + vup, 0.5 * (u + v - vup), 0.5 * (-u + v - vup), m)
    bt2 = get_buffer_t(-v + vup, 0.5 * (u + v + vup), 0.5 * (-u + v + vup), m)
    bu2 = get_buffer_u(u, v, vup, m)

    # cache vertex values for all lattice sites in temporary buffer
    get_Γ_avx!(r, bs1, bt1, bu1, a, temp, 1, ch_u = false)
    get_Γ_avx!(r, bs2, bt2, bu2, a, temp, 2)

    # compute contributions for all lattice sites
    @turbo unroll = 1 for i in eachindex(r.sites)
        # read cached values for site i
        v1s = temp[i, 1, 1]; v1d = temp[i, 2, 1]
        v2s = temp[i, 1, 2]; v2d = temp[i, 2, 2]

        # compute contribution at site i
        Γs = -p * (2.0 * v1s * v2s + v1s * v2d + v1d * v2s)
        Γd = -p * (3.0 * v1s * v2s + v1d * v2d)

        # parse result to output buffer
        buff[1, i] += dv * Γs
        buff[2, i] += dv * Γd
    end

    return nothing
end


text:
# POMDP model functions

abstract POMDP

abstract State
abstract Action
abstract Observation
typealias Reward Float64

abstract AbstractDistribution

# return the space sizes
@pomdp_func n_states(pomdp::POMDP)
@pomdp_func n_actions(pomdp::POMDP)
@pomdp_func n_observations(pomdp::POMDP)

# return the discount factor
@pomdp_func discount(pomdp::POMDP)

@pomdp_func transition(pomdp::POMDP, state::State, action::Action, distribution::AbstractDistribution=create_transition_distribution(pomdp))
@pomdp_func observation(pomdp::POMDP, state::State, action::Action, statep::State, distribution::AbstractDistribution=create_observation_distribution(pomdp))
@pomdp_func observation(pomdp::POMDP, state::State, action::Action, distribution::AbstractDistribution=create_observation_distribution(pomdp))
@pomdp_func reward(pomdp::POMDP, state::State, action::Action, statep::State)

@pomdp_func create_state(pomdp::POMDP)
@pomdp_func create_observation(pomdp::POMDP)

@pomdp_func isterminal(pomdp::POMDP, state::State) = false
@pomdp_func isterminal(pomdp::POMDP, observation::Observation) = false

@pomdp_func index(pomdp::POMDP, state::State)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
using LinearAlgebra

#
# KF
#

s = Statevector([0, 0, 0], Matrix(1.0I, 3, 3))
KalmanFilter.update!(s, [1, 1, 1], Matrix(1.0I, 3, 3), Matrix(1.0I, 3, 3))
@test s.xhat == [.5, .5, .5]
@test s.P == Matrix(0.5I, 3, 3)
KalmanFilter.predict!(s, [1 1 0;0 1 1;0 0 1], Matrix(0.5I, 3, 3))
@test s.xhat == [1, 1, .5]

#
# EKF
#

s = Statevector([0, 0, 0], Matrix(1.0I, 3, 3))
KalmanFilter.update!(s, [1, 1, 1], x -> [x[1], x[2], x[3]], x -> Matrix(1.0I, 3, 3), Matrix(1.0I, 3, 3))
@test s.xhat == [.5, .5, .5]
@test s.P == Matrix(0.5I, 3, 3)
KalmanFilter.predict!(s, x -> [1, 1, .5], x -> [1 1 0;0 1 1;0 0 1], Matrix(0.5I, 3, 3))
@test s.xhat == [1, 1, .5]

#
# EKF, numeric Jacobian
#

s = Statevector([0., 0., 0.], Matrix(1.0I, 3, 3))
KalmanFilter.update!(s, [1., 1., 1.], x -> [x[1], x[2], x[3]], Matrix(1.0I, 3, 3))
@test s.xhat == [.5, .5, .5]
@test s.P == Matrix(0.5I, 3, 3)
KalmanFilter.predict!(s, x -> [1, 1, .5], Matrix(0.5I, 3, 3))
@test s.xhat == [1, 1, .5]


text:
import FRS_Household_Getter: initialise, get_household
import Test: @testset, @test
import Model_Household: Household, Person, People_Dict, BUAllocation,
      PeopleArray, printpids,
      BenefitUnit, BenefitUnits, default_bu_allocation,
      get_benefit_units, get_head, get_spouse, num_people

start_year=2015

@testset "benefit unit allocations" begin
      rc = @timed begin
         num_households,total_num_people,nhh2 = initialise(
            household_name = "model_households_scotland",
            people_name    = "model_people_scotland",
            start_year = start_year )
      end
      println( "num_households=$num_households, num_people=$num_people")
      people_count = 0
      for hhno in 1:num_households
            hh = get_household( hhno )
            # println("people in HH")
            # printpids(hh.people)

            hhsize = size(collect(keys(hh.people)))[1]
            people_count += hhsize
            bus = get_benefit_units( hh )
            buallocation = default_bu_allocation( hh )
            nbus = size( bus )[1]
            @test 0 < nbus < 10
            by_bu_people_count = 0
            i = 0
            for bu in bus
                  # println( "from BU")
                  # printpids(bu)

                  i += 1
                  bu_people_count = 1
                  head = get_head( bu )
                  @test head.age >= 16
                  spouse = get_spouse( bu )
                  if spouse != nothing
                        bu_people_count += 1
                        @test spouse.age >= 16
                  end
                  for chno in bu.children
                        child = bu.people[chno]
                        @test child.age <= 19
                        bu_people_count += 1
                  end
                  by_bu_people_count += bu_people_count
                  @test bu_people_count == num_people( bu )
                  @test bu_people_count == size( buallocation[i])[1]
            end
            @test by_bu_people_count == hhsize

            bua_people_count = 0
            for bua in buallocation
                  bua_people_count += size( bua )[1]
            end
            # println("BU Allocation")
            # printpids(buallocation)
            @test bua_people_count == hhsize
      end
      @test people_count == total_num_people
end


text:
module IntertemporalChoiceHeuristics
    import DataFrames, Distributions, ForwardDiff, Optim, StatsBase

    include("models/base_model.jl")
    include("losses/base_loss.jl")

    include("links.jl")
    include("data_io.jl")
    include("weighting.jl")
    include("fitting.jl")
    include("average_loss.jl")
    include("simulate.jl")
    include("cov.jl")
    include("stderrs.jl")

    include("models/heuristics/baseline.jl")
    include("models/heuristics/ITCH.jl")
    include("models/heuristics/DRIFT.jl")
    include("models/heuristics/tradeoff.jl")

    include("models/exponential/standard.jl")
    include("models/exponential/homothetic.jl")
    include("models/exponential/intercept.jl")

    include("models/hyperbolic/standard.jl")
    include("models/hyperbolic/homothetic.jl")
    include("models/hyperbolic/intercept.jl")

    include("models/hyperboloid/standard.jl")
    include("models/hyperboloid/homothetic.jl")
    include("models/hyperboloid/intercept.jl")

    include("models/quasihyperbolic/standard.jl")
    include("models/quasihyperbolic/homothetic.jl")
    include("models/quasihyperbolic/intercept.jl")

    include("models/system2/standard.jl")
    include("models/system2/homothetic.jl")
    include("models/system2/intercept.jl")

    include("losses/absolute_error.jl")
    include("losses/extremity.jl")
    include("losses/log_loss.jl")
    include("losses/squared_error.jl")
    include("losses/zero_one_loss.jl")

    models = (
        Baseline,
        ITCH,
        DRIFT,
        Tradeoff,
        Exponential,
        HomotheticExponential,
        ExponentialIntercept,
        Hyperbolic,
        HomotheticHyperbolic,
        HyperbolicIntercept,
        Hyperboloid,
        HomotheticHyperboloid,
        HyperboloidIntercept,
        QuasiHyperbolic,
        HomotheticQuasiHyperbolic,
        QuasiHyperbolicIntercept,
        System2,
        HomotheticSystem2,
        System2Intercept,
    )

    training_losses = (
        SquaredError(),
        LogLoss(),
    )

    evaluation_losses = (
        AbsoluteError(),
        SquaredError(),
        LogLoss(),
        ZeroOneLoss(),
        Extremity(),
    )
end


text:
using Onda
using Documenter

makedocs(modules=[Onda],
         sitename="Onda",
         authors="Beacon Biosignals, Inc.",
         pages=["API Documentation" => "index.md",
                "Upgrading From Older Versions Of Onda" => "upgrading.md"])

deploydocs(repo="github.com/beacon-biosignals/Onda.jl.git", push_preview=true)


text:
mutable struct WirePanel
    points::Dict{Tuple{Int,Int},Char}
    distances::Dict{Tuple{Int,Int},Int}
    left::Int
    top::Int
    right::Int
    bottom::Int
end
WirePanel() = WirePanel(Dict{Tuple{Int,Int},Char}(), Dict{Tuple{Int,Int},Int}(), 0, 0, 0, 0)

function WirePanel(directions::AbstractString)
    wirepanel = WirePanel()
    pos = (0, 0)
    dist = 0
    for m in split(directions, ',')
        n = parse(Int, m[2:end])
        islr = occursin(m[1], "LR")
        r = if islr
            m[1] == 'L' ? (pos[1]-1:-1:pos[1]-n) : (pos[1]+1:pos[1]+n)
        else
            m[1] == 'U' ? (pos[2]-1:-1:pos[2]-n) : (pos[2]+1:pos[2]+n)
        end
        for k in r
            pos = islr ? (k, pos[2]) : (pos[1], k)
            wirepanel[pos] = if wirepanel[pos] != '·' || k == last(r)
                '╋'
            elseif islr
                '━'
            else
                '┃'
            end
            dist += 1
            wirepanel.distances[pos] = dist
        end
    end
    wirepanel
end

function Base.empty!(w::WirePanel)
    empty!(w.points)
    w.left = w.right = w.top = w.bottom = 0
    w
end

function Base.getindex(wirepanel::WirePanel, t::Tuple{Int,Int})
    if iszero(t[1]) && iszero(t[2])
        'o'
    else
        get(wirepanel.points, t, '·')
    end
end

function Base.setindex!(panel::WirePanel, c::AbstractChar, t::Tuple{Int,Int})
    panel.points[t] = c
    x, y = t
    panel.left = min(x, panel.left)
    panel.right = max(x, panel.right)
    panel.top = min(y, panel.top)
    panel.bottom = max(y, panel.bottom)
end

function Base.show(io::IO, w::WirePanel)
    println(io, abs(w.right - w.left) + 1, '×', abs(w.top - w.bottom) + 1, " WirePanel")
    spareline = "  " * ("·" ^ (w.right - w.left + 3))
    println(io, spareline)
    for y in w.top:w.bottom
        print(io, "  ·")
        for x in w.left:w.right
            print(io, w[(x, y)])
        end
        println(io, '·')
    end
    print(io, spareline)
end

loadpanels(filename::AbstractString) = WirePanel.(eachline(filename))

manhattandistance(t::Tuple{Int,Int}) = abs(t[1]) + abs(t[2])

"Returns the Manhattan distance to the closest short to the origin."
function closestshort(panels::Vararg{WirePanel})
    minimum(manhattandistance(k) for k in shorts(panels))
end

function shorts(panels::Vector{WirePanel})
    intersect((keys(panel.points) for panel in panels)...)
end

function part1(filename::AbstractString)
    closestshort(loadpanels(filename)...)
end

Base.in(t::Tuple, panel::WirePanel) = t in keys(panel.points)

function paneldistance(panels, point::Tuple{Int,Int})
    sum(panel.distances[point] for panel in panels)
end

function part2(filename::AbstractString)
    panels = loadpanels(filename)
    minimum(paneldistance(panels, short) for short in shorts(panels))
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# O(N^2) solution

n = 200
x = [parse(Int, readline()) for _ in 1:n]

comb = Set(x)
for a in x
    for b in x
        if 2020 - a - b ∈ comb
            println(a * b * (2020 - a - b))
            exit()
        end
    end
end


text:
function probOfBeingChosenAtLeastOnce(PP, freq, tps)
	#% PP(i,g) is the probability of type i touching page 'g'
	#%       note that sum(i,:) does not need to be one at all!
	#% tps(i) is the number of transactions of type i.
	#%
	#% T(g) is the probability that at least one of the tps transactions touch g

	effectiveTPS = freq;

	for i=1:size(effectiveTPS,1)
	    effectiveTPS[i,:] = effectiveTPS[i,:] * tps[i];
	end

	#%PP(1,:) = PP(1,:) / 50;
	PP = PP / 70;
	for i=1:length(PP)
		if PP[i] > 1
			PP[i] = 1
		end
	end

	#%P(i,p)= the probability of page 'p' staying CLEAN from a tran type 'i'
	#%after effectiveTPS recurrence of tran type 'i'
	P = (1-PP).^effectiveTPS;

	#%first let's calculate the probability that none of the tps transactions
	#%touch a page
	T = prod(P,1);

	#%now 1-T is the probability that at least one of the transactions touch them!
	T = 1-T;
	T
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
"""
    test trajectory optimization
"""
prob_traj = prob.prob

# objective
for i = 1:1
    z0 = rand(prob_traj.num_var)
    tmp_o(z) = eval_objective(prob_traj, z)
    ∇j = zeros(prob_traj.num_var)
    eval_objective_gradient!(∇j, z0, prob_traj)
    @assert norm(ForwardDiff.gradient(tmp_o, z0) - ∇j, Inf) < 1.0e-5
    println("i $i")
end

# constraints
c0 = zeros(prob_traj.num_con)
eval_constraint!(c0, z0, prob_traj)
tmp_c(c, z) = eval_constraint!(c, z, prob_traj)
∇c_fd = ForwardDiff.jacobian(tmp_c, c0, z0)
spar = sparsity_jacobian(prob_traj)
∇c_vec = zeros(length(spar))
∇c = zeros(prob_traj.num_con, prob_traj.num_var)
eval_constraint_jacobian!(∇c_vec, z0, prob_traj)
for (i,k) in enumerate(spar)
    ∇c[k[1],k[2]] = ∇c_vec[i]
end
@assert norm(vec(∇c) - vec(∇c_fd)) < 1.0e-10
@assert sum(∇c) - sum(∇c_fd) < 1.0e-10
norm(vec(∇c))# -


text:
# Autogenerated wrapper script for Pango_jll for x86_64-apple-darwin
export libpango, libpangocairo, libpangoft

using FriBidi_jll
using FreeType2_jll
using Glib_jll
using Fontconfig_jll
using HarfBuzz_jll
using Cairo_jll
JLLWrappers.@generate_wrapper_header("Pango")
JLLWrappers.@declare_library_product(libpango, "@rpath/libpango-1.0.0.dylib")
JLLWrappers.@declare_library_product(libpangocairo, "@rpath/libpangocairo-1.0.0.dylib")
JLLWrappers.@declare_library_product(libpangoft, "@rpath/libpangoft2-1.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(FriBidi_jll, FreeType2_jll, Glib_jll, Fontconfig_jll, HarfBuzz_jll, Cairo_jll)
    JLLWrappers.@init_library_product(
        libpango,
        "lib/libpango-1.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpangocairo,
        "lib/libpangocairo-1.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpangoft,
        "lib/libpangoft2-1.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Base: @pure

@pure function print_list(values)
  argument()  = PrintArgument()
  to_value(x) = PrintValue{Symbol(String(x))}()

  if length(values) == 1
    to_value(values[1])
  else
    if values[1] == ""
      PrintBoth(
        argument(),
        print_list(values[2:end])
      )
    else
      PrintBoth(
        to_value(values[1]),
        print_list([""; values[2:end]])
      )
    end
  end
end


text:
import ..UncertainDatasets:
    AbstractUncertainIndexValueDataset

import ..SamplingConstraints: SamplingConstraint
"""
    ConstrainedIndexValueResampling(constraints::NTuple{N_DATASETS, NTuple{N_VARIABLES, Union{SamplingConstraint, Vector{<:SamplingConstraint}}}}, n::Int)

Indicates that resampling should be performed with constraints on a set of uncertain index-value datasets.
See examples for usage.

## Fields

- `constraints`. The constraints for the datasets. The constraints are represented as a tuple of length
    `N_DATASETS`, where the `i`-th tuple element is itself a `N_VARIABLES`-length tuple containing the
    constraints for the `N_VARIABLES` different variables. See "Indexing" below for details.
    Constraints for each individual variable must be supplied as either a single sampling constraint,
    or as a vector of sampling constraints with length matching the length of the variable
    (`Union{SamplingConstraint, Vector{<:SamplingConstraint}}}`). For example, if the `j`-th variable
    for the `i`-th dataset contains 352 observations, then `constraints[i, j]` must be either a single
    sampling constraint (e.g. `TruncateStd(1.1)`) or a vector of 352 different sampling constraints
    (e.g. `[TruncateStd(1.0 + rand()) for i = 1:352]`).
- `n::Int`. The number of draws.

## Indexing

Assume `c` is an instance of `ConstrainedIndexValueResampling`. Then

- `c[i]` returns the `NTuple` of constraints for the `i`-th dataset, and
- `c[i, j]` returns the constraint(s) for the `j`-th variable of the `i`-th dataset.

## Example

### Defining `ConstrainedIndexValueResampling`s.

Assume we want to constraints three separate uncertain index-value datasets, with different
sampling constraints for the indices and the values for each of the datasets.

```julia
# (index constraints, value constraints) for the 1st, 2nd and 3rd datasets
c1 = (TruncateStd(1), TruncateStd(1.1))
c2 = (TruncateStd(0.5), TruncateQuantiles(0.1, 0.8))
c3 = (TruncateQuantiles(0.05, 0.95), TruncateQuantiles(0.33, 0.67))
c = ConstrainedIndexValueResampling(c1, c2, c3)
```

Now,

- `c[2]` returns the `NTuple` of constraints for the 2nd dataset, and
- `c[1, 2]` returns the constraint(s) for the 2nd variable of the 1st dataset.

### Controlling the number of draws

The number of draws defaults to 1 if not specified. To indicate that
more than one draw should be performed, just input the number of draws
before supplying the constraints to the constructor.

```
c1 = (TruncateStd(1), TruncateStd(1.1))
c2 = (TruncateStd(0.5), TruncateQuantiles(0.1, 0.8))

# A single draw
c_single = ConstrainedIndexValueResampling(c1, c2)

# Multiple (300) draws
c_multiple = ConstrainedIndexValueResampling(300, c1, c2)
```

### Detailed example

Let's say we have two uncertain index-value datasets `x` and `y`. We want to constrain
the furnishing distributions/population for both the time indices and values,
both for `x` and `y`. For `x`, truncate the indices at `0.8` times the standard
deviation around their mean, and for `y`, trucate the indices at `1.5` times the
standard deviation around their mean. Next, truncate `x`s values at roughly
(roughly) at their 20th percentile range, and truncate `y`s values at roughly
their 80th percentile range.

All this information can be combined in a  `ConstrainedIndexValueResampling` instance.
This instance can be passed on to any function that accepts uncertain index-value datasets,
to indicate that resampling should be performed on truncated versions of the
distributions/populations furnishing the datasets.

```julia
# some noise, so we don't truncate all furnishing distributions/population at
# exactly the same quantiles.
r = Uniform(0, 0.01)

constraints_x_inds = TruncateStd(0.8)
constraints_y_inds = TruncateStd(1.5)
constraints_x_vals = [TruncateQuantiles(0.4 + rand(r), 0.6 + rand(r)) for i = 1:length(x)];
constraints_y_vals = [TruncateQuantiles(0.1 + rand(r), 0.9 + rand(r)) for i = 1:length(x)];

cs_x = (constraints_x_inds, constraints_x_vals)
cs_y = (constraints_y_inds, constraints_y_vals)

resampling = ConstrainedIndexValueResampling(cs_x, cs_y)
```
"""
struct ConstrainedIndexValueResampling{N_VARIABLES, N_DATASETS} <: AbstractUncertainDataResampling where {N_VARIABLES, N_DATASETS}
    constraints::NTuple{N_DATASETS, NTuple{N_VARIABLES, Union{SamplingConstraint, Vector{<:SamplingConstraint}}}}
    n::Int
end

ConstrainedIndexValueResampling(constraints::NTuple{N_VARIABLES, Union{SamplingConstraint, Vector{<:SamplingConstraint}}}...) where N_VARIABLES =
    ConstrainedIndexValueResampling((constraints...,), 1)

ConstrainedIndexValueResampling(n::Int, constraints::NTuple{N_VARIABLES, Union{SamplingConstraint, Vector{<:SamplingConstraint}}}...) where N_VARIABLES =
    ConstrainedIndexValueResampling((constraints...,), n)

function Base.show(io::IO, constraints::ConstrainedIndexValueResampling{N_VARIABLES, N_DATASETS}) where {N_VARIABLES, N_DATASETS}
    s = "$(typeof(constraints)) for $N_DATASETS set(s) of index-value constraints (each a $N_VARIABLES-tuple) where n=$(constraints.n)"
    println(io, s)
end

Base.eachindex(c::ConstrainedIndexValueResampling) = Base.OneTo(length(c.constraints))
Base.length(c::ConstrainedIndexValueResampling) = length(c.constraints)
Base.firstindex(c::ConstrainedIndexValueResampling) = 1
Base.lastindex(c::ConstrainedIndexValueResampling) = length(c.constraints)
Base.getindex(c::ConstrainedIndexValueResampling, i) = c.constraints[i]
Base.getindex(c::ConstrainedIndexValueResampling, i, j) = c.constraints[i][j]
Base.iterate(c::ConstrainedIndexValueResampling, state = 1) = iterate(c.constraints, state)

export ConstrainedIndexValueResampling