text:
## LINEAR LEARNING NETWORKS (FOR INTERNAL USE ONLY)

# Here `model` is a function or `Unsupervised` model (possibly
# `Static`).  This function wraps `model` in an appropriate `Node`,
# `node`, and returns `(node, mode.machine)`.
function node_and_machine(model, X)
    if model isa Unsupervised
        if model isa Static
            mach = machine(model)
        else
            mach = machine(model, X)
        end
        return transform(mach, X), mach
    else
        n = node(model, X)
        return n, n.machine
    end
end

# `models_and_functions` can include both functions (static
# operatations) and bona fide models (including `Static`
# transformers). If `ys == nothing` the learning network is assumed to
# be unsupervised. Otherwise: If `target === nothing`, then no target
# transform is applied; if `target !== nothing` and `inverse ===
# nothing`, then corresponding `transform` and `inverse_transform` are
# applied; if neither `target` nor `inverse` are `nothing`, then both
# `target` and `inverse` are assumed to be StaticTransformations, to
# be applied respectively to `ys` and to the output of the
# `models_and_functions` pipeline. Note that target inversion is
# applied to the output of the *last* nodal machine in the pipeline,
# corresponding to the last element of `models_and_functions`.

# returns a learning network machine
function linear_learning_network_machine(super_type,
                                         Xs,
                                         ys,
                                         ws,
                                         target,
                                         inverse,
                                         invert_last,
                                         operation,
                                         models_and_functions...)

    n = length(models_and_functions)

    if ys !== nothing && target !== nothing
        yt, target_machine = node_and_machine(target, ys)
    else
        yt  = ys
    end

    function tip_after_inversion(tip)
        if target !== nothing
            if inverse === nothing
                return inverse_transform(target_machine, tip)
            else
                return node_and_machine(inverse, tip) |> first
            end
        end
        return tip
    end

    if ws !== nothing
        tail_args = (yt, ws)
    else
        tail_args = (yt,)
    end

    # initialize the tip (terminal node) of network:
    tip = Xs  # last node added to network

    for m in models_and_functions
        if m isa Supervised
            supervised_machine = machine(m, tip, tail_args...)
            tip = operation(supervised_machine, tip)
            invert_last ||
                (tip = tip_after_inversion(tip))
        else
            tip = node_and_machine(m, tip) |> first
        end
    end

    invert_last && (tip = tip_after_inversion(tip))

    if super_type <: Unsupervised
        return machine(super_type(), Xs; transform=tip)
    elseif ws == nothing
        return machine(super_type(), Xs, ys; predict=tip)
    else
        return machine(super_type(), Xs, ys, ws; predict=tip)
    end

end

## PREPROCESSING

# `M` is a model type and the return value a `Symbol`. The
# `existing_names` gets updated.
function generate_name!(M::DataType, existing_names)
    str = split(string(M), '{') |> first
    candidate = split(str, '.') |> last |> snakecase |> Symbol
    candidate in existing_names ||
        (push!(existing_names, candidate); return candidate)
    n = 2
    new_candidate = candidate
    while true
        new_candidate = string(candidate, n) |> Symbol
        new_candidate in existing_names || break
        n += 1
    end
    push!(existing_names, new_candidate)
    return new_candidate
end

generate_name!(model::Model, existing_names) =
    generate_name!(typeof(model), existing_names)

pipe_alert(message) = throw(ArgumentError("@pipeline error.\n"*
                                     string(message)))
pipe_alert(k::Int) = throw(ArgumentError("@pipeline error $k. "))

not_unsupervised_alert(v) =
    pipe_alert("The value, `$v`, for `target` is not allowed. "*
               "Expecting an "*
               "`Unsupervised` model, subtype of "*
               "`Unsupervised`, or function. ")
pipe_argument_error(v) =
    pipe_alert("Encountered `$v` where a "*
               "model instance, model type, function, "*
               "or key-word assignment was expected. ")

function super_type(prediction_type::Symbol)
    if prediction_type == :deterministic
        return Deterministic
    elseif prediction_type == :probabilistic
        return Probabilistic
    elseif prediction_type == :interval
        return Interval
    else
        return Unsupervised
    end
end

function super_type(M::Type{<:Model})
    options = [Deterministic, Probabilistic,
               Interval, Unsupervised]
    idx = findfirst(super -> M <: super, options)
    return options[idx]
end

super_type(model::Model) = super_type(typeof(model))

