text:
using PowerModels
using .PP2PM

function run_powermodels(json_path)
    pm = PP2PM.load_pm_from_json(json_path)
    model = PP2PM.get_model(pm["pm_model"])
    #
    solver = PP2PM.get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"], pm["pm_tol"])

    if haskey(pm["branch"]["1"],"c_rating_a")
        for (key, value) in pm["gen"]
           # value["pmin"] = 0
           value["pmax"] *= 0.01
           value["qmax"] *= 0.01
           value["qmin"] *= 0.01
           value["pg"] *= 0.01
           value["qg"] *= 0.01
           value["cost"] *= 100
        end

        for (key, value) in pm["branch"]
           value["c_rating_a"] *= 0.01
        end

        for (key, value) in pm["load"]
           value["pd"] *= 0.01
           value["qd"] *= 0.01
        end

        result = PowerModels._run_opf_cl(pm, model, solver,
                                        setting = Dict("output" => Dict("branch_flows" => true)))
    else
        result = PowerModels.run_opf(pm, model, solver,
                                        setting = Dict("output" => Dict("branch_flows" => true)))
    end

    result["termination_status"] = string(result["termination_status"])
    result["dual_status"] = string(result["dual_status"])
    result["primal_status"] = string(result["primal_status"])

    return result
end


text:
module EMLmixer
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("mixer/mixer2.jl")
	include("mixer/mixer3.jl")
	include("mixer/mixer_n.jl")
end

ext: jl


text:
module TokenizeMeta

using Tokenize

export Token
import Tokenize.Tokens.Token

include("Eval.jl")
include("Type.jl")

end


text:
using Pkg

dependencies = [
"JLD",
"LoopVectorization",
"Optim",
"NLsolve",
"StatsBase",
"Distributions",
"Plots",
"Revise",
"DataFrames",
"CSV"
]

Pkg.add(dependencies)

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
function Room1(;
    shape = RoomRectangle(6, 5),
    obstacles = [Rectangle((2.4, 1.9), 1.2, 0.5)],
    entrances = [
        HDoor((1.05, 5.0), 0.6),
        HDoor((2.7, 5.0), 0.6),
        HDoor((4.35, 5.0), 0.6),
    ],
    exits = [HDoor((2.7, 0.0), 0.6)],
    checkpoints = [
        Checkpoint((2.0, 4.0)),
        Checkpoint((3.0, 4.0)),
        Checkpoint((4.0, 4.0)),
        Checkpoint((2.0, 2.4)),
        Checkpoint((2.0, 1.9)),
        Checkpoint((4.0, 2.4)),
        Checkpoint((4.0, 1.9)),
    ],
)

    return Room(shape, obstacles, entrances, exits, checkpoints)
end

function Model1(;
    room = Room1(),
    time::TimeStrategy = Basic(),
    schedule::ScheduleStrategy = RandomScheduler(),
    target::TargetStrategy = Nearest(),
    move::MoveStrategy = JancaMove(),
    kwargs...
)

    return Model(room, time, schedule, target, move, kwargs...)
end

function Room2(;
    shape = RoomRectangle(6, 5),
    obstacles = Obstacle[],
    entrances = [
        HDoor((1.05, 5.0), 0.6),
        HDoor((2.7, 5.0), 0.6),
        HDoor((4.35, 5.0), 0.6),
    ],
    exits = [HDoor((2.7, 0.0), 0.6)],
    checkpoints = Target[],
)

    return Room(shape, obstacles, entrances, exits, checkpoints)
end

function Model2(;
    room = Room2(),
    time::TimeStrategy = Basic(; at = 40),
    schedule::ScheduleStrategy = RandomScheduler(),
    target::TargetStrategy = Nearest(),
    move::MoveStrategy = JancaMove(),
    kwargs...
)

    return Model(room, time, schedule, target, move, kwargs...)
end


text:
function countrecord(df)
    sort(combine(groupby(df, All()), nrow => :__count), :__count)
end


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

name = "LibVPX"
version = v"5.0.0"

# Collection of sources required to build LibVPX
sources = [
    "https://github.com/webmproject/libvpx/archive/v1.7.0/libvpx-1.7.0.tar.gz" =>
    "1fec931eb5c94279ad219a5b6e0202358e94a93a90cfb1603578c326abfc1238",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libvpx-1.7.0/
sed -i 's/cp -p/cp/' build/make/Makefile
mkdir libvpx-build
cd libvpx-build
CROSS=$target
apk add diffutils
../configure --prefix=$prefix --target=generic-gnu --enable-shared --disable-static
echo "SRC_PATH_BARE=.." >> config.mk
echo "target=libs" >> config.mk
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/jpsamaroo/YasmBuilder/releases/download/v1.3.0-pre/build_YasmBuilder.v1.3.0.jl"
]

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

max_stars_repo_name: ianshmean/LibVPXBuilder

max_issues_repo_name: ianshmean/LibVPXBuilder

max_forks_repo_name: ianshmean/LibVPXBuilder


text:
# FASTA Writer
# ============

"""
    FASTA.Writer(output::IO; width=70)

Create a data writer of the FASTA file format.

# Arguments
* `output`: data sink
* `width=70`: wrapping width of sequence characters
"""
struct Writer <: BioGenerics.IO.AbstractWriter
    output::IO
    # maximum sequence width (no limit when width ≤ 0)
    width::Int
end

function BioGenerics.IO.stream(writer::Writer)
    return writer.output
end

function Writer(output::IO; width::Integer = 70)
    return Writer(output, width)
end

function Base.write(writer::Writer, record::Record)
    checkfilled(record)
    output = writer.output
    n::Int = 0
    if writer.width ≤ 0
        n += write(output, record, '\n')
    else
        headerlen = hasdescription(record) ? last(record.description) : last(record.identifier)
        n += unsafe_write(output, pointer(record.data), headerlen)
        n += write(output, '\n')
        p = pointer(record.data, first(record.sequence))
        p_end = pointer(record.data, last(record.sequence))
        while p ≤ p_end
            w = min(writer.width, p_end - p + 1)
            n += unsafe_write(output, p, w)
            n += write(output, '\n')
            p += w
        end
    end
    return n
end


text:
using Lerche

grammar_path = joinpath(dirname(pathof(@__MODULE__)), "grammar.lark")

txt = read(open(grammar_path), String)
const parser = Lark(txt,
    parser="lalr", start="formula");

const parser_equation = Lark(txt,
    parser="lalr", start="equation");

const parser_assignment = Lark(txt,
    parser="lalr", start="assignment");

    const parser_equation_block = Lark(txt,
    parser="lalr", start="equation_block");

const parser_assignment_block = Lark(txt,
    parser="lalr", start="assignment_block");

LTree = Union{Tree, Token}

# sanitize:
# ensure that all variables have a date
# maybe we should convert numbers into actual numbers here ?

struct Sanitizer <: Transformer
    variables::Vector{String}
end

Sanitizer() = Sanitizer(String[])

@rule variable(sanitizer::Sanitizer, children) = length(children)==1 ?
    Tree("variable", Union{Tree,Token}[children[1], Tree("date", [Token("SIGNED_INT2", "0")])]) :
    Tree("variable", children)

@rule symbol(sanitizer::Sanitizer, children) = (children[1].value in sanitizer.variables) ?
    Tree("variable", Union{Tree,Token}[Tree("name", [Token("NAME", children[1].value )]), Tree("date", [Token("SIGNED_INT2", "0")])]) :
    Tree("symbol", children)

sanitize(expr::LTree; variables=AbstractString[])::LTree = Lerche.transform(Sanitizer(variables), expr)

struct TimeShifter <: Transformer
    Δt::Int64
end

@rule variable(ts::TimeShifter, children) = begin
    t0 = parse(Int, children[2].children[1].value)
    t1 = t0 + ts.Δt
    t1_s = string(t1)
    name = children[1].children[1].value
    Tree("variable", Union{Tree,Token}[Tree("name", [Token("NAME", name)]), Tree("date", [Token("SIGNED_INT2", t1_s)])])
end

time_shift(expr::LTree, Δt::Int64)::LTree = Lerche.transform(TimeShifter(Δt), expr)

struct SteadyState <: Transformer
end

@rule variable(t::SteadyState, children) = begin
    name = children[1].children[1].value
    Tree("variable", Union{Tree,Token}[Tree("name", [Token("NAME", name)]), Tree("date", [Token("SIGNED_INT2", "0")])])
end

steady_state(expr::LTree)::LTree = Lerche.transform(SteadyState(), expr)

import Dolang
import Dolang: stringify

## Converts an expression into a an expression, where all variables are replaced
## by their equation name

# struct Stringifier <: Transformer end

# stringify(variable::String, date::Int) = String(stringify(Symbol(variable), date))
# stringify(::typeof(Symbol), variable::String, date::Int) = (stringify(Symbol(variable), date))

# stringify(variable::String) = String(stringify(Symbol(variable)))
# stringify(::typeof(Symbol), variable::String) = (stringify(Symbol(variable)))

# @rule variable(t::Stringifier, tree) =  stringify(Symbol, tree[1].children[1].value, parse(Int,tree[2].children[1].value))
# @rule symbol(t::Stringifier, tree) =  stringify(Symbol,tree[1].value)

# stringify(expr::LTree) = Lerche.transform(Stringifier(), expr)

# converts to julia expression tree
struct Converter <: Transformer
    stringify::Bool
end

Converter() = Converter(false)

@rule variable(t::Converter, tree) = begin

    s = Symbol(tree[1].children[1].value)
    d = parse(Int,tree[2].children[1].value)
    if t.stringify
        stringify(s, d)
    else
        :($(s)[$(d)])
    end
end
@rule symbol(t::Converter, tree) = begin
    s = Symbol(tree[1].value)
    if s == :inf
        :(Inf)
    elseif t.stringify
        stringify(s)
    else
        s
    end
end

@rule call(t::Converter, children) =  Expr(:call, Symbol(children[1].value), children[2:end]...)
@rule add(t::Converter, children) =  Expr(:call, :+, children...)
@rule mul(t::Converter, children) =  Expr(:call, :*, children...)
@rule div(t::Converter, children) =  Expr(:call, :/, children...)
@rule pow(t::Converter, children) =  Expr(:call, :^, children...)
@rule neg(t::Converter, children) =  Expr(:call, :-, children...)
@rule sub(t::Converter, children) =  Expr(:call, :-, children...)
@rule number(t::Converter, children) =  parse(Float64, children[1].value)
@rule equality(t::Converter, children) =  :($(children[1]) == $(children[2]))

# this is not supposed to work for equation blocks.
convert(::typeof(Expr), expr::Tree; stringify=true) =  Lerche.transform(Converter(stringify), expr)
convert(::typeof(Expr), x::Number) = x
# convert(::typeof(Expr), s::Symbol) = s
# convert(::typeof(Expr), expr::Token) = expr.value

function parse_equation(s::String; variables=AbstractString[])::LTree
    expr = Lerche.parse(parser_equation, s)
    expr = sanitize(expr; variables=variables)
    return expr
end

function parse_assignment(s::String; variables=AbstractString[])::LTree
    expr = Lerche.parse(parser_assignment, s)
    expr = sanitize(expr; variables=variables)
    return expr
end

function parse_equation_block(s::String; variables=AbstractString[])::LTree
    expr = Lerche.parse(parser_equation_block, s)
    expr = sanitize(expr; variables=variables)
    return expr
end

function parse_assignment_block(s::String; variables=AbstractString[])::LTree
    expr = Lerche.parse(parser_assignment_block, s)
    expr = sanitize(expr; variables=variables)
    return expr
end

max_issues_repo_name: EconForge/Dolang


text:
if !isdefined(:runtests)
    addprocs(1)
end
srand(888)

description = """
Example showing off how to run GLVisualize in a different process
and visualize objects created on the main process.
"""

const workerid = workers()[]

using Images, GeometryTypes, GLVisualize, Reactive, GLWindow, Colors
using GLAbstraction

# Parallel helper function to launch our glvisualize operations on the process,
# on which we run GLVisualize
# I will move these functions to GLVisualize once this is tested better!
function p_setarg!(obj, key, value)
    @spawnat workerid begin
        GLAbstraction.set_arg!(fetch(obj), key, value)
        nothing
    end
end
function p_view(main, style=:default, cam=:perspective; kw_args...)
    @spawnat workerid begin
        obj = GLVisualize.visualize(main, style; kw_args...)
        GLVisualize._view(obj, camera=cam)
        obj
    end
end
function p_isopen()
    remotecall_fetch(workerid) do
        isempty(GLVisualize.get_screens()) && return false
        isopen(GLVisualize.current_screen())
    end
end
function p_empty!()
    remotecall(workerid) do
        empty!(GLVisualize.current_screen())
        nothing
    end
end

function solve_particles!(
        positions::AbstractVector{Point{N, T}}, s, dt::T = T(0.01)
    ) where {N, T}
    P = Point{N, T}
    # since Points are immutable, and we don't have comprehensions yet, there is
    # a map(index _> ..., ::Point) which we can use
    mask = map(i-> ifelse(isodd(i), -1, 1), P)
    @inbounds for i in eachindex(positions)
        vel = zero(P)
        posi = positions[i]
        for j in eachindex(positions)
            posj = positions[j]
            d = posj - posi
            distsq = dot(d, d) + T(1)
            vel = vel .+ (mask .* (s[j]*reverse(d)/distsq))
            any(x-> abs(x) > T(0.8), vel) && break # restrict velocity
        end
        positions[i] = posi .+ dt*vel
    end
    positions, s
end

"""
Generic implementation of initialising `N` dimensional points on an `N` dimensional
Sphere.
"""
function startpositions(N, radius::T, n) where T
    sphere = HyperSphere(Point{N, T}(0), T(radius))
    n = N == 3 ? floor(Int, sqrt(n)) : n # n must be n^2 for 3D Sphere
    decompose(Point{N, T}, sphere, n)
end

