text:
event_type(::Union{Val{XCB_KEY_PRESS}, Val{XCB_KEY_RELEASE}}) = xcb_key_press_event_t
event_type(::Union{Val{XCB_BUTTON_PRESS}, Val{XCB_BUTTON_RELEASE}}) = xcb_button_press_event_t
event_type(::Union{Val{XCB_ENTER_NOTIFY}, Val{XCB_LEAVE_NOTIFY}}) = xcb_enter_notify_event_t
event_type(::Val{XCB_MOTION_NOTIFY}) = xcb_motion_notify_event_t
event_type(::Val{XCB_EXPOSE}) = xcb_expose_event_t
event_type(::Val{XCB_CLIENT_MESSAGE}) = xcb_client_message_event_t
event_type(::Val{XCB_CONFIGURE_NOTIFY}) = xcb_configure_notify_event_t
event_type(::Val{85}) = xcb.xcb_xkb_state_notify_event_t # very hacky, but response type 85 is emitted instead of XCB_XKB_STATE_NOTIFY...
event_type(::Val{XCB_KEYMAP_NOTIFY}) = xcb_keymap_notify_event_t
event_type(rt) = nothing

response_type(data) = Int(data.response_type & 0x7f)

is_from_send_request(data) = Bool(data.response_type & 0x80)

function unsafe_load_event(xge_ptr; warn_unknown=false)
    xge = unsafe_load(xge_ptr)
    rt = response_type(xge)
    et = event_type(Val(rt))
    if isnothing(et)
        warn_unknown && @warn "Unknown event $(rt)"
        nothing
    else
        unsafe_load(convert(Ptr{et}, xge_ptr))
    end
end

window_id_field(event::xcb_button_press_event_t) = :event
window_id_field(event::xcb_key_press_event_t) = :event
window_id_field(event::xcb_enter_notify_event_t) = :event
window_id_field(event::xcb_motion_notify_event_t) = :event
window_id_field(event::xcb_expose_event_t) = :window
window_id_field(event::xcb_client_message_event_t) = :window
window_id_field(event::xcb_configure_notify_event_t) = :window

window_id(event) = getproperty(event, window_id_field(event))

location(event::xcb_button_press_event_t) = (event.event_x, event.event_y)
location(event::xcb_key_press_event_t) = (event.event_x, event.event_y)
location(event::xcb_enter_notify_event_t) = (event.event_x, event.event_y)
location(event::xcb_motion_notify_event_t) = (event.event_x, event.event_y)
location(event::xcb_expose_event_t) = (event.x, event.y)
location(event::xcb_configure_notify_event_t) = (event.x, event.y)

EventDetails(win::XCBWindow, data::EventData, event, t) =
    EventDetails(data, Int.(location(event)), t, win)
EventDetails(wh::XWindowHandler, win::XCBWindow, data::xcb_button_press_event_t, t) =
    EventDetails(win, MouseEvent(data), data, t)

function EventDetails(wh::XWindowHandler, win::XCBWindow, data::xcb_key_press_event_t, t)
    keycode_symbol = key_name(wh.keymap, data.detail)
    key_symbol = KeySymbol(wh.keymap, data.detail)
    input_char = Char(wh.keymap, data.detail)
    event_type = response_type(data) == XCB_KEY_PRESS ? KeyPressed() : KeyReleased()
    EventDetails(win, KeyEvent(keycode_symbol, key_symbol, input_char, KeyModifierState(data), event_type), data, t)
end

EventDetails(wh::XWindowHandler, win::XCBWindow, data::xcb_enter_notify_event_t, t) =
    EventDetails(win, response_type(data) == XCB_ENTER_NOTIFY ? PointerEntersWindowEvent() : PointerLeavesWindowEvent(), data, t)
EventDetails(wh::XWindowHandler, win::XCBWindow, data::xcb_motion_notify_event_t, t) =
    EventDetails(win, PointerMovesEvent(), data, t)
EventDetails(wh::XWindowHandler, win::XCBWindow, data::xcb_expose_event_t, t) =
    EventDetails(win, ExposeEvent((data.width, data.height), data.count), data, t)
EventDetails(wh::XWindowHandler, win::XCBWindow, data::xcb_configure_notify_event_t, t) =
    EventDetails(win, ResizeEvent((data.width, data.height)), data, t)

process_xevent(wh, xge::Nothing, t; warn_unknown=false, kwargs...) = nothing

function process_xevent(wh::XWindowHandler, xge, t; warn_unknown=false, kwargs...)
    event = unsafe_load_event(xge; warn_unknown)
    if !isnothing(event) # unknown event
        win = get_window(wh, event)
        if event isa xcb_client_message_event_t
            ed_8 = Int.(event.data.data8)
            event_data32_1 = ed_8[1] + ed_8[2] * 2^8 + ed_8[3] * 2^16 + ed_8[4] *2^24
            event_data32_1 == win.delete_request && throw(CloseWindow(win, ""))
        elseif event isa xcb_xkb_state_notify_event_t
            xkb_state_update_mask(wh.keymap.state, event.baseMods, event.latchedMods, event.lockedMods, event.baseGroup, event.latchedGroup, event.lockedGroup)
        elseif event isa xcb_keymap_notify_event_t
            wh.keymap = Keymap(wh.conn; setup_xkb=false)
        elseif !isnothing(win) # event happened on inexistant window
            details = EventDetails(wh, win, event, t)
            execute_callback(callbacks(wh, win), details; kwargs...)
        end
    end
end

function listen_for_events(wh::XWindowHandler, t0, next_event::Function; warn_unknown=false, on_iter_first=() -> nothing, on_iter_last=() -> nothing, kwargs...)
    while !isempty(wh.windows)
        try
            on_iter_first()
            xge = next_event(wh)
            t = time() - t0
            process_xevent(wh, xge, t; warn_unknown, kwargs...)
            on_iter_last()
        catch e
            if e isa WindowException
                win_callbacks = callbacks(wh, e.win)
                e isa CloseWindow && win_callbacks.on_close(wh, e)
                e isa InvalidWindow && win_callbacks.on_invalid(wh, e)
            else
                rethrow(e)
                break
            end
        end
    end
end

"""
Run an `EventLoop` attached to a `XWindowHandler` instance.
"""
function Base.run(wh::XWindowHandler, ::Synchronous; warn_unknown=false, poll=false, kwargs...)
    t0 = time()
    next_event = poll ? poll_for_event : wait_for_event
    listen_for_events(wh, t0, next_event; warn_unknown, kwargs...)
end

function Base.run(wh::XWindowHandler, ::Asynchronous; warn_unknown=false, kwargs...)
    t0 = time()
    @async listen_for_events(wh, t0, poll_for_event; warn_unknown, kwargs...)
end


text:
#=
   https://dtai.cs.kuleuven.be/problog/tutorial/basic/02_bayes.html
   """
   Since the random variables in the Bayesian network are all Boolean, we only need a single literal
   in the head of the rules. We can extend the Bayesian network to have a multi-valued variable by
   indicating the severity of the earthquake. The literal earthquake now has three possible values
   none, mild, heavy instead of previously two (no or yes).
   """

   Cf ~/blog/alarm_multi.blog
      ~/psi/alarm_multi.psi
      ~/webppl/alarm_multi.wppl

    cplint program:

    person(john).
    person(mary).

    0.7::burglary.
    0.01::earthquake(heavy); 0.19::earthquake(mild); 0.8::earthquake(none).

    0.90::alarm :-   burglary, earthquake(heavy).
    0.85::alarm :-   burglary, earthquake(mild).
    0.80::alarm :-   burglary, earthquake(none).
    0.10::alarm :- \+burglary, earthquake(mild).
    0.30::alarm :- \+burglary, earthquake(heavy).

    0.8::calls(X) :- alarm, person(X).
    0.1::calls(X) :- \+alarm, person(X).

    evidence(calls(john),true).
    evidence(calls(mary),true).

    query(burglary).
    query(earthquake(_)).

=#

#=
julia> @time include("alarm_multi.jl")
...

Summary Statistics
  parameters      mean       std   naive_se      mcse          ess      rhat
      Symbol   Float64   Float64    Float64   Float64      Float64   Float64

       alarm    1.0000    0.0000     0.0000    0.0000          NaN       NaN
    burglary    0.9885    0.1067     0.0005    0.0009   13653.3301    1.0001
    calls[1]    1.0000    0.0000     0.0000    0.0000          NaN       NaN
    calls[2]    1.0000    0.0000     0.0000    0.0000          NaN       NaN
  earthquake    2.7675    0.4532     0.0023    0.0037   14275.1136    1.0003

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       alarm    1.0000    1.0000    1.0000    1.0000    1.0000
    burglary    1.0000    1.0000    1.0000    1.0000    1.0000
    calls[1]    1.0000    1.0000    1.0000    1.0000    1.0000
    calls[2]    1.0000    1.0000    1.0000    1.0000    1.0000
  earthquake    2.0000    3.0000    3.0000    3.0000    3.0000

earthquake status: heavy=1 mild=2 none=3]
Distributions of variable (num:20)
3.00000 =>   31240  (0.781000)
2.00000 =>    8221  (0.205525)
1.00000 =>     539  (0.013475)
  5.424711 seconds (36.39 M allocations: 2.780 GiB, 23.75% gc time)

julia>
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function alarm_multi()

    heavy = 1
    mild  = 2
    none  = 3
    earthquake_status = [heavy,mild,none]

    john = 1
    mary = 2
    people = [john,mary]
    num_people = length(people)

    burglary ~ flip(0.7)
    earthquake ~ Categorical([0.01,0.19,0.8])

    alarm ~ if     burglary && earthquake  == heavy flip(0.9)
            elseif burglary && earthquake  == mild  flip(0.85)
            elseif burglary && earthquake  == none  flip(0.80)
            elseif !burglary && earthquake == heavy flip(0.30)
            elseif !burglary && earthquake == mild  flip(0.10)
            else flip(0.0)
            end

    calls = tzeros(num_people)
    for p in people
        if alarm
            calls[p] ~ flip(0.8)
        else
            calls[p] ~ flip(0.01)
        end
    end

    true ~ Dirac(calls[john] == true)
    true ~ Dirac(calls[mary] == true)

    # true ~ Dirac(calls[john] == true)
    # true ~ Dirac(calls[mary] == false)

    return earthquake
end

model = alarm_multi()

num_chains = 4

# chains = sample(model, Prior(), 10_000)
# chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
# chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(), 10_000)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)

chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 10_000)
# chains = sample(model, IS(), 10_000)
#
#chains = sample(model, Gibbs(HMC(0.1,5,:a,:b),PG(15,:p)), 10_000)
# chains = sample(model, Gibbs(NUTS(1000,0.65,:a,:b),PG(15,:p)), 10_000)
# chains = sample(model, Gibbs(HMC(0.1,5,:a,:b),SMC(1000,:p)), 10_000) # Nope

display(chains)
# display(plot(chains))

println("earthquake status: heavy=1 mild=2 none=3] ")
genq = generated_quantities(model, chains)
show_var_dist_pct(genq,20)

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
const intptr_t = UInt
const time_t = Int

const AV_NOPTS_VALUE = 0x8000000000000000 % Int64
AV_VERSION_INT(a, b, c) = (a << 16 | b << 8) | c
AV_VERSION(a,b,c) = nothing
const AVPROBE_SCORE_MAX = 100
MKTAG(a,b,c,d) = (UInt32(a) | (UInt32(b) << 8) | (UInt32(c) << 16) | UInt32(d) << 24)
FFERRTAG(a,b,c,d) = -MKTAG(a,b,c,d)
macro AV_PIX_FMT_NE(be, le)
    Symbol("AV_PIX_FMT_"*string(le))
end


text:
struct TunnelingModifier{T,F,V}
    f::F
    x::V
    fbest::Ref{T}
end

function tunneling(f,lb,ub,x0)
    N = length(ub)
    # Relevant configuration
    tolf=1e-8
    options = OptimizationOptions(g_abstol=5e-8)
    prob = ADScalarObjective(f,x0)
    # Minimisation phase
    prob = OptimizationProblem(obj=scalarobj,
    bounds=(lb,ub); inplace=true)
    opt_min.xtol_rel = 1e-8
    res =  NLSolvers.solve(prob, x0, ModifiedActiveBox(Newton();epsilon=tolbounds), options)

    best_f = res.info.minimum
    opt_x  = []
    best_x = res.info.minimizer
    append!(opt_x,[min_x])

    # Tunneling phase
    opt_tun = NLopt.Opt(:LD_MMA, length(ub))
    opt_tun.lower_bounds = lb
    opt_tun.upper_bounds = ub
    opt_tun.xtol_rel = 1e-8
    opt_tun.stopval = -1e-6

    for i in 1:10*N
        T0 = x -> (f(x)-f_best)*prod(exp(1e-1/sqrt(sum((x[i]-x_opt[j][i])^2 for i in 1:N))) for j in 1:k)
        T  = (x,g) -> NLopt_obj(T0,x,g)
        opt_tun.min_objective =  T

        r  = 2.0.*rand(Float64,(N)).-1.0
        ϵ1 = 2*(tolf)^(1/5)*(1+norm(best_x,2))
        x0 = r/norm(r,2).*ϵ1+best_x
        x0 = ub.*(x0.>=ub)+lb.*(x0.<=lb)+x0.*(ub.>x0.>lb)

        # Tunneling
        (new_f,new_x,status) = NLopt.optimize(opt_tun, x0)
        if minimum(res.info.∇fz) <=
            println(i)
            break
        end
        # Minimisation
        (min_f,min_x,status) = NLopt.optimize(opt_min, new_x)
        if min_f<best_f
            best_f = min_f
            best_x = min_x
            opt_x  = []
            append!(opt_x,[min_x])
        else
            append!(opt_x,[min_x])
        end
    end
    return (best_f,best_x)
end

function NLopt_obj(f,x,g)
        if length(g) > 0
            df = DiffResults.GradientResult(x)
            df = ForwardDiff.gradient!(df,f,x)
            g .= DiffResults.gradient(df)
            return DiffResults.value(df)
        else
            return f(x)
        end
end


text:
using MCMCChains, Test
using Statistics: std

@testset "Summarize to DataFrame tests" begin
    val = rand(1000, 8, 4)
    chns = Chains(val,
                ["a", "b", "c", "d", "e", "f", "g", "h"],
                Dict(:internals => ["c", "d", "e", "f", "g", "h"]),
                sorted=true)

    parm_df = summarize(chns, sections=[:parameters])

    @test 0.48 < parm_df[:a, :mean][1] < 0.52
    @test names(parm_df) == [:parameters, :mean, :std, :naive_se, :mcse, :ess, :r_hat]

    all_sections_df = summarize(chns, sections=[:parameters, :internals])
    @test all_sections_df[:,:parameters] == ["a", "b", "c", "d", "e", "f", "g", "h"]
    @test size(all_sections_df) == (8, 7)

    two_parms_two_funs_df = summarize(chns[[:a, :b]], mean, std)
    @test two_parms_two_funs_df[:, :parameters] == ["a", "b"]
    @test size(two_parms_two_funs_df) == (2, 3)

    three_parms_df = summarize(chns[[:a, :b, :c]], mean, std, sections=[:parameters, :internals])
    @test three_parms_df[:, :parameters] == ["a", "b", "c"]
    @test size(three_parms_df) == (3, 3)

    three_parms_df_2 = summarize(chns[[:a, :b, :g]], mean, std,
    sections=[:parameters, :internals], func_names=[:mean, :sd])
    @test three_parms_df_2[:, :parameters] == ["a", "b", "g"]
    @test size(three_parms_df_2) == (3, 3)
end


text: include("./sequencing/sequencingMain.jl")


text:
############################################################################################
## This file is part of RLinearAlgebra.jl
##
## Overview: abstractions and methods for stopping a randomized linear system solver
##
## Contents
## - Abstract Types
## - `check_stop_criterion` Docstring
## - Stopping Criteria
## - Export Statements
##
############################################################################################

# Dependencies:

#############################################
# Abstract Types
#############################################

"""
    LinSysStopCriterion

Abstract supertype for specifying stopping criteria for randomized linear solver.
"""
abstract type LinSysStopCriterion end

#############################################
# `check_stop_criterion` Docstring
#############################################

#############################################
# Stopping Criteria
#############################################
include("linear_solver_stops/stop_max_iter.jl")

#############################################
# Export Statements
#############################################
# See RLinearAlgebra.jl
#export LinSysStopCriterion
#export LSStopMaxIterations


text:
using .Gym
env = make("CartPole-v0", :human_pane)

actions = [sample(env._env.action_space) for i=1:1000]
i = 1
done = false
reset!(env)
while i <= length(actions)
    global i, done
    a, b, done, d = step!(env, actions[i])

    i += 1
end


text:
# =========================================================================
# Function related to terms, variable names and I/O

_diff(t::NTuple{N, T}) where {N, T} = ntuple(i->t[i + 1] - t[i], N - 1)
_diffn(t::NTuple{N, T}) where {N, T} = ntuple(i->t[i] - t[i + 1], N - 1)

"""
    coefnames(<model>, anova::Val{:anova})
    coefnames(<term>, anova::Val{:anova})

Customize coefnames for anova.
"""
StatsBase.coefnames(t::MatrixTerm, anova::Val{:anova}) = mapreduce(coefnames, vcat, t.terms, repeat([anova], length(t.terms)))
StatsBase.coefnames(t::FormulaTerm, ::Val{:anova}) = (coefnames(t.lhs), coefnames(t.rhs))
StatsBase.coefnames(::InterceptTerm{H}, ::Val{:anova}) where H = H ? "(Intercept)" : []
StatsBase.coefnames(t::ContinuousTerm, ::Val{:anova}) = string(t.sym)
StatsBase.coefnames(t::CategoricalTerm, ::Val{:anova}) = string(t.sym)
StatsBase.coefnames(t::FunctionTerm, ::Val{:anova}) = string(t.exorig)
StatsBase.coefnames(ts::StatsModels.TupleTerm, anova::Val{:anova}) = reduce(vcat, coefnames.(ts, anova))
StatsBase.coefnames(t::InteractionTerm, anova::Val{:anova}) = begin
    join(coefnames.(t.terms, anova), " & ")
end

StatsBase.coefnames(aov::AnovaResult) = coefnames(aov.model, Val(:anova))

# Base.show(io::IO, t::FunctionTerm) = print(io, "$(t.exorig)")

"""
    getterms(<term>)

Get an array of `Expr` or `Symbol` of terms.
"""
getterms(::AbstractTerm) = []
getterms(term::Union{Term, CategoricalTerm, ContinuousTerm}) = [term.sym]
getterms(::InterceptTerm{true}) = [Symbol(1)]
getterms(term::InteractionTerm) = mapreduce(getterms, union, term.terms)
getterms(term::FunctionTerm) = [term.exorig]
getterms(term::StatsModels.TupleTerm) = mapreduce(getterms, union, term)
getterms(term::MatrixTerm) = getterms(term.terms)
#=
"""
    getterm(term::AbstractTerm)
    getterm(term::FunctionTerm)

Get the `Symbol` of a single term.
"""=#
# getterm(term::AbstractTerm) = term.sym
# getterm(term::FunctionTerm) = term.exorig
# getterm(term::InterceptTerm) = Symbol(1)

# Determine selected terms for type 2 ss
"""
    isinteract(f::MatrixTerm, id1::Int, id2::Int)

Determine if `f.terms[id2]` is an interaction term of `f.terms[id1]` and other terms.
"""
isinteract(f::MatrixTerm, id1::Int, id2::Int) = issubset(getterms(f.terms[id1]), getterms(f.terms[id2]))

"""
    selectcoef(f::MatrixTerm, id::Int)
    selectcoef(f::MatrixTerm, ::Val{N})
    selectcoef(f::MatrixTerm, ::Val{1})

A set of index of `f.terms` which are interaction terms of `f.terms[id]` and other terms.
"""
selectcoef(f::MatrixTerm, id::Int) = Set([comp for comp in eachindex(f.terms) if isinteract(f, id, comp)])
selectcoef(f::MatrixTerm, ::Val{N}) where N = selectcoef(f, N)
selectcoef(f::MatrixTerm, ::Val{1}) = Set(eachindex(f.terms)) # For InterceptTerm

# Create sub-formula
"""
    subformula(lhs::AbstractTerm, rhs::MatrixTerm, id::Int; reschema::Bool = false)
    subformula(lhs::AbstractTerm, rhs::MatrixTerm, id; reschema::Bool = false)
    subformula(lhs::AbstractTerm, rhs::NTuple{N, AbstractTerm}, id::Int)

Create formula from existing `lhs` and `rhs` truncated to `1:id` or excluded collection `id`.
"""
subformula(lhs::AbstractTerm, rhs::MatrixTerm, id::Int; reschema::Bool = false) =
    id > 0 ? reschema_formula(lhs, rhs.terms[1:id], reschema) : reschema_formula(lhs, (InterceptTerm{false}(),), reschema)
