text:
mutable struct MainState
    src::CodeInfo
    pc::Int
    stmt::Any
end

mutable struct GateState
    src::CodeInfo
    pc::Int
    stmt::Any
    name::Symbol
    cargs::Dict{Int, Int}
    cargnames::Vector{SubString{String}}
end

function YaoCompiler.compile(target::OpenQASMTarget, f, tt::Type, options::HardwareFreeOptions)
    v"2.0" ≤ target.version < v"3.0" || error("only QASM 2.x is supported")
    interp = YaoInterpreter(;target, options)
    mi = method_instance(f, tt)
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)
    # we need to remove code coverage effect to pass validate
    src = CompilerPluginTools.rm_code_coverage_effect(src)
    validate(target, src)

    if target.toplevel
        st = MainState(src, 0, nothing)
    else
        op = tt.parameters[2]
        op.parameters[1] <: GenericRoutine || error("OpenQASM gate target does not support callable object")
        name = routine_name(op)
        cargnames = retrieve_cargnames(op)
        cargs = extract_cargs_map(src)
        st = GateState(src, 0, nothing, name, cargs, cargnames)
    end
    return emit_qasm(target, src, RegInfo(target, src), st)
end

function push_stmt!(list::Vector, stmt)
    isnothing(stmt) && return list
    if stmt isa Vector
        append!(list, stmt)
    else
        push!(list, stmt)
    end
end

function emit_qasm(target::OpenQASMTarget, ci::CodeInfo, ri::RegInfo, st::MainState)
    cregs, qregs = emit_registers(ri)
    prog = []
    if target.include_qelib1
        push!(prog, Include(qasm_str("qelib1.inc")))
    end
    append!(prog, cregs)
    append!(prog, qregs)

    st.pc = 1
    while st.pc <= length(ci.code)
        push_stmt!(prog, emit_stmt(target, ci.code[st.pc], ri, st))
        st.pc += 1
    end

    return MainProgram(target.version, prog)
end

function emit_qasm(target::OpenQASMTarget, ci::CodeInfo, ri::RegInfo, st::GateState)
    qargs = emit_qargs(ri)
    cargs = qasm_id.(st.cargnames)
    decl = GateDecl(qasm_id(st.name), cargs, qargs)
    prog = []
    st.pc = 1
    while st.pc <= length(ci.code)
        push_stmt!(prog, emit_stmt(target, ci.code[st.pc], ri, st))
        st.pc += 1
    end
    return Gate(decl, prog)
end

qreg_name(k::Int) = string("qreg_", k)

function emit_registers(ri::RegInfo)
    cregs = []
    for (name, size) in ri.creg_size
        push!(cregs, RegDecl(kw_creg, qasm_id(name), qasm_int(size)))
    end

    qregs = []
    for (qreg, locs) in ri.qreg_to_locs
        qreg_id = qasm_id(qreg_name(qreg))
        push!(qregs, RegDecl(kw_qreg, qreg_id, qasm_int(length(locs))))
    end
    return cregs, qregs
end

function emit_stmt(target::OpenQASMTarget, stmt, ri::RegInfo, st)
    @switch stmt begin
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, gate, QuoteNode(locs))
            emit_gate(target, gate, locs, ri, st)
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, gate, QuoteNode(locs), QuoteNode(ctrl))
            emit_ctrl(target, gate, locs, ctrl, ri, st)
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs))
            emit_measure(target, locs, ri, st)
        @case Expr(:invoke, _, GlobalRef(Intrinsics, :barrier), _, QuoteNode(locs))
            emit_barrier(target, locs, ri, st)
        @case GotoIfNot(SSAValue(id), dest)
            emit_goto_ifnot(target, id, dest, ri, st)
        @case _ # skip other statement
            nothing
    end
end

emit_qargs(ri::RegInfo) = [qasm_id(qreg_name(qreg)) for (qreg, _) in ri.qreg_to_locs]

function emit_exp(target::OpenQASMTarget, exp, ri::RegInfo, st)
    @switch exp begin
        @case ::Int
            qasm_int(exp)
        @case ::Irrational{:π}
            pi_token
        @case ::Real # convert other read number to f64
            qasm_f64(exp)
        @case SSAValue(id)
            emit_exp(target, st.src.code[id], ri, st)

        @case Expr(:invoke, mi, GlobalRef(_, :sin), arg) || Expr(:invoke, mi, &(QASM.sin), arg)
            Call(:sin, emit_exp(target, arg, ri, st))

        @case Expr(:invoke, mi, GlobalRef(_, :cos), arg) || Expr(:invoke, mi, &(QASM.cos), arg)
            Call(:cos, emit_exp(target, arg, ri, st))

        @case Expr(:invoke, mi, GlobalRef(_, :tan), arg) || Expr(:invoke, mi, &(QASM.tan), arg)
            Call(:tan, emit_exp(target, arg, ri, st))

        @case Expr(:invoke, mi, GlobalRef(_, :exp), arg) || Expr(:invoke, mi, &(QASM.exp), arg)
            Call(:exp, emit_exp(target, arg, ri, st))

        @case Expr(:invoke, mi, GlobalRef(_, :log), arg) || Expr(:invoke, mi, &(QASM.log), arg)
            Call(:ln, emit_exp(target, arg, ri, st))

        @case Expr(:invoke, mi, GlobalRef(_, :sqrt), arg) || Expr(:invoke, mi, &(QASM.sqrt), arg)
            Call(:sqrt, emit_exp(target, arg, ri, st))

        # unary
        @case Expr(:call, GlobalRef(mod, name), x)
            f = Core.Compiler.abstract_eval_global(mod, name)::Const
            sym = qasm_compatible_intrinsics[f.val]
            token = Token{:reserved}(sym)
            if sym == "-"
                Neg(emit_exp(target, x, ri, st))
            else
                Call(token, emit_exp(target, x, ri, st))
            end
        # binop
        @case Expr(:call, GlobalRef(mod, name), lhs, rhs)
            f = Core.Compiler.abstract_eval_global(mod, name)::Const
            token = Token{:reserved}(qasm_compatible_intrinsics[f.val])
            (emit_exp(target, lhs, ri, st), token, emit_exp(target, rhs, ri, st))
        # get k-th carg
        @case Expr(:call, GlobalRef(&Base, :getfield), &(SSAValue(1)), k::Int, true)
            Token{:id}(st.cargnames[k])
        @case _
            error("unknown expression: $exp")
    end
end

function emit_gate(target::OpenQASMTarget, @nospecialize(gate), locs::Locations, ri::RegInfo, st)
    qargs = [ri[loc] for loc in locs]
    if target.inline_intrinsics
        emit_gate_inline(target, gate, qargs, ri, st)
    else
        emit_gate_noinline(target, gate, qargs, ri, st)
    end
end

function emit_gate_noinline(target::OpenQASMTarget, @nospecialize(gate), qargs, ri::RegInfo, st)
    @switch gate begin
        @case QuoteNode(shift(θ))
            cargs = Any[emit_exp(target, θ, ri, st)]
            Instruction("shift", cargs, qargs)
        @case QuoteNode(Rx(θ))
            cargs = Any[emit_exp(target, θ, ri, st)]
            Instruction("rx", cargs, qargs)
        @case QuoteNode(Ry(θ))
            cargs = Any[emit_exp(target, θ, ri, st)]
            Instruction("ry", cargs, qargs)
        @case QuoteNode(Rz(θ))
            cargs = Any[emit_exp(target, θ, ri, st)]
            Instruction("rz", cargs, qargs)

        @case QuoteNode(val::IntrinsicRoutine)
            name = qasm_gate_name(typeof(val))
            Instruction(name, [], qargs)
        @case QuoteNode(val::Operation)
            name = string(routine_name(val))
            cargs = [emit_exp(target, x, ri, st) for x in val.args]
            Instruction(name, cargs, qargs)
        @case SSAValue(id) # parametric intrinsic gate
            emit_gate_parametric(target, id, qargs, ri, st)
        @case _
            error("unknown gate value: $gate")
    end
end

function emit_gate_inline(target::OpenQASMTarget, @nospecialize(gate), qargs, ri::RegInfo, st)
    @switch gate begin
        @case QuoteNode(shift(θ))
            intrinsic_qasm(typeof(shift), θ, qargs)
        @case QuoteNode(Rx(θ))
            intrinsic_qasm(typeof(Rx), θ, qargs)
        @case QuoteNode(Ry(θ))
            intrinsic_qasm(typeof(Ry), θ, qargs)
        @case QuoteNode(Rz(θ))
            intrinsic_qasm(typeof(Rz), θ, qargs)
        @case QuoteNode(val) # non parametric intrinsic
            intrinsic_qasm(typeof(val), qargs)
        @case SSAValue(id) # parametric intrinsic gate
            emit_gate_parametric(target, id, qargs, ri, st)
        @case _
            error("unknown gate value: $gate")
    end
end

function emit_gate_parametric(target::OpenQASMTarget, id, qargs, ri::RegInfo, st)
    target.toplevel && error("toplevel program cannot have non-constant gate parameter")
    stmt = st.src.code[id]

    if target.inline_intrinsics
        @switch stmt begin
            @case Expr(:new, gt::Type{<:IntrinsicRoutine}, theta)
                carg = emit_exp(target, theta, ri, st)
                return intrinsic_qasm(gt, carg, qargs)
            @case _
                error("expect gate definition, got $stmt")
        end
    else
        @switch stmt begin
            @case Expr(:new, gt::Type{<:IntrinsicRoutine}, theta)
                name = qasm_gate_name(gt)
            @case _
                error("expect gate definition, got $stmt")
        end
        carg = emit_exp(target, theta, ri, st)
        return Instruction(name, Any[carg], qargs)
    end
end

function emit_ctrl(::OpenQASMTarget, @nospecialize(gate), locs::Locations, ctrl::CtrlLocations, ri::RegInfo, st)
    @switch gate begin
        @case QuoteNode(&X)
            if length(ctrl) == 1 && length(locs) == 1
                return CXGate(ri[ctrl], ri[locs])
            elseif length(ctrl) == 2 && length(locs) == 1
                error(
                    "ccx is not valid intrinsic QASM instruction. posssible fix:\n",
                    "1. include a qasm stdlib e.g qelib<version>.inc and use ccx gate from it\n",
                    "2. define your own ccx gate via @device using CX and single qubit gate",
                )
            end
        @case _
            error("control statement is not compatible with QASM")
    end
end

function emit_measure(::OpenQASMTarget, locs::Locations, ri::RegInfo, st)
    cname = string(ri.ssa_creg_map[st.pc])
    return [Measure(ri[loc], Bit(cname, k-1)) for (k, loc) in enumerate(locs)]
end

function emit_barrier(::OpenQASMTarget, locs::Locations, ri::RegInfo, st)
    qargs = []; args = Dict{Int, Vector{Int}}()

    # scan location addrs to remove extra
    # index on qreg
    for loc in locs
        plain_loc = plain(loc)
        qreg = ri.locs_to_qreg[plain_loc]
        addr = ri.locs_to_addr[plain_loc]
        addrs = get!(args, qreg, Int[])
        push!(addrs, addr)
    end

    for (qreg, addrs) in args
        # do not index qreg explicitly if barrier size
        # is the same with register size
        if length(ri.qreg_to_locs[qreg]) == length(addrs)
            push!(qargs, Bit(qreg_name(qreg)))
        else
            for addr in addrs
                push!(qargs, Bit(qreg_name(qreg), addr))
            end
        end
    end
    return Barrier(qargs)
end

function emit_goto_ifnot(target::OpenQASMTarget, cond::Int, dest::Int, ri::RegInfo, st)
    @switch st.src.code[cond] begin
        @case Expr(:invoke, mi, GlobalRef(YaoCompiler, :measure_cmp), SSAValue(lhs), rhs::Int)
            creg, int = lhs, rhs
        @case Expr(:invoke, mi, GlobalRef(YaoCompiler, :measure_cmp), lhs::Int, SSAValue(rhs))
            int, creg = lhs, rhs
        @case _
            error("invalid statement $(st.src.code[cond])")
    end

    haskey(ri.ssa_creg_map, creg) || error("cannot have classical statement in toplevel QASM")
    cname = ri.ssa_creg_map[creg]
    stmt = move_to_next_quantum_statement(st)
    body = emit_stmt(target, stmt, ri, st)
    return IfStmt(qasm_id(cname), qasm_int(int), body)
end

function move_to_next_quantum_statement(st)
    ci = st.src
    pc′ = st.pc + 1
    local stmt
    while pc′ <= length(ci.code)
        stmt = ci.code[pc′]
        @switch stmt begin
            @case Expr(:invoke, _, GlobalRef(&Intrinsics, _), _...)
                break
            @case ::GotoNode
                pc′ = (stmt::Core.GotoNode).label
            @case _
                # NOTE:
                # we can allow some constant statements
                # but not control flows or other nodes
                error("unexpected statement inside control flow: $stmt")
        end
        pc′ += 1
    end

    st.pc = pc′
    st.stmt = stmt
    return stmt
end

function retrieve_cargnames(op)
    mi = method_instance(op.parameters[1].instance, op.parameters[2])
    return split(chop(mi.def.slot_syms), '\0')[2:end]