# does expression processing, syntax and semantic
# checks. is_probabilistic is `true`, `false` or `missing`.
#function pipeline_preprocess(modl, ex, is_probabilistic::Union{Missing,Bool})
function pipeline_preprocess(modl, exs...)

    length(exs) > 0 || pipe_alert("Cannot create an empty pipeline. ")

    options = [:deterministic, :probabilistic,
               :interval]
    options_pretty =
        join([string(':',opt) for opt in options], ", ")

    fieldnames_ = Symbol[]
    models_ = Symbol[]
    model_supertypes_ = []
    models_and_functions_ = []
    models_and_functions  = []
    trait_ex_given_name_ = Dict{Symbol,Any}()

    pipetype_ = nothing
    pred_type = nothing
    target = nothing
    target_ = :nothing
    inverse = nothing
    inverse_ = :nothing
    invert_last = false
    operation = predict

    function add_model(ex)
        model_, model = eval_and_reassign(modl, ex)
        fieldname_ = generate_name!(model, fieldnames_)
        push!(models_and_functions_, model_)
        push!(models_and_functions, model)
        push!(models_, model_)
        push!(model_supertypes_, super_type(model))
    end

    for ex in exs
        if ex isa Expr || ex isa Symbol
            if ex isa Expr && ex.head == :(=)

                # --- key-word assignments ---

                variable_ = ex.args[1]
                variable_ isa Symbol || pipe_alert(8)
                if variable_ == :name
                    value_ = ex.args[2]
                    value_ isa Symbol || value_ isa AbstractString ||
                        value_ isa QuoteNode ||
                    pipe_alert("You cannot use `$value_` as the name "*
                               "for a pipeline type. ")
                    if value_ isa String
                        value_ = Symbol(value_)
                    end
                    if value_ isa QuoteNode
                        value_.value isa Symbol ||
                            pipe_alert("You cannot use `$value_` as the name "*
                                       "for a pipeline type. ")
                        value_ = value_.value
                    end
                    isdefined(modl, value_) &&
                        pipe_alert("`$(ex.args[2])` is the name of "*
                                   "an existing object. ")
                    pipetype_ = value_
                else
                    value_, value = eval_and_reassign(modl, ex.args[2])
                    if variable_ == :invert_last
                        value isa Bool ||
                            pipe_alert("`invert_last` must be "*
                                       "`true` or `false`. ") # untested
                        invert_last = value
                    elseif variable_ == :operation
                        value isa Function ||
                            pipe_alert("`operation` needs to be a function "*
                                       "but a  $(typeof(value)) encountered.")
                        value in eval.(OPERATIONS) &&
                            !(value in [transform, inverse_transform]) ||
                            pipe_alert("`operation=$value` is "*
                                       "unsupported. ") # untested
                        operation = value
                    elseif variable_ == :target
                        if value isa Function
                            value_, value =
                                eval_and_reassign(
                                    modl,
                                    :(MLJBase.WrappedFunction($value)))
                        end
                        if value isa DataType
                            value <: Unsupervised ||
                                not_unsupervised_alert(value)
                            value_, value =
                                eval_and_reassign(modl, Expr(:call, value_))
                        elseif value isa Model
                            value isa Unsupervised ||
                                not_unsupervised_alert(value)
                        else
                            not_unsupervised_alert(value)
                        end
                        target_ = value_
                        target = value
                        push!(fieldnames_, :target)
                        push!(models_, value_)
                        push!(model_supertypes_, super_type(value))
                    elseif variable_ == :inverse
                        if value isa Function
                            value_, value =
                                eval_and_reassign(
                                    modl,
                                    :(MLJBase.WrappedFunction($value)))
                        else
                            pipe_alert("`inverse` must be a function. ")
                        end
                        inverse_ = value_
                        inverse = value
                        push!(fieldnames_, :inverse)
                        push!(models_, value_)
                        push!(model_supertypes_, super_type(value))
                    elseif variable_ == :prediction_type
                        if !(value in options)
                            pipe_alert("`prediction_type` can only "*
                                       "be one of these: "*
                                       "$options_pretty. ")
                        else
                            pred_type = value
                        end
                    else
                        pipe_alert("Unrecognized key-word `$variable_`.")
                    end
                end

            else

                # --- model instances, model types and functions ---

                thing = modl.eval(ex)
                if thing isa Function
                    push!(models_and_functions_, ex)
                    push!(models_and_functions, thing)
                elseif thing isa Model
                    add_model(ex)
                elseif thing isa DataType || thing isa UnionAll
                    thing <: Model ||
                        pipe_argument_error(ex)
                    add_model(Expr(:call, ex))
                else
                    pipe_argument_error(ex)
                end
            end
        else
            pipe_argument_error(ex)
        end
    end

    pipetype_ === nothing &&
        (pipetype_ = "Pipeline"*string(gensym())[3:end] |> Symbol)

    inverse !== nothing && target === nothing &&
        pipe_alert("You have specified `inverse=...` but no `target`.")

    supervised_components = filter(models_and_functions) do m
        !(m isa Unsupervised) && !(m isa Function)
    end

    length(supervised_components) < 2 ||
        pipe_alert("More than one component of the pipeline is a "*
                   "supervised model .")

    is_supervised  = (length(supervised_components) == 1)

    non_static_components = filter(models_and_functions) do m
        m isa Unsupervised &&
            !(m isa Static) ||
            m isa Supervised
    end

    is_static = isempty(non_static_components)

    # assign traits where this can be done reliably and does not
    # overwrite user-defined ones:
    if is_supervised
        supervised_model = supervised_components[1]
        supervised_is_last =
            supervised_model == last(models_and_functions) &&
            target == nothing
        if !haskey(trait_ex_given_name_, :supports_weights)
            trait_ex_given_name_[:supports_weights] =
                supports_weights(supervised_model)
        end
        if !haskey(trait_ex_given_name_, :target_scitype) &&
            supervised_is_last
            trait_ex_given_name_[:target_scitype] =
                :(MLJBase.target_scitype($supervised_model))
        end
    elseif !haskey(trait_ex_given_name_, :output_scitype)
        last_model = last(models_and_functions)
        if last_model isa Model
            trait_ex_given_name_[:output_scitype] =
                output_scitype(last_model)
        end
    end

    # check consistency of user-defined prediction_type and infer
    # final super type:
    if is_supervised
        if pred_type != nothing
            super = super_type(pred_type)
            supervised_is_last && !(supervised_model isa super) &&
                pipe_alert("The pipeline's last component model has type "*
                           "`$(typeof(supervised_model))`, which conflicts "*
                           "the declaration "*
                           "`prediction_type=$pred_type`. ")
        elseif supervised_is_last
            if operation != predict
                super = Deterministic
            else
                super = super_type(supervised_model)
            end
        else
            @info "Treating pipeline as a `Deterministic` predictor.\n"*
            "To override, specify `prediction_type=...` "*
            "(options: $options_pretty). "
            super = Deterministic
        end
    else
        if pred_type != nothing
            @warn "Pipeline appears to have no supervised "*
            "component models. Ignoring declaration "*
            "`prediction_type=$(pred_type)`. "
        end
        super = is_static ? Static : Unsupervised
    end

    target isa MLJBase.WrappedFunction && inverse == nothing &&
        pipe_alert("It appears `target` is a function. "*
                   "You must therefore specify `inverse=...` .")

    target == nothing || is_supervised ||
        pipe_alert("`target=...` has been specified but no "*
                   "supervised components have been specified. ")

    # check if attempting a target inversion immediately after a
    # probablisitic prediction:
    if target != nothing  &&
        supervised_model isa Probabilistic &&
        operation == predict && (supervised_is_last || !invert_last)
        @warn "Pipeline is applying a target inversion "*
        "immediately after a probabilistic target prediction. "*
        "Perhaps you want to set "*
        "`operation=predict_mean` or `operation=predict_mode`. "*
        "Alternatively, you may want to "*
        "insert a function/model after the "*
        "supervised model and set `invert_last=true`. "
    end

    return (pipetype_, fieldnames_, models_, model_supertypes_,
            models_and_functions_, target_, inverse_,
            invert_last, operation,
            super, trait_ex_given_name_)