# For nestedmodels & GLM

function subformula(lhs::AbstractTerm, rhs::MatrixTerm, id; reschema::Bool = false)
    terms = rhs.terms[setdiff(eachindex(rhs.terms), id)]
    1 in id && (terms = (InterceptTerm{false}(), terms...))
    reschema_formula(lhs, terms, reschema)
end
# For type 2/3 & GLM

subformula(lhs::AbstractTerm, rhs::NTuple{N, AbstractTerm}, id::Int) where N =
    id > 0 ? FormulaTerm(lhs, (collect_matrix_terms(first(rhs).terms[1:id]), rhs[2:end]...)) : FormulaTerm(lhs, (collect_matrix_terms((InterceptTerm{false}(),)), rhs[2:end]...))
# For nestedmodels & MixedModels

"""
    extract_contrasts(f::FormulaTerm)

Extract a dictionary of contrasts.
"""
extract_contrasts(f::FormulaTerm) =
    Dict{Symbol, Any}(t.sym => t.contrasts.contrasts for t in f.rhs.terms if isa(t, CategoricalTerm))

"""
    clearschema(<terms with schema>) = <terms without schema>

Clear any applied schema on terms.
"""
clearschema(::InterceptTerm{true}) = ConstantTerm(1)
clearschema(::InterceptTerm{false}) = ConstantTerm(0)
clearschema(t::FunctionTerm) = t
clearschema(t::Union{CategoricalTerm, ContinuousTerm}) = Term(t.sym)
clearschema(t::InteractionTerm) = InteractionTerm(clearschema.(t.terms))

# reschema only happen when using TupleTerm rather than MatrixTerm
reschema_formula(lhs::AbstractTerm, ts::StatsModels.TupleTerm, reschema::Bool) =
    reschema ? FormulaTerm(clearschema(lhs), clearschema.(ts)) : FormulaTerm(lhs, collect_matrix_terms(ts))

# test name
tname(::Type{FTest}) = "F test"
tname(::Type{LRT}) = "Likelihood-ratio test"
#tname(M::AnovaStatsRao) = "Rao score test"
#tname(M::AnovaStatsCp) = "Mallow's Cp"

# ============================================================================================================================
# AnovaTable mostly from CoefTable
mutable struct AnovaTable
    cols::Vector
    colnms::Vector
    rownms::Vector
    pvalcol::Int
    teststatcol::Int
    function AnovaTable(cols::Vector, colnms::Vector, rownms::Vector,
                       pvalcol::Int=0, teststatcol::Int = 0)
        nc = length(cols)
        nrs = map(length, cols)
        nr = nrs[1]
        length(colnms) in [0, nc] || throw(ArgumentError("colnms should have length 0 or $nc"))
        length(rownms) in [0, nr] || throw(ArgumentError("rownms should have length 0 or $nr"))
        all(nrs .== nr) || throw(ArgumentError("Elements of cols should have equal lengths, but got $nrs"))
        pvalcol in 0:nc || throw(ArgumentError("pvalcol should be between 0 and $nc"))
        teststatcol in 0:nc || throw(ArgumentError("teststatcol should be between 0 and $nc"))
        new(cols, colnms, rownms, pvalcol, teststatcol)
    end

    function AnovaTable(mat::Matrix, colnms::Vector, rownms::Vector,
                       pvalcol::Int=0, teststatcol::Int=0)
        nc = size(mat, 2)
        cols = Any[mat[:, i] for i in 1:nc]
        AnovaTable(cols, colnms, rownms, pvalcol, teststatcol)
    end
end

"""
Show a p-value using 6 characters, either using the standard 0.XXXX
representation or as <Xe-YY.
"""
struct PValue
    v::Real
    function PValue(v::Real)
        0 <= v <= 1 || isnan(v) || error("p-values must be in [0; 1]")
        new(v)
    end
end

function show(io::IO, pv::PValue)
    v = pv.v
    if isnan(v)
        print(io,"")
    elseif v >= 1e-4
        @printf(io,"%.4f", v)
    else
        @printf(io,"<1e%2.2d", ceil(Integer, max(nextfloat(log10(v)), -99)))
    end
end

"""Show a test statistic using 2 decimal digits"""
struct TestStat <: Real
    v::Real
end

show(io::IO, x::TestStat) = isnan(x.v) ? print(io,"") : @printf(io, "%.4f", x.v)

"""Wrap a string so that show omits quotes"""
struct NoQuote
    s::String
end

show(io::IO, n::NoQuote) = print(io, n.s)

"""Filter NaN"""
struct OtherStat <: Real
    v::Real
end

function show(io::IO, x::OtherStat)
    v = x.v
    if isnan(v)
        print(io, "")
    elseif floor(v) == v
        print(io, Int(v))
    elseif v < 1e-4 && v > 0
        @printf(io,"<1e%2.2d", ceil(Integer, max(nextfloat(log10(v)), -99)))
    elseif abs(v) < 10
        @printf(io,"%.4f", v)
    else
        @printf(io,"%.2f", v)
    end
end

function show(io::IO, at::AnovaTable)
    cols = at.cols; rownms = at.rownms; colnms = at.colnms;
    nc = length(cols)
    nr = length(cols[1])
    if length(rownms) == 0
        rownms = [lpad("[$i]",floor(Integer, log10(nr)) + 3) for i in 1:nr]
    end
    mat = [j == 1 ? NoQuote(rownms[i]) :
           j - 1 == at.pvalcol ? PValue(cols[j - 1][i]) :
           j - 1 in at.teststatcol ? TestStat(cols[j - 1][i]) :
           cols[j-1][i] isa AbstractString ? NoQuote(cols[j - 1][i]) : OtherStat(cols[j - 1][i])
           for i in 1:nr, j in 1:nc + 1]
    # Code inspired by print_matrix in Base
    io = IOContext(io, :compact=>true, :limit=>false)
    A = Base.alignment(io, mat, 1:size(mat, 1), 1:size(mat, 2),
                       typemax(Int), typemax(Int), 3)
    nmswidths = pushfirst!(length.(colnms), 0)
    A = [nmswidths[i] > sum(A[i]) ? (A[i][1] + nmswidths[i] - sum(A[i]), A[i][2]) : A[i]
         for i in 1:length(A)]
    totwidth = sum(sum.(A)) + 2 * (length(A) - 1)
    println(io, repeat('─', totwidth))
    print(io, repeat(' ', sum(A[1])))
    for j in 1:length(colnms)
        print(io, "  ", lpad(colnms[j], sum(A[j + 1])))
    end
    println(io, '\n', repeat('─', totwidth))
    for i in 1:size(mat, 1)
        Base.print_matrix_row(io, mat, A, i, 1:size(mat, 2), "  ")
        i != size(mat, 1) && println(io)
    end
    print(io, '\n', repeat('─', totwidth))
    nothing
end

# ====================================================================================================================================
# AnovaTable api
function anovatable(aov::AnovaResult{<: RegressionModel}; kwargs...)
    at = _anovatable(aov; kwargs...)
    cfnames = coefnames(aov)
    # when first factor is null
    length(at.rownms) == length(cfnames) - 1 && popfirst!(cfnames)
    at.rownms = cfnames
    at
end

# check first and last modeltype
anovatable(aov::AnovaResult{<: Tuple}; kwargs...) =
    _anovatable(aov, typeof(first(aov.model)), typeof(last(aov.model)); kwargs...)

# default anovatable api for comparing multiple models
_anovatable(aov::AnovaResult{<: Tuple}, modeltype1, modeltype2; kwargs...) = _anovatable(aov; kwargs...)

function _anovatable(aov::AnovaResult{<: Tuple, FTest}; kwargs...)
    AnovaTable(hcat(vectorize.((
                    dof(aov),
                    [NaN, _diff(dof(aov))...],
                    dof_residual(aov),
                    deviance(aov),
                    [NaN, _diffn(deviance(aov))...],
                    teststat(aov),
                    pval(aov)
                    ))...),
              ["DOF", "ΔDOF", "Res.DOF", "Deviance", "ΔDeviance", "F value", "Pr(>|F|)"],
              ["$i" for i in eachindex(pval(aov))], 7, 6)
end

function _anovatable(aov::AnovaResult{<: Tuple, LRT}; kwargs...)
    AnovaTable(hcat(vectorize.((
                    dof(aov),
                    [NaN, _diff(dof(aov))...],
                    dof_residual(aov),
                    deviance(aov),
                    teststat(aov),
                    pval(aov)
                    ))...),
              ["DOF", "ΔDOF", "Res.DOF", "Deviance", "χ²", "Pr(>|χ²|)"],
              ["$i" for i in eachindex(pval(aov))], 6, 5)
end

# Show function that delegates to anovatable
function show(io::IO, aov::AnovaResult{<: RegressionModel, T}) where {T <: GoodnessOfFit}
    at = anovatable(aov)
    println(io, "Analysis of Variance")
    println(io)
    println(io, "Type $(aov.type) test / $(tname(T))")
    println(io)
    println(io, formula(aov.model))
    println(io)
    println(io, "Table:")
    show(io, at)
end

function show(io::IO, aov::AnovaResult{<: Tuple, T}) where {T <: GoodnessOfFit}
    at = anovatable(aov)
    println(io,"Analysis of Variance")
    println(io)
    println(io, "Type $(aov.type) test / $(tname(T))")
    println(io)
    for(id, m) in enumerate(aov.model)
        println(io,"Model $id: ", formula(m))
    end
    println(io)
    println(io,"Table:")
    show(io, at)
end