end

function extract_cargs_map(ci::CodeInfo)
    vars = findfirst(ci.code) do stmt
        @match stmt begin
            Expr(:call, GlobalRef(&Base, :getfield), Argument(3), QuoteNode(:args)) => true
            _ => false
        end
    end

    ssa_cargs = Dict{Int, Int}()
    isnothing(vars) && return ssa_cargs

    for (v, stmt) in enumerate(ci.code)
        @switch stmt begin
            @case Expr(:call, getfield, SSAValue(&vars), k, true)
                ssa_cargs[v] = k
            @case _
                continue
        end
    end
    return ssa_cargs
end


text:
const DATASETPARAMS = [:is_sparse, :max_bin, :data_random_seed, :categorical_feature]

const BOOSTERPARAMS = [:application, :learning_rate, :num_leaves, :max_depth, :tree_learner,
                       :num_threads, :histogram_pool_size, :min_data_in_leaf,
                       :min_sum_hessian_in_leaf, :lambda_l1, :lambda_l2, :min_gain_to_split,
                       :feature_fraction, :feature_fraction_seed, :bagging_fraction,
                       :bagging_freq, :bagging_seed, :early_stopping_round, :sigmoid,
                       :is_unbalance, :metric, :is_training_metric, :ndcg_at, :num_machines,
                       :local_listen_port, :time_out, :machine_list_file, :num_class]

const INDEXPARAMS = [:categorical_feature]

const MAXIMIZE_METRICS = ["auc", "ndcg"]

"""
    savemodel(estimator, filename; [num_iteration = -1])

Save the fitted model in `estimator` as `filename`.

# Arguments
* `estimator::LGBMEstimator`: the estimator to use in the prediction.
* `filename::String`: the name of the file to save the model in.
* `num_iteration::Integer`: keyword argument that sets the number of iterations of the model that
    should be saved. `< 0` for all iterations.
"""
function savemodel(estimator::LGBMEstimator, filename::String; num_iteration::Integer = -1)
    @assert(estimator.booster.handle != C_NULL, "Estimator does not contain a fitted model.")
    LGBM_BoosterSaveModel(estimator.booster, num_iteration, filename)
    return nothing
end

"""
    loadmodel(estimator, filename)

Load the fitted model `filename` into `estimator`. Note that this only loads the fitted model—not
the parameters or data of the estimator whose model was saved as `filename`.

# Arguments
* `estimator::LGBMEstimator`: the estimator to use in the prediction.
* `filename::String`: the name of the file that contains the model.
"""
function loadmodel(estimator::LGBMEstimator, filename::String)
    estimator.booster = LGBM_BoosterCreateFromModelfile(filename)
    return nothing
end

function log_fatal(verbosity::Integer, msg...)
    warn(msg...)
end

function log_warning(verbosity::Integer, msg...)
    verbosity >= 0 && warn(msg...)
end

function log_info(verbosity::Integer, msg...)
    verbosity >= 1 && print(msg...)
end

function log_debug(verbosity::Integer, msg...)
    verbosity >= 2 && print(msg...)
end

function shrinkresults!(results, last_retained_iter::Integer)
    for test_key in keys(results)
        test = results[test_key]
        for metric_key in keys(test)
            test[metric_key] = test[metric_key][1:last_retained_iter]
        end
    end
    return nothing
end

ext: jl


text:
#!/usr/bin/env julia

# Positive: number of ticks. Negative: Distance in hours between ticks
const NUMBER_OF_TICKS = -3

using Colors
# For color preview, look here: https://www.december.com/html/spec/colorsvg.html

include("src/structures.jl")

# if no size is given, this size is used for the plot
const DEFAULT_SIZE = (1400, 430)

# make nights slightly blue
const DRAW_NIGHT_BACKGROUND = true

# List of all commands and keywords to generate graphs for
global QUERY = Array{QWidget,1}([

    QPlot(
        "Network",
        days = 3,
        data = [
            QGroup(
                "kB",
                log = true,
                min = 1,
                data = [QStack([Sysstat("txpck/s"), Sysstat("rxpck/s")], ctxs = [])],
            ),
            QGroup("%", min = 0, max = 100, ctxs = [], data = [Sysstat("%ifutil")]),
        ],
    ),

    QPlot(
        "I/O",
        days = 3,
        data = [
            QGroup(
                "packets/s",
                log = true,
                data = [
                    QStack([
                        QStyle(Sysstat("wtps"), Dict("color" => :red)),
                        QStyle(Sysstat("rtps"), Dict("color" => :royalblue)),
                    ]),
                    # QStyle(Sysstat("tps"), Dict("color" => :navy))
                ],
            ),
            QGroup(
                "blocks/s",
                data = [QStyle(
                    QStack([
                        QStyle(Sysstat("bread/s"), Dict("color" => :navy)),
                        QStyle(Sysstat("bwrtn/s"), Dict("color" => :maroon)),
                    ]),
                    Dict("fillalpha" => 0.4),
                )],
            ),
        ],
    ),

    QPlot(
        "CPU Long",
        days = 30,
        data = [QGroup(
            "%",
            min = 0,
            max = 100,
            data = [
                QStack([
                    Sysstat("%steal"),
                    Sysstat("%nice"),
                    Sysstat("%iowait"),
                    Sysstat("%system"),
                    Sysstat("%user"),
                ]),
                QStyle(Sysstat("%memused"), Dict("color" => :red)),
            ],
        )],
    ),

    QPlot(
        "CPU",
        days = 3,
        data = [QGroup(
            "%",
            min = 0,
            max = 100,
            data = [
                QStack([
                    Sysstat("%steal"),
                    Sysstat("%nice"),
                    Sysstat("%iowait"),
                    Sysstat("%system"),
                    Sysstat("%user"),
                ]),
                QStyle(Sysstat("%memused"), Dict("color" => :red)),
            ],
        )],
    ),

    QPlot(
        "RAM",
        days = 3,
        data = [
            QGroup(
                "kB",
                data = [
                    Sysstat("kbdirty"),
                    QStack([
                        #Sysstat("kbswpused"), # TODO: I want this!
                        Sysstat("kbinact"),
                        Sysstat("kbactive"),
                        QStyle(
                            Sysstat("kbavail"),
                            Dict("color" => :lightblue, "fillrange" => nothing),
                        ),
                    ]),
                ],
            ),
            QGroup(
                "%",
                min = 0,
                data = [
                #QStyle(Sysstat("%memused"), Dict("color"=>:darkmagenta)),
                QStyle(Sysstat("%commit"), Dict("color" => :red))],
            ),
        ],
    ),

    QPlot(
        "Handles",
        days = 3,
        data = [
            QGroup(
                "Number",
                min = 0,
                max = 50,
                data = [QStack([QStyle(
                    Sysstat("pty-nr"),
                    Dict("color" => :darkmagenta, "fillalpha" => 0.3),
                )])],
            ),
            QGroup(
                "Number",
                min = 0,
                data = [
                    QStack([Sysstat("file-nr"), Sysstat("inode-nr")]),
                    Sysstat("dentunusd"),
                ],
            ),
        ],
    ),

    QPlot(
        "Task Creation & Context Switches",
        days = 3,
        data = [
            QGroup("Number / s", min = 0, data = [QStack([Sysstat("cswch/s")])]),
            QGroup("Number / s", min = 0, data = [QStack([Sysstat("proc/s")])]),
        ],
    ),

    QPlot(
        "Network failures",
        days = 3,
        data = [
            QGroup(
                "Number / s",
                min = 0,
                ctxs = [],
                data = [
                    QStack([
                        Sysstat("coll/s"),
                        Sysstat("rxerr/s"),
                        Sysstat("txerr/s"),
                        Sysstat("rxdrop/s"),
                        Sysstat("txdrop/s"),
                        Sysstat("txcarr/s"),
                        Sysstat("rxfram/s"),
                        Sysstat("rxfifo/s"),
                        Sysstat("txfifo/s"),
                    ]),
                ],
            ),
        ],
    ),

    QPlot(
        "Sockets",
        days = 3,
        data = [
            QGroup(
                "Number",
                min = 0,
                data = [
                    Sysstat("totsck"),
                    QStack([Sysstat("tcpsck"), Sysstat("udpsck"), Sysstat("rawsck")]),
                    Sysstat("tcp-tw"),
                ],
            ),
            QGroup("Number", min = 0, data = [QStack([Sysstat("ip-frag")])]),
        ],
    ),

    QPlot(
        "Disk I/O",
        days = 3,
        data = [
            QGroup(
                "kB",
                log = true,
                data = [QStack([Sysstat("wkB/s"), Sysstat("rkB/s")], ctxs = [])],
            ),
            QGroup("%", min = 0, max = 100, ctxs = [], data = [QStack([Sysstat("%util")])]),
        ],
    ),

    #=
    QPlot(
        "AWS UsageQuantity",
        days = 30,
        data = [QGroup("X",
            min = 0,
            data = [
                QStack([
                    Aws("UsageQuantity AWS Lambda"),
                    Aws("UsageQuantity AWS Key Management Service"),
                    Aws("UsageQuantity AWS Cost Explorer"),
                    Aws("UsageQuantity AWS Budgets"),
                    Aws("UsageQuantity Tax"),
                    Aws("UsageQuantity AmazonCloudWatch"),
                    Aws("UsageQuantity Amazon Simple Storage Service"),
                    Aws("UsageQuantity Amazon Simple Notification Service"),
                    Aws("UsageQuantity Amazon DynamoDB"),
                ]),
            ],
        )],
    ),

    QPlot(
        "AWS BlendedCost",
        days = 30,
        data = [QGroup("X",
            min = 0,
            max = 0.2, # TODO: "minmax" - scale it down if there isn't much so if the value increases it stands out immediately!
            data = [
                QStack([
                    Aws("BlendedCost AWS Lambda"),
                    Aws("BlendedCost AWS Key Management Service"),
                    Aws("BlendedCost AWS Cost Explorer"),
                    Aws("BlendedCost AWS Budgets"),
                    Aws("BlendedCost Tax"),
                    Aws("BlendedCost AmazonCloudWatch"),
                    Aws("BlendedCost Amazon Simple Storage Service"),
                    Aws("BlendedCost Amazon Simple Notification Service"),
                    Aws("BlendedCost Amazon DynamoDB"),
                ]),
            ],
        )],
    ),

    =#
])

# Color Theme
# See https://github.com/JuliaPlots/PlotThemes.jl
global PLOT_THEME = :default # :wong

# Show the last n days in the plot.
const SHOW_DAYS = 3

const MAKE_OVERVIEW = true
const OVERVIEW_WIDTH = DEFAULT_SIZE[1] ÷ 3
const OVERVIEW_HEIGHT = DEFAULT_SIZE[2] ÷ 3

# Sysstat Corruption detection
const SYSSTAT_CORRUPTION_THRESHOLD = 1_000_000_000_000
const SYSSTAT_CORRUPTION_MARGIN = 4

const USE_PNGQUANT = true

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
const sym2loader = Dict{Symbol,Vector{Symbol}}()
const sym2saver  = Dict{Symbol,Vector{Symbol}}()

function is_installed(pkg::Symbol)
    isdefined(pkg) && isa(getfield(Main, pkg), Module) && return true # is a module defined in Main scope
    path = Base.find_in_path(string(pkg)) # hacky way to determine if a Package is installed
    path == nothing && return false
    return isfile(path)
end

function checked_import(pkg::Symbol)
    !is_installed(pkg)      && throw(NotInstalledError(pkg, ""))
    !isdefined(Main, pkg)   && eval(Main, Expr(:import, pkg))
    return getfield(Main, pkg)
end

for (applicable_, add_, dict_) in (
        (:applicable_loaders, :add_loader, :sym2loader),
        (:applicable_savers,  :add_saver,  :sym2saver))
    @eval begin
        $applicable_{sym}(::Union{Type{DataFormat{sym}}, Formatted{DataFormat{sym}}}) = get($dict_, sym, [:FileIO]) # if no loader is declared, fallback to FileIO
        function $add_{sym}(::Type{DataFormat{sym}}, pkg::Symbol)
            list = get($dict_, sym, Symbol[])
            $dict_[sym] = push!(list, pkg)
        end
    end
end

"`add_loader(fmt, :Package)` triggers `using Package` before loading format `fmt`"
add_loader
"`add_saver(fmt, :Package)` triggers `using Package` before saving format `fmt`"
add_saver

"""
- `load(filename)` loads the contents of a formatted file, trying to infer
the format from `filename` and/or magic bytes in the file.
- `load(strm)` loads from an `IOStream` or similar object. In this case,
the magic bytes are essential.
- `load(File(format"PNG",filename))` specifies the format directly, and bypasses inference.
- `load(f; options...)` passes keyword arguments on to the loader.
"""
load(s::Union{AbstractString,IO}, args...; options...) =
    load(query(s), args...; options...)

"""
- `save(filename, data...)` saves the contents of a formatted file,
trying to infer the format from `filename`.
- `save(Stream(format"PNG",io), data...)` specifies the format directly, and bypasses inference.
- `save(f, data...; options...)` passes keyword arguments on to the saver.
"""
save(s::Union{AbstractString,IO}, data...; options...) =
    save(query(s), data...; options...)