end

function pipeline_(modl, exs...)

    (pipetype_, fieldnames_, models_, model_supertypes_,
     models_and_functions_,
     target_, inverse_, invert_last, operation,
     super_type, trait_ex_given_name_) =
         pipeline_preprocess(modl, exs...)

    if prediction_type === :unknown # unsupervised case
        ys_ = :nothing
    else
        ys_ = :(source())
    end

    if haskey(trait_ex_given_name_, :supports_weights) &&
        trait_ex_given_name_[:supports_weights]
        ws_ = :(source())
    else
        ws_ = nothing
    end

    mach_ = quote
        $(isdefined(modl, :MLJ) ? :(import MLJ.MLJBase) : :(import MLJBase))
        MLJBase.linear_learning_network_machine($super_type,
                                                source(nothing), $ys_, $ws_,
                                                $target_, $inverse_,
                                                $invert_last,
                                                $operation,
                                                $(models_and_functions_...))
    end

    mach_, mach = eval_and_reassign(modl, mach_)

    field_declarations =
        [:($(fieldnames_[j])::$(model_supertypes_[j])=$(models_[j]))
                                for j in eachindex(models_)]

    supertype_ = _exported_type(mach.model)

    struct_ = :(mutable struct $pipetype_ <: $supertype_
                        $(field_declarations...)
                    end)

    no_fields = length(models_) == 0

    from_network_(modl,
                  mach_,
                  pipetype_,
                  struct_,
                  no_fields,
                  trait_ex_given_name_)

    return pipetype_

end