text:
function assemble_bilinear_form!(
    sysmatrix::SystemMatrix,
    cutmeshbfs::BilinearForms,
    cutmesh,
)

    dofspernode = dimension(cutmesh)
    ncells = number_of_cells(cutmesh)
    cellsign = cell_sign(cutmesh)

    uniformvals1 = vec(cutmeshbfs[+1])
    uniformvals2 = vec(cutmeshbfs[-1])

    for cellid = 1:ncells
        s = cellsign[cellid]
        if s == +1
            nodeids = nodal_connectivity(cutmesh, +1, cellid)
            assemble_cell_matrix!(sysmatrix, nodeids, dofspernode, uniformvals1)
        elseif s == -1
            nodeids = nodal_connectivity(cutmesh, -1, cellid)
            assemble_cell_matrix!(sysmatrix, nodeids, dofspernode, uniformvals2)
        elseif s == 0
            nodeids1 = nodal_connectivity(cutmesh, +1, cellid)
            vals1 = vec(cutmeshbfs[+1, cellid])
            assemble_cell_matrix!(sysmatrix, nodeids1, dofspernode, vals1)

            nodeids2 = nodal_connectivity(cutmesh, -1, cellid)
            vals2 = vec(cutmeshbfs[-1, cellid])
            assemble_cell_matrix!(sysmatrix, nodeids2, dofspernode, vals2)
        else
            error("Expected cellsign ∈ {+1,0,-1}, got cellsign = $s")
        end
    end
end

function assemble_interface_condition!(
    sysmatrix::SystemMatrix,
    interfacecondition::InterfaceCondition,
    cutmesh;
    eta = 1,
)

    @assert eta == 1 || eta == 0 || eta == -1
    dofspernode = dimension(cutmesh)
    cellsign = cell_sign(cutmesh)

    cellids = findall(cellsign .== 0)
    for cellid in cellids
        negativenodeids = nodal_connectivity(cutmesh, -1, cellid)
        positivenodeids = nodal_connectivity(cutmesh, +1, cellid)

        pptop = -0.5 * vec(traction_operator(interfacecondition, +1, +1, cellid))
        pntop = -0.5 * vec(traction_operator(interfacecondition, +1, -1, cellid))
        nptop = -0.5 * vec(traction_operator(interfacecondition, -1, +1, cellid))
        nntop = -0.5 * vec(traction_operator(interfacecondition, -1, -1, cellid))

        pptopT =
            -0.5 *
            eta *
            vec(transpose(traction_operator(interfacecondition, +1, +1, cellid)))
        pntopT =
            -0.5 *
            eta *
            vec(transpose(traction_operator(interfacecondition, +1, -1, cellid)))
        nptopT =
            -0.5 *
            eta *
            vec(transpose(traction_operator(interfacecondition, -1, +1, cellid)))
        nntopT =
            -0.5 *
            eta *
            vec(transpose(traction_operator(interfacecondition, -1, -1, cellid)))

        assemble_cell_matrix!(sysmatrix, positivenodeids, dofspernode, -pptop)
        assemble_cell_matrix!(sysmatrix, positivenodeids, dofspernode, -pptopT)

        assemble_couple_cell_matrix!(
            sysmatrix,
            positivenodeids,
            negativenodeids,
            dofspernode,
            -pntop,
        )
        assemble_couple_cell_matrix!(
            sysmatrix,
            positivenodeids,
            negativenodeids,
            dofspernode,
            nptopT,
        )

        assemble_cell_matrix!(sysmatrix, negativenodeids, dofspernode, nntop)
        assemble_cell_matrix!(sysmatrix, negativenodeids, dofspernode, nntopT)

        assemble_couple_cell_matrix!(
            sysmatrix,
            negativenodeids,
            positivenodeids,
            dofspernode,
            nptop,
        )
        assemble_couple_cell_matrix!(
            sysmatrix,
            negativenodeids,
            positivenodeids,
            dofspernode,
            -pntopT,
        )

        ppmass = vec(mass_operator(interfacecondition, +1, +1, cellid))
        pnmass = vec(mass_operator(interfacecondition, +1, -1, cellid))
        npmass = vec(mass_operator(interfacecondition, -1, +1, cellid))
        nnmass = vec(mass_operator(interfacecondition, -1, -1, cellid))

        assemble_cell_matrix!(sysmatrix, negativenodeids, dofspernode, nnmass)
        assemble_cell_matrix!(sysmatrix, positivenodeids, dofspernode, ppmass)
        assemble_couple_cell_matrix!(
            sysmatrix,
            negativenodeids,
            positivenodeids,
            dofspernode,
            -npmass,
        )
        assemble_couple_cell_matrix!(
            sysmatrix,
            positivenodeids,
            negativenodeids,
            dofspernode,
            -pnmass,
        )
    end
end

function assemble_body_force_linear_form!(systemrhs, rhsfunc, basis, cellquads, cutmesh)

    ncells = number_of_cells(cutmesh)
    cellsign = cell_sign(cutmesh)
    dofspernode = dimension(cutmesh)

    for cellid = 1:ncells
        s = cellsign[cellid]
        @assert s == -1 || s == 0 || s == 1
        if s == +1 || s == 0
            cellmap = cell_map(cutmesh, +1, cellid)
            pquad = cellquads[+1, cellid]
            rhs = linear_form(rhsfunc, basis, pquad, cellmap)
            nodeids = nodal_connectivity(cutmesh, +1, cellid)
            edofs = element_dofs(nodeids, dofspernode)
            assemble!(systemrhs, edofs, rhs)
        end
        if s == -1 || s == 0
            cellmap = cell_map(cutmesh, -1, cellid)
            nquad = cellquads[-1, cellid]
            rhs = linear_form(rhsfunc, basis, nquad, cellmap)
            nodeids = nodal_connectivity(cutmesh, -1, cellid)
            edofs = element_dofs(nodeids, dofspernode)
            assemble!(systemrhs, edofs, rhs)
        end
    end
end

function assemble_traction_force_linear_form!(
    systemrhs,
    tractionfunc,
    basis,
    facequads,
    cutmesh,
    onboundary,
)

    dim = dimension(cutmesh)
    facemidpoints = reference_face_midpoints()
    cellconnectivity = cell_connectivity(cutmesh)
    cellids = findall(is_boundary_cell(cutmesh))
    nfaces = length(facemidpoints)
    facedetjac = face_determinant_jacobian(cutmesh)

    refnormals = reference_face_normals()

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cell_sign(cutmesh, cellid)
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 0 || s == +1
                        rhs = linear_form(
                            tractionfunc,
                            basis,
                            facequads[+1, faceid, cellid],
                            cellmap,
                            facedetjac[faceid],
                        )
                        nodeids = nodal_connectivity(cutmesh, +1, cellid)
                        edofs = element_dofs(nodeids, dim)
                        assemble!(systemrhs, edofs, rhs)
                    end
                    if s == 0 || s == -1
                        rhs = linear_form(
                            tractionfunc,
                            basis,
                            facequads[-1, faceid, cellid],
                            cellmap,
                            facedetjac[faceid],
                        )
                        nodeids = nodal_connectivity(cutmesh, -1, cellid)
                        edofs = element_dofs(nodeids, dim)
                        assemble!(systemrhs, edofs, rhs)
                    end
                end
            end
        end
    end
end

function assemble_traction_force_component_linear_form!(
    systemrhs,
    tractionfunc,
    basis,
    facequads,
    cutmesh,
    onboundary,
    component,
)

    dim = dimension(cutmesh)
    facemidpoints = reference_face_midpoints()
    cellconnectivity = cell_connectivity(cutmesh)
    cellids = findall(is_boundary_cell(cutmesh))
    nfaces = length(facemidpoints)
    facedetjac = face_determinant_jacobian(cutmesh)

    refnormals = reference_face_normals()

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cell_sign(cutmesh, cellid)
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 0 || s == +1
                        rhs = component_linear_form(
                            tractionfunc,
                            basis,
                            facequads[+1, faceid, cellid],
                            component,
                            cellmap,
                            facedetjac[faceid],
                        )
                        nodeids = nodal_connectivity(cutmesh, +1, cellid)
                        edofs = element_dofs(nodeids, dim)
                        assemble!(systemrhs, edofs, rhs)
                    end
                    if s == 0 || s == -1
                        rhs = component_linear_form(
                            tractionfunc,
                            basis,
                            facequads[-1, faceid, cellid],
                            component,
                            cellmap,
                            facedetjac[faceid],
                        )
                        nodeids = nodal_connectivity(cutmesh, -1, cellid)
                        edofs = element_dofs(nodeids, dim)
                        assemble!(systemrhs, edofs, rhs)
                    end
                end
            end
        end
    end
end

function assemble_penalty_displacement_bc!(sysmatrix, sysrhs, dispcondition, cutmesh)

    ncells = number_of_cells(cutmesh)
    dim = dimension(cutmesh)

    for cellid = 1:ncells
        for faceid = 1:4
            for s in [-1, +1]
                if has_operator(dispcondition, s, faceid, cellid)
                    nodeids = nodal_connectivity(cutmesh, s, cellid)

                    mop = mass_operator(dispcondition, s, faceid, cellid)
                    top = traction_operator(dispcondition, s, faceid, cellid)
                    op = mop - top
                    assemble_cell_matrix!(sysmatrix, nodeids, dim, vec(op))

                    rhs = displacement_rhs(dispcondition, s, faceid, cellid)
                    assemble_cell_rhs!(sysrhs, nodeids, dim, rhs)
                end
            end
        end
    end
end

function assemble_bulk_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    cellquads,
    cutmesh,
)

    ncells = number_of_cells(cutmesh)
    cellsign = cell_sign(cutmesh)
    dofspernode = dimension(cutmesh)
    jac = jacobian(cutmesh)

    for cellid = 1:ncells
        s = cellsign[cellid]
        @assert s == -1 || s == 0 || s == +1
        if s == +1 || s == 0
            pquad = cellquads[+1, cellid]
            rhs = bulk_transformation_rhs(basis, pquad, transfstress, jac)
            nodeids = nodal_connectivity(cutmesh, +1, cellid)
            edofs = element_dofs(nodeids, dofspernode)
            assemble!(systemrhs, edofs, rhs)
        end
    end
end

function assemble_coherent_interface_transformation_rhs!(
    systemrhs,
    transfstress,
    basis,
    interfacequads,
    cutmesh,
)

    cellsign = cell_sign(cutmesh)
    dofspernode = dimension(cutmesh)

    cellids = findall(cellsign .== 0)
    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        normals = interface_normals(interfacequads, cellid)

        for s in [+1, -1]
            quad = interfacequads[s, cellid]
            rhs =
                s *
                0.5 *
                interface_transformation_rhs(basis, quad, normals, transfstress, cellmap)
            nodeids = nodal_connectivity(cutmesh, s, cellid)
            assemble_cell_rhs!(systemrhs, nodeids, dofspernode, rhs)
        end

    end
end