# Forced format
function save{sym}(df::Type{DataFormat{sym}}, f::AbstractString, data...; options...)
    libraries = applicable_savers(df)
    checked_import(libraries[1])
    eval(Main, :($save($File($(DataFormat{sym}), $f),
                       $data...; $options...)))
end

function save{sym}(df::Type{DataFormat{sym}}, s::IO, data...; options...)
    libraries = applicable_savers(df)
    checked_import(libraries[1])
    eval(Main, :($save($Stream($(DataFormat{sym}), $s),
                       $data...; $options...)))
end

# Fallbacks
function load{F}(q::Formatted{F}, args...; options...)
    if unknown(q)
        isfile(filename(q)) || open(filename(q))  # force systemerror
        throw(UnknownFormat(q))
    end
    libraries = applicable_loaders(q)
    failures  = Any[]
    for library in libraries
        try
            Library = checked_import(library)
            if !has_method_from(methods(Library.load), Library)
                throw(LoaderError(string(library), "load not defined"))
            end
            return eval(Main, :($(Library.load)($q, $args...; $options...)))
        catch e
            push!(failures, (e, q))
        end
    end
    handle_exceptions(failures, "loading \"$(filename(q))\"")
end
function save{F}(q::Formatted{F}, data...; options...)
    unknown(q) && throw(UnknownFormat(q))
    libraries = applicable_savers(q)
    failures  = Any[]
    for library in libraries
        try
            Library = checked_import(library)
            if !has_method_from(methods(Library.save), Library)
                throw(WriterError(string(library), "save not defined"))
            end
            return eval(Main, :($(Library.save)($q, $data...; $options...)))
        catch e
            push!(failures, (e, q))
        end
    end
    handle_exceptions(failures, "saving \"$(filename(q))\"")
end

function has_method_from(mt, Library)
    for m in mt
        if getmodule(m) == Library
            return true
        end
    end
    false
end

if VERSION < v"0.5.0-dev+3543"
    getmodule(m) = m.func.code.module
else
    getmodule(m) = m.module
end


text:
"""
	euler_characteristic(V::Lar.Points, copEV::Lar.ChainOp,	copFE::Lar.ChainOp)

Return the Euler Characteristic number of a 2D- or 3D-model.
`V` is given by rows.
"""
function euler_characteristic(
		V::Lar.Points, # by rows
		copEV::Lar.ChainOp,
		copFE::Lar.ChainOp)

	# χ = size(V,1) - size(copEV,1) + size(copFE,1)
	# return χ
	return euler_characteristic_general(V, copEV, copFE)
end

"""
	euler_characteristic_general(V::Lar.Points, cc::Lar.ChainOp...)

Return the Euler Characteristic number of a n-dimensional model.
`V` is given by rows.
"""
function euler_characteristic_general(V::Lar.Points, cc::Lar.ChainOp...)
	χ = size(V,1)
	for k in 1:length(cc)
		χ = χ + (-1)^k * size(cc[k],1)
	end
	return χ
end


text:
const DOC_PLOT_PARAMS = """
- **`title`** : Text to display on the top of the plot.

- **`xlabel`** : Text to display on the x axis of the plot

- **`ylabel`** : Text to display on the y axis of the plot

- **`labels`** : Boolean. Can be used to hide the labels by
  setting `labels = false`.

- **`border`** : The style of the bounding box of the plot.
  Supports `:corners`, `:solid`, `:bold`, `:dashed`, `:dotted`,
  `:ascii`, and `:none`.

- **`margin`** : Number of empty characters to the left of the
  whole plot.

- **`padding`** : Space of the left and right of the plot between
  the labels and the canvas.

- **`color`** : Color of the drawing.
  Can be any of `:green`, `:blue`, `:red`, `:yellow`, `:cyan`,
  `:magenta`, `:white`, `:normal`, an integer in the range `0`-`255` for Color256 palette,
  or a tuple of three integers as RGB components.

- **`width`** : Number of characters per row that should be used
  for plotting.
"""

transform_name(::typeof(identity), basename = "") = basename
function transform_name(f, basename = "")
    fname = string(nameof(f))
    fname = occursin("#", fname) ? "custom" : fname
    string(basename, " [", fname, "]")
end

roundable(num::Number) = isinteger(num) & (typemin(Int) <= num < typemax(Int))
compact_repr(num::Number) = repr(num, context=:compact => true)

ceil_neg_log10(x) = roundable(-log10(x)) ? ceil(Integer, -log10(x)) : floor(Integer, -log10(x))
round_neg_log10(x) = roundable(-log10(x)) ? round(Integer, -log10(x), RoundNearestTiesUp) : floor(Integer, -log10(x))
round_up_tick(x,m) = x == 0. ? 0. : (x > 0 ? ceil(x, digits=ceil_neg_log10(m)) : -floor(-x, digits=ceil_neg_log10(m)))
round_down_tick(x,m) = x == 0. ? 0. : (x > 0 ? floor(x, digits=ceil_neg_log10(m)) : -ceil(-x, digits=ceil_neg_log10(m)))
round_up_subtick(x,m) = x == 0. ? 0. : (x > 0 ? ceil(x, digits=ceil_neg_log10(m)+1) : -floor(-x, digits=ceil_neg_log10(m)+1))
round_down_subtick(x,m) = x == 0. ? 0. : (x > 0 ? floor(x, digits=ceil_neg_log10(m)+1) : -ceil(-x, digits=ceil_neg_log10(m)+1))
float_round_log10(x::F,m) where {F<:AbstractFloat} = x == 0. ? F(0) : (x > 0 ? round(x, digits=ceil_neg_log10(m)+1)::F : -round(-x, digits=ceil_neg_log10(m)+1)::F)
float_round_log10(x::Integer,m) = float_round_log10(float(x), m)
float_round_log10(x) = x > 0 ? float_round_log10(x,x) : float_round_log10(x,-x)

function plotting_range(xmin, xmax)
    diff = xmax - xmin
    xmax = round_up_tick(xmax, diff)
    xmin = round_down_tick(xmin, diff)
    Float64(xmin), Float64(xmax)
end

function plotting_range_narrow(xmin, xmax)
    diff = xmax - xmin
    xmax = round_up_subtick(xmax, diff)
    xmin = round_down_subtick(xmin, diff)
    Float64(xmin), Float64(xmax)
end

function extend_limits(vec, limits)
    mi, ma = map(Float64, extrema(limits))
    if mi == 0. && ma == 0.
        mi, ma = map(Float64, extrema(vec))
    end
    diff = ma - mi
    if diff == 0
        ma = mi + 1
        mi = mi - 1
    end
    (limits == (0.,0.) || limits == [0.,0.]) ? plotting_range_narrow(mi, ma) : (mi, ma)
end

sort_by_keys(dict::Dict) = sort!(collect(dict), by=x->x[1])

function sorted_keys_values(dict::Dict; k2s=true)
    if k2s  # check and force key type to be of AbstractString type if necessary
        kt, vt = eltype(dict).types
        if !(kt <: AbstractString)
            dict = Dict(string(k) => v  for (k, v) in pairs(dict))
        end
    end
    keys_vals = sort_by_keys(dict)
    first.(keys_vals), last.(keys_vals)
end

const bordermap = Dict{Symbol,Dict{Symbol,String}}()
const border_solid  = Dict{Symbol,String}()
const border_corners = Dict{Symbol,String}()
const border_barplot = Dict{Symbol,String}()
const border_bold   = Dict{Symbol,String}()
const border_none   = Dict{Symbol,String}()
const border_dashed = Dict{Symbol,String}()
const border_dotted = Dict{Symbol,String}()
const border_ascii  = Dict{Symbol,String}()
border_solid[:tl] = "┌"
border_solid[:tr] = "┐"
border_solid[:bl] = "└"
border_solid[:br] = "┘"
border_solid[:t] = "─"
border_solid[:l] = "│"
border_solid[:b] = "─"
border_solid[:r] = "│"
border_corners[:tl] = "┌"
border_corners[:tr] = "┐"
border_corners[:bl] = "└"
border_corners[:br] = "┘"
border_corners[:t] = " "
border_corners[:l] = " "
border_corners[:b] = " "
border_corners[:r] = " "
border_barplot[:tl] = "┌"
border_barplot[:tr] = "┐"
border_barplot[:bl] = "└"
border_barplot[:br] = "┘"
border_barplot[:t] = " "
border_barplot[:l] = "┤"
border_barplot[:b] = " "
border_barplot[:r] = " "
border_bold[:tl] = "┏"
border_bold[:tr] = "┓"
border_bold[:bl] = "┗"
border_bold[:br] = "┛"
border_bold[:t] = "━"
border_bold[:l] = "┃"
border_bold[:b] = "━"
border_bold[:r] = "┃"
border_none[:tl] = " "
border_none[:tr] = " "
border_none[:bl] = " "
border_none[:br] = " "
border_none[:t] = " "
border_none[:l] = " "
border_none[:b] = " "
border_none[:r] = " "
border_dashed[:tl] = "┌"
border_dashed[:tr] = "┐"
border_dashed[:bl] = "└"
border_dashed[:br] = "┘"
border_dashed[:t] = "╌"
border_dashed[:l] = "┊"
border_dashed[:b] = "╌"
border_dashed[:r] = "┊"
border_dotted[:tl] = "⡤"
border_dotted[:tr] = "⢤"
border_dotted[:bl] = "⠓"
border_dotted[:br] = "⠚"
border_dotted[:t] = "⠤"
border_dotted[:l] = "⡇"
border_dotted[:b] = "⠒"
border_dotted[:r] = "⢸"
border_ascii[:tl] = "+"
border_ascii[:tr] = "+"
border_ascii[:bl] = "+"
border_ascii[:br] = "+"
border_ascii[:t] = "-"
border_ascii[:l] = "|"
border_ascii[:b] = "-"
border_ascii[:r] = "|"
bordermap[:solid]  = border_solid
bordermap[:corners] = border_corners
bordermap[:barplot] = border_barplot
bordermap[:bold]   = border_bold
bordermap[:none]   = border_none
bordermap[:dashed] = border_dashed
bordermap[:dotted] = border_dotted
bordermap[:ascii]  = border_ascii

const UserColorType = Union{Integer,Symbol,NTuple{3,Integer},Nothing}  # allowed color type
const JuliaColorType = Union{Symbol,Int}  # color type for printstyled (defined in base/util.jl)
const ColorType = Union{Nothing,UInt8}  # internal UnicodePlots color type

const color_cycle = [:green, :blue, :red, :magenta, :yellow, :cyan]

function print_color(color::UserColorType, io::IO, args...)
    printstyled(io, string(args...); color = julia_color(color))
end

function crayon_256_color(color::UserColorType)::ColorType
    color in (:normal, :default, :nothing, nothing) && return nothing
    ansicolor = Crayons._parse_color(color)
    if ansicolor.style == Crayons.COLORS_16
        return Crayons.val(ansicolor) % 60
    elseif ansicolor.style == Crayons.COLORS_24BIT
        return Crayons.val(Crayons.to_256_colors(ansicolor))
    end
    Crayons.val(ansicolor)
end

function julia_color(color::UserColorType)::JuliaColorType
    if color isa Nothing
        :normal
    elseif color isa Symbol
        color
    elseif color isa Integer
        Int(color)
    else
        julia_color(crayon_256_color(color))
    end
end

@inline function set_color!(colors::Array{ColorType,2}, x::Int, y::Int, color::ColorType; force::Bool=false)
    if color === nothing || colors[x, y] === nothing || force
        colors[x, y] = color
    else
        colors[x, y] |= color
    end
    nothing
end


text:
#=
Main program for two phase flow inversion
=#

include("args.jl")

function sw_p_to_lambda_den(sw, p)
    sw = tf.reshape(sw, (1, m, n, 1))
    p = tf.reshape(p, (1, m, n, 1))
    sw = tf.image.resize_bilinear(sw, (nz, nx))
    p = tf.image.resize_bilinear(p, (nz, nx))
    sw = cast(sw, Float64)
    p = cast(p, Float64)
    sw = squeeze(sw)
    p = squeeze(p)
    # tran_lambda, tran_den = Gassman(sw)
    # tran_lambda, tran_den = RockLinear(sw) # test linear relationship
    # tran_lambda, tran_den = Patchy(sw)
    tran_lambda, tran_den = Brie(sw)
    tran_lambda_pad =  tf.pad(tran_lambda, [nPml (nPml+nPad); nPml nPml], constant_values=3500.0^2*2200.0/3.0) /1e6
    tran_den_pad = tf.pad(tran_den, [nPml (nPml+nPad); nPml nPml], constant_values=2200.0)
    return tran_lambda_pad, tran_den_pad
end