"""
    @pipeline model1 model2 ... modelk

Create an instance of an automatically generated composite model type,
in which the specified models are composed in order. This means
`model1` receives inputs, whose output is passed to `model2`, and so
forth. Model types or instances may be specified.

**Important.** By default a new model *type* name is automatically
generated. To specify a different name add a keyword argument such as
`name=MyPipeType`. This is necessary if serializing the pipeline; see
[`MLJ.save`](@ref).

At most one of the models may be a supervised model, but this model
can appear in any position.

The `@pipeline` macro accepts several key-word arguments discussed
further below.

Static (unlearned) transformations - that is, ordinary functions - may
also be inserted in the pipeline as shown in the following example:

    @pipeline X->coerce(X, :age=>Continuous) OneHotEncoder ConstantClassifier

### Target transformation and inverse transformation

A learned target transformation (such as standardization) can also be
specified, using the key-word `target`, provided the transformer
provides an `inverse_transform` method:

    @pipeline OneHotEncoder KNNRegressor target=UnivariateTransformer

A static transformation can be specified instead, but then
an `inverse` must also be given:

    @pipeline(OneHotEncoder, KNNRegressor,
              target = v -> log.(v),
              inverse = v -> exp.(v))

*Important.* By default, the target inversion is applied *immediately
 following* the (unique) supervised model in the pipeline. To apply
 at the end of the pipeline, specify `invert_last=true`.

### Optional key-word arguments

- `target=...`  -  any `Unsupervised` model or `Function`

- `inverse=...`  -  any `Function` (unspecified if `target` is `Unsupervised`)

- `invert_last`  -  set to `true` to delay target
  inversion to end of pipeline (default=`true`)

- `prediction_type`  -
  prediction type of the pipeline; possible values: `:deterministic`,
  `:probabilistic`, `:interval` (default=`:deterministic` if not inferable)

- `operation` - operation applied to the supervised component model,
  when present; possible values: `predict`, `predict_mean`,
  `predict_median`, `predict_mode` (default=`predict`)

- `name` - new composite model type name; can be any name not already
  in current global namespace (autogenerated by default(

See also: [`@from_network`](@ref)

"""
macro pipeline(exs...)

    pipetype_ = pipeline_(__module__, exs...)

    esc(quote
        $pipetype_()
        end)

end

ext: jl


text:
# load pkgs
using Ipopt, HarmonicPowerModels, PowerModels

# pkg const
const _PMs = PowerModels
const _HPM = HarmonicPowerModels

# path to the data
path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_Yd11_simplified.m")
# path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_Yd11_simplified_noload.m")

# transformer excitation data
xfmr = Dict("voltage_harmonics" => [1,3],
            "current_harmonics" => [1,3], #this one has to match the data in the matpower file
            "N" => 50,
            "current_type" => :rectangular,
            "excitation_type" => :sigmoid,
            "inom" => 0.13,
            "ψmax" => 1,
            "voltage_type" => :rectangular,
            "dv" => [0.1,0.1],
            "vmin" => -[1.1,1.1],
            "vmax" => [1.1,1.1],
            "dθ" => [π/5,π/5],
            "θmin" => [0.0,0.0],
            "θmax" => [2π,2π])

# load data
data  = _PMs.parse_file(path)

vmmin = 0.8
for (b,branch) in data["branch"]
    branch["c_rating"] = branch["rate_a"]/vmmin
end
for (d,load) in data["load"]
    load["c_rating"] = abs(load["pd"] + im* load["qd"])/vmmin
end
for (g,gen) in data["gen"]
    gen["c_rating"] = abs(gen["pmax"] + im* gen["qmax"])/vmmin
end

# hdata = _HPM.replicate(data)
hdata = _HPM.replicate(data, xfmr_exc=xfmr)

# set the solver
solver = Ipopt.Optimizer

# solve the hopf
# run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
pm = _PMs.instantiate_model(hdata, _PMs.IVRPowerModel, _HPM.build_hopf_iv; ref_extensions=[_HPM.ref_add_xfmr!]);
result = optimize_model!(pm, optimizer=solver, solution_processors=[ _HPM.sol_data_model!])
_HPM.append_indicators!(result, hdata)

##
# print(pm.model)

println("Harmonic 3")
_PMs.print_summary(result["solution"]["nw"]["2"])
println("Harmonic 1")
_PMs.print_summary(result["solution"]["nw"]["1"])
result["objective"]
result["termination_status"]

@show result["solution"]["nw"]["1"]["xfmr"]["1"]["pexc"]
##
@show target = 77/1000/100


text:
n = NamedArray(randn(5, 10))

@test array(n) == n.array
@test convert(Array, n) == n.array
@test convert(NamedArray, n) == n

for T in [Float32, Float16]
	for m in (map(T, n), convert(NamedArray{T}, n))
		@test m.array == map(T, n.array)
		@test allnames(m) == allnames(n)
		@test dimnames(m) == dimnames(n)
	end
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CovidSurvey
using CSV
using DataFrames
using Dates
# using RCall
using Turing
using StatsPlots
using Plots.PlotMeasures
using Random