function assemble_incoherent_interface_transformation_rhs!(
    systemrhs,
    transfstress,
    basis,
    interfacequads,
    cutmesh,
)

    cellsign = cell_sign(cutmesh)
    dofspernode = dimension(cutmesh)

    cellids = findall(cellsign .== 0)
    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        normals = interface_normals(interfacequads, cellid)
        components = normals

        for s in [+1, -1]
            quad = interfacequads[s, cellid]
            rhs =
                s *
                0.5 *
                interface_transformation_component_rhs(
                    basis,
                    quad,
                    components,
                    normals,
                    transfstress,
                    cellmap,
                )
            nodeids = nodal_connectivity(cutmesh, s, cellid)
            assemble_cell_rhs!(systemrhs, nodeids, dofspernode, rhs)
        end

    end
end

function assemble_stress_mass_matrix!(sysmatrix, basis, cellquads, mesh)
    ncells = number_of_cells(mesh)
    detjac = determinant_jacobian(mesh)
    dim = dimension(basis)
    sdim = number_of_symmetric_degrees_of_freedom(dim)

    uniformquad = uniform_cell_quadrature(cellquads)
    uniformcellvals = vec(stress_cell_mass_matrix(basis, uniformquad, detjac))

    for cellid = 1:ncells
        s = cell_sign(mesh, cellid)
        @assert s == -1 || s == 0 || s == +1
        if s == -1 || s == +1
            nodeids = nodal_connectivity(mesh, s, cellid)
            assemble_cell_matrix!(sysmatrix, nodeids, sdim, uniformcellvals)
        else
            pquad = cellquads[+1, cellid]
            pvals = vec(stress_cell_mass_matrix(basis, pquad, detjac))
            pnodeids = nodal_connectivity(mesh, +1, cellid)
            assemble_cell_matrix!(sysmatrix, pnodeids, sdim, pvals)

            nquad = cellquads[-1, cellid]
            nvals = vec(stress_cell_mass_matrix(basis, nquad, detjac))
            nnodeids = nodal_connectivity(mesh, -1, cellid)
            assemble_cell_matrix!(sysmatrix, nnodeids, sdim, nvals)
        end
    end
end

function assemble_stress_linear_form!(
    sysrhs,
    basis,
    cellquads::CellQuadratures,
    stiffness,
    nodaldisplacement,
    mesh,
)

    dim = dimension(basis)
    sdim = number_of_symmetric_degrees_of_freedom(dim)
    ncells = number_of_cells(mesh)
    jac = jacobian(mesh)
    uniformquad = uniform_cell_quadrature(cellquads)

    uniformop =
        [stress_cell_rhs_operator(basis, uniformquad, stiffness[s], jac) for s in [+1, -1]]

    for cellid = 1:ncells
        s = cell_sign(mesh, cellid)
        @assert s == -1 || s == 0 || s == +1
        if s == -1 || s == +1
            nodeids = nodal_connectivity(mesh, s, cellid)
            dispdofs = element_dofs(nodeids, dim)

            idx = cell_sign_to_row(s)
            op = uniformop[idx]
            vals = op * nodaldisplacement[dispdofs]

            stressdofs = element_dofs(nodeids, sdim)
            assemble!(sysrhs, stressdofs, vals)
        else
            pnodeids = nodal_connectivity(mesh, +1, cellid)
            pdispdofs = element_dofs(pnodeids, dim)
            pdisp = nodaldisplacement[pdispdofs]
            pquad = cellquads[+1, cellid]
            pvals = stress_cell_rhs(basis, pquad, stiffness[+1], pdisp, jac)
            pstressdofs = element_dofs(pnodeids, sdim)
            assemble!(sysrhs, pstressdofs, pvals)

            nnodeids = nodal_connectivity(mesh, -1, cellid)
            ndispdofs = element_dofs(nnodeids, dim)
            ndisp = nodaldisplacement[ndispdofs]
            nquad = cellquads[-1, cellid]
            nvals = stress_cell_rhs(basis, nquad, stiffness[-1], ndisp, jac)
            nstressdofs = element_dofs(nnodeids, sdim)
            assemble!(sysrhs, nstressdofs, nvals)
        end
    end
end

function assemble_transformation_stress_linear_form!(
    sysrhs,
    transfstress,
    basis,
    cellquads::CellQuadratures,
    mesh,
)

    dim = dimension(basis)
    sdim = number_of_symmetric_degrees_of_freedom(dim)
    ncells = number_of_cells(mesh)
    detjac = determinant_jacobian(mesh)

    uniformquad = uniform_cell_quadrature(cellquads)
    uniformtransfrhs = constant_linear_form(transfstress, basis, uniformquad, detjac)

    for cellid = 1:ncells
        s = cell_sign(mesh, cellid)
        @assert s == -1 || s == 0 || s == +1
        if s == +1
            nodeids = nodal_connectivity(mesh, +1, cellid)
            stressdofs = element_dofs(nodeids, sdim)

            assemble!(sysrhs, stressdofs, -uniformtransfrhs)
        elseif s == 0
            pquad = cellquads[+1, cellid]
            transfrhs = constant_linear_form(transfstress, basis, pquad, detjac)
            nodeids = nodal_connectivity(mesh, +1, cellid)
            stressdofs = element_dofs(nodeids, sdim)

            assemble!(sysrhs, stressdofs, -transfrhs)
        end
    end
end

function assemble_penalty_displacement_transformation_rhs!(
    sysrhs,
    transfstress,
    basis,
    facequads,
    cutmesh,
    onboundary,
)

    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    facedetjac = face_determinant_jacobian(cutmesh)

    dim = dimension(cutmesh)
    refnormals = reference_face_normals()
    refquads = reference_face_quadrature(facequads)

    uniformrhs = [
        face_traction_transformation_rhs(basis, q, n, transfstress, s)
        for (q, n, s) in zip(refquads, refnormals, facedetjac)
    ]

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cell_sign(cutmesh, cellid)
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 1
                        rhs = uniformrhs[faceid]
                        nodeids = nodal_connectivity(cutmesh, +1, cellid)
                        assemble_cell_rhs!(sysrhs, nodeids, dim, -rhs)
                    elseif s == 0
                        pquad = facequads[+1, faceid, cellid]
                        rhs = face_traction_transformation_rhs(
                            basis,
                            pquad,
                            refnormals[faceid],
                            transfstress,
                            facedetjac[faceid],
                        )
                        nodeids = nodal_connectivity(cutmesh, +1, cellid)
                        assemble_cell_rhs!(sysrhs, nodeids, dim, -rhs)
                    end
                end
            end
        end
    end
end

function assemble_penalty_displacement_component_transformation_rhs!(
    sysrhs,
    transfstress,
    basis,
    facequads,
    cutmesh,
    onboundary,
    component,
)

    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    facedetjac = face_determinant_jacobian(cutmesh)

    dim = dimension(cutmesh)
    refnormals = reference_face_normals()
    refquads = reference_face_quadrature(facequads)

    uniformrhs = [
        face_traction_component_transformation_rhs(basis, q, component, n, transfstress, s) for (q, n, s) in zip(refquads, refnormals, facedetjac)
    ]

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cell_sign(cutmesh, cellid)
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 1
                        rhs = uniformrhs[faceid]
                        nodeids = nodal_connectivity(cutmesh, +1, cellid)
                        assemble_cell_rhs!(sysrhs, nodeids, dim, -rhs)
                    elseif s == 0
                        pquad = facequads[+1, faceid, cellid]
                        rhs = face_traction_component_transformation_rhs(
                            basis,
                            pquad,
                            component,
                            refnormals[faceid],
                            transfstress,
                            facedetjac[faceid],
                        )
                        nodeids = nodal_connectivity(cutmesh, +1, cellid)
                        assemble_cell_rhs!(sysrhs, nodeids, dim, -rhs)
                    end
                end
            end
        end
    end
end