# NOTE Generate Data
if args["generate_data"]
    println("Generate Test Data...")
    K = 20.0 .* ones(m,n) # millidarcy
    K[8:10,:] .= 120.0
    # K[17:21,:] .= 100.0
    # for i = 1:m
    #     for j = 1:n
    #         if i <= (14 - 24)/(30 - 1)*(j-1) + 24 && i >= (12 - 18)/(30 - 1)*(j-1) + 18
    #             K[i,j] = 100.0
    #         end
    #     end
    # end
    tfCtxTrue = tfCtxGen(m,n,h,NT,Δt,Z,X,ρw,ρo,μw,μo,K,g,ϕ,qw,qo, sw0, true)
    out_sw_true, out_p_true = imseq(tfCtxTrue)
    lambdas = Array{PyObject}(undef, n_survey)
    dens = Array{PyObject}(undef, n_survey)
    for i = 1:n_survey
        sw = out_sw_true[survey_indices[i]]
        p = out_p_true[survey_indices[i]]
        lambdas[i], dens[i] = sw_p_to_lambda_den(sw, p)
    end

    misfit = Array{PyObject}(undef, n_survey)
    for i = 1:n_survey
        if !isdir("./$(args["version"])/Data$i")
            mkdir("./$(args["version"])/Data$i")
        end
        para_fname = "./$(args["version"])/para_file$i.json"
        survey_fname = "./$(args["version"])/survey_file$i.json"
        paraGen(nz_pad, nx_pad, dz, dx, nSteps, dt, f0, nPml, nPad, para_fname, survey_fname, "./$(args["version"])/Data$i/")
        # shot_inds = collect(1:3:length(z_src)) .+ mod(i-1,3) # 5src rotation
        # shot_inds = i # 1src rotation
        shot_inds = collect(1:length(z_src)) # all sources
        surveyGen(z_src[shot_inds], x_src[shot_inds], z_rec, x_rec, survey_fname)
        tf_shot_ids0 = constant(collect(0:length(shot_inds)-1), dtype=Int32)
        misfit[i] = fwi_obs_op(lambdas[i], tf_shear_pad, dens[i], tf_stf, tf_gpu_id0, tf_shot_ids0, para_fname)
    end
    config = tf.ConfigProto()
    config.intra_op_parallelism_threads = 24
    config.inter_op_parallelism_threads = 24
    sess = Session(config=config); init(sess);
    run(sess, misfit)
    error("Generate Data: Stop")
end

tfCtxInit = tfCtxGen(m,n,h,NT,Δt,Z,X,ρw,ρo,μw,μo,K_init,g,ϕ,qw,qo, sw0, false)
out_sw_init, out_p_init = imseq(tfCtxInit)
lambdas = Array{PyObject}(undef, n_survey)
dens = Array{PyObject}(undef, n_survey)
for i = 1:n_survey
    sw = out_sw_init[survey_indices[i]]
    p = out_p_init[survey_indices[i]]
    lambdas[i], dens[i] = sw_p_to_lambda_den(sw, p)
end

# NOTE Compute FWI loss
loss = constant(0.0)
for i = 1:n_survey
    global loss
    para_fname = "./$(args["version"])/para_file$i.json"
    survey_fname = "./$(args["version"])/survey_file$i.json"
    # shot_inds = collect(1:3:length(z_src)) .+ mod(i-1,3)
    # shot_inds = i
    shot_inds = collect(1:length(z_src)) # all sources
    tf_shot_ids0 = constant(collect(0:length(shot_inds)-1), dtype=Int32)
    loss += fwi_op(lambdas[i], tf_shear_pad, dens[i], tf_stf, tf_gpu_id_array[mod(i,nGpus)], tf_shot_ids0, para_fname) # mod(i,2)
end
gradK = gradients(loss, tfCtxInit.K)

if args["verbose"]
    sess = Session(); init(sess)
    println("Initial loss = ", run(sess, loss))
    g = gradients(loss, tfCtxInit.K)
    G = run(sess, g)
    pcolormesh(G); savefig("test.png"); close("all")
end

# Optimization
__cnt = 0
# invK = zeros(m,n)
function print_loss(l, K, gradK, brie_coef)
    global __cnt, __l, __K, __gradK, __brie_coef
    if mod(__cnt,1)==0
        println("\niter=$__iter, eval=$__cnt, current loss=",l)
        # println("a=$a, b1=$b1, b2=$b2")
    end
    __cnt += 1
    __l = l
    __K = K
    __gradK = gradK
    __brie_coef = brie_coef
end

__iter = 0
function print_iter(rk)
    global __iter, __l
    if mod(__iter,1)==0
        println("\n************* ITER=$__iter *************\n")
    end
    __iter += 1
    open("./$(args["version"])/loss.txt", "a") do io
        writedlm(io, Any[__iter __l])
    end
    open("./$(args["version"])/K$__iter.txt", "w") do io
        writedlm(io, __K)
    end
    open("./$(args["version"])/gradK$__iter.txt", "w") do io
        writedlm(io, __gradK)
    end
    open("./$(args["version"])/brie_coef.txt", "a") do io
        writedlm(io, Any[__iter __brie_coef])
    end
end

config = tf.ConfigProto()
config.intra_op_parallelism_threads = 24
config.inter_op_parallelism_threads = 24
sess = Session(config=config); init(sess);
# opt = ScipyOptimizerInterface(loss, var_list=[tfCtxInit.K], var_to_bounds=Dict(tfCtxInit.K=> (10.0, 130.0)), method="L-BFGS-B",
#     options=Dict("maxiter"=> 100, "ftol"=>1e-6, "gtol"=>1e-6))
opt = ScipyOptimizerInterface(loss, var_list=[tfCtxInit.K, tf_brie_coef], var_to_bounds=Dict(tfCtxInit.K=> (10.0, 130.0), tf_brie_coef=>(1.0,100.0)), method="L-BFGS-B",
    options=Dict("maxiter"=> 100, "ftol"=>1e-6, "gtol"=>1e-6))
@info "Optimization Starts..."
# ScipyOptimizerMinimize(sess, opt, loss_callback=print_loss, step_callback=print_iter, fetches=[loss,tfCtxInit.K,gradK])
ScipyOptimizerMinimize(sess, opt, loss_callback=print_loss, step_callback=print_iter, fetches=[loss,tfCtxInit.K,gradK, tf_brie_coef])


text:
# Naming convention:
# ------------------
# {type}[S]_{dims}[_d][_n]
# type: C: Continuous, D: Discrete
# scalar: add "S" if any input is scalar
# dims: "nxnuny"
# feedthrough: append "_d" if `D` is present
# names: append "_n" if some inputs/outputs/states are named
SS = HeteroStateSpace
@testset "test_statespace" begin
    # SCALARS
    for SS in (StateSpace, HeteroStateSpace)
        a_2 = [-5 -3; 2 -9]
        CS_111 = SS(-5, 2, 3, [0])
        CS_111_d = SS([3], 2, 1, 1)
        CS_211 = SS(a_2, [1; 2], [1 0], 0)
        CS_221 = SS(a_2, [1 0; 0 2], [1 0], 0)
        CS_222 = SS(a_2, [1 0; 0 2], eye_(2), 0)

        # CONTINUOUS
        a_1 = [-5]
        C_111 = SS(a_1, [2], [3], [0])
        C_211 = SS(a_2, [1; 2], [1 0], [0])
        C_212 = SS(a_2, [1; 2], eye_(2), [0; 0])
        C_221 = SS(a_2, [1 0; 0 2], [1 0], [0 0])
        C_222 = SS(a_2, [1 0; 0 2], eye_(2), zeros(Int,2,2))
        C_222_d = SS(a_2, [1 0; 0 2], eye_(2), eye_(2))
        C_022 = SS(4.0*eye_(2))

        # DISCRETE
        da_1 = [-0.5]
        da_2 = [0.2 -0.8; -0.8 0.07]
        D_111 = SS(da_1, [2], [3], [0], 0.005)
        D_211 = SS(da_2, [1; 2], [1 0], [0], 0.005)
        D_221 = SS(da_2, [1 0; 0 2], [1 0], [0 0], 0.005)
        D_222 = SS(da_2, [1 0; 0 2], eye_(2), zeros(2,2), 0.005)
        D_222_d = SS(da_2, [1 0; 0 2], eye_(2), eye_(2), 0.005)
        D_022 = SS(4.0*eye_(2), 0.005)

        # Definition of input, output and state names
        C_222_d_n = SS(a_2, [1 0; 0 2], eye_(2), eye_(2))

        # TESTS
        # Contstuct with scalars
        @test CS_111 == C_111
        @test CS_111_d == SS([3],[2],[1],[1])
        @test CS_211 == C_211
        @test CS_221 == C_221
        @test CS_222 == C_222

        # Addition
        @test C_111 + C_111 == SS([-5 0; 0 -5],[2; 2],[3 3],[0])
        @test C_222 + C_222 == SS([-5 -3 0 0; 2 -9 0 0; 0 0 -5 -3;
        0 0 2 -9],[1 0; 0 2; 1 0; 0 2], [1 0 1 0; 0 1 0 1],[0 0; 0 0])
        @test C_222 + 1 == SS([-5 -3; 2 -9],[1 0; 0 2],[1 0; 0 1],[1 1; 1 1])
        @test D_111 + D_111 == SS([-0.5 0; 0 -0.5],[2; 2],[3 3],[0], 0.005)

        # Subtraction
        @test C_111 - C_211 == SS([-5 0 0; 0 -5 -3; 0 2 -9],[2; 1; 2],[3 -1 -0],[0])
        @test 1 - C_222 == SS([-5 -3; 2 -9],[1 0; 0 2],[-1 -0; -0 -1],[1 1; 1 1])
        @test D_111 - D_211 == SS([-0.5 0 0; 0 0.2 -0.8; 0 -0.8 0.07],[2; 1; 2],
        [3 -1 -0],[0], 0.005)

        # Multiplication
        @test C_111 * C_221 == SS([-5 2 0; 0 -5 -3; 0 2 -9],
        [0 0; 1 0; 0 2],[3 0 0],[0 0])
        @test C_212 * C_111 == SS([-5 -3 3; 2 -9 6; 0 0 -5],
        [0; 0; 2],[1 0 0; 0 1 0],[0; 0])
        @test 4*C_222 == SS([-5 -3; 2 -9],[1 0; 0 2],[4 0; 0 4],[0 0; 0 0])
        @test D_111 * D_221 == SS([-0.5 2 0; 0 0.2 -0.8; 0 -0.8 0.07],
        [0 0; 1 0; 0 2],[3 0 0],[0 0],0.005)

        # Division
        @test 1/C_222_d == SS([-6 -3; 2 -11],[1 0; 0 2],[-1 0; -0 -1],[1 -0; 0 1])
        @test C_221/C_222_d == SS([-5 -3 -1 0; 2 -9 -0 -2; 0 0 -6 -3;
        0 0 2 -11],[1 0; 0 2; 1 0; 0 2],[1 0 0 0],[0 0])
        @test 1/D_222_d == SS([-0.8 -0.8; -0.8 -1.93],[1 0; 0 2],[-1 0; -0 -1],
        [1 -0; 0 1],0.005)

        # Indexing
        @test size(C_222) == (2, 2)
        @test size(C_212) == (2, 1)
        @test C_222[1,1] == SS([-5 -3; 2 -9],[1; 0],[1 0],[0])
        @test C_222[1:1,1] == SS([-5 -3; 2 -9],[1; 0],[1 0],[0])
        @test C_222[1,1:2] == C_221
        @test size(C_222[1,[]]) == (1,0)

        A = [-1.0 -2.0; 0.0 -1.0]
        B = [0.0; -2.0]
        C = [1.0 1.0]
        D = 1.0
        sys = SS(A, B, C, D)

        @test sys + 1.0 == SS(A, B, C, D + 1.0)
        @test 2.0 + sys == SS(A, B, C, D + 2.0)

        @test -sys == SS(A, B, -C, -D)

        # Printing
        if SS <: StateSpace
            res = ("StateSpace{Int64,Array{Int64,2}}\nA = \n -5  -3\n  2  -9\nB = \n 1  0\n 0  2\nC = \n 1  0\n 0  1\nD = \n 0  0\n 0  0\n\nContinuous-time state-space model")
            @test sprint(show, C_222) == res
            res = ("StateSpace{Float64,Array{Float64,2}}\nA = \n  0.2  -0.8 \n -0.8   0.07\nB = \n 1.0  0.0\n 0.0  2.0\nC = \n 1.0  0.0\n 0.0  1.0\nD = \n 0.0  0.0\n 0.0  0.0\n\nSample Time: 0.005 (seconds)\nDiscrete-time state-space model")
            @test sprint(show, D_222) == res
            res = ("StateSpace{Float64,Array{Float64,2}}\nD = \n 4.0  0.0\n 0.0  4.0\n\nStatic gain")
            @test sprint(show, C_022) == res
            res = "StateSpace{Float64,Array{Float64,2}}\nD = \n 4.0  0.0\n 0.0  4.0\n\nSample Time: 0.005 (seconds)\nStatic gain"
            @test sprint(show, D_022) == res
        end

    # Errors
        @test_throws ErrorException C_111 + C_222             # Dimension mismatch
        @test_throws ErrorException C_111 - C_222             # Dimension mismatch
        @test_throws ErrorException C_111 * C_222             # Dimension mismatch
        @test_throws ErrorException D_111 + C_111             # Sampling time mismatch
        @test_throws ErrorException D_111 - C_111             # Sampling time mismatch
        @test_throws ErrorException D_111 * C_111             # Sampling time mismatch
        D_diffTs = SS([1], [2], [3], [4], 0.1)
        @test_throws ErrorException D_111 + D_diffTs            # Sampling time mismatch
        @test_throws ErrorException D_111 - D_diffTs            # Sampling time mismatch
        @test_throws ErrorException D_111 * D_diffTs            # Sampling time mismatch
        @test_throws ErrorException 1/C_222                     # Not invertible
        @test_throws ErrorException 1/C_212                     # Not invertible
        @test_throws ErrorException SS([1 2], [1], [2], [3])      # Not square A
        @test_throws ErrorException SS([1], [2 0], [1], [2])      # I/0 dim mismatch
        @test_throws ErrorException SS([1], [2], [3 4], [1])      # I/0 dim mismatch
        @test_throws ErrorException SS([1], [2], [3], [4], -0.1)  # Negative samping time
        @test_throws ErrorException SS(eye_(2), eye_(2), eye_(2), [0]) # Dimension mismatch

    end