## ===========================================================================
# load old model without predictors and use posterior samples as new priors
# ============================================================================
fdir = "/home/and/.julia/dev/CovidSurvey/out/"
fname = "chains_cases=2020-06-01_chains=1_hospit=true_model=v1_observ=2021-03-25_seed=10_sero=true_steps=800_warmup=300.jls"

# ps = CovidSurvey.parseparams(fdir*fname)
# data = National.load_data(ps.observ, nothing, ps.cases, ps.hospit, ps.sero)
# turing_data = data.turing_data;
# model = name2model[ps.model]
chain = chainscat(
    read(fdir*fname, Chains),
    read("/home/and/.julia/dev/CovidSurvey/out/chains_cases=2020-06-01_chains=1_hospit=true_model=v1_observ=2021-03-25_seed=11_sero=true_steps=800_warmup=300.jls", Chains),
    read("/home/and/.julia/dev/CovidSurvey/out/chains_cases=2020-06-01_chains=1_hospit=true_model=v1_observ=2021-03-25_seed=12_sero=true_steps=800_warmup=300.jls", Chains),
)
chain = chain[201:end,:,:]
#------------------------------------------------------------------------------
# use posterior samples as priors
elementary_priors = (
    :y => Gamma{Float64},
    :ϕ => Gamma{Float64},
    :ϕ2 => Normal{Float64},
    :ϕ3 => Normal{Float64},
    :R0 => Normal{Float64},
    :R1 => Normal{Float64},
    :ifr => Normal{Float64},
    :ihr => Normal{Float64},
    :iar0 => Gamma{Float64},
    :σ_rt => Gamma{Float64},
    :σ_iar => Gamma{Float64},
);
array_priors = (
    :logRt => Normal{Float64},
    :seroprev_σ => Normal{Float64},
    :logit_iar => Normal{Float64},
);

function fitsamples(chain)
    dfits = Dict()
    for (n, d) in elementary_priors
        dfits[n] = fit(d, vec(chain[n]) )
    end
    for (n, d) in array_priors
        c = group(chain, n)
        ns = names(c)
        @assert issorted( [parse(Int, match(r"\d+" , s).match) for s in String.(ns)] )
        dfits[n] = [fit(d, c[n] |> vec) for n in ns]
    end
    dfits[:logRt] = [Normal{Float64}(mean(d), std(d)*2) for d in dfits[:logRt]]
    dfits[:logRt] = dfits[:logRt][1:272] #DEBUG
    dfits[:logit_iar] = dfits[:logit_iar][1:29] #DEBUG
    (; dfits...)
end

priors = fitsamples(chain);

# histogram(chain["logRt[25]"]|>vec, normalize=true)
# plot!(priors[:logRt][1])

## ===========================================================================
# setup new model
# ============================================================================

name2model = Dict(
    "v1" => National.model_v1,
    "gp" => National.model_gp,
    "contacts" => National.model_contacts,
    "fitted" => ,
)

ps = (
    warmup = 200,
    steps  = 600,
    seed   = 20000,
    observ = "2020-12-16",#"2021-03-25"
    cases  = "2020-06-01",
    model  = "fitted",
    preds = "CC,CS,CF,CR",
    hospit = true,
    sero = true,
)

fname = savename("chains", ps, "")
Random.seed!(ps.seed);
data = National.load_data(ps.observ, ps.preds |> parse_predictors, ps.cases, ps.hospit, ps.sero; update=false)
turing_data = data.turing_data;
#----------------------------------------------------------------------------
# sample model
m = National.model_contacts_fitted(turing_data..., priors, false)
@time chain = sample(m, NUTS(ps.warmup, 0.95), MCMCThreads(), ps.steps + ps.warmup, 4)
chain = chain[ps.warmup+1:end,:,:]
@info "Saving at: $(projectdir("out", fname))"
safesave(projectdir("out/tmp", fname*".jls"), chain)
#-----------------------------------------------------------------------------
# diagnostics
diagnostics = gelmandiag(chain)
fname_diagnostics = projectdir("out/tmp", fname*"_GELMANDIAG.csv")
safesave(fname_diagnostics, diagnostics)
run(`firefox $(fname_diagnostics)`, wait=false)
#-----------------------------------------------------------------------------
# make predictions
m_pred = model(turing_data..., priors, true)
gq = Turing.generated_quantities(m_pred, chain)
generated_posterior = vectup2tupvec( reshape(gq, length(gq)) );
observed_daily_cases, expected_daily_deaths, expected_daily_hospit, Rt, cumulative_cases, iar = generated_posterior;
#---------------------------------------------------------------------------
# plot posterior predictions
plotlyjs()
p = National.plot_results(data, generated_posterior...)
savefig(p, "/home/and/tmp/figures/"*fname*".html")
# savefig(p, projectdir("figures/tmp", fname*".png") )
run(`firefox $("/home/and/tmp/figures/"*fname*".html")`, wait=false)
#---------------------------------------------------------------------------
# plot chain statistics
plotlyjs()
n = filter( x->!occursin(r"\[", x), String.(names(chain)))
p = plot(chain[n]);
savefig(p, "/home/and/tmp/figures/"*fname*"_CHAINS.html" )
run(`firefox  $("/home/and/tmp/figures/"*fname*"_CHAINS.html")`, wait=false)
#---------------------------------------------------------------------------
# plot more chain statistics (meanplot)
using MCMCChains
plotlyjs()
p = meanplot(chain[n]);
savefig(p, projectdir("/home/and/tmp/figures/", fname*"_MEANPLOT.html") )
run(`firefox $(projectdir("/home/and/tmp/figures/", fname*"_MEANPLOT.html"))`, wait=false)
## ===========================================================================
# evaluation
# ============================================================================
pgfplotsx()
using LaTeXStrings
## posterior predictions
# scalefontsizes(2.)
p = let
    ocases  = hcat(observed_daily_cases...)#hcat([observed_daily_cases[t] for t = 1:data.num_total]...)
    ehospit  = hcat(expected_daily_hospit...)

    tickdates = data.dates[1]:Month(1):data.dates[end]
    p3 = plot(; legend = :outertopright, ylims=(0,7000), bottom_margin=3mm)
    National.plot_cases_posterior!(p3, data, ocases)
    p4 = plot(; legend = :outertopright, ylims=(0,500), top_margin=3mm)
    National.plot_hospit_posterior!(p4, data, ehospit)
    plot(p3, p4, layout=(2,1), sharex=true, link=:x, size=(1000, 700)) #, size=(1500, 1000)