text:
## Spectral clustering
using DelimitedFiles
using SparseArrays
using LinearAlgebra
using Random
ei,ej = copy.(eachcol(Int.(readdlm("undir-graph-14.edges"))))
xy = readdlm("undir-graph-14.xy")'
n = max(maximum(ei),maximum(ej))
A = sparse(ei,ej,1,n,n)
A = A+A'
m = length(ei) # number of edges
B = sparse(1:m,ei,1,m,n) - sparse(1:m,ej,1,m,n)
L = B'*B
d = vec(sum(A,dims=2))
D = Diagonal(d)
Di = inv(D) # Julia does this efficiently...
Pt = Di*A
P = A*Di # A is symmetric ...
## Make a vector orthogonal to degrees
Random.seed!(1)
x = randn(n)
#x .-= ((d'*x)/(d'*d)).*d
x .-= (d'*x)/(sum(d))
x ./= x'*D*x

##

##
Random.seed!(1)
function simple_spectral_eigenvector(A,nsteps=500)
  @assert issymmetric(A) # expensive, but useful...
  n = size(A,1)
  d = vec(sum(A,dims=1))
  nd2 = d'*d
  x = randn(n)
  # project x orthogonal to d
  x .-= ((d'*x)/nd2).*d
  x ./= x'*(d.*x) # normalize
  for i=1:nsteps
    x = (A*x)./d + x       # matrix-vector
    x .-= ((d'*x)/nd2).*d
    x ./= sqrt(x'*(d.*x)) # normalize
    # make sure the first component is positive
    x .*= sign(x[1])
  end

  return x
end
f = simple_spectral_eigenvector(A)
##
using MatrixNetworks
f2 = fiedler_vector(A)[1]
# make sure the first component is positive
f2 .*= sign(f2[1])
##
norm(f-f2)

##
using Plots
using GraphRecipes
pyplot(size=(300,300))
p = graphplot(ei, ej, x =xy[1,:], y=xy[2,:],
  markercolor=:white, markerstrokecolor=:black,
  markersize=13, linecolor=3, linealpha=1.0, linewidth=1.0,
  marker_z = f,
  axis_buffer=0.5, background=:transparent,dpi=300)
xlims!(0,1)
ylims!(-0.05, 1.05)
for i=1:size(A,1)
  if f[i] >= 0
    annotate!(xy[1,i],xy[2,i],text("$i",9,:black))
  else
    annotate!(xy[1,i],xy[2,i],text("$i",9,:white))
  end
end
plot!()
savefig("spectral-clustering-undir.pdf")


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite  -  HS109
# This file is JuMP implementation of the model described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model description, at
# http://www.ai7.uni-bayreuth.de/downloads.htm
#
# Nonconvex objective with cubics, and constraints contain trignometry and
# products of variables.
#############################################################################

using JuMP
using Base.Test

let

a  = 50.176
b1 = 0.25
b = sin(b1)
c = cos(b1)

L = [0.0, 0.0, -0.55, -0.55, 196, 196, 196, -400, -400]
U = [Inf, Inf,  0.55,  0.55, 252, 252, 252,  800,  800]

m = Model()
@variable(m, L[i] <= x[i=1:9] <= U[i], start = 0.0)

@NLobjective(m, Min, 3 * x[1] + 1e-6 * x[1]^3 + 2 * x[2] + .522074e-6 * x[2]^3)

@NLconstraint(m, x[4] - x[3] + 0.55 >= 0)
@NLconstraint(m, x[3] - x[4] + 0.55 >= 0)
@NLconstraint(m, 2250000 - x[1]^2 - x[8]^2 >= 0)
@NLconstraint(m, 2250000 - x[2]^2 - x[9]^2 >= 0)
@NLconstraint(m,
    x[5] * x[6] * sin(-x[3] - .25) + x[5] * x[7] * sin(-x[4] - .25) +
    2 * b * x[5]^2 - a * x[1] + 400 * a == 0)
@NLconstraint(m,
    x[5] * x[6] * sin(x[3] - .25) + x[6] * x[7] * sin(x[3] - x[4] - .25) +
    2 * b * x[6]^2 - a * x[2] + 400 * a == 0)
@NLconstraint(m,
    x[5] * x[7] * sin(x[4] - .25) + x[6] * x[7] * sin(x[4] - x[3] - .25) +
    2 * b * x[7]^2 + 881.779 * a == 0)
@NLconstraint(m,
    a * x[8] + x[5] * x[6] * cos(-x[3] - .25) +
    x[5] * x[7] * cos(-x[4] - .25) - 200 * a - 2 * c * x[5]^2 +
    .7533e-3 * a * x[5]^2 == 0)
@NLconstraint(m,
    a * x[9] + x[5] * x[6] * cos(x[3] - .25) +
    x[6] * x[7] * cos(x[3] - x[4] - .25) - 2 * c * x[6]^2 +
    .7533e-3 * a * x[6]^2 - 200 * a == 0)
@NLconstraint(m,
    x[5] * x[7] * cos(x[4] - .25) + x[6] * x[7] * cos(x[4] - x[3] - .25) -
    2 * c * x[7]^2 + 22.938 * a + .7533e-3 * a * x[7] ^2 == 0)

solve(m)

@test_approx_eq_eps getobjectivevalue(m) 5326.851310161077 1e-5

end


text:
__precompile__()

module GeometricTransforms

### Dependencies
using Compat
using Reexport

@reexport using StaticArrays

### Imports
import Base: in

### Implementation
path = dirname(realpath(@__FILE__)) # works even for symlinks
include(joinpath(path, "types.jl"))
include(joinpath(path, "utils.jl"))
include(joinpath(path, "in.jl"))
include(joinpath(path, "transforms.jl"))

### Exports
export Shape, Cube, Cylinder, Ellipsoid, EllipticCylinder, HollowCylinder,
       Parallelepiped, Point, RectangularPyramid, Ring, Sphere, SphericalCap,
       SquarePyramid, TSP, TriangularToroid, TruncatedSquarePyramid, Torus,
       Vec, domain, halflengths, ftransform, ptransform

end # module


text:
module NeRCA

using LinearAlgebra
using Statistics
using Printf
using Sockets

using DocStringExtensions
using Parameters
using StaticArrays
using UnROOT
using PGFPlotsX
using Colors
using RecipesBase
using JuMP
using Ipopt
using IntervalArithmetic
using LandauDistribution
using Distributions: pdf, Normal

using CSV
using DataFrames
using HTTP

import Base: read, +, -, *, getindex, length, eltype

export
    OnlineFile, OfflineFile,
    Position, Direction,
    MCEventInfo, MCTrack, RecoTrack, NoRecoTrack,
    AbstractDAQHit, AbstractMCHit,
    Hit, SnapshotHit, CalibratedHit, MCHit, TimesliceHit,
    DAQEvent,
    ismxshower, is3dmuon, is3dshower, isnb,
    Calibration, calibrate, OMKey, omkey2domid, floordist,
    duhits, domhits, nfoldhits, triggered,
    CHClient, CHTag, CHPrefix, CHMessage, subscribe,
    ztplot,
    savefigs,
    most_frequent, categorize,
    SingleDUParams, SingleDURecoParams,
    @ip_str,
    initdb, streamds, detx  # db.jl

include("definitions/daqdatatypes.jl")
include("definitions/fitparameters.jl")
include("definitions/reconstruction.jl")
include("definitions/trigger.jl")
include("definitions/w2list_genhen.jl")
include("definitions/w2list_gseagen.jl")

include("constants.jl")
include("types.jl")
include("tools.jl")
include("io/detx.jl")
include("io/root.jl")
include("calibration.jl")
include("math.jl")
include("hits.jl")
include("mc.jl")
include("fit.jl")
include("controlhost.jl")
include("plot.jl")
include("db.jl")

end


text:
# 1. Create example point data
x = 8 .* rand(100) .- 4
y = 8 .* rand(100) .- 4
z = sin.(x) + cos.(y)
# Contour plot on the left without color (lines labelled by default)
subplot(1, 2, 1)
contour(x, y, z, colorbar = false)
# 2. Create example grid data with a callable
x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
f(x, y) = sin(x) + cos(y)
# Contour plot on the right with color
# and explicit labels every three lines
subplot(1, 2, 2)
contour(x, y, f, majorlevels = 3)


text:
# cubic L2 norm (times a constant

export CubeNormL2

"""
**Cubic Euclidean norm (weighted)**

    CubeNormL2(λ=1.0)

With a nonnegative scalar `λ`, returns the function
```math
f(x) = λ\\|x\\|^3.
```
"""
struct CubeNormL2{R <: Real} <: ProximableFunction
    lambda::R
    function CubeNormL2{R}(lambda::R) where R
        if lambda < 0
            error("coefficient λ must be nonnegative")
        else
            new(lambda)
        end
    end
end

is_convex(f::CubeNormL2) = true
is_smooth(f::CubeNormL2) = true

CubeNormL2(lambda::R=1.0) where {R <: Real} = CubeNormL2{R}(lambda)

function (f::CubeNormL2{R})(x::AbstractArray{T}) where {R, T <: RealOrComplex{R}}
    return f.lambda * norm(x)^3
end

function gradient!(y::AbstractArray{T}, f::CubeNormL2{R}, x::AbstractArray{T}) where {R, T <: RealOrComplex{R}}
    norm_x = norm(x)
    y .= (3 * f.lambda * norm_x) .* x
    return f.lambda * norm_x^3
end

function prox!(y::AbstractArray{T}, f::CubeNormL2{R}, x::AbstractArray{T}, gamma::R=R(1)) where {R, T <: RealOrComplex{R}}
    norm_x = norm(x)
    scale = 2 / (1 + sqrt(1 + 12 * gamma * f.lambda * norm_x))
    y .= scale .* x
    return f.lambda * (scale * norm_x)^3
end

function prox_naive(f::CubeNormL2{R}, x::AbstractArray{T}, gamma=R(1)) where {R, T <: RealOrComplex{R}}
    y = 2 / (1 + sqrt(1 + 12 * gamma * f.lambda * norm(x))) * x
    return y, f.lambda * norm(y)^3
end


text:
module Events

import ..Dagger
import ..Dagger: Event, Chunk
import ..Dagger: init_similar

"""
    CoreMetrics

Tracks the timestamp, category, and kind of the `Event` object generated by log
events.
"""
struct CoreMetrics end

(::CoreMetrics)(ev::Event{T}) where T =
    (;timestamp=ev.timestamp,
      category=ev.category,
      kind=T)

"""
    IDMetrics

Tracks the ID of `Event` objects generated by log events.
"""
struct IDMetrics end

(::IDMetrics)(ev::Event{T}) where T = ev.id

"""
    TimelineMetrics

Tracks the timeline of `Event` objects generated by log events.
"""
struct TimelineMetrics end

(::TimelineMetrics)(ev::Event{T}) where T = ev.timeline

"""
    FullMetrics

Tracks the full `Event` object generated by log events.
"""
struct FullMetrics end

(::FullMetrics)(ev) = ev

"""
    BytesAllocd

Tracks memory allocated for `Chunk`s.
"""
mutable struct BytesAllocd
    allocd::Int
end
BytesAllocd() = BytesAllocd(0)
init_similar(::BytesAllocd) = BytesAllocd()

function (ba::BytesAllocd)(ev::Event{:start})
    if ev.category in (:move, :evict) && ev.timeline.data isa Chunk
        sz = Int(ev.timeline.data.handle.size)
        if ev.category == :move && !haskey(Dagger.Sch.CHUNK_CACHE, ev.timeline.data)
            ba.allocd += sz
        elseif ev.category == :evict && haskey(Dagger.Sch.CHUNK_CACHE, ev.timeline.data)
            ba.allocd -= sz
        end
    end
    ba.allocd
end
(ba::BytesAllocd)(ev::Event{:finish}) = ba.allocd

"""
    CPULoadAverages

Monitors the CPU load averages.
"""
struct CPULoadAverages end

(::CPULoadAverages)(ev::Event) = Sys.loadavg()[1]

"""
    MemoryFree

Monitors the percentage of free system memory.
"""
struct MemoryFree end

(::MemoryFree)(ev::Event) = Sys.free_memory() / Sys.total_memory()

"""
    EventSaturation

Tracks the compute saturation (running tasks) per-processor.
"""
mutable struct EventSaturation
    saturation::Dict{Symbol,Int}
end
EventSaturation() = EventSaturation(Dict{Symbol,Int}())
init_similar(::EventSaturation) = EventSaturation()

function (es::EventSaturation)(ev::Event{:start})
    old = get(es.saturation, ev.category, 0)
    es.saturation[ev.category] = old + 1
    NamedTuple(filter(x->x[2]>0, es.saturation))
end
function (es::EventSaturation)(ev::Event{:finish})
    old = get(es.saturation, ev.category, 0)
    es.saturation[ev.category] = old - 1
    NamedTuple(filter(x->x[2]>0, es.saturation))
end

"""
    WorkerSaturation

Tracks the compute saturation (running tasks).
"""
mutable struct WorkerSaturation
    saturation::Int
end
WorkerSaturation() = WorkerSaturation(0)
init_similar(::WorkerSaturation) = WorkerSaturation()

function (ws::WorkerSaturation)(ev::Event{:start})
    if ev.category == :compute
        ws.saturation += 1
    end
    ws.saturation
end
function (ws::WorkerSaturation)(ev::Event{:finish})
    if ev.category == :compute
        ws.saturation -= 1
    end
    ws.saturation
end

"""
    ProcessorSaturation

Tracks the compute saturation (running tasks) per-processor.
"""
mutable struct ProcessorSaturation
    saturation::Dict{Dagger.Processor,Int}
end
ProcessorSaturation() = ProcessorSaturation(Dict{Dagger.Processor,Int}())
init_similar(::ProcessorSaturation) = ProcessorSaturation()

function (ps::ProcessorSaturation)(ev::Event{:start})
    if ev.category == :compute
        proc = ev.timeline.to_proc
        old = get(ps.saturation, proc, 0)
        ps.saturation[proc] = old + 1
    end
    filter(x->x[2]>0, ps.saturation)
end
function (ps::ProcessorSaturation)(ev::Event{:finish})
    if ev.category == :compute
        proc = ev.timeline.to_proc
        old = get(ps.saturation, proc, 0)
        ps.saturation[proc] = old - 1
    end
    filter(x->x[2]>0, ps.saturation)
end

"""
    LogWindow

Aggregator that prunes events to within a given time window.
"""
mutable struct LogWindow
    window_length::UInt64
    core_name::Symbol
    creation_handlers::Vector{Any}
    deletion_handlers::Vector{Any}
end
LogWindow(window_length, core_name) = LogWindow(window_length, core_name, [], [])

function (lw::LogWindow)(logs::Dict)
    core_log = logs[lw.core_name]

    # Inform creation hooks
    if length(core_log) > 0
        log = Dict{Symbol,Any}()
        for key in keys(logs)
            log[key] = last(logs[key])
        end
        for obj in lw.creation_handlers
            creation_hook(obj, log)
        end
    end

    # Find entries outside the window
    window_start = time_ns() - lw.window_length
    idx = something(findfirst(ev->ev.timestamp >= window_start, core_log),
                    length(core_log))

    # Return if all events are in the window
    if length(core_log) == 0 || (core_log[1].timestamp >= window_start)
        return
    end

    # Inform deletion hooks
    for obj in lw.deletion_handlers
        deletion_hook(obj, idx)
    end

    # Remove entries outside the window
    for name in keys(logs)
        deleteat!(logs[name], 1:idx)
    end
end

creation_hook() = nothing
deletion_hook() = nothing

end # module Events


text:
""" Julia program to Check if a number is an Armstrong Number or not.
An armstrong number is a number such that the sum of its digits each
raised to the power of the number of digits.
"""

function check_armstrong(num)
    temp = num
    sum = 0
    cnt = floor(log10(num) + 1)
    while (temp > 0)
        rem = temp % 10
        sum = sum + (rem ^ cnt)
        temp = temp ÷ 10;
    end
    if(num == sum)
        return true
    else
        return false
    end
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
if(num < 0)
    println("The given number $num is not a Armstrong Number.")
    exit()
end
res = check_armstrong(num)
if res
    println("The given number $num is a Armstrong Number.")
else
    println("The given number $num is not a Armstrong Number.")
end

"""
Time Complexity: O(log(num)), where 'num' is the given number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1
Enter the number: 371
The given number 371 is a Armstrong Number.

SAMPLE 2
Enter the number: 1234
The given number 1234 is not a Armstrong Number.

"""


text:
using Test
using TimerOutputs
using CUDA

using Athesis

@testset "Athesis tests" begin
    include("unit.jl")
    include("validation.jl")
    include("benchmark.jl")
end


text:
using Quaternions
using Documenter

DocMeta.setdocmeta!(Quaternions, :DocTestSetup, :(using Quaternions); recursive=true)

makedocs(;
    modules=[Quaternions],
    repo="https://github.com/JuliaGeometry/Quaternions.jl/blob/{commit}{path}#{line}",
    sitename="Quaternions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaGeometry.github.io/Quaternions.jl",
        assets = ["assets/custom.css", "assets/favicon.ico"],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaGeometry/Quaternions.jl",
)


text:
using DifferencesBase
using SparseArrays
using Test

@testset "DifferencesBase" begin
    include("vector.jl")
    include("matrix.jl")
    include("dict.jl")
    include("namedtuple.jl")
    include("set.jl")
    include("show.jl")
end


text:
using FdeSolver
using Plots

## inputs
I0 = 0.001             # intial value of infected
tSpan = [0, 100]       # [intial time, final time]
y0 = [1 - I0, I0, 0]   # initial values [S0,I0,R0]
α = [1, 1, 1]          # order of derivatives
h = 0.1                # step size of computation (default=0.01)

## ODE model
par = [0.4, 0.04] # parameters [β, recovery rate]

function F(t, y, par)

    # parameters
    β = par[1]    # infection rate
    γ = par[2]    # recovery rate

    S = y[1]   # Susceptible
    I = y[2]   # Infectious
    R = y[3]   # Recovered

    # System equation
    dSdt = - β .* S .* I
    dIdt = β .* S .* I .- γ .* I
    dRdt = γ .* I

    return [dSdt, dIdt, dRdt]

end

## Jacobian of ODE system
function JacobF(t, y, par)

    # parameters
    β = par[1]     # infection rate
    γ = par[2]     # recovery rate

    S = y[1]    # Susceptible
    I = y[2]    # Infectious
    R = y[3]    # Recovered

    # System equation
    J11 = - β * I
    J12 = - β * S
    J13 =  0
    J21 =  β * I
    J22 =  β * S - γ
    J23 =  0
    J31 =  0
    J32 =  γ
    J33 =  0

    J = [J11 J12 J13
         J21 J22 J23
         J31 J32 J33]

    return J

end
## Solution
t, Yapp = FDEsolver(F, tSpan, y0, α, par, h = h)
t1, Yapp1 = FDEsolver(F, tSpan, y0, α, par, JF = JacobF, h = h)

## plotting
plot(t, Yapp, linewidth = 5, title = "Numerical solution of SIR model",
     xaxis = "Time (t)", yaxis = "SIR populations", label=["Susceptible" "Infectious" "Recovered"])
plot!(t1, Yapp1, linewidth = 5, ls = :dash, label=["S_Jacob" "I_Jacob" "R_Jacob"])


text:
using ABCDeGraphGenerator
using Random

@info "Usage: julia deg_sampler.jl filename τ₁ d_min d_max n max_iter [seed]"
@info "Example: julia deg_sampler.jl degrees.dat 3 5 50 10000 1000 42"

filename = ARGS[1]
τ₁ = parse(Float64, ARGS[2])
d_min = parse(Int, ARGS[3])
d_max = parse(Int, ARGS[4])
n = parse(Int, ARGS[5])
max_iter = parse(Int, ARGS[6])
length(ARGS) == 7 && Random.seed!(parse(Int, ARGS[7]))

@info "Expected value of degree: $(ABCDeGraphGenerator.get_ev(τ₁, d_min, d_max))"

degs = ABCDeGraphGenerator.sample_degrees(τ₁, d_min, d_max, n, max_iter)

open(io -> foreach(d -> println(io, d), degs), filename, "w")


text:
# render a block of a RegressionTable
function render(io::IO, block::Array{String, 2}, colWidths::Vector{Int64}, align::String, settings::RenderSettings = asciiSettings; isHeader::Bool = false)

    c = size(block,2)

    if length(colWidths) != c
        error("colWidths has invalid length.")
    end
    if length(align) != c
        error("align string has invalid length.")
    end

    # print the whole thing
    for row = 1:size(block,1)
        s = settings.linestart
        for col = 1:c
            # if the string is too long, truncate it
            # (this sometimes happens with column headers)
            # if length(block[row,col]) > colWidths[col]
            #     printstring = block[row,col][1:colWidths[col]]
            # else
            # end
            printstring = block[row,col]

            if align[col] == 'l'
                s = s * rpad(printstring,colWidths[col])
            elseif align[col] == 'r'
                s = s * lpad(printstring,colWidths[col])
            elseif align[col] == 'c'
                l = iseven(colWidths[col]-length(printstring)) ? Int64((colWidths[col]-length(printstring))/2) : Int64((colWidths[col]-length(printstring)+1)/2)
                r = iseven(colWidths[col]-length(printstring)) ? Int64((colWidths[col]-length(printstring))/2) : Int64((colWidths[col]-length(printstring)-1)/2)
                # if the printstring is too long, so be it
                l = max(l,0)
                r = max(r,0)
                s = s * (" " ^ l) * printstring * (" " ^ r)
            end
            if col < c
                s = s * (isHeader ? settings.headercolsep : settings.colsep)
            end
        end
        s = s * settings.linebreak
        println(io, s)
    end

end

# render a whole table
function render(io::IO, tab::RegressionTable, align::String, settings::RenderSettings)

    c = columns(tab)

    if length(align) != c
        error("align string has invalid length.")
    end

    # construct column width, first from the maximum of the bodies' column widths.
    colWidths = zeros(Int64, c)
    for colIndex = 1:c
        if (align[colIndex] == 'l') || (align[colIndex] == 'r') || (align[colIndex] == 'c')
            colWidths[colIndex] = maximum([length(b[r,colIndex]) for b in tab.bodies for r=1:size(b,1) ])
        else
            error("Invalid character in align string. Only 'l', 'r', 'c' are allowed.")
        end
    end

    # construct, but not print, the header
    # header
    headerLabels = Vector{String}(undef,0)
    headerWidths = Vector{Int64}(undef,0)
    headerCellStartEnd = Vector{Vector{Int64}}(undef,0)
    # first column is empty (top left)
    push!(headerLabels, "")
    push!(headerWidths, colWidths[1])
    push!(headerCellStartEnd, [1,1])
    # first regression result
    push!(headerLabels, tab.header[1,2])
    push!(headerWidths, colWidths[2])
    push!(headerCellStartEnd, [2,2])
    if columns(tab)>2
        for rIndex = 3:size(tab.header,2)
            if tab.header[1,rIndex] == tab.header[1,rIndex-1]
                headerWidths[end] += length(settings.headercolsep) + colWidths[rIndex]
                headerCellStartEnd[end][2] += 1
            else
                push!(headerLabels, tab.header[1,rIndex])
                push!(headerWidths, colWidths[rIndex])
                push!(headerCellStartEnd, [rIndex, rIndex])
            end
        end
    end
    # second line
    headerArray = Array{String}(undef,1,length(headerLabels))
    headerArray[1,1] = ""
    for i = 2:size(headerArray,2)
        headerArray[1,i] = settings.encapsulateRegressand(headerLabels[i],headerCellStartEnd[i][1],headerCellStartEnd[i][2] )
    end
    # now it could be that a columns of the header is wider than the colWidth of the bodies.
    # in that case, make the last column of the bodies wider.
    for i = 2:size(headerArray,2)
        totalWidth = sum([colWidths[cind] for cind in headerCellStartEnd[i][1] : headerCellStartEnd[i][2]]) + length(settings.colsep)*(headerCellStartEnd[i][2] - headerCellStartEnd[i][1])
        if length(headerArray[1,i])>totalWidth
            # extend width of cells
            colWidths[headerCellStartEnd[i][2]] +=  (length(headerArray[1,i])-totalWidth)
            headerWidths[i] += (length(headerArray[1,i])-totalWidth)
        end
    end
    # second line
    # distinguish two cases:
    #   if headerrule gives a string of length one, put into table, and repeat the string
    hr = settings.headerrule(headerCellStartEnd)
    if length(hr) == 1
        secondRow = Array{String}(undef,1,length(headerLabels))
        secondRow[1,1] = ""
        for i = 2:size(secondRow,2)
            secondRow[1,i] = (length(settings.midrule) == 1 ? settings.midrule ^ headerWidths[i] : settings.midrule)
        end
        headerArray = [headerArray; secondRow]
        print_headerrule_separately = false
    else
        print_headerrule_separately = true
    end

    # START RENDERING

    # header
    println(io, settings.header(c,align))

    # headerString
    # if tab.headerString == ""
    #     # don't print anything
    # elseif length(tab.headerString) == 1
    #     println(io, tab.headerString ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep)))
    # else
    #     println(io, tab.headerString)
    # end

    # print toprule
    if length(settings.toprule)==1
        # one character, extend over the whole line
        println(io, settings.toprule ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep))  )
    else
        println(io, settings.toprule)
    end

    # header
    render(io, headerArray, headerWidths, ("c" ^ size(headerArray,2)), settings, isHeader=true)
    if print_headerrule_separately
        println(io, hr)
    end

    # bodies
    for b = 1:size(tab.bodies,1)

        render(io, tab.bodies[b], colWidths, align, settings)

        # if we're not at the last block, print the midrule
        if b < size(tab.bodies,1)
            if length(settings.midrule)==1
                # one character, extend over the whole line
                println(io, settings.midrule ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep))  )
            else
                println(io, settings.midrule)
            end
        end
    end

    # print bottomrule
    if length(settings.bottomrule)==1
        # one character, extend over the whole line
        println(io, settings.bottomrule ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep))  )
    else
        println(io, settings.bottomrule)
    end

    # footerString
    # if tab.footerString == ""
    #     # don't print anything
    # elseif length(tab.footerString) == 1
    #     println(io, tab.footerString ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep)))
    # else
    #     println(io, tab.footerString)
    # end

    # footer
    println(io, settings.footer(c,align))