"""
The main simulation method, setting everything up and starting the simulation.
`T` controles what number type to use and defaults to Float32, but
code is written to support T == Float64 as well. So you can quickly try out both.
`N` changes the dimensionality between 2D and 3D
"""
function main(n, T=Float32, N=3)
    p_empty!() # empty window, so that we can call main multiple times

    positions = startpositions(N, T(0.5), n)
    # random factor for turbulence
    s = (-1 + 2 * rand(T, length(positions)))

    # boundingbox is an expensive operation, so if you don't need it
    # you can simply set it to nothing.
    # The signature for all kind of particles is:
    # visualize((primitive, positions), keyword_arguments...)
    # Note that we use p_view here, to visualize it on the worker process
    pointobj = p_view(
        (Circle(Point2f0(0), 0.01f0), positions),
        color = RGBA(1f0, 1f0, 1f0, 0.6f0),
        stroke_color = RGBA(1f0, 1f0, 1f0, 1f0),
        stroke_width = 0.001f0,
        boundingbox = nothing,
    )
    # maximum number of line segments we use to trace particle trajectory
    max_history = 80
    # To visualize a large number of lines with the same length, Matrix{Point} is
    # the way to go.
    lines = fill(Point{N, T}(NaN), max_history, length(positions))
    lines_color = fill(RGBA{Float32}(0, 0, 0, 0), max_history, length(positions))

    lines[1, :] = positions
    lines_color[1, :] = RGBA{Float32}(1, 1, 1, 0.2)

    linesobj = p_view(
        lines, :lines,
        color = vec(lines_color), # needs to be 1D right now, to lessen the amount of automatic conversions
        boundingbox = nothing,
        # this is an internal detail, need to figure out a better API
        # this makes the lines ignore NaN to seperate multiple lines
        startend = nothing,
        thickness = 0.5f0
    )
    t = 0.0
    dt = T(0.1)
    wait(linesobj) # wait until actually visualized
    while p_isopen()
        # lets print the times, to observe the slowing down when the velocity
        # cuttoff doesn't get reached as quickly anymore. Lucky for us that
        # GLVisualize runs in another process, so that the camera etc is still smooth
        tic()
        @time solve_particles!(positions, s, dt)

        # this seems a lot faster than circshift(lines, (1, 0))
        c = RGBA{Float32}(1, 1, (sin(t)+1.)/2., 0.2)
        for (i, p) in enumerate(positions)
            lines[:, i] = circshift(view(lines, :, i), 1)
            lines[1, i] = p
            lines_color[:, i] = circshift(view(lines_color, :, i), 1)
            lines_color[1, i] = c
        end

        p_setarg!(linesobj, :vertex, map(Point{N, Float32}, vec(lines)))
        p_setarg!(linesobj, :color, vec(lines_color))
        p_setarg!(pointobj, :position, map(Point{N, Float32}, positions))
        t += 0.05
        dt = T(toq())
    end
end

# start GLVisualize worker process
@spawnat workerid begin
    window = GLVisualize.glscreen(color = RGBA(0f0, 0f0, 0f0, 1f0))
    @async GLWindow.renderloop(window)
    nothing
end
# start main
main(10_000)

# clean up process
remotecall_fetch(workerid) do
    GLVisualize.cleanup()
    return
end
nothing


text: # stopped at https://juliadocs.github.io/Documenter.jl/stable/man/guide/#Building-an-empty-document-1


text:
# optmizie.jl - optimize ExGraph

const OPT_PHS = [:X, :Y, :Z, :V, :W,
                 :i, :j, :k, :l, :m, :n, :p, :q, :r, :s, :t]