end
# savefig(p, projectdir("figures/", "posterior_check.png") )
#---------------------------------------------------------------------------
# plot posterior effect size distribution vs. prior
pgfplotsx()
p = let
    xs = chain["α[1]"]|>vec
    p1 = histogram(xs, normalize=true, title="contacts to colleagues", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.015), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = chain["α[2]"]|>vec
    p2 = histogram(xs, normalize=true, title="contacts to strangers", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.015), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = chain["α[3]"]|>vec
    p3 = histogram(xs, normalize=true, title="contacts to family", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.015), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = chain["α[4]"]|>vec
    p4 = histogram(xs, normalize=true, title="contacts to friends", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.015), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    plot(p1,p2,p3,p4, layout=(2,2))
end
# savefig(p, projectdir("figures/", "total_effect_analysis.png") )
#-----------------------------------------------------------------------------
# evaluate Rt
ib = findfirst(==(Date("2020-05-15")), data.dates)
ie = findfirst(==(Date(ps.observ)), data.dates)
Rt_array = hcat(Rt...)[ib:ie,:]

qs = [quantile(v, [0.025, 0.25, 0.5, 0.75, 0.975]) for v in eachrow(Rt_array)]
llq, lq, mq, uq, uuq = (eachrow(hcat(qs...))..., )

t = data.dates[ib:ie]
plot(t, mq, ribbon=(mq - llq, uuq - mq), linewidth=0, label="95% quantiles")
plot!(t, mq, ribbon=(mq - lq, uq - mq), linewidth=2, label="50% quantiles", hover=mq)
# save(projectdir("out", "Rt_corrected_with_contacts.csv"), DataFrame((;llq, lq, mq, uq, uuq)))
#------------------------------------------------------------------------------
# Rt_with_vs_without_predictors
pgfplotsx()
Rt_original = load(projectdir("out", "Rt.csv")) |> DataFrame

p = let
    p = plot(legend=:topright)
    plot!(t, mq, ribbon=(mq - llq, uuq - mq), linewidth=2, label="with predictors (contact data)")
    plot!(Rt_original.date, Rt_original.mq, ribbon=(Rt_original.mq - Rt_original.llq, Rt_original.uuq - Rt_original.mq), linewidth=2, label="without predictors")
    xlims!((t[1],t[end]))
end
#------------------------------------------------------------------------------
# Rt_vs_predictors
pgfplotsx()
cov = turing_data.covariates
p = let
    p = plot(legend=:topright)
    plot!(Rt_original.date, Rt_original.mq, ribbon=(Rt_original.mq - Rt_original.llq, Rt_original.uuq - Rt_original.mq), linewidth=2, label="Rt without predictors", c=:black)
    p2 = twinx()
    for (i, n) in enumerate(["colleagues","strangers","family","friends"])
        plot!(p2, t, cov[1:length(t),i], lw=2, label=n, grid=:off, xticks=:none)
    end
    xlims!((t[1],t[end]))