end


text:
########################################################################################
#
# jFlow: a 3-D saturated-unsaturated groundwater flow model
#
# (1) Methodology is the integral finite difference method (IFDM)
# (2) Pressure head in a cell is defined at z of node point (i.e., if P = 0, S is typically = 0.5)
# (3) Unsaturated properties defined with Van Genuchten model, or as simple unconfined flow
#
########################################################################################

using SparseArrays

### data types

mutable struct Material
    index::Int64
    Kh::Float64                         # saturated hydraulic conductivity components
    Kz::Float64
    Ss::Float64                         # specific storage
    phi::Float64                        # porosity
    alpha::Float64                      # Van Genuchten model parameters
    N::Float64
    m::Float64
end

mutable struct Cell
    index::Int64                        # index tracking number for cell (for initial struct array setup)
    x::Float64                          # grid cell center coordinates
    y::Float64
    z::Float64                          # node point elevation
    vol::Float64                        # cell volume
    b::Float64                          # cell vertical thickness (==0 for VG unsat'd properties; >0 for Sy formulation)
    P::Float64                          # pressure head
    S::Float64                          # saturation
    mat::Material                       # corresponding material
    connects::Array{Int64, 1}           # list of associated connection indices
    neighbors::Array{Int64, 1}          # list of connected cells
end

mutable struct Connection
    cell1::Int64                        # connected cells (index numbers)
    cell2::Int64
    del1::Float64                       # respective distances from cell centerpoints to interface
    del2::Float64
    A::Float64                          # interfacial area of the connection
    zConnect::Bool                      # connection is designated as vertically-oriented
	KMean::Float64 						# mean saturated hydraulic conductivity across connection
end

mutable struct Params
    endTime::Float64                    # end of simulation time
    dt0::Float64                        # initial (and minimum) time step
    dtMax::Float64                      # maximum time step
    dhMax::Float64                      # maximum head change in any model cell, per time step
    f::Float64                          # time step multiplier
    c::Float64                          # Crank-Nicolson implicit solution weighting factor
	mesh::Bool 							# read externally generated mesh (true) or generate from input file (false)
end

include("jFlowIOModule.jl")             # various input and output functions included in separate module

### hydraulic parameterization functions (Van Genuchten or linear models for unsaturated conditions)

function kr(cell::Cell)::Float64
    # relative permeability as a function of water saturation
    if cell.b == 0
        kRel = cell.S^0.5 * (1.0 - (1.0 - cell.S^(1.0/cell.mat.m))^cell.mat.m)^2    # Van Genuchten formulation
    else
        kRel = cell.S           # for an unconfined aquifer condition
    end
    return kRel
end

function Sat(cell::Cell)::Float64
    # (effective) saturation as a function of pressure head
    if cell.b == 0.0
        # Van Genuchten formulation
        if cell.P < 0
            S = (1. + abs(cell.mat.alpha * cell.P)^cell.mat.N)^(-cell.mat.m)
        else
            S = 1.0
        end
    else
        # unconfined aquifer formulation
        if cell.P < cell.b/2
            S = maximum([cell.P/cell.b + 0.5, 0.0])
        else
            S = 1.0
        end
    end
    return S
end

function dSdP(cell::Cell)::Float64
    # storage term a function of change in effective saturation per change in capillary pressure head
    if cell.b == 0
        dS = -(cell.mat.alpha*abs(cell.P))^cell.mat.N * ((cell.mat.alpha*abs(cell.P))^cell.mat.N + 1.0)^(-cell.mat.m-1.0) * cell.mat.m * cell.mat.N/cell.P
    else
        dS = 1.0/cell.b
    end
    return dS * cell.mat.phi
end

function Storage(cell::Cell)::Float64
    # return appropriate storage term, depending on saturation state of cell
    if Sat(cell) < 1.0
        stor = dSdP(cell) * cell.vol
    else
        stor = cell.mat.Ss * cell.vol
    end
    return stor
end

### utility functions

function SortCell(tempCell::Array{Cell, 1})::Array{Cell, 1}
    # insertion sort approach to create ordered array of cell structs
    cell = [tempCell[1]]
    seq = [tempCell[1].index]
    for i = 2:length(tempCell)
        index = length(findall(seq.<tempCell[i].index)) + 1
        insert!(seq, index, tempCell[i].index)
        insert!(cell, index, tempCell[i])
    end
    println("Sorted cell struct array.")
    return cell
end

function MatchConnects(cell::Array{Cell, 1}, connect::Array{Connection, 1})::Array{Cell, 1}
    # populate cell connection lists
    for (i, cn) in enumerate(connect)
        push!(cell[cn.cell1].connects, i)
        push!(cell[cn.cell2].connects, i)
        push!(cell[cn.cell1].neighbors, cn.cell2)
        push!(cell[cn.cell2].neighbors, cn.cell1)
    end
    return cell
end

function FindStressPeriod(t::Float64, stressPeriods::Array{Float64, 1})::Int64
    # index number of stress period to which "t" belongs
    return sum((t.>=stressPeriods)*1)
end

function BlendedK(connect::Array{Connection, 1}, cell::Array{Cell, 1})::Array{Connection, 1}
    # weighted harmonic mean for saturated hydraulic conductivity heterogeneous cell interfaces
	for cn in connect
		if cn.zConnect == false
			cn.KMean = (cn.del1 + cn.del2)/(cn.del1/cell[cn.cell1].mat.Kh + cn.del2/cell[cn.cell2].mat.Kh)
		else
			cn.KMean = (cn.del1 + cn.del2)/(cn.del1/cell[cn.cell1].mat.Kz + cn.del2/cell[cn.cell2].mat.Kz)
		end
	end
    return connect
end

function EffectiveK(cn::Connection, cell::Array{Cell, 1})::Float64
    # combo weighted harmonic mean & arithmetic means for unsaturated cell interfaces
    krMean = (cn.del1*kr(cell[cn.cell1]) + cn.del2*kr(cell[cn.cell2]))/(cn.del1 + cn.del2)
    return krMean * cn.KMean
end

function TimeStep(t::Float64, dt::Float64, params::Params, tStops::Array{Float64, 1})
    # revise time step
    finFlag = false
    dtNew = minimum([dt*params.f, params.dtMax])
    for tStop in tStops
        if (t < tStop) && (t+dtNew > tStop)
            finFlag = (t+dtNew > tStops[end])       # indicates end-of-simulation-time has been reached
            dtNew = tStop - t
            break
        end
    end
    return dtNew, finFlag
end

function Vel(cell::Array{Cell, 1}, connect::Array{Connection, 1})
    # compute pore velocity vector arrays at cell node points
    vx = zeros(Float64, length(cell))
    vy = zeros(Float64, length(cell))
    vz = zeros(Float64, length(cell))
    xCompCell = zeros(Float64, length(cell))    # connection component projections along respective axes, per cell
    yCompCell = zeros(Float64, length(cell))
    zCompCell = zeros(Float64, length(cell))
    for cn in connect
        conduct = EffectiveK(cn, cell)        # effective hydraulic conductivity
        d = cn.del1+cn.del2
        xComp = (cell[cn.cell2].x - cell[cn.cell1].x)/d^2
        yComp = (cell[cn.cell2].y - cell[cn.cell1].y)/d^2
        zComp = (cell[cn.cell2].z - cell[cn.cell1].z)/d^2
        q = conduct * ((cell[cn.cell1].P+cell[cn.cell1].z)-(cell[cn.cell2].P+cell[cn.cell2].z)) / d
        vx[cn.cell1] += xComp * q/(cell[cn.cell1].mat.phi*cell[cn.cell1].S)
        vx[cn.cell2] += xComp * q/(cell[cn.cell2].mat.phi*cell[cn.cell2].S)
        vy[cn.cell1] += yComp * q/(cell[cn.cell1].mat.phi*cell[cn.cell1].S)
        vy[cn.cell2] += yComp * q/(cell[cn.cell2].mat.phi*cell[cn.cell2].S)
        vz[cn.cell1] += zComp * q/(cell[cn.cell1].mat.phi*cell[cn.cell1].S)
        vz[cn.cell2] += zComp * q/(cell[cn.cell2].mat.phi*cell[cn.cell2].S)
        xCompCell[cn.cell1] += abs(xComp)
        xCompCell[cn.cell2] += abs(xComp)
        yCompCell[cn.cell1] += abs(yComp)
        yCompCell[cn.cell2] += abs(yComp)
        zCompCell[cn.cell1] += abs(zComp)
        zCompCell[cn.cell2] += abs(zComp)
    end
    # normalize velocity components at each cell node by number of connections
    for i = 1:length(cell)
        vx[i] /= (xCompCell[i] + 1e-10)                 # 1e-10 factor added to prevent divide-by-zero errors
        vy[i] /= (yCompCell[i] + 1e-10)
        vz[i] /= (zCompCell[i] + 1e-10)
    end
    println("Computed pore velocity vectors.")
    return vx, vy, vz
end

### matrix operations functions

function AssembleMatrix(cell::Array{Cell, 1}, connect::Array{Connection, 1}, params::Params, dt::Float64, t::Float64, stressPeriods::Array{Float64, 1}, sourceCells::Array{Int64, 1}, source::Array{Float64, 2})
    # fill out the LHS of the equation matrix
    row_index = Int64[]                     # indexing system for sparse matrix
    col_index = Int64[]
    data = Float64[]
    b = Float64[]
    for (i, ce) in enumerate(cell)          # for each row
        diag = Storage(ce)/dt
        bSum = 0.
        for (j, nghbr) in enumerate(ce.neighbors)           # for each neighboring cell
            conduct = EffectiveK(connect[ce.connects[j]], cell) *
                connect[ce.connects[j]].A / (connect[ce.connects[j]].del1 + connect[ce.connects[j]].del2)
            push!(row_index, i)                                             # left-hand-side matrix
            push!(col_index, nghbr)
            push!(data, -params.c * conduct)
            diag += params.c * conduct
            bSum += conduct * ((cell[nghbr].P + cell[nghbr].z) - (ce.P + ce.z))         # right-hand-side (explicit) vector
        end
        push!(row_index, i)         # set diagonal term
        push!(col_index, i)
        push!(data, diag)
        push!(b, bSum)              # set explicit vector
    end
    # add source/sink terms to b[], as warranted
    j = FindStressPeriod(t, stressPeriods)
    for (i, cellIndex) in enumerate(sourceCells)
        b[cellIndex] += source[i, j]
    end
    return data, row_index, col_index, b
end

### main script

function jFlow(mode::Int64)

    # read and process input files
    params = ReadParams()                               # read model properties
    material = ReadMaterials()                          # read materials properties
    if params.mesh == true
        cell = ReadCells1(material)                     # read cell network from external grid generator
        connect = ReadGridConnects1()                   # read grid connections from externally-generated file
    else
        tempCell = ReadCells2(material)                 # read and assign cells (possibly out of sequence) from raw constraints
        cell = SortCell(tempCell)                       # use insertion sort approach to sequence cell struct array
        #cell = ReadCellMods(cell, material)            # modify specific cell properties, as specified
        connect = ReadGridConnects2()                        # read internal grid connections
        #connect = ReadSpecConnects(connect)                # append unique or boundary connections
    end
    cell = MatchConnects(cell, connect)                 # assign cell connection lists
	connect = BlendedK(connect, cell) 					# pre-calculated connection saturated K's
    stressPeriods, sourceCells, source = ReadSources()  # tabulate stress periods and source/sink cells
	drainCell, drainP = ReadDrains() 					# read cells designated as drains
    monitor = ReadMonitors()                            # note monitor locations (results tabulated for every time step)

    # create container arrays for monitoring points
    timeMon = Float64[]
    timeSeries = []
    for m in monitor
        push!(timeSeries, Float64[])
    end

    # summarize model setup
    WriteCells(cell, "CellSummary.csv", 1, drainCell, drainP)
    WriteConnections(connect)
    WriteSources(stressPeriods, sourceCells, source)

    if mode != 0            # if mode == 0, do not run model (just print out model setup and exit)

        # initialize model
        t = 0.0
        dt = params.dt0                     # first time step size
        tStops = copy(stressPeriods)        # times for which dt must be adjusted for "t" to match ...
        push!(tStops, params.endTime)
        finFlag = false
        convgFlag = false
        dP = zeros(Float64, length(cell))

        while finFlag==false

            while convgFlag == false

                # assemble matrix
                data, row_index, col_index, b = AssembleMatrix(cell, connect, params, dt, t, stressPeriods, sourceCells, source)
                A = sparse(row_index, col_index, data, length(cell), length(cell))

                # solve system of equations
                dP = \(A, b)

                # adjust cells with drains
				for (i, indx) in enumerate(drainCell)
					dP[indx] = minimum([dP[indx], drainP[i]-cell[indx].P])
				end

                # check convergence (maximum head change); adjust dt if needed
                if maximum(abs.(dP)) >= params.dhMax
                    dt /= params.f
                    @assert(dt > params.dt0)             # stop execution if dt falls below minimum
                else
                    convgFlag = true
                end

            end

            # update results
            t += dt
            for (i, ce) in enumerate(cell)
                ce.P += dP[i]
                ce.S = Sat(ce)
            end
            println("\t", "t = ", string(t))

            # write to monitor file(s)
            push!(timeMon, t)
            for (i, monCell) in enumerate(monitor)
                push!(timeSeries[i], cell[monCell].P)
            end

            # determine next time step
            dt, finFlag = TimeStep(t, dt, params, tStops)

            convgFlag = false           # reset head-change flag for next time step

        end

    end

    # write output
    WriteCells(cell, "FinalState.csv", 0, drainCell, drainP)	# write cell output file
    vx, vy, vz = Vel(cell, connect)                 			# compute and write velocity vectors to file
    WriteVectors(cell, vx, vy, vz)
    WriteTimeSeries(monitor, timeMon, timeSeries)   			# write monitor point time series to file

    println("Finished.")