end


text: #This script will handle all pressure conversion.


text:
@discord_object struct ActivitySecrets
    join::String
    spectate::String
    watch::String
end


text:
__precompile__()

module Clustering
    using Compat
    using Distances
    using NearestNeighbors
    using StatsBase

    import Base: show
    import StatsBase: IntegerVector, RealVector, RealMatrix, counts
    import Compat: view, issymmetric

    export
    # reexport from StatsBase
    sample, sample!,

    # common
    ClusteringResult,
    nclusters, counts, assignments,

    # seeding
    SeedingAlgorithm,
    RandSeedAlg, KmppAlg, KmCentralityAlg,
    copyseeds, copyseeds!,
    initseeds, initseeds!,
    initseeds_by_costs, initseeds_by_costs!,
    kmpp, kmpp_by_costs,

    # kmeans
    kmeans, kmeans!, KmeansResult, kmeans_opts,

    # kmedoids
    kmedoids, kmedoids!, KmedoidsResult,

    # affprop
    AffinityPropResult, affinityprop,

    # dbscan
    DbscanResult, DbscanCluster, dbscan,

    # silhouette
    silhouettes,

    # varinfo
    varinfo,

    # randindex
    randindex,

    # hclust
    Hclust, hclust, cutree,

    # MCL
    mcl, MCLResult

    ## source files

    include("utils.jl")
    include("seeding.jl")

    include("kmeans.jl")
    include("kmedoids.jl")
    include("affprop.jl")
    include("dbscan.jl")
    include("mcl.jl")

    include("silhouette.jl")
    include("randindex.jl")
    include("varinfo.jl")

    include("hclust.jl")

    include("deprecate.jl")