const OPT_VEC_RULES = [
    :(Z = X * transpose(Y)) => :(Z = X * Y'),
    :(Z = X * Y') => :(Z = X * Y'),
    :(Z = transpose(X) * Y) => :(Z = X' * Y),
    :(Z = X' * Y) => :(Z = X' * Y),
]

function remove_unused(g::AbstractExGraph, outvars::Vector{Symbol})
    deps = collect_deps(g, outvars)
    push!(deps, outvars...)
    gr = reset_tape(g)
    for nd in g.tape
        if in(varname(nd), deps)
            push!(gr, nd)
        end
    end
    return gr
end

remove_unused(g::AbstractExGraph, outvar::Symbol) = remove_unused(g, [outvar])
remove_unused(g::AbstractExGraph) = remove_unused(g, [varname(g[end])])

"""
Removes unused variables from multiline expressions, e.g. in:

    x = u * v
    y = x + 1
    z = 2x

`y` isn't used to compute output variable `z`, so it's removed:

    x = u * v
    z = 2x

"""
function remove_unused(ex::Expr; output_vars=nothing)
    ex.head == :block || return ex  # nothing to remove
    g = ExGraph(ex; fuse=false)
    output_vars = output_vars != nothing ? output_vars : [varname(g[end])]
    deps = collect_deps(g, output_vars)
    for vname in output_vars
        push!(deps, vname)
    end
    res = quote end
    for subex in ex.args
        if subex.head == :(=)
            vname = split_indexed(subex.args[1])[1]
            if in(vname, deps)
                push!(res.args, subex)
            end
        else
            push!(res.args, subex)
        end
    end
    return res
end

function reset_tape(g::G) where {G <: AbstractExGraph}
    # new_g = deepcopy(g)
    new_g = G()
    new_g.tape = []
    new_g.idx = Dict()
    new_g.ctx = g.ctx
    return new_g
end

function tryoptimize(ex::Expr)
    for (pat, subs_ex) in OPT_VEC_RULES
        new_ex = tryrewrite(ex, pat, subs_ex; phs=OPT_PHS, allow_ex=false)
        if new_ex != nothing
            genname_patterns = unique(findex(:(genname__(_n)), new_ex))
            if !isempty(genname_patterns)
                new_names = gennames(length(genname_patterns))
                st = Dict(zip(genname_patterns, new_names))
                return subs(new_ex, st)
            else
                return new_ex
            end
        end
    end
    return nothing
end

function expand_fixed_sequences(g::ExGraph, nd::ExNode)
    ex = to_expr(nd)
    changed = false
    for dep in dependencies(nd)
        expanded = subs(ex, Dict(dep => getexpr(g[dep])))
        new_ex = tryoptimize(expanded)
        if new_ex != nothing
            ex = new_ex
            changed = true
        end
    end
    return changed ? copy(nd; category=:opaque, var=ex.args[1], ex=ex.args[2]) : nd
end

"""
Look at each node's dependencies and, if there are known pattern sequences,
rewrite them in a more optimal way.
"""
function expand_fixed_sequences(g::ExGraph)
    new_g = reset_tape(g)
    for nd in g.tape
        if isa(nd, ExNode{:input})
            push!(new_g, nd)
        else
            new_nd = expand_fixed_sequences(g, nd)
            push!(new_g, new_nd)
        end
    end
    new_g = remove_unused(new_g,  varname(new_g[end]))
    # new_g = fuse_assigned(new_g)
    return new_g
end

# common subexpression elimination

safe_size(a::AbstractArray) = size(a)
safe_size(a) = nothing

function eliminate_common(g::AbstractExGraph)
    g = reindex_from_beginning(g)
    new_g = reset_tape(g)
    existing = Dict()             # index of existing expressions
    st = Dict{Symbol,Symbol}()    # later var -> previous var with the same expression
    for nd in g.tape
        new_full_ex = subs(to_expr_kw(nd), st)
        vname, vidxs = split_indexed(new_full_ex.args[1])
        key = string((vidxs, new_full_ex.args[2]))
        if haskey(existing, key) &&
            (g[vname] |> getvalue |> safe_size) == (g[existing[key]] |> getvalue |> safe_size)
            st[vname] = existing[key]
        else
            # C = getcategory(nd)
            new_nd = copy(nd)
            setexpr_kw!(new_nd, new_full_ex.args[2])
            push!(new_g, new_nd)
            existing[key] = vname
        end
    end
    rename!(new_g, st)
    return new_g
end

# fuse broadcasting

function is_bcast_vec(nd::ExNode{:opaque})
    # convert to expression, re-parse and check every node
    return all(is_bcast_vec(nd) for nd in ExGraph(to_expr(nd)).tape)
end

const OLD_BCAST_OPS = Set([:.+, :.-, :.*, :./, :.^])

is_bcast_vec(nd::ExNode{:call}) = getexpr(nd).args[1] in OLD_BCAST_OPS
is_bcast_vec(nd::ExNode{:bcast}) = true
is_bcast_vec(nd::ExNode) = false

function fuse_broadcasting_node(g::ExGraph, new_g::ExGraph, nd::ExNode)
    dep_nds = [new_g[dep] for dep in dependencies(nd)]
    if any(is_bcast_vec(dep_nd) for dep_nd in dep_nds)
        # at least one dependency is broadcasting
        # transform node to :opaque and expand dep expression
        st = Dict{Any, Any}()
        for dep_nd in dep_nds
            if is_bcast_vec(dep_nd)
                # if dependency is broadcasting, replace its name with its expression
                dep_ex = getexpr(dep_nd)
                # idx_st = Dict(zip(varidxs(dep_nd),
                #                   get_indices_in_expr(getexpr(nd), varname(dep_nd))))
                # new_dep_ex = subs(dep_ex, idx_st)
                new_dep_ex = dep_ex
                st[getvar(dep_nd)] = new_dep_ex
            end
        end
        new_ex = subs(getexpr(nd), st)
        return copy(nd; ex=new_ex, category=:opaque)
    else
        return nd
    end
end

function fuse_broadcasting(g::ExGraph)
    new_g = reset_tape(g)
    for nd in g.tape
        if is_bcast_vec(nd)
            push!(new_g, fuse_broadcasting_node(g, new_g, nd))
        else
            push!(new_g, nd)
        end
    end
    return remove_unused(new_g, varname(g[end]))
end


text:
function analysis_functions(c::Coupler)
    # Extract number of branches
    n = length(c.a)

    # Form ABCD matrix
    M = _ABCD_open_stub(c.a[1])
    for i in eachindex(c.b)
        M *= (_ABCD_line(c.b[i]) * _ABCD_open_stub(c.a[i+1]))
    end
    A_ = M[1,1]
    B_ = M[1,2]
    C_ = M[2,1]
    D_ = M[2,2]

    # Γₑ_N
    Γₑ_N = (A_ - D_) + (B_ - C_)

    # Γₑ_D
    Γₑ_D = (A_ + D_ + B_ + C_)

    # Tₑ_N
    Tₑ_N = 2*(Polynomial([1, 0, -1], :t)^(n-1))

    # Create functions
    # Reflection coefficient
    Γ(t) = (1/2) * (Γₑ_N(t)/Γₑ_D(t) + Γₑ_N(1/t)/Γₑ_D(1/t))

    # Directivity
    D(t) = (1/2) * (Γₑ_N(t)/Γₑ_D(t) - Γₑ_N(1/t)/Γₑ_D(1/t))

    # Insertion Loss
    IL(t)  = (1/2) * (Tₑ_N(t)/Γₑ_D(t) + (-1)^(n) * Tₑ_N(1/t)/Γₑ_D(1/t))

    # Coupling
    C(t) = (1/2) * (Tₑ_N(t)/Γₑ_D(t) - (-1)^(n) * Tₑ_N(1/t)/Γₑ_D(1/t))

    return (Γ = Γ, D = D, IL = IL, C = C)
end

db(x) = 10*log10(abs2(x))

function _ABCD_open_stub(a)
    A = Polynomial(Float64[1], :t)
    B = Polynomial(Float64[0], :t)
    C = Polynomial(Float64[0, a], :t)
    D = Polynomial(Float64[1], :t)
    return [A B;C D]
end

function _ABCD_line(b)
    A = Polynomial(Float64[1, 0, 1], :t)
    B = Polynomial(Float64[0, 2/b], :t)
    C = Polynomial(Float64[0, 2*b], :t)
    D = Polynomial(Float64[1, 0, 1], :t)
    return [A B;C D]
end


text:
# Combined 1D/2D plots across multiple frames, multi-process version.
#
# To run on a single node,
# julia -p $ncores demo_1d2d_parallel_pyplot.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Distributed, ParallelDataTransfer, Glob
@everywhere using Vlasiator, PyPlot, Printf, LaTeXStrings
@everywhere using Vlasiator: set_args, prep2d, set_colorbar, set_plot

@assert matplotlib.__version__ ≥ "3.4" "Require Matplotlib version 3.4+ to use subfigure!"

@everywhere function init_figure(x1, x2)
   fig = plt.figure(myid(), constrained_layout=true, figsize=(12, 12))
   subfigs = fig.subfigures(1, 2, wspace=0.05)

   axsL = subfigs[1].subplots(4, 1, sharex=true)
   axsR = subfigs[2].subplots(2, 1, sharex=true, sharey=true)

   # Set line plots' axes
   axsL[end].set_xlim(x1, x2)

   axsL[1].set_ylim(ρmin, ρmax)
   axsL[2].set_ylim(vmin, vmax)
   axsL[3].set_ylim(pmin, pmax)
   axsL[4].set_ylim(bmin, bmax)
   for ax in axsL
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.grid(true)
   end

   axsL[end].set_xlabel(L"x [$R_E$]"; fontsize)
   axsL[1].set_ylabel("Density [amu/cc]";    fontsize)
   axsL[2].set_ylabel("Velocity [km/s]";     fontsize)
   axsL[3].set_ylabel("Pressure [nPa]";      fontsize)
   axsL[4].set_ylabel("Magnetic field [nT]"; fontsize)

   fakeline = loc
   l1 = axsL[1].plot(loc, fakeline, label="Proton density", color="#1f77b4")
   l2 = axsL[2].plot(loc, fakeline, label="Vx",             color="#1f77b4")
   l3 = axsL[3].plot(loc, fakeline, label="Dynamic",        color="#1f77b4")
   l4 = axsL[3].plot(loc, fakeline, label="Thermal",        color="#ff7f0e")
   l5 = axsL[4].plot(loc, fakeline, label="Bz",             color="#1f77b4")

   ls = (l1, l2, l3, l4, l5)

   axsL[2].legend(;loc="upper right", fontsize)
   axsL[3].legend(;loc="lower right", fontsize)
   axsL[4].legend(;loc="upper right", fontsize)

   vl1 = axsL[1].vlines(loc[imagnetopause_], ρmin, ρmax;
      colors="r", linestyle="dashed", alpha=0.5)

   vl2 = axsL[2].vlines(loc[imagnetopause_], vmin, vmax;
      colors="r", linestyle="dashed", alpha=0.5)

   vl3 = axsL[3].vlines(loc[imagnetopause_], pmin, pmax;
      colors="r", linestyle="dashed", alpha=0.5)

   vl4 = axsL[4].vlines(loc[imagnetopause_], bmin, bmax;
       colors="r", linestyle="dashed", alpha=0.5)
   hl4 = axsL[4].hlines(0.0, loc[1], loc[end]; colors="k", linestyle="dashed", alpha=0.2)

   vlines = (vl1, vl2, vl3, vl4)

   for ax in axsR
      ax.set_aspect("equal")

      # Set border line widths
      for loc in ("left", "bottom", "right", "top")
         edge = get(ax.spines, loc, nothing)
         edge.set_linewidth(2.0)
      end
      ax.xaxis.set_tick_params(width=2.0, length=3)
      ax.yaxis.set_tick_params(width=2.0, length=3)
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())

      ax.set_ylabel(L"Y [$R_E$]"; fontsize)
   end

   axsR[2].set_xlabel(L"X [$R_E$]"; fontsize)

   axsR[1].set_title("Alfven speed", fontsize="x-large")
   axsR[2].set_title("Sound speed", fontsize="x-large")

   x, y = Vlasiator.get_axis(pArgs1)
   fakedata = zeros(Float32, length(y), length(x))

   c1 = axsR[1].pcolormesh(x, y, fakedata, norm=cnorm1, cmap=cmap, shading="nearest")
   c2 = axsR[2].pcolormesh(x, y, fakedata, norm=cnorm2, cmap=cmap, shading="nearest")

   cb1 = colorbar(c1; ax=axsR[1], ticks=cticks1, fraction=0.046, pad=0.04)
   cb1.ax.set_ylabel("[km/s]"; fontsize)
   cb1.outline.set_linewidth(1.0)

   cb2 = colorbar(c2; ax=axsR[2], ticks=cticks2, fraction=0.046, pad=0.04)
   cb2.ax.set_ylabel("[km/s]"; fontsize)
   cb2.outline.set_linewidth(1.0)

   fig.suptitle("Density Pulse Run", fontsize="xx-large")

   cs = (c1, c2)

   return fig, subfigs, ls, vlines, cs
end

function update_vline(h, x)
   seg_old = h.get_segments()
   ymin = seg_old[1][1, 2]
   ymax = seg_old[1][2, 2]

   seg_new = [[[x, ymin] [x, ymax]]]

   h.set_segments(seg_new)
end

@everywhere function process(subfigs, ls, vlines, cs, file, cellids)
   isfile("../out/"*file[end-8:end-5]*".png") && return

   println("file = $(basename(file))")
   meta = load(file)

   p_extract = readvariable(meta, "vg_pressure", cellids) .* 1e9 |> vec # [nPa]
   rho_extract = readvariable(meta, "proton/vg_rho", cellids) |> vec
   v_extract = readvariable(meta, "proton/vg_v", cellids)
   vmag2_extract = sum(v_extract.^2, dims=1) |> vec
   pdyn_extract = rho_extract .* Vlasiator.mᵢ .* vmag2_extract .* 1e9 # [nPa]

   bz = readvariable(meta, "vg_b_vol", cellids)[3,:] .* 1e9 #[nT]

   ls[1][1].set_ydata(rho_extract ./ 1e6)
   ls[2][1].set_ydata(v_extract[1,:] ./ 1e3)
   ls[3][1].set_ydata(pdyn_extract)
   ls[4][1].set_ydata(p_extract)
   ls[5][1].set_ydata(bz)

   imagnetopause_ = findfirst(<(0.0), bz)
   for vline in vlines
      update_vline(vline, loc[imagnetopause_])
   end

   str_title = @sprintf "Sun-Earth line, t= %4.1fs" meta.time
   subfigs[1].suptitle(str_title, fontsize="x-large")

   data = prep2d(meta, "VA", :z)'
   cs[1].set_array(data ./ 1e3)

   data = prep2d(meta, "VS", :z)'
   cs[2].set_array(data ./ 1e3)

   savefig("../out/"*file[end-8:end-5]*".png", bbox_inches="tight")
   return
end

function make_jobs(files)
   for f in files
       put!(jobs, f)
   end
end

@everywhere function do_work(jobs, status)
   fig, subfigs, ls, vlines, cs = init_figure(x1, x2)
   while true
      file = take!(jobs)
      process(subfigs, ls, vlines, cs, file, cellids)
      put!(status, true)
   end
   close(fig)
end

############################################################################################
files = glob("bulk*.vlsv", "../run_rho2_bz-5_timevarying_startfrom300s")
const nfile = length(files)
@passobj 1 workers() files

@broadcast begin
   # Set contour plots' axes and colorbars
   const cmap = matplotlib.cm.turbo
   colorscale = Linear
   axisunit = RE
   # Upper/lower limits for each variable
   const ρmin, ρmax = 0.0, 10.0     # [amu/cc]
   const vmin, vmax = -640.0, 0.0   # [km/s]
   const pmin, pmax = 0.0, 1.82     # [nPa]
   const bmin, bmax = -25.0, 60.0   # [nT]
   const vamin, vamax = 0.0, 250.0  # [km/s]
   const vsmin, vsmax = 30.0, 350.0 # [km/s]

   meta = load(files[1])

   const pArgs1 = set_args(meta, "VA", axisunit; normal=:none)
   const cnorm1, cticks1 = set_colorbar(colorscale, vamin, vamax)
   const cnorm2, cticks2 = set_colorbar(colorscale, vsmin, vsmax)

   const fontsize = 14
end

const jobs   = RemoteChannel(()->Channel{String}(nfile))
const status = RemoteChannel(()->Channel{Bool}(nworkers()))

Re = Vlasiator.Re # Earth radii
x1, x2 = 8.0, 29.0
point1 = [x1, 0, 0] .* Re
point2 = [x2, 0, 0] .* Re

meta = load(files[1])
cellids, _, _ = getcellinline(meta, point1, point2)

passobj(1, workers(), [:x1, :x2, :cellids])
@broadcast const loc = range(x1, x2, length=length(cellids))

println("Total number of files: $nfile")
println("Running with $(nworkers()) workers...")

@async make_jobs(files) # Feed the jobs channel with all files to process.

@sync for p in workers()
   @async remote_do(do_work, p, jobs, status)
end

let n = nfile
   t = @elapsed while n > 0 # wait for all jobs to complete
      take!(status)
      n -= 1
   end
   println("Finished in $(round(t, digits=2))s.")
end


text:
"""
    latticebasis(A::AbstractVector...)

Given a tuple of basis vectors `A = (a1,  a2, ...)` returns a SMatrix with the vector basis as columns.
"""
function latticebasis(A::AbstractVector...)
    dim = length(A)
    for ai in A
        @assert length(ai) == dim "number of basis vectors, must coincide with their dimension"
    end

    return SMatrix{dim, dim}(hcat(A...))
end

"""
    reciprocalbasis(A::AbstractMatrix)

Given a basis `A = [a1 a2 ...]` (with the basis vectors forming the matrix columns), returns the reciprocal basis `B = [b1 b2 ...]` (with the reciprocal basis vectors forming the matrix columns) such that:
dot(aᵢ, bⱼ) = 2πδᵢⱼ
"""
reciprocalbasis(A::AbstractMatrix) = 2π*transpose(inv(A))

"""
    spanlattice(B::SMatrix, cutoff::Float64)

Span the lattice points generated by the lattice basis vectors B[:, i], with length shorter than `cutoff`.
"""
function spanlattice(B::SMatrix, cutoff::Float64)

    Ns = [ceil(Int, cutoff/norm(B[:, i])) for i=1:size(B, 2)]
    dim = size(B, 1)

    itr = Iterators.product([-Ni:Ni for Ni in Ns]...)

    list=SVector{dim, Float64}[]

    for I in itr

        pt = B*SVector{dim}(I...)
        if norm(pt) <= cutoff
            push!(list, pt)
        end
    end

    return list
end

"""
    addlattice(lats::Vector{Vector{StaticVectors}}, cutoff::Float64)

Given n spanned lattices [lats[1], ..., lats[N]], form tuples of vector of the form (R₁, ..., Rₙ) with Rᵢ belonging to lats[i],  such that norm(R₁ + ... + Rₙ) < cutoff.
"""
function addlattice(lats::Vector{Vector{SVector{D, T}}}, cutoff::Float64) where {D, T}

    itr = Iterators.product(lats...)
    Nlat = length(lats)
    list = NTuple{Nlat, SVector{D, T}}[]

    for G in itr
        if norm(sum(G)) < cutoff
            push!(list, G)
        end
    end

    return list
end


text:
"""
    propagate(f, g, aggr; xi, xj, e)  ->  m̄

Performs message passing on graph `g`. Takes care of materializing the node features on each edge,
applying the message function, and returning an aggregated message ``\\bar{\\mathbf{m}}``
(depending on the return value of `f`, an array or a named tuple of
arrays with last dimension's size `g.num_nodes`).

It can be decomposed in two steps:

```julia
m = apply_edges(f, g, xi, xj, e)
m̄ = aggregate_neighbors(g, aggr, m)
```

GNN layers typically call `propagate` in their forward pass,
providing as input `f` a closure.

# Arguments

- `g`: A `GNNGraph`.
- `xi`: An array or a named tuple containing arrays whose last dimension's size
        is `g.num_nodes`. It will be appropriately materialized on the
        target node of each edge (see also [`edge_index`](@ref)).
- `xj`: As `xj`, but to be materialized on edges' sources.
- `e`: An array or a named tuple containing arrays whose last dimension's size is `g.num_edges`.
- `f`: A generic function that will be passed over to [`apply_edges`](@ref).
      Has to take as inputs the edge-materialized `xi`, `xj`, and `e`
      (arrays or named tuples of arrays whose last dimension' size is the size of
      a batch of edges). Its output has to be an array or a named tuple of arrays
      with the same batch size.
- `aggr`: Neighborhood aggregation operator. Use `+`, `mean`, `max`, or `min`.

# Usage Examples

```julia
using GraphNeuralNetworks, Flux

struct GNNConv <: GNNLayer
    W
    b
    σ
end

Flux.@functor GNNConv

function GNNConv(ch::Pair{Int,Int}, σ=identity)
    in, out = ch
    W = Flux.glorot_uniform(out, in)
    b = zeros(Float32, out)
    GNNConv(W, b, σ)
end

function (l::GNNConv)(g::GNNGraph, x::AbstractMatrix)
    message(xi, xj, e) = l.W * xj
    m̄ = propagate(message, g, +, xj=x)
    return l.σ.(m̄ .+ l.bias)
end

l = GNNConv(10 => 20)
l(g, x)
```

See also [`apply_edges`](@ref).
"""
function propagate end

propagate(l, g::GNNGraph, aggr; xi=nothing, xj=nothing, e=nothing) =
    propagate(l, g, aggr, xi, xj, e)

function propagate(l, g::GNNGraph, aggr, xi, xj, e)
    m = apply_edges(l, g, xi, xj, e)
    m̄ = aggregate_neighbors(g, aggr, m)
    return m̄
end

## APPLY EDGES

"""
    apply_edges(f, g, xi, xj, e)
    apply_edges(f, g; [xi, xj, e])

Returns the message from node `j` to node `i` .
In the message-passing scheme, the incoming messages
from the neighborhood of `i` will later be aggregated
in order to update the features of node `i`.

The function operates on batches of edges, therefore
`xi`, `xj`, and `e` are tensors whose last dimension
is the batch size, or can be named tuples of
such tensors.

# Arguments

- `g`: A `GNNGraph`.
- `xi`: An array or a named tuple containing arrays whose last dimension's size
        is `g.num_nodes`. It will be appropriately materialized on the
        target node of each edge (see also [`edge_index`](@ref)).
- `xj`: As `xj`, but to be materialized on edges' sources.
- `e`: An array or a named tuple containing arrays whose last dimension's size is `g.num_edges`.
- `f`: A function that takes as inputs the edge-materialized `xi`, `xj`, and `e`.
       These are arrays (or named tuples of arrays) whose last dimension' size is the size of
       a batch of edges. The output of `f` has to be an array (or a named tuple of arrays)
       with the same batch size.

See also [`propagate`](@ref).
"""
function apply_edges end

apply_edges(l, g::GNNGraph; xi=nothing, xj=nothing, e=nothing) =
    apply_edges(l, g, xi, xj, e)

function apply_edges(f, g::GNNGraph, xi, xj, e)
    s, t = edge_index(g)
    xi = _gather(xi, t)   # size: (D, num_nodes) -> (D, num_edges)
    xj = _gather(xj, s)
    m = f(xi, xj, e)
    return m
end

_gather(x::NamedTuple, i) = map(x -> _gather(x, i), x)
_gather(x::Tuple, i) = map(x -> _gather(x, i), x)
_gather(x::AbstractArray, i) = NNlib.gather(x, i)
_gather(x::Nothing, i) = nothing

##  AGGREGATE NEIGHBORS

function aggregate_neighbors(g::GNNGraph, aggr, m)
    s, t = edge_index(g)
    return _scatter(aggr, m, t)
end

_scatter(aggr, m::NamedTuple, t) = map(m -> _scatter(aggr, m, t), m)
_scatter(aggr, m::Tuple, t) = map(m -> _scatter(aggr, m, t), m)
_scatter(aggr, m::AbstractArray, t) = NNlib.scatter(aggr, m, t)

### SPECIALIZATIONS OF PROPAGATE ###
"""
    copy_xj(xi, xj, e) = xj
"""
copy_xj(xi, xj, e) = xj

"""
    copy_xi(xi, xj, e) = xi
"""
copy_xi(xi, xj, e) = xi

"""
    xi_dot_xj(xi, xj, e) = sum(xi .* xj, dims=1)
"""
xi_dot_xj(xi, xj, e) = sum(xi .* xj, dims=1)

"""
    e_mul_xj(xi, xj, e) = reshape(e, (...)) .* xj

Reshape `e` into broadcast compatible shape with `xj`
(by prepending singleton dimensions) then perform
broadcasted multiplication.
"""
function e_mul_xj(xi, xj::AbstractArray{Tj,Nj}, e::AbstractArray{Te,Ne}) where {Tj,Te, Nj, Ne}
    @assert Ne <= Nj
    e = reshape(e, ntuple(_ -> 1, Nj-Ne)..., size(e)...)
    return e .* xj
end

function propagate(::typeof(copy_xj), g::GNNGraph, ::typeof(+), xi, xj::AbstractMatrix, e)
    A = adjacency_matrix(g)
    return xj * A
end

## avoid the fast path on gpu until we have better cuda support
function propagate(::typeof(copy_xj), g::GNNGraph{<:Union{COO_T,SPARSE_T}}, ::typeof(+), xi, xj::AnyCuMatrix, e)
    propagate((xi,xj,e)->copy_xj(xi,xj,e), g, +, xi, xj, e)
end

# function propagate(::typeof(copy_xj), g::GNNGraph, ::typeof(mean), xi, xj::AbstractMatrix, e)
#     A = adjacency_matrix(g)
#     D = compute_degree(A)
#     return xj * A * D
# end

# # Zygote bug. Error with sparse matrix without nograd
# compute_degree(A) = Diagonal(1f0 ./ vec(sum(A; dims=2)))

# Flux.Zygote.@nograd compute_degree


text:
"""
    struct Token <: Any
Hold a pointer to a `clang::Token` object.
"""
struct Token
    ptr::CXToken_
end

Base.unsafe_convert(::Type{CXToken_}, x::Token) = x.ptr
Base.cconvert(::Type{CXToken_}, x::Token) = x

"""
    struct AnnotationValue <: Any
Hold a pointer to an "AnnotationValue".
"""
struct AnnotationValue
    ptr::CXAnnotationValue
end

Base.unsafe_convert(::Type{CXAnnotationValue}, x::AnnotationValue) = x.ptr
Base.cconvert(::Type{CXAnnotationValue}, x::AnnotationValue) = x


text:
@testset "classify matrices" begin
    A = [0..2 1..1;-1.. -1 0..2]
    @test !is_H_matrix(A)
    @test !is_strongly_regular(A)

    B = [-2.. -2 1..1; 5..6 -2.. -2]
    @test is_strongly_regular(B)
    @test !is_Z_matrix(B)
    @test !is_M_matrix(B)

    C = [2..2 1..1; 0..2 2..2]
    @test is_H_matrix(C)
    @test !is_strictly_diagonally_dominant(C)

    D = [2..2 -1..0; -1..0 2..2]
    @test is_strictly_diagonally_dominant(D)
    @test is_Z_matrix(D)
    @test is_M_matrix(D)

    E = [2..4 -2..1;-1..2 2..4]
    @test !is_Z_matrix(E)
    @test !is_M_matrix(E)
    @test !is_H_matrix(E)
    @test is_strongly_regular(E)
end


text:
import Dates
println("Last updated: ", Dates.now(), " (PT)")

# Load environment
import Pkg; Pkg.activate("../../../")

# Import Libraries
using Turing
using Turing: Variational
import Turing.RandomMeasures.DirichletProcess
import Turing.RandomMeasures.ChineseRestaurantProcess
using Distributions
using JSON3
using PyPlot
using StatsFuns
import Random
import StatsBase.countmap
include(joinpath(@__DIR__, "../util/BnpUtil.jl"));

# Directory where all simulation data are stored.
data_dir = joinpath(@__DIR__, "../data")
path_to_data = joinpath(data_dir, "gmm-data-n50.json")

# Load data in JSON format.
data = let
    x = open(f -> read(f, String), path_to_data)
    JSON3.read(x, Dict{Symbol, Vector{Any}})
end

# Convert data to vector of floats
y = Float64.(data[:y]);

# Visualize data
plt.hist(y, bins=50, density=true)
plt.xlabel("y")
plt.ylabel("density")
plt.title("Histogram of data");

function counts_plot(x; density=false, color="C0", lw=4)
    cm_x = countmap(x)
    number = collect(keys(cm_x))
    count = let
        c = collect(values(cm_x))
        density ? c / sum(c) : c
    end
    plt.vlines(number, 0, count, color=color, lw=lw)
end

exp_num_clus(a, n) = a * log1p(n / a)

alpha_prior = Uniform(0, 0.5)
a = rand(alpha_prior, 10000)
nclus = rand.(Poisson.(exp_num_clus.(a, length(y))))

counts_plot(nclus, density=true)

plt.xlabel("counts")
plt.ylabel("prior probability")
plt.title("prior distribution for number of clusters");

# a = [0.1, 0.5, 1.0, 1.5, 2.0]
# plt.scatter(a, exp_num_clus.(a, 50))

@model function dp_gmm_crp_delayed(y, Kmax, ::Type{T}=Vector{Float64}) where T
    nobs = length(y)

    alpha ~ Gamma(1, 1/10)  # mean = a*b
    # alpha ~ Uniform(0, 1.0)
    rpm = DirichletProcess(alpha)

    # Base measure.
    H_mu = Normal(0, 3)
    H_sigma = Uniform(0, 0.3)
    # H_sigma = LogNormal(-1, 0.5)

    # Latent assignment.
    z = tzeros(Int, nobs)
    n = tzeros(Int, Kmax)
    K = 0

    mu = T(undef, Kmax)
    sigma = T(undef, Kmax)

    for i in 1:nobs
        z[i] ~ ChineseRestaurantProcess(rpm,  n)
        n[z[i]] += 1

        if z[i] > K
            mu[z[i]] ~ H_mu
            sigma[z[i]] ~ H_sigma
            K += 1
        end

        y[i] ~ Normal(mu[z[i]], sigma[z[i]])
    end
end;

# Set random seed for reproducibility
Random.seed!(0);

# Compile time approx. 32s.
# Run time approx. 70s

@time chain = begin
    burn = 3000  # NOTE: The burn in is also returned. Discard manually.
    n_samples = 1000
    iterations = burn + n_samples

    sample(dp_gmm_crp_delayed(y, 20), SMC(), iterations)

    # sample(dp_gmm_crp_delayed(y, 30), PG(10), iterations)
    # sample(dp_gmm_crp_delayed(y, 15), SMC(), iterations)
    # sample(dp_gmm_crp_delayed(y, 30), Gibbs(PG(20, :z), PG(10, :alpha, :mu, :sigma)), iterations)
    # sample(dp_gmm_crp_delayed(y, 15), Gibbs(PG(20, :z), PG(10, :mu, :sigma)), iterations)
    # sample(dp_gmm_crp_delayed(y, 30), Gibbs(MH(0.1, :alpha), PG(50, :z), MH(0.01, :mu, :sigma)), iterations)

    # sample(dp_gmm_crp_delayed(y, 30), Gibbs(MH(0.1, :alpha), PG(10, :z), HMC(0.01, 100, :mu, :sigma)), iterations)

    # NOTE: Works very poorly...
    # sample(dp_gmm_crp_delayed(y, 30), Gibbs(PG(5, :z), HMC(0.01, 10, :alpha, :mu, :sigma)), iterations)

    # sample(dp_gmm_crp(y), SMC(), iterations)
    # sample(dp_gmm_crp(y), IS(), iterations)
    # sample(dp_gmm_crp(y), Gibbs(PG(5, :z), PG(5, :mu_sigma)), iterations)
    # sample(dp_gmm_crp(y), PG(10), iterations)
end;

# TODO: Maybe try initializing at MAP?

burn=0
chain.value

z = chain.value.data[burn+1:end, (end-length(y)+1):end, 1]
nclus = [length(unique(z[i, :])) for i in 1:size(z, 1)];

plt.figure(figsize=(12,4))
plt.subplot(1, 3, 1)
plt.plot(nclus);
plt.xlabel("iteration")
plt.ylabel("number of clusters")

plt.subplot(1, 3, 2)
counts_plot(nclus, density=true);
plt.xlabel("number of clusters")
plt.ylabel("probability")

plt.subplot(1, 3, 3)
i = 1000
counts_plot(z[i, :], density=true);
plt.xlabel("cluster label for sample $(i)")
plt.ylabel("probability")

plt.tight_layout();

countmap(vec(Float64.(chain[:alpha].data)))

function extract_vec(chain, sym)
    syms = filter(x -> x != nothing, match.(Regex("$(sym)\\[\\d+\\]"), String.(chain.name_map[1])))
    syms = map(x -> Symbol(x.match), syms)
    return chain[syms]
end

plt.figure(figsize=(12, 4))

plt.subplot(1, 3, 1)
mu = coalesce.(extract_vec(chain, :mu).value.data[burn+1:end, :, 1], 0.0)
plt.boxplot(mu);
foreach(line -> plt.axhline(line, ls=":"), data[:mu])
plt.ylabel("mu")
plt.xlabel("mixture component")

plt.subplot(1, 3, 2)
sigma = coalesce.(extract_vec(chain, :sigma).value.data[burn+1:end, :, 1], 0.0)
plt.boxplot(coalesce.(sigma, 0.0));
foreach(line -> plt.axhline(line, ls=":"), data[:sig])
plt.ylabel("sigma")
plt.xlabel("mixture component")

plt.subplot(1, 3, 3)
plt.hist(Float64.(chain[:alpha].data[:, 1, 1]), density=true)
plt.xlabel("alpha")
plt.ylabel("density")

plt.tight_layout();

plt.plot(mu)
foreach(line -> plt.axhline(line, ls=":"), data[:mu])

jitter(x; eps=std(x)/30) = x + randn(length(x)) * eps

plt.scatter(jitter(vec(sigma)), jitter(vec(mu)), label="posterior samples (jittered)", s=1);
plt.scatter(data[:sig], data[:mu], marker="X", s=60, label="truth")
plt.xlabel("sigma")
plt.ylabel("mu")
# foreach(line -> plt.axhline(line, ls=":"), data[:mu])
# foreach(line -> plt.axvline(line, ls=":"), data[:sig])
plt.legend(loc="upper left");

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


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

"""
    JMatrix{T}

Square matrix of ones. See also [`J`](@ref).
"""
struct JMatrix{T}
  λ::T
end

(J::JMatrix{T})(d::Integer) where {T} = J.λ * Ones{T}(d, d)

*(λ::Number, J::JMatrix) = JMatrix(J.λ * λ)

+(J::JMatrix, A) = J(size(A, 1)) + A
+(A, J::JMatrix) = J + A

-(J::JMatrix, A) = J + (-A)
-(A, J::JMatrix) = -(J - A)

*(J::JMatrix, A) = repeat(J.λ * sum(A, dims=1), size(A, 1))
*(A, J::JMatrix) = (J*A')'

adjoint(J::JMatrix) = J

"""
    FMatrix{T}

`F` matrix, as defined by Aitchison 1986. See also [`F`](@ref).
"""
struct FMatrix{T}
  λ::T
end

(F::FMatrix{T})(d::Integer) where {T} = F.λ * [I(d) -Ones{T}(d)]

*(λ::Number, F::FMatrix) = FMatrix(F.λ * λ)

+(F::FMatrix, A) = F(size(A, 1)) + A
+(A, F::FMatrix) = F + A

-(F::FMatrix, A) = F + (-A)
-(A, F::FMatrix) = -(F - A)

*(F::FMatrix, A) = begin
  R = F.λ * (A[begin:end-1,:] .- A[end:end,:])
  ndims(A) == 1 ? vec(R) : R
end
*(A, F::FMatrix) = F.λ * [A -sum(A, dims=2)]

adjoint(F::FMatrix) = F′Matrix(F.λ)

"""
    F′Matrix{T}

Lazy adjoint of `FMatrix{T}`.
"""
struct F′Matrix{T}
  λ::T
end

(F′::F′Matrix{T})(d::Integer) where {T} = F′.λ * [I(d); -Ones{T}(1,d)]

*(λ::Number, F′::F′Matrix) = F′Matrix(F′.λ * λ)

+(F′::F′Matrix, A) = F′(size(A, 2)) + A
+(A, F′::F′Matrix) = F′ + A

-(F′::F′Matrix, A) = F′ + (-A)
-(A, F′::F′Matrix) = -(F′ - A)

*(F′::F′Matrix, A) = F′.λ * [A; -sum(A, dims=1)]
*(A, F′::F′Matrix) = F′.λ * (A[:,begin:end-1] .- A[:,end:end])

adjoint(F′::F′Matrix) = FMatrix(F′.λ)

"""
    GMatrix{T}

`G` matrix, as defined by Aitchison 1986. See also [`G`](@ref).
"""
struct GMatrix{T}
  λ::T
end

(G::GMatrix{T})(D::Integer) where {T} = G.λ * I(D) - (G.λ / D) * J(D)

*(λ::Number, G::GMatrix) = GMatrix(G.λ * λ)

+(G::GMatrix, A) = G(size(A, 1)) + A
+(A, G::GMatrix) = G + A

-(G::GMatrix, A) = G + (-A)
-(A, G::GMatrix) = -(G - A)

*(G::GMatrix, A) = G.λ * (A .- sum(A, dims=1) / size(A, 1))
*(A, G::GMatrix) = (G*A')'

adjoint(G::GMatrix) = G

"""
    HMatrix{T}

`H` matrix, as defined by Aitchison 1986. See also [`H`](@ref).
"""
struct HMatrix{T}
  λ::T
end

(H::HMatrix{T})(d::Integer) where {T} = H.λ * I(d) +  H.λ * J(d)

*(λ::Number, H::HMatrix) = HMatrix(H.λ * λ)

+(H::HMatrix, A) = H(size(A, 1)) + A
+(A, H::HMatrix) = H + A

-(H::HMatrix, A) = H + (-A)
-(A, H::HMatrix) = -(H - A)

*(H::HMatrix, A) = H.λ * (A .+ sum(A, dims=1))
*(A, H::HMatrix) = H.λ * (A .+ sum(A, dims=2))

inv(H::HMatrix) = H⁻¹Matrix(inv(H.λ))

"""
    H⁻¹Matrix{T}

Lazy inverse of `HMatrix{T}`.
"""
struct H⁻¹Matrix{T}
  λ::T
end

(H⁻¹::H⁻¹Matrix{T})(d::Integer) where {T} = H⁻¹.λ * I(d) - (H⁻¹.λ / (d+1)) * J(d)

*(λ::Number, H⁻¹::H⁻¹Matrix) = H⁻¹Matrix(H⁻¹.λ * λ)

+(H⁻¹::H⁻¹Matrix, A) = H⁻¹(size(A, 1)) + A
+(A, H⁻¹::H⁻¹Matrix) = H⁻¹ + A

-(H⁻¹::H⁻¹Matrix, A) = H⁻¹ + (-A)
-(A, H⁻¹::H⁻¹Matrix) = -(H⁻¹ - A)

*(H⁻¹::H⁻¹Matrix, A) = H⁻¹.λ * (A .- sum(A, dims=1) / (size(A, 1) + 1))
*(A, H⁻¹::H⁻¹Matrix) = H⁻¹.λ * (A .- sum(A, dims=2) / (size(A, 2) + 1))

inv(H⁻¹::H⁻¹Matrix) = HMatrix(inv(H⁻¹.λ))

"""
    J
    J(d)

User interface for [`JMatrix`](@ref), a `d` by `d` matrix of ones.

## Examples

```
julia> J(3)
julia> J*v
julia> v'*J
```
"""
const J = JMatrix{Bool}(true)

"""
    F
    F(d)

User interface for [`FMatrix`](@ref), as defined by Aitchison.

`F` is a `d` by `D` matrix that can be defined as

`F[i, j] = 1`, if `i==j`

`F[i, j] = -1`, if `j==D`

`F[i, j] = 0`, otherwise

## Examples

```
julia> F(3)
julia> F*v
julia> v'*F
```
"""
const F = FMatrix{Int}(1)

"""
    G
    G(D)

User interface for [`GMatrix`](@ref), as defined by Aitchison.

`G` is an `D` by `D` matrix that can be defined as

`G[i, j] = I[i, j] - J[i, j] / D`

## Examples

```
julia> G(3)
julia> G*v
julia> v'*G
```
"""
const G = GMatrix{Float64}(1.0)

"""
    H
    H(d)

User interface for [`HMatrix`](@ref), as defined by Aitchison.

`H` is a `d` by `d` matrix that can be defined as

`H[i, j] = I[i, j] + J[i, j]`

## Examples

```
julia> H(3)
julia> H*v
julia> v'*H
```
"""
const H = HMatrix{Int}(1)


text:
using Utilities

@time for i=1:1000
    mat = rand(10,10).>0.5
     mat = putZeroDiag(mat)
     f = transitivityDirBin_funErgm(mat)
     i,j = (2,1)
     delta = transitivityDirBin_funErgm_delta(mat,i,j)
     mat[i,j] = false ; f1 =  transitivityDirBin_funErgm(mat);
     mat[i,j] = true ; f2 =  transitivityDirBin_funErgm(mat);
     (f2 - f1) != delta ? error():()
end

function pseudoLogLikelihood(θ::Real,A::BitArray{2},statF::statErgmFun)
     N1,N2 = size(A)
     statVal = statF.f(A) #the value of the statistics on the real networks
     # compute a matrix of changes in the statistics when each element is toggled from 0 to 1
     # come lo faccio in modo efficiente??
     deltaFun(l,m) = statF.delta(A,l,m)
     logPsLike = 0
     for i=1:N1,j=1:N2
         if i!=j
            π_ij = 1/(1+ exp(- θ * deltaFun(i,j) ))
            if A[i,j] #is non zero
                logPsLike = logPsLike + log(π_ij)
            else
                logPsLike = logPsLike + log( 1 - π_ij )
            end
        end
     end
     return logPsLike
 end

testStat = statErgmFun(transitivityDirBin_funErgm,transitivityDirBin_funErgm_delta)

mat = rand(10,10).>0.5
pseudoLogLikelihood(0.0,mat,testStat)


text:
using Documenter
using LeapSeconds

makedocs(
    sitename = "LeapSeconds",
    format = Documenter.HTML(),
    modules = [LeapSeconds]
)

deploydocs(
    repo = "github.com/JuliaTime/LeapSeconds.jl.git"
)


text:
module ElementarySymmetricFunctions
    import SpecialFunctions.logabsbinomial

    include("util.jl")
    include("esf.jl")
    include("poisbin.jl")

    function esf(x::AbstractVector{T}; scale::Symbol=:identity, method=:recursive, k::D=2) where {T <: Real, D <: Integer}
        if scale == :identity
            if method == :recursive
                return esf_sum(x)
            elseif method == :DCFFT
                return esf_dc_fft(x,k)
            elseif method == :FFT
                return poisbin_fft(x ./ (1 .+ x)) ./ prod(1 ./(1 .+ x))
            else
                error("method $method not implemented with scaling $scale")
            end
        elseif scale == :norm
            p = x ./ (1 .+ x)
            if method == :recursive
                return poisbin_sum_taub(p)
            elseif method == :DCFFT
                return poisbin_dc_fft(p,k)
            elseif method == :FFT
                return poisbin_fft(p)
            else
                error("method $method not implemented with scaling $scale")
            end
        elseif scale == :log
            if method == :recursive
        end
    end

    export esf_sum, esf_sum_reg, esf_dc_fft, esf_dc_fft_reg, esf_sum_log
    export esf_sum_dervs_1, esf_sum_dervs_2, esf_sum_dervs_1_reg, esf_sum_dervs_2_reg
    export esf_diff, esf_diff_reg, esf_diff_updown
    export poisbin_sum_taub, poisbin_fft, poisbin_chen, poisbin_dc_fft, poisbin_fft_cf, poisbin_sum_taub_log
end


text:
write_test_output = false
path = dirname(@__FILE__)

if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

###################################################################
# Test: scalar_reshape()
###################################################################
s1 = SMC.scalar_reshape(1, 0)
s2 = SMC.scalar_reshape([1, 0], [1, 0])
###################################################################
@testset "Scalar reshape" begin
    @test s1 == Array{Float64,1}[[1.0], [0.0]]
    @test s2 == Array{Float64,1}[[1.0, 0.0], [1.0, 0.0]]
end

###################################################################
# Test: vector_reshape()
###################################################################
v1 = SMC.vector_reshape(1.0, 2.0)
v2 = SMC.vector_reshape([1, 0], 1, 0)
v3 = SMC.vector_reshape([1.0, 0.0], 1.0, 0.0)

if write_test_output
    JLD2.jldopen(string("reference/vector_reshape_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["v1"] = v1
        file["v2"] = v2
        file["v3"] = v3
    end
end

test_v1, test_v2, test_v3 = JLD2.jldopen(string("reference/vector_reshape_version=", ver, ".jld2"), "r") do file
    file["v1"], file["v2"], file["v3"]
end
###################################################################
@testset "Vector reshape" begin
    @test v1 == test_v1
    @test v2 == test_v2
    @test v3 == test_v3
end

###################################################################
# Test: scalar_reduce()
###################################################################
s1_r = SMC.scalar_reduce([s1]...)
s2_r = SMC.scalar_reduce([s2 for i in 1:5]...)

if write_test_output
    JLD2.jldopen(string("reference/scalar_reduce_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["s1_r"] = s1_r
        file["s2_r"] = s2_r
    end
end

test_s1_r, test_s2_r = JLD2.jldopen(string("reference/scalar_reduce_version=", ver, ".jld2"), "r") do file
    file["s1_r"], file["s2_r"]
end
###################################################################
@testset "Scalar reduce" begin
    @test s1_r == test_s1_r
    @test s2_r == test_s2_r
end

###################################################################
# Test: vector_reduce()
###################################################################
v1_r = SMC.vector_reduce([v1]...)
v2_r = SMC.vector_reduce([v2 for i in 1:5]...)
v3_r = SMC.vector_reduce([v3 for i in 1:5]...)

if write_test_output
    JLD2.jldopen(string("reference/vector_reduce_version=", ver, ".jld2"), true, true, true, IOStream) do file
        file["v1_r"] = v1_r
        file["v2_r"] = v2_r
        file["v3_r"] = v3_r
    end
end

test_v1_r, test_v2_r, test_v3_r = JLD2.jldopen(string("reference/vector_reduce_version=", ver, ".jld2"), "r") do file
    file["v1_r"], file["v2_r"], file["v3_r"]
end

###################################################################
@testset "Vector reduce" begin
    @test v1_r == test_v1_r
    @test v2_r == test_v2_r
    @test v3_r == test_v3_r
end

###################################################################
# Test: <, min, max
###################################################################
@testset "<, min, max" begin
    @test 1 < 2
    @test 1 < Complex(2, 1)
    @test Complex(0.95, 0.5) < 1
    @test Complex(0.5, 20) < Complex(0.99, 0.5)

    @test min(1, 2) == 1
    @test min(Complex(0.99, 0.5), 1) == 0.99
    @test min(1, Complex(20, 0.5)) == 1
    @test min(Complex(0.99, 0.5), Complex(20, 0.5)) == 0.99

    @test max(1, 2) == 2
    @test max(Complex(0.99, 0.5), 1) == 1
    @test max(1, Complex(20, 0.5)) == 20
    @test max(Complex(0.99, 0.5), Complex(20, 0.5)) == 20
end


text:
using ModelingToolkit
using Test

# Define some variables
@Param t σ ρ β
@Unknown x(t) y(t) z(t)
@Deriv D'~t # Default of first derivative, Derivative(t,1)
@Const c=0

# Define a differential equation
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
de = DiffEqSystem(eqs,[t],[x,y,z],[σ,ρ,β])
ModelingToolkit.generate_ode_function(de)
ModelingToolkit.generate_ode_function(de;version=ModelingToolkit.SArrayFunction)
jac_expr = ModelingToolkit.generate_ode_jacobian(de)
jac = ModelingToolkit.calculate_jacobian(de)
f = ODEFunction(de)
ModelingToolkit.generate_ode_iW(de)

# Differential equation with automatic extraction of variables
de2 = DiffEqSystem(eqs, [t])

function test_vars_extraction(de, de2)
    for el in (:ivs, :dvs, :ps)
        names2 = sort(collect(var.name for var in getfield(de2,el)))
        names = sort(collect(var.name for var in getfield(de,el)))
        @test names2 == names
    end
end
test_vars_extraction(de, de2)

# Conversion to first-order ODEs #17
@Deriv D3'''~t
@Deriv D2''~t
@Unknown u(t) u_tt(t) u_t(t) x_t(t)
eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
       D2(x) ~ D(x) + 2]
de = DiffEqSystem(eqs, [t])
de1 = ode_order_lowering(de)
lowered_eqs = [D(u_tt) ~ 2u_tt + u_t + x_t + 1
               D(x_t)  ~ x_t + 2
               D(u_t)  ~ u_tt
               D(u)    ~ u_t
               D(x)    ~ x_t]
function test_eqs(eqs1, eqs2)
    length(eqs1) == length(eqs2) || return false
    eq = true
    for (eq1, eq2) in zip(eqs1, eqs2)
        lhs1, lhs2 = eq1.lhs, eq2.lhs
        typeof(lhs1) === typeof(lhs2) || return false
        for f in fieldnames(typeof(lhs1))
            eq = eq & isequal(getfield(lhs1, f), getfield(lhs2, f))
        end
        eq = eq & isequal(eq1.rhs, eq2.rhs)
    end
    eq
end
@test_broken test_eqs(de1.eqs, lowered_eqs)

# Internal calculations
a = y - x
eqs = [D(x) ~ σ*a,
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
de = DiffEqSystem(eqs,[t],[x,y,z],[σ,ρ,β])
ModelingToolkit.generate_ode_function(de)
jac = ModelingToolkit.calculate_jacobian(de)
f = ODEFunction(de)

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs,[x,y,z],[t,σ,ρ,β])
ns2 = NonlinearSystem(eqs)
for el in (:vs, :ps)
    names2 = sort(collect(var.name for var in getfield(ns2,el)))
    names = sort(collect(var.name for var in getfield(ns,el)))
    @test names2 == names
end

ModelingToolkit.generate_nlsys_function(ns)

@Deriv D'~t
@Param A B C
_x = y / C
eqs = [D(x) ~ -A*x,
       D(y) ~ A*x - B*_x]
de = DiffEqSystem(eqs,[t],[x,y],[A,B,C])
test_vars_extraction(de, DiffEqSystem(eqs,[t]))
test_vars_extraction(de, DiffEqSystem(eqs))
@test eval(ModelingToolkit.generate_ode_function(de))([0.0,0.0],[1.0,2.0],[1,2,3],0.0) ≈ -1/3

# Now nonlinear system with only variables
@Unknown x y z
@Param σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs)
nlsys_func = ModelingToolkit.generate_nlsys_function(ns)
jac = ModelingToolkit.generate_nlsys_jacobian(ns)
f = @eval eval(nlsys_func)

# Intermediate calculations
# Define a nonlinear system
eqs = [a ~ y-x,
       0 ~ σ*a,
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])
nlsys_func = ModelingToolkit.generate_nlsys_function(ns)
jac = ModelingToolkit.calculate_jacobian(ns)
jac = ModelingToolkit.generate_nlsys_jacobian(ns)


text:
export BrickletOLED128x64DisplayConfiguration
struct BrickletOLED128x64DisplayConfiguration
    contrast::Integer
    invert::Bool
end

export BrickletOLED128x64Identity
struct BrickletOLED128x64Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletOLED128x64
"""
3.3cm (1.3") OLED display with 128x64 pixels
"""
mutable struct BrickletOLED128x64 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletOLED128x64(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_oled_128x64")
        deviceInternal = package.BrickletOLED128x64(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export write
"""
    $(SIGNATURES)

Appends 64 byte of data to the window as set by :func:`New Window`.

Each row has a height of 8 pixels which corresponds to one byte of data.

Example: if you call :func:`New Window` with column from 0 to 127 and row
from 0 to 7 (the whole display) each call of :func:`Write` (red arrow) will
write half of a row.

.. image:: /Images/Bricklets/bricklet_oled_128x64_display.png
   :scale: 100 %
   :alt: Display pixel order
   :align: center
   :target: ../../_images/Bricklets/bricklet_oled_128x64_display.png

The LSB (D0) of each data byte is at the top and the MSB (D7) is at the
bottom of the row.

The next call of :func:`Write` will write the second half of the row
and the next two the second row and so on. To fill the whole display
you need to call :func:`Write` 16 times.
"""
function write(device::BrickletOLED128x64, data)
    device.deviceInternal.write(data)
end

export new_window
"""
    $(SIGNATURES)

Sets the window in which you can write with :func:`Write`. One row
has a height of 8 pixels.
"""
function new_window(device::BrickletOLED128x64, column_from, column_to, row_from, row_to)
    device.deviceInternal.new_window(column_from, column_to, row_from, row_to)
end

export clear_display
"""
    $(SIGNATURES)

Clears the current content of the window as set by :func:`New Window`.
"""
function clear_display(device::BrickletOLED128x64)
    device.deviceInternal.clear_display()
end

export set_display_configuration
"""
    $(SIGNATURES)

Sets the configuration of the display.

You can set a contrast value from 0 to 255 and you can invert the color
(black/white) of the display.
"""
function set_display_configuration(device::BrickletOLED128x64, contrast, invert)
    device.deviceInternal.set_display_configuration(contrast, invert)
end

export get_display_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Display Configuration`.
"""
function get_display_configuration(device::BrickletOLED128x64)
    return device.deviceInternal.get_display_configuration()
end

export write_line
"""
    $(SIGNATURES)

Writes text to a specific line with a specific position.
The text can have a maximum of 26 characters.

For example: (1, 10, "Hello") will write *Hello* in the middle of the
second line of the display.

You can draw to the display with :func:`Write` and then add text to it
afterwards.

The display uses a special 5x7 pixel charset. You can view the characters
of the charset in Brick Viewer.
"""
function write_line(device::BrickletOLED128x64, line, position, text)
    device.deviceInternal.write_line(line, position, text)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletOLED128x64)
    return device.deviceInternal.get_identity()
end


text:
# 10/24/2020 © Shuowen Chen
# This script simulates CAPM data using Tauchen & Hussey (1991) method
# as in Wright (2003) and explores sensitivity of Euler equation estimation
# to number of discretization points.

# load packages
using Random, Distributions, LinearAlgebra, Statistics, FastGaussQuadrature, Optim
# set seed
Random.seed!(123)

#======= Defining Functions =======#
# univariate tauchen
function tauchen(m, ρ, N, σ²)
    # INPUTS:
    # 1. m: multiple of unconditional standard deviation
    # 2. ρ: AR(1) coefficient
    # 3. N: number of grid points
    # 4. σ²: variance of the white noise

    # OUTPUTS:
    # 1. Transition probability matrix
    # 2. Sequence of grid points

    # unconditional standard deviation of the AR(1)
    std_un = sqrt(σ²/(1 - ρ^2));
    # max and min of the grid
    ymax = m*std_un;
    ymin = -ymax;
    # length of adjacent grid point
    w = (ymax - ymin)/(N - 1);
    # lay out the grids
    ygrid = collect(ymin:w:ymax);
    # Transition probability matrix (N by N)
    P = zeros(N, N);
    P[:, 1] = cdf.(Normal(0, 1), ((ymin .- ρ*ygrid .+ w/2)./sqrt(σ²)));
    P[:, end] = 1 .- cdf.(Normal(0, 1), ((ymax .- ρ*ygrid .- w/2)./sqrt(σ²)));
    for k = 2:(N - 1)
        P[:, k] = cdf.(Normal(0, 1), (ygrid[k] .- ρ*ygrid .+ w/2)./sqrt(σ²)) -
        cdf.(Normal(0, 1), (ygrid[k] .- ρ*ygrid .- w/2)./sqrt(σ²));
    end
    return(transition = P, grid = ygrid);
end

# Univariate Tauchen and Hussey (1991)
function tauchenhussey(ρ, a, σ², N)
    # PURPOSE:
    # Discretize the following AR(1) process:
    # z_{t+1} = a + ρ*z_{t} + epsilon_{t+1}
    # INPUTS:
    # 1. ρ: AR(1) coefficient
    # 2. a: intercept
    # 3. σ²: variance of epsilon
    # 4. N: number of grid points

    # OUTPUTS:
    # 1. P: transition probability matrix
    # 2. grid: vector of grid points

    # A function from FastGaussQuadrature
    (x, w) = gausshermite(N);
    # Gaussian nodes and weights for normal distribution
    # Acknowledgement: I benefit from reading slides by Karen Kopecky
    # http://www.karenkopecky.net/Teaching/eco613614/Notes_DiscretizingAR1s.pdf

    # unconditional mean of the process
    μ = a / (1 - ρ);
    # Establish the grid using Gaussian-Hermite node
    grid = x .* sqrt(2*σ²) .+ μ;
    # Rescale the Gaussian-Hermite weight
    w = w ./ sqrt(pi);
    # Compute transition probability matrix
    P = zeros(N, N);
    for i = 1:N
        for j = 1:N
            EZprime = a + ρ*grid[i];
            # Gauchen-Hussey weight
            P[i, j] = w[j] * pdf(Normal(EZprime, sqrt(σ²)), grid[j]) /
            pdf(Normal(μ, sqrt(σ²)), grid[j]);
        end
    end
    # Normalization
    for i = 1:N
        P[i, :] = P[i, :] ./ sum(P[i, :])
    end
    return(P, grid);
end

# An auxiliary function to be called in tauchenhusseyvar
function productrule(N, EZ, μ, Σ, grid_j, weight_j, k)
    # Purpose:
    # Compute individual transition prob from state i to j and aggregate
    # by Gauss product rule (tensor product)

    # Inputs:
    # 1. N: number of grid points
    # 2. EZ: expected value in state i (k by 1)
    # 3. μ: unconditional mean of VAR process (k by 1)
    # 4. Σ: var-cov matrix of VAR(1) process
    # 5. grid_j: grid points for the next state
    # 6. weight_j: Gauss-Hermite weight for the next state
    # 7. k: number of variables in the VAR

    # Outputs:
    # 1. prob: aggregate transition probability from state i to j

    den_ind = zeros(1, k);
    for i = 1:k
        den_ind[i] = weight_j/sqrt(pi) * pdf(Normal(EZ[i], diag(sqrt(Σ))[i]),
        grid_j[i]) / pdf(Normal(μ[i], diag(sqrt(Σ))[i]), grid_j[i])
    end
    den_agg = prod(den_ind);
    return(den_agg);
end

# Bivariate Tauchen-Hussey (1991)
function tauchenhusseyvar(k, A, B, Σ, N)
    # PURPOSE: Discretize the following VAR(1) process
    # Y_t = A + B * Y_t-1 + e_t,
    # where e_t ∼ N(0, Σ)
    # Acknowledgement: I benefit from reading the quadrature code by Tauchen,
    # which is available in his website.

    # INPUTS:
    # 1. k: number of variables in the VAR
    # 2. A: constant of VAR(1) (k by 1)
    # 3. B: VAR(1) coefficient matrix (k by k)
    # 4. Σ: iid var-cov matrix
    # 5. N: number of discrete points for each variable k in the VAR

    # OUTPUTS:
    # 1. P: transition matrix (Nᵏ by Nᵏ)
    # 2. Grid: grid points (k by Nᵏ).
    # Note: number of states in the system is equal to number of grids * number of lags,
    # since we consider VAR(1), number of states = number of grid points

    # unconditional mean of the process
    μ = inv(I(k) - B) * A;
    # obtain Gaussian-Hermite nodes and weights
    (x, w) = gausshermite(N);
    # Grid for each variable
    grid_ind = repeat(x, outer = [1, k])' .* sqrt.(2 * diag(Σ)) .+ μ;
    # total grid (k by N^k) and GH weights (N^k) using product rule
    grid = zeros(k, N^k);
    w_tot = zeros(k, N^k);
    for i = 1:k
        grid[i, :] = repeat(grid_ind[i, :], inner = [N^(k-i), 1],
                            outer = [N^(i-1), 1]);
        w_tot[i, :] = repeat(w, inner = [N^(k-i), 1], outer = [N^(i-1), 1]);
    end
    w_tot = prod(w_tot, dims = 1);
    # Compute transition probability matrix
    P = zeros(N^k, N^k);
    for i = 1:N^k
        for j = 1:N^k
            # expected value in state i at time t-1
            EZprime = A + B*grid[:, i]; # k by 1
            # Gauchen-Hussey
            P[i, j] = productrule(N, EZprime, μ, Σ, grid[:, j], w_tot[j], k);
        end
    end
    # Normalization
    for i = 1:N^k
        P[i, :] = P[i, :] ./ sum(P[i, :])
    end
    return(grid, P);
end

# Simulate data using Markov transition probability and grid points
function simulationMC(Pr_mat, grid, T, sim)
    # INPUTS:
    # 1. Pr_mat: transition matrix
    # 2. grid: grid points
    # 3. T: number of observations in each simulation
    # 4. sim: number of simulations
    # 5. k: number of variables

    # OUTPUT:
    # 1. Z: simulated states using MC (k by T for # sim times)
    # Economic meaning: log consumption and dividend growth data
    # 2. state: state index for stock return simulation

    # first obs of the MC (middle state)
    initial_state = convert(Int, floor((length(grid)/size(grid)[1]+1)/2));
    # placeholder for simulated states
    Z = zeros(k, T, sim);
    # initial state
    Z[:, 1, :] .= grid[:, initial_state]; # need to broadcast

    # placeholder for state index
    state = zeros(T, sim);
    # initial state
    state[1, :] .= initial_state;

    # This follows from Adda and Cooper (2003)
    # standard uniformly distributed rv for
    # determining the new state
    u = rand(Uniform(), (T-1, sim));
    # store ∑_{l=1}^{j}π_{i,l} for j = 1,...,N
    Pr_sum = cumsum(Pr_mat, dims = 2);
    old_state = initial_state;
    for i = 1:sim
        for t = 2:T
            new_state = findfirst(u[t-1, i] .<= Pr_sum[old_state, :])
            Z[:, t, i] = grid[:, new_state];
            # store the state index
            state[t, i] = new_state;
            # move to the next obs
            old_state = new_state;
        end
    end
    state = convert.(Int, state);
    return(Z, state);
end

# a function that directly simulates the VAR
function simulationVAR(T, sim, A, B, Σ, k)
    # INPUTS:
    # 1. T: number of observations in each simulation
    # 2. sim: number of simulations
    # 3. A: constant of VAR(1) (k by 1)
    # 4. B: VAR(1) coefficient matrix (k by k)
    # 5. Σ: iid var-cov matrix
    # 6. K: number of variables

    # OUTPUTS:
    # 1. TS: simulated time series (k by 1)

    # normal errors e_t ∼ N(0, Σ) for t = 2,...,T and each simulation
    err = zeros(k, T-1, sim);
    for i = 1:sim
        err[:, :, i] = rand(MvNormal([0; 0], Σ), T-1);
    end
    # placeholder for simulated VAR
    TS = zeros(k, T, sim);
    # initial observation (steady state)
    for i = 1:sim
        TS[:, 1, i] = inv(I(k) - B) * A;
    end
    # Then use the relationship Y_t = A + B*Y_t-1 + e_t
    for i = 1:sim
        for j = 2:T
            # Note the index for err is due to specification of err
            TS[:, j, i] = A + B*TS[:, j-1, i] + err[:, j-1, i];
        end
    end
    return(TS)
end

# Simulate stock return
function simulateReturn(β, γ, P, CG, DG, N, stateindex)
    # INPUTS:
    # 1. β: discount factor
    # 2. γ: CRRA preference
    # 3. P: probability transition matrix
    # 4. CG: consumption growth
    # 5. DG: dividend growth
    # 6. N: number of grid points for each variable

    # OUTPUTS:
    # 1. smreturn: simulated return data

    # Define transition matrix for return
    B = β*P;
    for j = 1:N^2
        B[:, j] = B[:, j] .* CG[j].^(-γ) .* DG[j];
    end
    # State contingent price-dividend ratio (N^2 by 1)
    psm = inv(I(N^2) - B) * B * ones(N^2);
    # Realized stock return from state i to j (N^2 by N^2)
    Rmat = zeros(N^2, N^2);
    for i = 1:N^2
        for j = 1:N^2
            Rmat[i, j] = DG[j]*(1+psm[j])/psm[i];
        end
    end
    # simulate the return data
    T = size(stateindex)[1] - 1;
    smreturn = zeros(T, 1);
    for t = 1:T
        smreturn[t] = Rmat[stateindex[t], stateindex[t+1]];
    end
    return(smreturn);
end

# A function to conduct 2 step GMM
function GMM2step(moments, gradients, k, initial, samplesize)
    # INPUTS:
    # 1. moments: an anonymous function that maps theta to moment conditions
    # 2. gradients: an anoymous function that mapes theta to gradient conditions
    # 3. k: number of moments
    # 4. initial: initial values for estimation (must be floating)
    # 5. samplesize: sample size of the moments
    # OUTPUTS:
    # 1. est: parameter estimate
    # 2. se: standard errors

    # sample average moments
    sample_m = theta -> vec(mean(moments(theta), dims = 1));
    # first step weighting matrix
    W = I(k);
    # objective function in the first step
    obj = theta -> sample_m(theta)'*W*sample_m(theta);
    # minimization using limited BFGS
    firststep = Optim.optimize(obj, initial, LBFGS());
    # second step weighting matrix
    W2 = inv(cov(moments(firststep.minimizer)));
    # objective function in the second step
    obj2 = theta -> sample_m(theta)'*W2*sample_m(theta);
    secondstep = Optim.optimize(obj2, initial, LBFGS());
    # compute the gradient
    # dimension of D: k (# of pars) by p (# of moments)
    D = reshape(mean(gradients(secondstep.minimizer), dims = 1), k, length(initial));
    var = inv(D'*inv(cov(moments(secondstep.minimizer)))*D);
    # compute standard errors
    se = sqrt.(diag(var)./samplesize);
    return(firststep.minimizer, secondstep.minimizer, se);
end

#======= Main Program Starts Here =======#
# Setting parameters
β = 1.139; # utility discount factor
γ = 13.7; # CRRA parameter
k = 2; # number of variables
A = [0.021; 0.04]; # Calibrated values as in Wright (2003)
B = [-0.161 0.017; 0.414 0.117];
Σ = [0.0012 0.00177; 0.00177 0.014];
N = 4; # number of grids for each variable
T = 400; # number of time periods
sim = 400; # number of simulations

# Generate grid points and transition probability
(grid, P) = tauchenhusseyvar(k, A, B, Σ, N);
# Simulate log consumption and dividend growth
(timeseries, stateindex) = simulationMC(P, grid, T, sim);

# Consumption and dividend growth
CG_tot = exp.(timeseries[1, :, :]);
DG_tot = exp.(timeseries[2, :, :]);
# use state index to generate return stock market return
R = zeros(T-1, sim);
for i = 1:sim
    R[:, i] = simulateReturn(β, γ, P, CG_tot[:, i], DG_tot[:, i], N, stateindex[:, i]);
end

# conduct estimation
container = zeros(sim, 2);
for i = 1:sim
    # intercept, lagged return and lagged consumption growth
    ins = [ones(T-2) R[1:end-1, i] CG_tot[1:end-2, i]];
    # sample size of the moments
    t_mom = size(instruments)[1];
    # Construct Euler equation
    euler = theta -> (theta[1] .* R[2:end, i] .* CG_tot[2:end-1, i] .^ (-theta[2]) .-1) .* ins;

    grad = theta -> hcat(R[2:end, i] .* CG_tot[2:end-1, i] .^ (-theta[2]) .* ins,
                             -(theta[1] .* R[2:end, i] .* CG_tot[2:end-1, i] .^ (-theta[2]) .* log.(CG_tot[2:end-1, i])) .* ins);

    gest = GMM2step(euler, grad, 3, [0.0, 13.0], t_mom);
    container[i,:] = gest[2];
end

mean(container, dims =1)

max_stars_repo_name: shuowencs/Robust-Discretization

max_issues_repo_name: shuowencs/Robust-Discretization

max_forks_repo_name: shuowencs/Robust-Discretization


text:
###############################################################################
#
#  From relative to absolute
#
###############################################################################

function _relative_to_absoluteQQ(L::NfRelNS{nf_elem}, auts::Vector{NfRelNSToNfRelNSMor_nf_elem})
  K, gK = number_field(NfAbsNS, L)
  Ks, mKs = simplified_simple_extension(K, isabelian = true, cached = false)
  #Now, I have to translate the automorphisms.
  #First, to automorphisms of K.
  autsK = Vector{NfAbsNSToNfAbsNS}(undef, length(auts))
  Qxy = parent(K.pol[1])
  for i = 1:length(auts)
    embs = Vector{NfAbsNSElem}(undef, ngens(K))
    imgs = image_generators(auts[i])
    for j = 1:length(imgs)
      embs[j] = K(map_coefficients(FlintQQ, imgs[j].data, parent = Qxy))
    end
    autsK[i] = hom(K, K, embs, check = false)
  end
  #Final step: translate the automorphisms to Ks
  Hecke._assure_has_inverse_data(mKs)
  autsKs = Vector{NfToNfMor}(undef, length(autsK))
  for i = 1:length(autsK)
    autsKs[i] = hom(Ks, Ks, mKs\(image_primitive_element(mKs*autsK[i])), check = false)
  end
  return Ks, autsKs
end

function _relative_to_absolute(L::NfRelNS{nf_elem}, auts::Vector{NfRelNSToNfRelNSMor_nf_elem})
  Ks, mKs = simplified_absolute_field(L)
  Hecke._assure_has_inverse_data(mKs)
  #Now, I have to translate the automorphisms.
  #First, to automorphisms of K.
  autsKs = Vector{NfToNfMor}(undef, length(auts))
  imggen = mKs(gen(Ks))
  for i = 1:length(auts)
    autsKs[i] = hom(Ks, Ks, mKs\(auts[i](imggen)))
  end
  #finally, the embedding of the base field of L into Ks
  k = base_field(L)
  img = mKs\(L(gen(k)))
  embed = hom(k, Ks, img, check = false)
  return Ks, autsKs, embed
end


text:
using Test
using JMPS
include("test_function.jl")
# include("test_JMPS.jl")
# include("test_CUDA.jl")

max_stars_repo_name: qiyang-ustc/JMPS

max_issues_repo_name: qiyang-ustc/JMPS

max_forks_repo_name: qiyang-ustc/JMPS


text:
using Random, Distributions, StatsPlots, Plots; pyplot()
Random.seed!(0)

mu = 20
d1, d2 = Normal(mu,mu), Exponential(mu)

n = 100
data1 = rand(d1,n)
data2 = rand(d2,n)

qqnorm(data1, c=:blue, ms=3, msw=0, label="Normal Data")
qqnorm!(data2, c=:red, ms=3, msw=0, label="Exponential Data",
        xlabel="Normal Theoretical Quantiles",
        ylabel="Data Quantiles", legend=true)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
## Pure saturation conditions solver
function x0_bubble_pressure(model,T,x)

    #strategy: supose virial max pressure as the saturation pressure
    #suppose ideal gas and calculate y0

    #TODO
    #on sufficiently large temps,
    #the joule-thompson inversion occurs
    #making the virial coeff positive
    #on those cases, use an strategy that supposes pure gas on that side
    #Pbi = inf
    #xi = 0

    #check each T with T_scale, if treshold is over, replace Pi with inf

    T_threshold = T_scales(model,x)
    @show T_threshold
    replaceP = ifelse.(T_threshold .< T,true,false)
    @show replaceP

    eachx = eachcol(Diagonal(ones(eltype(x),length(x))))
    Bi = second_virial_coefficient.(model,T,eachx)
    #using P_B(2B) as a sat aproximation
    #z = 1 + B/v
    #P_B = RT/v(1+B/v)
    #P_B(2B) = -RT/2B(1-B/2B)
    #P_B(2B) = -0.25*RT/B

    P_Bi = @. -0.25*R̄*T/Bi
    #=xP0 = yP
    #dot(x,P0) = P
    P = dot(x,P0)
    =#
    P = zero(T)
    for i in 1:length(x)
        if !replaceP[i]
            P+=x[i]*P_Bi[i]
        end
    end

    #@show P_Bi
    #P = dot(x,P_Bi)
    y = @. x*P_Bi/P
    ysum = 1/∑(y)
    return y .* ysum
    #return y

    #=
    lb_v = lb_volume(model,x)
    k0 = 10.0
    y0    = k0 .*x./(1 .+x.*(k0 .- 1))
    y0    = y0 ./sum(y0)
    v0    = [log10(lb_v/0.45),
    log10(lb_v/1e-4)]
    return append!(v0,y0[1:end-1])
    =#
end

## Mixture saturation solver
function bubble_pressure(model, T, x; v0 =nothing)
    TYPE = promote_type(eltype(T),eltype(x))
    lb_v = lb_volume(model,x)
    ts = T_scales(model,x)
    pmix = p_scale(model,x)
    if v0 === nothing
        v0 = x0_bubble_pressure(model,T,x)
    end
    f! = (F,z) -> Obj_bubble_pressure(model, F, T, exp10(z[1]), exp10(z[2]), x,@view(z[3:end]),ts,pmix)
    r  =Solvers.nlsolve(f!,v0,TrustRegion(Newton(),NTR()))
    v_l = exp10(r.info.zero[1])
    v_v = exp10(r.info.zero[2])
    y = FractionVector(r.info.zero[3:end])
    P_sat = pressure(model,v_l,T,x)
    return (P_sat, v_l, v_v, y)
end

function Obj_bubble_pressure(model, F, T, v_l, v_v, x, y,ts,ps)
    y = FractionVector(y) #julia magic, check misc.jl
    μ_l = vt_chemical_potential(model,v_l,T,x)
    μ_v = vt_chemical_potential(model,v_v,T,y)
    p_l = pressure(model,v_l,T,x)
    p_v = pressure(model,v_v,T,y)
    for i in 1:length(x)
        F[i] = (μ_l[i]-μ_v[i])/(R̄*ts[i])
    end
    F[end] = (p_l-p_v)/ps
end

    #j! = (J,z) -> Jac_bubble_pressure(model, J, T, exp10(z[1]), exp10(z[2]), x[i,:], z[3:end])
    #=
    _y0 = collect(FractionVector(v0[3:end]))
    @show _y0
    @show _pl0 = pressure(model,exp10(v0[1]),T,x)
    @show _pv0 = pressure(model,exp10(v0[2]),T,_y0)
    _P = 0.5*(_pl0+_pv0)
    for _ in 1:3
        _y0,_P = rr_bubble_pressure_refine(model,x,_y0,_P,T)
    end
    =#
function bubble_pressure_rr(model, T, x; P = 40000)
    sol0 = x0_bubble_pressure(model,T,x)
    vl0 = exp10(sol0[1])
    vv0 = exp10(sol0[2])

    @show y0 = collect(FractionVector(sol0[3:end]))
    @show vl0 = volume(model,P,T,x,phase=:l)
    @show vv0 = volume(model,P,T,y0,phase=:v)
    @show μ_l = vt_chemical_potential(model,vl0,T,x)
    @show μ_v = vt_chemical_potential(model,vv0,T,y0)
    y1 =  μ_l ./ μ_v .* x
    @show y1 = y1 ./ sum(y1)
    @show pl0 = pressure(model,vl0,T,x)
    @show pv0 = pressure(model,vv0,T,y1)

    @show P = (pl0 - pv0)/(log(pl0) - log(pv0))

    @show vl = volume(model,P,T,x,phase=:l)
    @show vv = volume(model,P,T,y1,phase=:v)
    #=
    μ_l = vt_chemical_potential(model,vl,T,x)
    μ_v = vt_chemical_potential(model,vv,T,y)
    K = log.(μ_v) ./ log.(μ_l)
    y = K .* x
    y = y./sum(y)
    pl = pressure(model,vl,T,x)
    pv = pressure(model,vv,T,y)
    P = (pl+pv)/2
    @show vl,vv,y
    return y,P
    =#
end

export bubble_pressure


text:
# reusing data loaded by test/hurdle.jl
X=Matrix{Float64}(bioChemists[!,[:femWomen,:marMarried,:kid5,:phd,:ment]])
Xwconst=[ones(size(X,1)) X]
y=convert(Array{Float64,1},bioChemists.art)
p = size(Xwconst,2)

###########################################################
# Xpos == Xzero
###########################################################
@testset "inrep with Xpos == Xzero" begin

coefsR1=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_coefsR1.csv"), DataFrame)))
yhatR1=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_yhatR1.csv"), DataFrame)))
yhatR1partial=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_yhatR1partial.csv"), DataFrame)))

# simple inrep with GLM underlying
increpfit = fit(InclusionRepetition,GeneralizedLinearModel,Xwconst,y)
showres = IOBuffer()
show(showres, increpfit)
showstr = String(take!(copy(showres)))
@test occursin("InclusionRepetition regression",showstr)
@test occursin("Positive part coefficients",showstr)
@test occursin("Zero part coefficients",showstr)

coefsJ=vcat(coef(increpfit)...)
@test coefsJ[1+p:end] ≈ coefsR1[1+p:end] rtol=1e-6
yhatJ=predict(increpfit, Xwconst)
@test yhatJ ≈ yhatR1 rtol=1e-2
yhatJpartial=predict(increpfit, Xwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR1partial rtol=1e-2

# same simple inrep through UNREGULATED lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,X,y;λ=[0.0, 0.01])
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+p:end] ≈ coefsR1[1+p:end] rtol=1e-4
# rdist(coefsJ,coefsR1)
yhatJ = predict(increpglrfit, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=1e-2
yhatJpartial=predict(increpglrfit, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=1e-2
yhatJ = predict(increpglrfit, X; select=AllSeg())
@test size(yhatJ) == (size(X,1),2)

# regulated lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,X,y; γ=0.0)
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+p:end] ≈ coefsR1[1+p:end] rtol=0.4
# rdist(coefsJ[1+p:end],coefsR1[1+p:end])
yhatJ = predict(increpglrfit, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=0.05
yhatJpartial=predict(increpglrfit, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=0.05

coefsJpos, coefsJzero = coef(increpglrfit;select=AllSeg())
@test size(coefsJpos,1) == 6
@test size(coefsJzero,1) == 6

showres = IOBuffer()
show(showres, increpglrfit)
showstr = String(take!(copy(showres)))
@test occursin("InclusionRepetition regression",showstr)
@test occursin("Positive part regularization path",showstr)
@test occursin("Zero part regularization path",showstr)

# this one throws an error because we did not specify the same λ vector for both submodels so they have different lengths
@test_throws AssertionError predict(increpglrfit, X; select=AllSeg())

Random.seed!(1)
coefsJCVmin=vcat(coef(increpglrfit, MinCVKfold{MinCVmse}(5))...)
@test coefsJCVmin[1+p:end] ≈ coefsR1[1+p:end] rtol=0.30
# rdist(coefsJCVmin,coefsR1)

# try with SharedArray
Xs = convert(SharedArray,X)
increpglrfitShared = fit(InclusionRepetition,GammaLassoPath,Xs,y; γ=0.0)
coefsJShared=vcat(coef(increpglrfitShared;select=MinAICc())...)
@test coefsJShared ≈ coefsJ
yhatJ = predict(increpglrfitShared, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=0.05
yhatJpartial=predict(increpglrfitShared, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=0.05

# regulated gamma lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,X,y; γ=10.0)
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+p:end] ≈ coefsR1[1+p:end] rtol=0.2
# rdist(coefsJ[1+p:end],coefsR1[1+p:end])
yhatJ = predict(increpglrfit, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=0.05
yhatJpartial=predict(increpglrfit, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=0.05

# using ProfileView
# Profile.init(delay=0.001)
# Profile.clear()
# @profile increpglrfit2 = fit2(InclusionRepetition,GammaLassoPath,X,y; γ=10.0);
# ProfileView.view()
# Profile.print()

# using DataFrames formula interface
increpfit = fit(InclusionRepetition,GeneralizedLinearModel,
    @formula(art ~ femWomen + marMarried + kid5 + phd + ment), bioChemists)
coefsJ=vcat(coef(increpfit)...)
@test coefsJ[1+p:end] ≈ coefsR1[1+p:end] rtol=1e-6
yhatJ=predict(increpfit, Xwconst)
@test yhatJ ≈ yhatR1 rtol=1e-2
yhatJpartial=predict(increpfit, Xwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR1partial rtol=1e-2

end

###########################################################
# Xpos ≠ Xzero
###########################################################
@testset "inrep with Xpos ≠ Xzero" begin

coefsR2=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_coefsR2.csv"), DataFrame)))
yhatR2=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_yhatR2.csv"), DataFrame)))
yhatR2partial=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_yhatR2partial.csv"), DataFrame)))

Xpos = X[:,1:3]
Xzero = X[:,4:5]
Xzerowconst=[ones(size(X,1)) Xzero]
Xposwconst=[ones(size(X,1)) Xpos]
ppos = size(Xposwconst,2)
# ixpos = y.>0
# ypos = y[ixpos]
# countmap(ypos)

# simple inrep with GLM underlying
increpfit = fit(InclusionRepetition,GeneralizedLinearModel,Xzerowconst,y; Xpos=Xposwconst)
coefsJ=vcat(coef(increpfit)...)
@test coefsJ[1+ppos:end] ≈ coefsR2[1+ppos:end] rtol=1e-5
yhatJ=predict(increpfit, Xzerowconst; Xpos=Xposwconst)
@test yhatJ ≈ yhatR2 rtol=1e-2
yhatJpartial=predict(increpfit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR2partial rtol=1e-2

# same simple inrep through UNREGULATED lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,Xzero,y; Xpos=Xpos, λ=[0.0,0.0001])
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+ppos:end] ≈ coefsR2[1+ppos:end] rtol=1e-4
# rdist(coefsJ,coefsR2)
yhatJ=predict(increpglrfit, Xzero; Xpos=Xpos, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=1e-2
yhatJpartial=predict(increpglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=1e-2
yhatJ = predict(increpglrfit, Xzero; Xpos=Xpos, select=AllSeg())
@test size(yhatJ) == (size(X,1),2)

# regulated lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,Xzero,y; Xpos=Xpos, γ=0.0)
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+ppos:end] ≈ coefsR2[1+ppos:end] rtol=0.3
# rdist(coefsJ,coefsR2)
yhatJ = predict(increpglrfit, Xzero; Xpos=Xpos, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=0.05
yhatJpartial=predict(increpglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=0.05

coefsJpos, coefsJzero = coef(increpglrfit;select=AllSeg())
@test size(coefsJpos,1) == 4
@test size(coefsJzero,1) == 3

# try with SharedArray
Xzeros = convert(SharedArray,Xzero)
Xposs = convert(SharedArray,Xpos)
increpglrfitShared = fit(InclusionRepetition,GammaLassoPath,Xzeros,y; Xpos=Xposs, γ=0.0)
coefsJShared=vcat(coef(increpglrfitShared;select=MinAICc())...)
@test coefsJShared ≈ coefsJ
yhatJ = predict(increpglrfitShared, Xzeros; Xpos=Xposs, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=0.05
yhatJpartial=predict(increpglrfitShared, Xzeros[ixpartial,:]; Xpos=Xposs[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=0.05

# regulated gamma lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,Xzero,y; Xpos=Xpos, γ=10.0)
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+ppos:end] ≈ coefsR2[1+ppos:end] rtol=0.3
# rdist(coefsJ,coefsR2)
yhatJ = predict(increpglrfit, Xzero; Xpos=Xpos, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=0.05
yhatJpartial=predict(increpglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=0.05

# using DataFrames formula interface
increpfit = fit(InclusionRepetition,GeneralizedLinearModel,@formula(art ~ phd + ment), bioChemists; fpos = @formula(art ~ femWomen + marMarried + kid5))
coefsJ=vcat(coef(increpfit)...)
@test coefsJ[1+ppos:end] ≈ coefsR2[1+ppos:end] rtol=1e-5
# rdist(coefsJ,coefsR2)
yhatJ=predict(increpfit, Xzerowconst; Xpos=Xposwconst)
@test yhatJ ≈ yhatR2 rtol=1e-3
yhatJpartial=predict(increpfit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR2partial rtol=1e-3

end

###########################################################
# Xpos ≠ Xzero
###########################################################
@testset "inrep with Xpos ≠ Xzero AND offset specified" begin

coefsR3=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_coefsR3.csv"), DataFrame)))
yhatR3=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_yhatR3.csv"), DataFrame)))
yhatR3partial=vec(Matrix{Float64}(CSV.read(joinpath(testdir,"data","hurdle_yhatR3partial.csv"), DataFrame)))

offpos = convert(Vector{Float64},bioChemists.offpos)
offzero = convert(Vector{Float64},bioChemists.offzero)
Xpos = X[:,1:3]
Xzero = X[:,4:5]
Xzerowconst=[ones(size(X,1)) Xzero]
Xposwconst=[ones(size(X,1)) Xpos]
ppos = size(Xposwconst, 2)
# ixpos = y.>0
# ypos = y[ixpos]
# countmap(ypos)

# simple inrep with GLM underlying
increpfit = fit(InclusionRepetition,GeneralizedLinearModel,Xzerowconst,y; Xpos=Xposwconst, offsetzero=offzero, offsetpos=offpos)
coefsJ=vcat(coef(increpfit)...)
@test coefsJ[1+ppos:end] ≈ coefsR3[1+ppos:end] rtol=1e-5
yhatJ=predict(increpfit, Xzerowconst; Xpos=Xposwconst, offsetzero=offzero, offsetpos=offpos)
@test yhatJ ≈ yhatR3 rtol=0.1
yhatJpartial=predict(increpfit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial])
@test yhatJpartial ≈ yhatR3partial rtol=0.2

# same simple inrep through UNREGULATED lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,Xzero,y; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, λ=[0.0,0.0001])
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+ppos:end] ≈ coefsR3[1+ppos:end] rtol=1e-4
# rdist(coefsJ,coefsR3)
yhatJ=predict(increpglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=0.1
yhatJpartial=predict(increpglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=0.2
yhatJ = predict(increpglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=AllSeg())
@test size(yhatJ) == (size(X,1),2)

# regulated lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,Xzero,y; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, γ=0.0)
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+ppos:end] ≈ coefsR3[1+ppos:end] rtol=0.10
# rdist(coefsJ,coefsR3)
yhatJ = predict(increpglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=0.1
yhatJpartial=predict(increpglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=0.2

coefsJpos, coefsJzero = coef(increpglrfit;select=AllSeg())
@test size(coefsJpos,1) == 4
@test size(coefsJzero,1) == 3

# try with SharedArray
Xzeros = convert(SharedArray,Xzero)
Xposs = convert(SharedArray,Xpos)
offzeros = convert(SharedArray,offzero)
offposs = convert(SharedArray,offpos)
increpglrfitShared = fit(InclusionRepetition,GammaLassoPath,Xzeros,y; Xpos=Xposs, offsetzero=offzeros, offsetpos=offposs, γ=0.0)
coefsJShared=vcat(coef(increpglrfitShared;select=MinAICc())...)
@test coefsJShared ≈ coefsJ
yhatJShared = predict(increpglrfitShared, Xzeros; Xpos=Xposs, offsetzero=offzeros, offsetpos=offposs, select=MinAICc())
@test yhatJShared ≈ yhatJ
yhatJSharedpartial=predict(increpglrfitShared, Xzeros[ixpartial,:]; Xpos=Xposs[ixpartial,:], offsetzero=offzeros[ixpartial], offsetpos=offposs[ixpartial], select=MinAICc())
@test yhatJSharedpartial ≈ yhatJpartial

# regulated gamma lasso path
increpglrfit = fit(InclusionRepetition,GammaLassoPath,Xzero,y; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, γ=10.0)
coefsJ=vcat(coef(increpglrfit;select=MinAICc())...)
@test coefsJ[1+ppos:end] ≈ coefsR3[1+ppos:end] rtol=0.10
# rdist(coefsJ,coefsR3)
yhatJ = predict(increpglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=0.1
yhatJpartial=predict(increpglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=0.2

# using DataFrames formula interface
increpfit = fit(InclusionRepetition,GeneralizedLinearModel,@formula(art ~ phd + ment), bioChemists; fpos = @formula(art ~ femWomen + marMarried + kid5), offsetzero=offzero, offsetpos=offpos)
coefsJ=vcat(coef(increpfit)...)
@test coefsJ[1+ppos:end] ≈ coefsR3[1+ppos:end] rtol=1e-5
# rdist(coefsJ,coefsR3)
yhatJ=predict(increpfit, Xzerowconst; Xpos=Xposwconst, offsetzero=offzero, offsetpos=offpos)
@test yhatJ ≈ yhatR3 rtol=0.1
yhatJpartial=predict(increpfit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial])
@test yhatJpartial ≈ yhatR3partial rtol=0.2

end

###########################################################
# degenrate cases
###########################################################
@testset "inrep degenerate cases" begin

# degenerate positive counts data case 1
include(joinpath(testdir,"data","degenerate_hurdle_1.jl"))
inrep = fit(InclusionRepetition,GammaLassoPath,X,y; showwarnings=true)
coefsJ=vcat(coef(inrep;select=MinAICc())...)
@test vec(coefsJ) ≈ [-4.8941, 0.0, -6.04112, 0.675767] rtol=1e-4
coefsJpos, coefsJzero = coef(inrep;select=AllSeg())
@test size(coefsJpos,1) == 2
@test size(coefsJzero,1) == 2

ydeg = zero(y)
ydeg[1] = 2
ydeg[2] = 3
increpglm = @test_logs (:warn, r"failed to fit truncated counts model to positive") fit(InclusionRepetition,GeneralizedLinearModel,[ones(size(X,1)) X],ydeg; showwarnings=true)

# degenerate positive counts data case 1 without >1
y0or1 = deepcopy(y)
y0or1[y.>1] .= 1
inrep = @test_logs (:warn, r"ypos has no elements larger than 0") fit(InclusionRepetition,GammaLassoPath,X,y0or1; showwarnings=true)
coefs0or1=vcat(coef(inrep;select=MinAICc())...)
@test coefs0or1 ≈ [0.0, 0.0, -6.04112, 0.675767] rtol=1e-4
coefs0or1pos, coefs0or1zero = coef(inrep;select=AllSeg())
@test all(iszero,coefs0or1pos)
@test coefs0or1zero == coefsJzero

# degenerate positive counts data case 1 without zeros
y0or1 = deepcopy(y)
y0or1[y.==0] .= 1
y0or1[1] = 3.0

inrep = @test_logs (:warn, r"I\(y\) is all ones") fit(InclusionRepetition,GammaLassoPath,X,y0or1; verbose=true, showwarnings=true)
# coefs0or1=vcat(coef(inrep;select=MinAICc())...)
coefs0or1pos, coefs0or1zero = coef(inrep;select=MinAICc())
@test coefs0or1pos ≈ [-10.2903, 0.820934] rtol=0.1
@test all(iszero,coefs0or1zero)
coefs0or1pos, coefs0or1zero = coef(inrep;select=AllSeg())
@test all(iszero,coefs0or1zero)

increpglm = @test_logs (:warn, r"I\(y\) is all ones") fit(InclusionRepetition,GeneralizedLinearModel,[ones(size(X,1)) X],y0or1; verbose=true, showwarnings=true)

# degenerate positive counts data case 1 with only zeros
y0or1 = zero(y)
@test_throws ErrorException fit(InclusionRepetition,GammaLassoPath,X,y0or1)
@test_logs (:warn, r"I\(y\) is all zeros") @test_throws ErrorException fit(InclusionRepetition,GammaLassoPath,X,y0or1; verbose=true, showwarnings=true)

# degenerate positive counts data case 2
include(joinpath(testdir,"data","degenerate_hurdle_2.jl"))
inrep = fit(InclusionRepetition,GammaLassoPath,X,y; showwarnings=true)
coefsJ=vcat(coef(inrep;select=MinAICc())...)
@test vec(coefsJ) ≈ [-3.26059, -0.220401, -5.3013, 0.185418] rtol=1e-4
coefsJpos, coefsJzero = coef(inrep;select=AllSeg())
@test size(coefsJpos,1) == 2
@test size(coefsJzero,1) == 2

# degenerate positive counts data case 3
include(joinpath(testdir,"data","degenerate_hurdle_3.jl"))
inrep = @test_logs (:warn, r"ypos has no elements larger than 0") fit(InclusionRepetition,GammaLassoPath,X,y; showwarnings=true)
coefsJ=vcat(coef(inrep;select=MinAICc())...)
# @test vec(coefsJ) ≈ [0.0, 0.0, -4.54363, 0.000458273] rtol=1e-4
@test length(coefsJ) == 4
coefsJpos, coefsJzero = coef(inrep;select=AllSeg())
Matrix( coefsJzero)
@test size(coefsJpos,1) == 2
@test coefsJpos == zero(coefsJpos)
@test size(coefsJzero,1) == 2

# this test case used to give numerical headaches to devresid(PositivePoisson(),...)
include(joinpath(testdir,"data","degenerate_hurdle_5.jl"))
inrep = fit(InclusionRepetition,GammaLassoPath,X,y; Xpos=Xpos, offset=offset)
@test !(ismissing(inrep.mpos))
@test !(ismissing(inrep.mzero))
@test typeof(inrep.mpos) <: GammaLassoPath
@test typeof(inrep.mzero) <: GammaLassoPath

end