end
# savefig(p, projectdir("figures/", "Rt_vs_predictors.png") )
#------------------------------------------------------------------------------
# compare to SSI data to predictors
Rt_SSI = National.readsurvey("/home/and/data/covidsurvey/contacts.csv")
Rt_SSI = select( Rt_SSI, [:Rt,:uncertainty_lower,:uncertainty_upper, :date])
p = let
    p = plot(legend=:topright)
    plot!(Rt_SSI.date, Rt_SSI.Rt, ribbon=(Rt_SSI.Rt - Rt_SSI.uncertainty_lower, Rt_SSI.uncertainty_upper - Rt_SSI.Rt), linewidth=2, label="Rt (SSI) without predictors", c=:black)
    p2 = twinx()
    for (i, n) in enumerate(["colleagues","strangers","family","friends"])
        plot!(p2, t, cov[1:length(t),i], lw=2, label=n, grid=:off, xticks=:none)
    end
    xlims!((t[1],t[end]))
end
# savefig(p, projectdir("figures/", "Rt_from_SSI_vs_predictors.png") )
#------------------------------------------------------------------------------
# compare to SSI data to inferred Rt
p = let
    p = plot(legend=:topright)
    plot!(Rt_original.date, Rt_original.mq, ribbon=(Rt_original.mq - Rt_original.llq, Rt_original.uuq - Rt_original.mq), linewidth=2, label="Rt without predictors")
    plot!(Rt_SSI.date, Rt_SSI.Rt, ribbon=(Rt_SSI.Rt - Rt_SSI.uncertainty_lower, Rt_SSI.uncertainty_upper - Rt_SSI.Rt), linewidth=2, label="Rt (SSI)", c=:black)
    xlims!((t[1],t[end]))
end
# savefig(p, projectdir("figures/", "Rt_without_predictors_vs_Rt_from_SSI.png") )
#------------------------------------------------------------------------------
# compare to Rt from SSI or model to predictors
survey = National.readsurveys_raw()
p = let
    ps = []
    ts = survey.date
    for (i, n) in [(3,:family), (1,:colleagues), (4,:friends), (2,:strangers)] #"CC,CS,CF,CR"
        p = plot()
        plot!(
         #Rt_SSI.date, Rt_SSI.Rt, ribbon=(Rt_SSI.Rt - Rt_SSI.uncertainty_lower, Rt_SSI.uncertainty_upper - Rt_SSI.Rt),
         Rt_original.date, Rt_original.mq, ribbon=(Rt_original.mq - Rt_original.llq, Rt_original.uuq - Rt_original.mq),
         linewidth=2, label=nothing, c=:black, α=.5, ylabel="\$R_t\$",
         right_margin=20mm, bottom_margin=10mm)
        p2 = twinx()
        xs = survey[:, n]
        x = cov[1:length(t),i]
        National.surveynorm!(xs)
        plot!(p2, ts, xs / std(xs) * std(x), lw=2,
            label="$n (raw)", grid=:off, xticks=:none, c=:darkred,
            ylabel="contacts (normalized)",
            right_margin=20mm, bottom_margin=10mm,
            extra_kwargs=:subplot, legend=:outertop, legend_columns=-1)
        plot!(p2, t, x, lw=4, label="$n (smoothed)",
            grid=:off, xticks=:none, c=:darkorange,
            right_margin=20mm, bottom_margin=10mm,
            extra_kwargs=:subplot, legend=:outertop, legend_columns=-1)
        xlims!((t[1],t[end]))
        push!(ps, p)
    end
    #l = @layout [a{0.2w} b{0.2w}; c d]
    plot(ps..., layout=4, size=(1300, 1000))
end
# savefig(p, projectdir("figures/", "Rt_vs_predictors_smoothed_and_raw.png") )

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
using Base.Test
using Lens

clear_all_listeners!()
@test length(values(Lens.lens_to_listeners)) == 0

function simpletest1()
  a = 10
  b = 10
  lens(:after_rands1;a=a,b=b)
end

gfunc(m) = @test sum(m[:a]) + sum(m[:b]) == 20
register!(:after_rands1, Listener(:gfunc, gfunc, true, true))
@test nlisteners() == 1
enable_all_listeners!()
simpletest1()

clear_all_listeners!()
function simpletest2()
  a = 10
  b = 10
  lens(:after_rands2,a,b)
end

register!(:gunit, :after_rands2, false) do x,y
  @test x + y == 20
end

simpletest2()
enable_all_listeners!()
disable_all_listeners!()

function f()
  x = 20
  y = 30
  z = 50
  lens(:x,x)
  lens(:y,y)
  lens(:z,z)
end

value, res = capture(f,[:x,:y,:z])


text:
module adapters_postgresql_execute_result_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo Schema Raw

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

Repo.debug_sql()

include("options.jl")

Repo.connect(;
    adapter = Octo.Adapters.PostgreSQL,
    Options.for_postgresql...
)

result = Repo.execute([DROP TABLE IF EXISTS Employee])
@test result === nothing