end

### run script
jFlow(1)


text:
## Lifted operators

importall Base.Operators
import Base: promote_op, abs, abs2, sqrt, cbrt, scalarmin, scalarmax, isless
using Compat: @functorize

if isdefined(Base, :fieldname) && Base.fieldname(Nullable, 1) == :hasvalue # Julia 0.6
    _Nullable(R, x, hasvalue::Bool) = Nullable{R}(x, hasvalue)
else
    _Nullable(R, x, hasvalue::Bool) = Nullable{R}(x, !hasvalue)
end

# Methods adapted from Base Julia 0.5
if VERSION < v"0.5.0-dev+5096"
    promote_op(::Any, T) = T
    promote_op{T}(::Type{T}, ::Any) = T

    promote_op{R<:Number}(op, ::Type{R}) = typeof(op(one(R)))
    promote_op{R<:Number,S<:Number}(op, ::Type{R}, ::Type{S}) = typeof(op(one(R), one(S)))
end

"""
    null_safe_op(f::Any, ::Type)::Bool
    null_safe_op(f::Any, ::Type, ::Type)::Bool

Returns whether an operation `f` can safely be applied to any value of the passed type(s).
Returns `false` by default.

Custom types should implement methods for some or all operations `f` when applicable:
returning `true` means that the operation may be called on any bit pattern without
throwing an error (though returning invalid or nonsensical results is not a problem).
In particular, this means that the operation can be applied on the whole domain of the
type *and on uninitialized objects*. As a general rule, these proporties are only true for
safe operations on `isbits` types.

Types declared as safe can benefit from higher performance for operations on nullable: by
always computing the result even for null values, a branch is avoided, which helps
vectorization.
"""
null_safe_op(f::Any, ::Type) = false
null_safe_op(f::Any, ::Type, ::Type) = false

typealias SafeSignedInts Union{Int128,Int16,Int32,Int64,Int8}
typealias SafeUnsignedInts Union{Bool,UInt128,UInt16,UInt32,UInt64,UInt8}
typealias SafeInts  Union{SafeSignedInts,SafeUnsignedInts}
typealias SafeFloats Union{Float16,Float32,Float64}

# Float types appear in both since they promote to themselves,
# and therefore can't fail due to conversion of negative numbers
typealias SafeSigned Union{SafeSignedInts,SafeFloats}
typealias SafeUnsigned Union{SafeUnsignedInts,SafeFloats}
typealias SafeTypes Union{SafeInts,SafeFloats}

# Unary operators

# Note this list does not include sqrt since it can raise an error,
# nor cbrt (for which there is no functor on Julia 0.4)
for op in (:+, :-, :abs, :abs2)
    @eval begin
        null_safe_op{T<:SafeTypes}(::typeof(@functorize($op)), ::Type{T}) = true
    end
end

# No functors for these methods on 0.4: use the slow path
if VERSION >= v"0.5.0-dev"
    null_safe_op{T<:SafeInts}(::typeof(~), ::Type{T}) = true
    null_safe_op{T<:SafeTypes}(::typeof(cbrt), ::Type{T}) = true
    null_safe_op(::typeof(!), ::Type{Bool}) = true
end

for op in (:+, :-, :!, :~, :abs, :abs2, :sqrt, :cbrt)
    @eval begin
        @inline function $op{S}(x::Nullable{S})
            R = promote_op($op, S)
            if null_safe_op(@functorize($op), S)
                _Nullable(R, $op(x.value), !isnull(x))
            else
                isnull(x) ? Nullable{R}() :
                            Nullable{R}($op(x.value))
            end
        end
        $op(x::Nullable{Union{}}) = Nullable()
    end
end

# Binary operators

# Note this list does not include ^, ÷ and %
# Operations between signed and unsigned types are not safe: promotion to unsigned
# gives an InexactError for negative numbers
if VERSION >= v"0.5.0-dev"
    for op in (:+, :-, :*, :/, :&, :|, :<<, :>>, :(>>>),
               :(==), :<, :>, :<=, :>=,
               :scalarmin, :scalarmax,
               :isless)
        @eval begin
            # to fix ambiguities
            null_safe_op{S<:SafeFloats,
                         T<:SafeFloats}(::typeof($op), ::Type{S}, ::Type{T}) = true
            null_safe_op{S<:SafeSigned,
                         T<:SafeSigned}(::typeof($op), ::Type{S}, ::Type{T}) = true
            null_safe_op{S<:SafeUnsigned,
                         T<:SafeUnsigned}(::typeof($op), ::Type{S}, ::Type{T}) = true
        end
    end
else # No functors for all methods on 0.4: use the slow path for missing ones
    for op in (:+, :-, :*, :/, :&, :|,
               :<, :>,
               :scalarmin, :scalarmax)
        @eval begin
            # to fix ambiguities
            null_safe_op{S<:SafeFloats,
                         T<:SafeFloats}(::typeof(@functorize($op)), ::Type{S}, ::Type{T}) = true
            null_safe_op{S<:SafeSigned,
                         T<:SafeSigned}(::typeof(@functorize($op)), ::Type{S}, ::Type{T}) = true
            null_safe_op{S<:SafeUnsigned,
                         T<:SafeUnsigned}(::typeof(@functorize($op)), ::Type{S}, ::Type{T}) = true
        end
    end
end

for op in (:+, :-, :*, :/, :%, :÷, :&, :|, :^, :<<, :>>, :(>>>),
           :(==), :<, :>, :<=, :>=,
           :scalarmin, :scalarmax)
    @eval begin
        @inline function $op{S,T}(x::Nullable{S}, y::Nullable{T})
            R = promote_op(@functorize($op), S, T)
            if null_safe_op(@functorize($op), S, T)
                _Nullable(R, $op(x.value, y.value), !(isnull(x) | isnull(y)))
            else
                (isnull(x) | isnull(y)) ? Nullable{R}() :
                                          Nullable{R}($op(x.value, y.value))
            end
        end
        $op(x::Nullable{Union{}}, y::Nullable{Union{}}) = Nullable()
        $op{S}(x::Nullable{Union{}}, y::Nullable{S}) = Nullable{S}()
        $op{S}(x::Nullable{S}, y::Nullable{Union{}}) = Nullable{S}()
    end
end

if !method_exists(isless, Tuple{Nullable{Int}, Nullable{Int}})
    function isless{S,T}(x::Nullable{S}, y::Nullable{T})
        # NULL values are sorted last
        if null_safe_op(@functorize(isless), S, T)
            (!isnull(x) & isnull(y)) |
            (!isnull(x) & !isnull(y) & isless(x.value, y.value))
        else
            if isnull(x)
                return false
            elseif isnull(y)
                return true
            else
                return isless(x.value, y.value)
            end
        end
    end
    isless(x::Nullable{Union{}}, y::Nullable{Union{}}) = false
    isless(x::Nullable{Union{}}, y::Nullable) = false
    isless(x::Nullable, y::Nullable{Union{}}) = !isnull(x)
end


text:
mutable struct KalmanSmoothed{T}
	predicted::Array{T}
	filtered::Array{T}
	smoothed::Array{T}
	error_cov::Array{T}
	model::StateSpaceModel{T}
	y::Array{T}
  u::Array{T}
	loglik::T
end

mutable struct KalmanSmoothedMinimal{T}
	y::Array{T}
	smoothed::Array{T}
	error_cov::Union{Vector{Matrix{T}}, Vector{SparseMatrixCSC{T, Int}}}
  u::Array{T}
	model::StateSpaceModel{T}
	loglik::T
end

function show(io::IO, smoothed::KalmanSmoothed{T}) where T
    n = size(smoothed.y, 1)
    dx, dy = smoothed.model.nx, smoothed.model.ny
    println("KalmanSmoothed{$T}")
    println("$n observations, $dx-D process x $dy-D observations")
    println("Negative log-likelihood: $(smoothed.loglik)")
end

function show(io::IO, smoothed::KalmanSmoothedMinimal{T}) where T
    n = size(smoothed.y, 1)
    dx, dy = smoothed.model.nx, smoothed.model.ny
    println("KalmanSmoothedMinimal{$T}")
    println("$n observations, $dx-D process x $dy-D observations")
    println("Negative log-likelihood: $(smoothed.loglik)")
end