end


text:
# Wrap the functions defined in C++
module CppHello
  include(joinpath(@__DIR__, "testcommon.jl"))
  @wrapmodule CxxWrap.CxxWrapCore.libhello()

  function __init__()
    @initcxx
  end
end

using Test

# Output:
@show CppHello.greet()

@testset "$(basename(@__FILE__)[1:end-3])" begin

# Test the result
@test CppHello.greet() == "hello, world"

end


text:
using GitWorkers

## ---------------------------------------------------------------
# cmd for run the server
# julia 'SERVER_SCRIPT_PATH'

## ---------------------------------------------------------------
# run to reset all
# gw_create_devland(;
#     sys_root = "SYSTEM_ROOT",
#     clear_repos = true,
#     clear_scripts = false,
#     verbose = false,
#     deb = false,
# )

## ---------------------------------------------------------------
gw_setup(;
    sys_root = "CLIENT_ROOT",
    url = "REMOTE_URL",
)

## ---------------------------------------------------------------
# Make your stuff here

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
println("==============flag! bb_type_correct.jl======================================")

"""
    is_type_correct(x, var_type, atol)

Check whether a variable x has the correct type
"""
function is_type_correct(x, var_type, atol)
    if var_type != :Cont
        if !isapprox(round(x)-x, 0; atol=atol)
            return false
        end
    end
    return true
end

"""
    are_type_correct(sol, types, disc2var_idx, atol)

Check whether all variables have the correct type
"""
function are_type_correct(sol, types, disc2var_idx, atol)
    for i in disc2var_idx
        if !isapprox(round(sol[i])-sol[i], 0; atol=atol)
            return false
        end
    end
    return true
end

"""
    all_reasonable_type_correct(sol, disc2var_idx, reasonable_idx, atol)

Check whether all reasonable variable have are discrete alreadyy
"""
function all_reasonable_type_correct(sol, disc2var_idx, reasonable_idx, atol)
    for i in reasonable_idx
        idx = disc2var_idx[i]
        if !isapprox(round(sol[idx])-sol[idx], 0; atol=atol)
            return false
        end
    end
    return true
end


text:
__precompile__()
module PlotSeis
using Printf, SeisIO, PyPlot
using Statistics: mean
using DSP: spectrogram, time, freq, power
import SeisIO: t_expand, μs, t_win
export plotseis,
       uptimes,
       logspec

# Time series plots
include("TimeSeries/ts_internals.jl")
include("TimeSeries/plotseis.jl")

# Spectral plots
include("Spectral/logspec.jl")

# Other plots
include("Other/uptimes.jl")

end


text:
function solve_rosenbrock()
    function get_rosenbrock()
        function eval_obj(x::AbstractArray{T}) where T
            f = (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2

            f
        end

        function eval_cons(x::AbstractArray{T}) where T
            g = vcat(-5. .- x, x .- 5.)

            g
        end

        return generate_autodiff_solver_fn(eval_obj,eval_cons,[],1:4)
    end

    solver_fn = get_rosenbrock()

    x0 = [0.0; 0.0]
    options = Dict{String, Any}()
    options["Derivative option"] = 1
    options["Verify level"] = 1
    xopt, info = Bilevel.snopt(solver_fn, 0, 4, x0, options)

    function get_barnesgrad()
        a1 = 75.196
        a3 = 0.12694
        a5 = 1.0345e-5
        a7 = 0.030234
        a9 = 3.5256e-5
        a11 = 0.25645
        a13 = 1.3514e-5
        a15 = -5.2375e-6
        a17 = 7.0e-10
        a19 = -1.6638e-6
        a21 = 0.0005
        a2 = -3.8112
        a4 = -2.0567e-3
        a6 = -6.8306
        a8 = -1.28134e-3
        a10 = -2.266e-7
        a12 = -3.4604e-3
        a14 = -28.106
        a16 = -6.3e-8
        a18 = 3.4054e-4
        a20 = -2.8673

        function eval_obj(x::AbstractArray{T}) where T
            x1 = x[1]
            x2 = x[2]
            y1 = x1*x2
            y2 = y1*x1
            y3 = x2^2
            y4 = x1^2

            f = a1 + a2*x1 + a3*y4 + a4*y4*x1 + a5*y4^2 +
                a6*x2 + a7*y1 + a8*x1*y1 + a9*y1*y4 + a10*y2*y4 +
                a11*y3 + a12*x2*y3 + a13*y3^2 + a14/(x2+1) +
                a15*y3*y4 + a16*y1*y4*x2 + a17*y1*y3*y4 + a18*x1*y3 +
                a19*y1*y3 + a20*exp(a21*y1)

            f
        end

        function eval_cons(x::AbstractArray{T}) where T
            x1 = x[1]
            x2 = x[2]
            y1 = x1*x2
            y2 = y1*x1
            y3 = x2^2
            y4 = x1^2

            g = zeros(T, 3)
            g[1] = 1 - y1/700.0
            g[2] = y4/25.0^2 - x2/5.0
            g[3] = (x1/500.0- 0.11) - (x2/50.0-1)^2

            g
        end

        fres = DiffResults.HessianResult(zeros(2))
        fcfg = ForwardDiff.HessianConfig(eval_obj, fres, zeros(2))
        gres = DiffResults.JacobianResult(zeros(3), zeros(2))
        gcfg = ForwardDiff.JacobianConfig(eval_cons, zeros(2))

        return Bilevel.generate_autodiff_solver_fn(eval_obj,fres,fcfg,eval_cons,gres,gcfg,[],1:3)
    end

    solver_fn = get_barnesgrad()

    x0 = [10., 10.]
    options = Dict{String, Any}()
    options["Derivative option"] = 1
    options["Verify level"] = 1
    xopt, info = Bilevel.snopt(solver_fn, 0, 3, x0, options)

    xopt
end