result = Repo.execute(Raw("""
    CREATE TABLE Employee (
        ID SERIAL,
        Name VARCHAR(255),
        Salary FLOAT(8),
        PRIMARY KEY (ID)
    )"""))
@test result === nothing

changes = (Name="John", Salary=10000.50)
@test Repo.insert!(Employee, changes).id == 1

changes = (Salary=3000, Name="Mike")
@test Repo.insert!(Employee, changes).id == 2

struct Test1
end
Schema.model(Test1, table_name="test1", primary_key=nothing)
result = Repo.execute([DROP TABLE IF EXISTS Test1])
@test result === nothing

result = Repo.execute(Raw("""
CREATE TABLE IF NOT EXISTS test1 (a boolean, b text)
"""))
@test result === nothing

result = Repo.insert!(Test1, (a=true, b="sic est"))
@test result == (num_affected_rows = 1,)
result = Repo.insert!(Test1, (a=true, b="sic est"); returning=nothing)
@test result == (num_affected_rows = 1,)

struct Test2
end
Schema.model(Test2, table_name="test2", primary_key=("a", "b"))
Repo.execute([DROP TABLE IF EXISTS Test2])
Repo.execute(Raw("""
CREATE TABLE IF NOT EXISTS test2 (a boolean, b text, PRIMARY KEY (a, b))
"""))

result = Repo.insert!(Test2, (a=true, b="sic est1"))
@test result == (a=true, b="sic est1", num_affected_rows=1)
result = Repo.insert!(Test2, (a=true, b="sic est2"); returning=nothing)
@test result == (num_affected_rows = 1,)
result = Repo.insert!(Test2, (a=true, b="sic est3"); returning=[:a])
@test result == (a=true, num_affected_rows=1)
result = Repo.insert!(Test2, (a=true, b="sic est4"); returning=[:a, :b])
@test result == (a=true, b="sic est4", num_affected_rows=1)
result = Repo.insert!(Test2, (a=true, b="sic est5"); returning=[:b, :a])
@test result == (b="sic est5", a=true, num_affected_rows=1)

Repo.disconnect()

end # module adapters_postgresql_execute_result_test


text:
using Images
using Printf
using FileIO
using STC
using STC.SUniward
using ArgParse
using Random
using DataFrames

function parsearguments(args)
    s = ArgParseSettings("experiment_1.jl")

    @add_arg_table s begin
        "-p"
            default = 0.3
            arg_type=Float64
            help = "payload to embed (default 0.3bpp"
        "--input"
            default = "../test/1.pgm"
            arg_type=String
            help = "path to the input image (in pgm format)"
        "--output"
            default = "/tmp/1.pgm"
            arg_type=String
            help = "path to the output image (in pgm format)"
        "--height"
            default = 3
            arg_type=Int
            help = "height of the embedding matrix"
        "--method"
            default = "dynamic"
            help = "a methods to use during embedding"
    end

    parse_args(args, s,as_symbols=true)
end

function saveandshow(cover,stego,d1,ofname)
    δ = Int.(cover) - Int.(stego.image)
    d2 = SUniward.suniwardcosts(cover,stego.image);
    cp = sum(δ .> 0);
    cm = sum(δ .< 0);
    println(@sprintf("stc cost = %g true cost = %g pixel increased / decreased %d / %d",d1,d2,cp,cm));
    !isdir(dirname(ofname)) && mkpath(dirname(ofname))
    save(ofname,Gray.(stego.image/255))
    open(ofname[1:end-3]*"meta","w") do fid
        println(fid,@sprintf("%g %g %d %d",d1,d2,cp,cm))
    end
    true
end

settings=parsearguments(ARGS)
h = settings[:height]
cover=rawview(channelview(load(settings[:input])));
numofcols=Int(round(1/(0.66*settings[:p])));
hhat=rand(1:3^h-1,numofcols);
numofblocks=Int(floor(length(cover)/numofcols));
message=rand(0:2,numofblocks);
embpath = randperm(length(cover));

if settings[:method] == "dynamic"
    stegos = SUniward.IncrementalSUniward(cover,3^h;sigma=1.0,T=Float16);
    (stego,d1)=STC.variablecoding(cover,stegos,message,hhat,h,embpath,3);
    saveandshow(cover,stego,d1,settings[:output])
elseif settings[:method] == "additive"
    ρ = SUniward.suniwardcosts(Float64.(cover));
    stegos = [SUniward.SUniwardAdd(cover, ρ) for i in 1:3^h];
    (stego, d1) = STC.variablecoding(cover,stegos,message,hhat,h,embpath,3);
    saveandshow(cover,stego,d1,settings[:output])
else
    error("unknown embedding method ",settings[:method])
end


text:
using Test

@testset "Sampling" begin
  include("Sampling.jl")
end

@testset "SchemaBase" begin
  include("SchemaBase.jl")
end

@testset "SchemaEvents" begin
  include("SchemaEvents.jl")
end