# Durbin-Koopman Smoothing
function kalman_smooth(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=zeros(size(y,1), model.nu)) where T

    Ksparse = issparse(model.A(1)) && issparse(model.V(1)) && issparse(model.C(1))

    y_orig = copy(y)
    y = y'
    u_orig = copy(u)
    u = u'
    n = size(y, 2)

    @assert !any(isnan, u)
    # @assert size(y, 1) == model.ny
    if size(y, 1) != model.ny
        @show size(y)
        @show model.ny
        error("")
    end
    @assert size(u, 1) == model.nu
    @assert size(u, 2) == n

    I0ny = zeros(model.ny, model.ny)

    y_notnan = (!).(isnan.(y))
    y = y .* y_notnan

    x_pred_t, P_pred_t  = model.x1, model.P1
    x, P                = zeros(model.nx, n), Array{typeof(P_pred_t)}(0)

    innov_t             = zeros(model.ny)
    innov_cov_t         = zeros(model.ny, model.ny)
    innov_cov_inv_t     = copy(innov_cov_t)
    Kt                  = Ksparse ? spzeros(model.nx, model.ny) : zeros(model.nx, model.ny)
    KCt                 = Ksparse ? spzeros(model.nx, model.nx) : zeros(model.nx, model.nx)

    innov               = zeros(model.ny, n)
    innov_cov_inv       = zeros(model.ny, model.ny, n)
    KC                  = Array{typeof(KCt)}(0)

    ut, Ct, Wt     = zeros(model.nu), model.C(1), model.W(1)

    log_likelihood = n*model.ny*log(2pi)/2
    marginal_likelihood(innov::Vector, innov_cov::Matrix,
          innov_cov_inv::Matrix) =
              (dot(innov, innov_cov_inv * innov) + logdet(innov_cov))/2

    # Prediction and filtering
    for t = 1:n

        # Predict using last iteration's values
        if t > 1

            At, But, Vt = model.A(t-1), model.B(t-1) * ut, model.V(t-1)
            AP_pred_t   = At * P_pred_t
            Kt          = AP_pred_t * Ct' * innov_cov_inv_t
            Ksparse && (Kt = sparse(Kt))
            KCt         = Kt * Ct
            push!(KC, KCt)

            x_pred_t = At * x_pred_t + But + Kt * innov_t
            P_pred_t = any(Wt .!= 0) ? AP_pred_t * (At - KCt)' + Vt : Vt

        end #if
        x[:, t]  = x_pred_t
        push!(P, (P_pred_t + P_pred_t')/2)

        # Assign new values
        ut  = u[:, t]
        Ct  = model.C(t)
        Dut = model.D(t) * ut
        Wt  = model.W(t)

        # Check for and handle missing observation values
        missing_obs = !all(y_notnan[:, t])
        if missing_obs
            ynnt = y_notnan[:, t]
            I1, I2 = spdiagm(ynnt), spdiagm(.!ynnt)
            Ct, Dut = I1 * Ct, I1 * Dut
            Wt = I1 * Wt * I1 + I2 * Wt * I2
        end #if

        # Compute nessecary filtering quantities
        innov_t           = y[:, t] - Ct * x_pred_t - Dut
        innov_cov_t       = Ct * P_pred_t * Ct' + Wt |> full
        nonzero_innov_cov = all(diag(innov_cov_t) .!= 0)
        innov_cov_inv_t   = nonzero_innov_cov ? inv(innov_cov_t) : I0ny

        innov[:, t]             = innov_t
        innov_cov_inv[:, :, t]  = innov_cov_inv_t

        nonzero_innov_cov && (log_likelihood += marginal_likelihood(innov_t, innov_cov_t, innov_cov_inv_t))

    end #for
    push!(KC, zeros(KC[1]))

    # Smoothing
    r = zeros(model.nx)
    N = zeros(model.P1)

    for t = n:-1:1

        Ct = !all(y_notnan[:, t]) ? spdiagm(y_notnan[:, t]) * model.C(t) : model.C(t)

        CF = Ksparse ? sparse(Ct' * innov_cov_inv[:, :, t]) : Ct' * innov_cov_inv[:, :, t]
        P_pred_t = P[t]
        L = model.A(t) - KC[t]

        r = CF * innov[:, t] + L' * r
        N = CF * Ct  + L' * N * L

        x[:, t] = x[:, t] + P_pred_t * r
        P_smooth_t = P_pred_t  - P_pred_t * N * P_pred_t
        P[t] = (P_smooth_t + P_smooth_t')/2

    end #for

    return KalmanSmoothedMinimal(y_orig, x', P, u_orig, model, log_likelihood)

end #smooth

function lag1_smooth(y::Array{T}, u::Array{T}, m::StateSpaceModel{T}) where T

    A_0, A_I  = zeros(m.A(1)), eye(m.A(1))
    B_0, V_0, C_0 = zeros(m.B(1)), zeros(m.V(1)), zeros(m.C(1))
    x1_0, P1_0  = zeros(m.x1), zeros(m.P1)

    A_stack(t) = [m.A(t) A_0; A_I A_0]
    B_stack(t) = [m.B(t); B_0]
    V_stack(t) = [m.V(t) V_0; V_0 V_0]
    C_stack(t) = [m.C(t) C_0]
    x1_stack   = [m.x1; x1_0]
    P1_stack   = [m.P1 P1_0; P1_0 P1_0]
    stack_model = StateSpaceModel(A_stack, V_stack, C_stack, m.W,
                                      x1_stack, P1_stack, B=B_stack, D=m.D)
    stack_smoothed = kalman_smooth(y, stack_model, u=u)

    x     = stack_smoothed.smoothed[:, 1:m.nx]'

    err_cov_type = typeof(m.P1)
    P     = err_cov_type[P[1:m.nx, 1:m.nx] for P in stack_smoothed.error_cov]
    Plag1 = err_cov_type[P[1:m.nx, (m.nx+1):end] for P in stack_smoothed.error_cov]

    return x, P, Plag1, stack_smoothed.loglik

end #function

# Rauch-Tung-Striebel Smoothing
function kalman_smooth(filt::KalmanFiltered)

    n = size(filt.y, 1)
    model = filt.model
    x_pred = filt.predicted'
    x_filt = filt.filtered'
    x_smooth = zeros(x_filt)
    P_pred = filt.pred_error_cov
    P_filt = filt.error_cov
    P_smoov = zeros(P_filt)

    x_smooth[:, n] = x_filt[:, n]
    P_smoov[:, :, n] = P_filt[:, :, n]
    for i = (n-1):-1:1
        J = P_filt[:, :, i] * model.A(i)' * inv(P_pred[:,:,i+1])
        x_smooth[:, i] = x_filt[:, i] + J * (x_smooth[:, i+1] - x_pred[:, i+1])
        P_smoov[:, :, i] = P_filt[:, :, i] + J * (P_smoov[:, :, i+1] - P_pred[:,:,i+1]) * J'
    end

    return KalmanSmoothed(x_pred', x_filt', x_smooth', P_smoov, model, filt.y, filt.u, filt.loglik)
end


text:
using Test, LazyWAVFiles, WAV, BenchmarkTools

@testset "LazyWAVFiles" begin
    @info "Testing LazyWAVFiles"

    d = mktempdir()
    a,b,c = randn(Float32,10), randn(Float32,10), randn(Float32,10)
    WAV.wavwrite(a, joinpath(d,"f1.wav"), Fs=8000)
    WAV.wavwrite(b, joinpath(d,"f2.wav"), Fs=8000)
    p1 = joinpath(d,"f1.wav")
    p2 = joinpath(d,"f2.wav")

    @testset "One file, one channel" begin
        @info "Testing One file, one channel"

        lf = LazyWAVFile(p1)
        @test lf.fs == 8000
        @test size(lf) == (10,)
        @test path(lf) == p1
        @test lf[1] == a[1]
        @test lf[1:2] == a[1:2]
        @test lf[1:10] == a
        @test lf[:] == a

        dst = zeros(Float32, 10)
        copyto!(dst, a)
        @test dst == a

        dst = zeros(Float32, 10)
        @views copyto!(dst[1:10], a[1:10])
        @test dst == a

        dst = zeros(Float32, 10)
        copyto!(dst,1, a, 1, 10)
        @test dst == a

        lf2 = LazyWAVFile(p2)
        df = [lf; lf2]
        @test df isa DistributedWAVFile
        @test size(df) == (20,)

    end

    @testset "Distributed file" begin
        @info "Testing Distributed file"

        df = DistributedWAVFile(d)
        @show df
        @test df.fs == 8000
        @test df[1] == a[1]
        @test df[1:2] == a[1:2]
        @test df[1:10] == a
        @test df[:] == [a;b]
        @test df[9:11] == [a[9:end];b[1]]
        @test df[[1,3,5]] == a[[1,3,5]]
        @test df[[1,3,5,12]] == [a[[1,3,5]];b[2]]
        @inferred df[1]
        @inferred df[1:2]
        @inferred df[1:10]
        @inferred df[:]
        @inferred df[9:11]
        @inferred df[[1,3,5]]
        @inferred df[[1,3,5,12]]

        @test size(df) == (20,)
        @test length(df) == 20
        @test length(df.files[1]) == 10
        @test length([df; df]) == 40
        @test_nowarn display(df)
        @test_nowarn display(df.files[1])
        @test ndims(df.files[1]) == 1

    end

    @testset "Multiple channels" begin
        @info "Testing Multiple channels"

        WAV.wavwrite([a b], joinpath(d,"f3.wav"), Fs=8000)
        p3 = joinpath(d,"f3.wav")
        lf = LazyWAVFile(p3)
        @test lf.fs == 8000
        @test size(lf) == (10,2)
        @test path(lf) == p3
        @test lf[1] == a[1]
        @test lf[1,2] == b[1]
        @test lf[:,1] == a
        @test lf[:,2] == b
        @test lf[1:10] == [a b]
        @test lf[:,:] == [a b]

        WAV.wavwrite([a b c], joinpath(d,"f4.wav"), Fs=8000)
        p4 = joinpath(d,"f4.wav")
        lf = LazyWAVFile(p4)
        @test lf.fs == 8000
        @test size(lf) == (10,3)
        @test path(lf) == p4
        @test @inferred(lf[1]) == a[1]
        @test @inferred(lf[1,2]) == b[1]
        @test @inferred(lf[1,3]) == c[1]
        @test @inferred(lf[:,1]) == a
        @test @inferred(lf[:,2]) == b
        @test @inferred(lf[:,3]) == c
        @test @inferred(lf[1:10]) == [a b c]
        @test @inferred(lf[:,:]) == [a b c]
        @test @inferred(lf[1,:]) == [a[1], b[1], c[1]]

        d2 = mktempdir()
        WAV.wavwrite([a b], joinpath(d2,"f5.wav"), Fs=8000)
        WAV.wavwrite([a b], joinpath(d2,"f6.wav"), Fs=8000)
        df = DistributedWAVFile(d2)
        @test df.fs == 8000
        @test size(df) == (20,2)
        @test @inferred(df[1]) == a[1]
        @test @inferred(df[1,2]) == b[1]
        @test @inferred(df[:,1]) == vcat(a, a)
        @test @inferred(df[:,2]) == vcat(b, b)
        @test @inferred(df[:,:]) == vcat([a b], [a b])
        @inferred(df[10:12,1]) == vcat(a[10:10], a[1:2])

    end

    @testset "Misc" begin
        @info "Testing Misc"

        @test_throws MethodError DistributedWAVFile(d)

        d = mktempdir()
        WAV.wavwrite(a, joinpath(d,"fs800.wav"), Fs=800)
        WAV.wavwrite(a, joinpath(d,"fs8000.wav"), Fs=8000)
        @test_throws ErrorException DistributedWAVFile(d)

        d1 = mktempdir()
        WAV.wavwrite(a, joinpath(d1,"fs800.wav"), Fs=800)
        d2 = mktempdir()
        WAV.wavwrite(b, joinpath(d2,"fs800.wav"), Fs=800)

        wavpaths = vcat(joinpath.(d1, readdir(d1)), joinpath.(d2, readdir(d2)))
        df = DistributedWAVFile(wavpaths)
        @test df[1:10] == a
        @test df[11:end] == b
    end

    @testset "Benchmarks" begin
        @info "Testing Benchmarks"

        path = mktempdir()
        y = sin.((0:99999999)/48000*2pi*440);
        wavwrite(y, joinpath(path, "test1.wav"), Fs=48000)

        dfile = DistributedWAVFile(path)

        store = zeros(10000)
        t = @belapsed copyto!($store, 1, $dfile.files[1], 2, 10000) evals=3 samples=3
        @test store == dfile[2:10001]
        @test t < 0.01

        indices1 = 1:96000
        t = @belapsed $dfile[$indices1] evals=3 samples=3
        @test t < 0.1
        t = @belapsed @view($dfile[$indices1]) evals=3 samples=3
        @test t < 0.1
        indices2 = 96000:96000*2
        t = @belapsed $dfile[$indices2] evals=3 samples=3
        @test t < 0.1
        t = @belapsed @view($dfile[$indices2]) evals=3 samples=3
        @test t < 0.1
        indices3 = 96000*10:96000*11
        t = @belapsed $dfile[$indices3] evals=3 samples=3
        @test t < 0.1
        t = @belapsed @view($dfile[$indices3]) evals=3 samples=3
        @test t < 0.1

        filepath = joinpath(path, "test1.wav")
        lfile = LazyWAVFile(filepath)

        indices1 = 1:96000
        t = @belapsed $lfile[$indices1] evals=3 samples=3
        @test t < 0.01
        t = @belapsed @view($lfile[$indices1]) evals=3 samples=3
        @test t < 0.01
        indices2 = 96000:96000*2
        t = @belapsed $lfile[$indices2] evals=3 samples=3
        @test t < 0.01
        t = @belapsed @view($lfile[$indices2]) evals=3 samples=3
        @test t < 0.01
        indices3 = 96000*10:96000*11
        t = @belapsed $lfile[$indices3] evals=3 samples=3
        @test t < 0.01
        t = @belapsed @view($lfile[$indices3]) evals=3 samples=3
        @test t < 0.01

    end

end


text:
using NavAbilitySDK
using Documenter

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

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

deploydocs(;
    repo="github.com/NavAbility/NavAbilitySDK.jl",
    devbranch="main",
)


text:
abstract type AbstractTreeWalker end

"""
A `DFSWalker` does DFS walk on a tree to compute some property of the tree.  It should implement at least
one of the two functions.

    visitfirst(walker::AbstractTreeWalker, degree_sequence, node_index, parent_index)
    visitsecond(walker::AbstractTreeWalker, degree_sequence, node_index, parent_index)

The first is called when the DFS walk first enters a node. The second is call when the walk finally
leaves the node.
"""
abstract type DFSWalker <: AbstractTreeWalker end

"""
A DFS walker that only needed to be called when the walk enters a node.
"""
abstract type FirstOnlyDFSWalker <: DFSWalker end

"""
    visitfirst(walker::AbstractTreeWalker, degree_sequence, node_index, parent_index)

The function is called when the `walker` first enters a node.
"""
function visitfirst(walker::AbstractTreeWalker, degree_sequence, node_index, parent_index)
end

"""
    visitsecond(walker::AbstractTreeWalker, degree_sequence, node_index, parent_index)

The function is called when the `walker` finally leaves a node.
"""
function visitsecond(walker::AbstractTreeWalker, degree_sequence, node_index, parent_index)
end

"""
    result(walker::AbstractTreeWalker)

Return of the result of the walk on the tree.
"""
function result(walker::AbstractTreeWalker)
    error("Please implement this function")
end

# DFS walk

function walk(degree_sequence, walker::AbstractTreeWalker, first_only=false::Bool)
    # initialize the stack
    cur_index::Int = 1
    deg, state = iterate(degree_sequence)

    stack = Stack{Tuple{Int, Int, Int}}()
    push!(stack, (deg, cur_index, 0))
    visitfirst(walker, degree_sequence, cur_index, 0)

    # while the stack is not empty
    while length(stack)> 0
        #print stack
        child_left, node_index, parent_index = pop!(stack)

        if child_left > 0
            # If we go down wards
            # put this back node in stack

            # no need to put in the stack if we do not need it anymore
            if (child_left != 1) || (!first_only)
                push!(stack, (child_left - 1, node_index, parent_index))
            end

            # put the child node in the stack
            cur_index += 1
            deg, state = iterate(degree_sequence, state)
            push!(stack, (deg, cur_index, node_index))
            visitfirst(walker, degree_sequence, cur_index, node_index)
        else
            if !first_only
                visitsecond(walker, degree_sequence, node_index, parent_index)
            end
        end
    end
    degree_sequence
end

function walk(tree::FiniteTree, walker::FirstOnlyDFSWalker)
    degree_sequence = degrees(tree)
    walk(degree_sequence, walker, true)
end

"""
    walk(tree::FiniteTree, walker::DFSWalker)

Make `walker` do a DFS walk on `tree`.
"""
function walk(tree::FiniteTree, walker::DFSWalker)
    degree_sequence = degrees(tree)
    walk(degree_sequence, walker)
end

# DepthWalker
"A walker computes the depths of all nodes in a tree."
struct DepthWalker <: FirstOnlyDFSWalker
    depth_seq::Vector{Int}
end

DepthWalker(tree::FiniteTree) = DepthWalker(zeros(size(tree)))

result(walker::DepthWalker) = walker.depth_seq

function visitfirst(walker::DepthWalker, degree_sequence, node_index, parent_index)
    if (node_index > 1)
        walker.depth_seq[node_index] = walker.depth_seq[parent_index] + 1
    end
end

# SubtreeSizeWalker

"A walker computes the size of all subtrees of a tree."
struct SubtreeSizeWalker <: DFSWalker
    subtree_size_seq::Vector{Int}
end

SubtreeSizeWalker(tree::FiniteTree) = SubtreeSizeWalker(ones(size(tree)))

result(walker::SubtreeSizeWalker) = walker.subtree_size_seq

function visitsecond(walker::SubtreeSizeWalker, degree_sequence, node_index, parent_index)
    if parent_index >= 1
        walker.subtree_size_seq[parent_index] += walker.subtree_size_seq[node_index]
    end
end

# KcutWalker

"A walker computes the (random) k-cut number of a tree."
mutable struct KcutWalker <: FirstOnlyDFSWalker
    k::Int
    record_num::Vector{Int}
    record::Float64
end
"""
    KcutWalker(k::Integer)

Construct a `KcutWalker` to compute the `k`-cut number of a tree.
"""
KcutWalker(k::Integer) = KcutWalker(k, zeros(Int, k), Inf)

result(walker::KcutWalker) = walker.record_num

function visitfirst(walker::KcutWalker, degree_sequence, node_index, parent_index)
    t::Float64 = 0.0
    for i in 1:walker.k
        t += randexp()
        if t < walker.record
            walker.record_num[i]+=1.0
        end
    end
    if t < walker.record
        walker.record = t
    end
end

# GraphWalker

"A walker that converts a tree represented in DFS degree sequence to a `FixedTreeGraph`."
mutable struct GraphWalker <: FirstOnlyDFSWalker
    tree_digraph::FixedTreeGraph
end

result(walker::GraphWalker) = walker.tree_digraph

GraphWalker(size::Int) = GraphWalker(FixedTreeGraph(size))
GraphWalker(tree::FiniteTree) = GraphWalker(size(tree))

function visitfirst(walker::GraphWalker, degree_sequence, node_index, parent_index)
    if parent_index >= 1
        addedge!(walker.tree_digraph, parent_index, node_index)
    end
end


text:
# QuantumOptics.jl functions
export Basis                  , basis                   , logarithmic_negativity  , potentialoperator       , spre,
       Bra                    , basisstate              , manybody                , printing                , states,
       CompositeBasis         , bosonstates             , manybodyoperator        , projector               , steadystate,
       DenseOperator          , coherentstate           , metrics                 , ptrace                  , stochastic,
       DenseSuperOperator     , create                  , momentum                , ptranspose              , subspace,
       FockBasis              , dagger                  , negativity              , qfunc                   , superoperators,
       GenericBasis           , destroy                 , nlevel                  , random                  , tensor,
       Ket                    , diagonaljumps           , nlevelstate             , randoperator            , timecorrelations,
       LazyProduct            , diagonaloperator        , normalize               , randstate               , timeevolution,
       LazySum                , displace                , normalize!              , samplepoints            , tracedistance,
       LazyTensor             , dm                      , number                  , semiclassical           , tracedistance_h,
       ManyBodyBasis          , eigenenergies           , occupation              , sigmam                  , tracedistance_nh,
       MomentumBasis          , eigenstates             , onebodyexpect           , sigmap                  , tracenorm,
       NLevelBasis            , embed                   , operators               , sigmax                  , tracenorm_h,
       Operator               , entropy_vn              , operators_dense         , sigmay                  , tracenorm_nh,
       PPT                    , #eval                    ,
       operators_lazyproduct  , sigmaz                  , transform,
       PositionBasis          , expect                  , operators_lazysum       , simdiag                 , transformations,
       SparseOperator         , fermionstates           , operators_lazytensor    , sortedindices           , transition,
       SparseSuperOperator    , fidelity                , operators_sparse        , sparsematrix            , variance,
       SpinBasis              , fock                    , particle                , spectralanalysis        , wigner,
       StateVector            , fockstate               , permutesystems          , spin                    , ⊗,
       SubspaceBasis          , gaussianstate           , phasespace              , spindown,
       SuperOperator          , identityoperator        , polynomials             , spinup,
       bases                  , liouvillian             , position                , spost

# basic
export outer,
       inner,
       qubit,
       qudit,
       quxit,
       sparse_spinallup_dm,
       sparse_spinalldown_dm,
       ghz,
       measure

# gates
export hadamard,
       sgate,
       sdggate,
       tgate,
       tdggate,
       phaseshift,
       singleunitary,
       controlsgate,
       cnot,
       toffoli,
       swap

# printing
export tex


text:
##############################################################################
#
# Simple Python program to benchmark several Python Excel writing modules.
#
# python bench_excel_writers.jl [num_rows] [num_cols]
#
# Original code
# Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
#

using XlsxWriter

row_max = 10000
col_max = 50

function time_xlsxwriter(const_mem=false)
	rm = Int64(floor(row_max)/2)

    wb = const_mem ? Workbook("xlsxwriter_const_mem.xlsx", Dict("constant_memory"=>true)) : Workbook("xlsxwriter.xlsx")
    ws = add_worksheet!(wb)

    for row in 1:rm
        for col in 1:col_max
            write_string!(ws, row * 2, col, "Row: $row Col: $col")
		end
        for col in 1:col_max
            write_number!(ws, row * 2 + 1, col, row + col)
		end
	end
    close(wb)
end

time_xlsxwriter(false)

println("False")
@time time_xlsxwriter(false)
println("True")
@time time_xlsxwriter(true)


text:
### Reference:
### Gareth O. Roberts and Jeffrey S. Rosenthal
### Examples of Adaptive MCMC
### Journal of Computational and Graphical Statistics, 2009, 18 (2), pp 349-367

### AMWG state

abstract type AMWGState{F<:VariateForm} <: MHSamplerState{F} end

mutable struct UnvAMWGState <: AMWGState{Univariate}
  proposal::Union{Distribution{Univariate, Continuous}, Void}
  pstate::ParameterState{Continuous, Univariate}
  tune::UnvRobertsRosenthalMCTune
  ratio::Real
  accept::Bool
  diagnosticindices::Dict{Symbol, Integer}

  function UnvAMWGState(
    proposal::Union{Distribution{Univariate, Continuous}, Void},
    pstate::ParameterState{Continuous, Univariate},
    tune::UnvRobertsRosenthalMCTune,
    ratio::Real,
    accept::Bool,
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    new(proposal, pstate, tune, ratio, accept, diagnosticindices)
  end
end

UnvAMWGState(
  proposal::Union{Distribution{Univariate, Continuous}, Void},
  pstate::ParameterState{Continuous, Univariate},
  tune::UnvRobertsRosenthalMCTune,
) =
  UnvAMWGState(proposal, pstate, tune, NaN, true, Dict{Symbol, Integer}())

mutable struct MuvAMWGState <: AMWGState{Multivariate}
  proposal::Union{Distribution{Univariate, Continuous}, Void}
  pstate::ParameterState{Continuous, Multivariate}
  tune::MuvRobertsRosenthalMCTune
  ratio::Real
  accept::Vector{Bool}
  diagnosticindices::Dict{Symbol, Integer}

  function MuvAMWGState(
    proposal::Union{Distribution{Univariate, Continuous}, Void},
    pstate::ParameterState{Continuous, Multivariate},
    tune::MuvRobertsRosenthalMCTune,
    ratio::Real,
    accept::Vector{Bool},
    diagnosticindices::Dict{Symbol, Integer}
  )
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    new(proposal, pstate, tune, ratio, accept, diagnosticindices)
  end
end

MuvAMWGState(
  proposal::Union{Distribution{Univariate, Continuous}, Void},
  pstate::ParameterState{Continuous, Multivariate},
  tune::MuvRobertsRosenthalMCTune
) =
  MuvAMWGState(proposal, pstate, tune, NaN, Array{Bool}(pstate.size), Dict{Symbol, Integer}())

### Adaptive Metropolis-within-Gibbs (AMWG) sampler

abstract type AMWG <: MHSampler end

struct UnvAMWG <: AMWG
  lower::Real
  upper::Real
  logσ0::Real
end

UnvAMWG(; lower::Real=-Inf, upper::Real=Inf, σ::Real=1.) = UnvAMWG(lower, upper, log(σ))

function initialize!(
  pstate::ParameterState{Continuous, Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::UnvAMWG,
  outopts::Dict
)
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array{Any}(length(pstate.diagnostickeys))
  end
end

tuner_state(parameter::Parameter{Continuous, Univariate}, sampler::UnvAMWG, tuner::RobertsRosenthalMCTuner) =
  UnvRobertsRosenthalMCTune(logσ=sampler.logσ0, accepted=0, proposed=0, totproposed=tuner.period)

function sampler_state(
  parameter::Parameter{Continuous, Univariate},
  sampler::UnvAMWG,
  tuner::RobertsRosenthalMCTuner,
  pstate::ParameterState{Continuous, Univariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sstate = UnvAMWGState(nothing, generate_empty(pstate), tuner_state(parameter, sampler, tuner))
  set_diagnosticindices!(sstate, [:accept, :logσ], diagnostickeys)
  sstate
end

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::UnvAMWG
)
  pstate.value = x
  parameter.logtarget!(pstate)
end

function reset!(tune::UnvRobertsRosenthalMCTune, sampler::UnvAMWG, tuner::RobertsRosenthalMCTuner)
  tune.logσ = sampler.logσ0
  tune.δ = NaN
  tune.batch = 0
  (tune.accepted, tune.proposed, tune.totproposed, tune.rate) = (0, 0, tuner.period, NaN)
end

reset!(
  sstate::UnvAMWGState,
  pstate::ParameterState{Continuous, Univariate},
  parameter::Parameter{Continuous, Univariate},
  sampler::UnvAMWG,
  tuner::RobertsRosenthalMCTuner
) =
  reset!(sstate.tune, sampler, tuner)

struct MuvAMWG <: AMWG
  lower::RealVector
  upper::RealVector
  logσ0::RealVector
end

function MuvAMWG(σ::RealVector)
  d = length(σ)
  MuvAMWG(fill(-Inf, d), fill(Inf, d), map(log, σ))
end

MuvAMWG(d::Integer; lower::RealVector=fill(-Inf, d), upper::RealVector=fill(Inf, d), σ::Real=1.) =
  MuvAMWG(lower, upper, fill(log(σ), d))

function initialize!(
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::MuvAMWG,
  outopts::Dict
)
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array{Any}(length(pstate.diagnostickeys))
  end
end

tuner_state(parameter::Parameter{Continuous, Multivariate}, sampler::MuvAMWG, tuner::RobertsRosenthalMCTuner) =
  MuvRobertsRosenthalMCTune(sampler.logσ0, accepted=fill(0, length(sampler.logσ0)), proposed=0, totproposed=tuner.period)

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::MuvAMWG,
  tuner::RobertsRosenthalMCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
)
  sstate = MuvAMWGState(nothing, generate_empty(pstate), tuner_state(parameter, sampler, tuner))
  set_diagnosticindices!(sstate, [:accept, :logσ], diagnostickeys)
  sstate
end

function initialize_diagnosticvalues!(pstate::ParameterState{Continuous, Multivariate}, sstate::MuvAMWGState)
  if haskey(sstate.diagnosticindices, :accept)
    pstate.diagnosticvalues[sstate.diagnosticindices[:accept]] = Array{Bool}(pstate.size)
  end

  if haskey(sstate.diagnosticindices, :logσ)
    pstate.diagnosticvalues[sstate.diagnosticindices[:logσ]] = Array{Real}(pstate.size)
  end
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::MuvAMWG
)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

function reset!(tune::MuvRobertsRosenthalMCTune, sampler::MuvAMWG, tuner::RobertsRosenthalMCTuner)
  d = length(sampler.logσ0)
  tune.logσ = copy(sampler.logσ0)
  tune.δ = NaN
  tune.batch = 0
  tune.accepted = fill(0, d)
  tune.proposed = 0
  tune.totproposed = 0
  tune.rate = fill(NaN, d)
end

reset!(
  sstate::MuvAMWGState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::MuvAMWG,
  tuner::RobertsRosenthalMCTuner
) =
  reset!(sstate.tune, sampler, tuner)