text:
using DrWatson
using FileIO, BSON, DataFrames
using PrettyTables
using Statistics
using StatsBase
import PGFPlots
include("./utils/pgf_boxplot.jl")
include("./utils/ranks.jl")

using GenerativeAD.Evaluation: MODEL_MERGE, MODEL_ALIAS, DATASET_ALIAS, MODEL_TYPE, apply_aliases!
using GenerativeAD.Evaluation: _prefix_symbol, aggregate_stats_mean_max, aggregate_stats_max_mean
using GenerativeAD.Evaluation: PAT_METRICS, PATN_METRICS, PAC_METRICS, BASE_METRICS, TRAIN_EVAL_TIMES
using GenerativeAD.Evaluation: rank_table, print_rank_table, latex_booktabs, convert_anomaly_class

const PAT_METRICS_NAMES = ["\$PR@\\%0.01\$","\$PR@\\%0.1\$","\$PR@\\%1\$","\$PR@\\%5\$","\$PR@\\%10\$","\$PR@\\%20\$"]
const PAC_METRICS_NAMES = ["\$AUC@\\#5\$","\$AUC@\\#10\$","\$AUC@\\#50\$","\$AUC@\\#100\$","\$AUC@\\#500\$","\$AUC@\\#1000\$"]
const PATN_METRICS_NAMES = ["\$PR@\\#5\$","\$PR@\\#10\$","\$PR@\\#50\$","\$PR@\\#100\$","\$PR@\\#500\$","\$PR@\\#1000\$"]

AE_MERGE = Dict(
    "aae_full" => "aae",
    "aae_vamp" => "aae",
    "wae_full" => "wae",
    "wae_vamp" => "wae",
    "vae_full" => "vae")

DOWNSAMPLE = Dict(
    zip(["aae", "vae", "wae", "MAF", "RealNVP"], [100, 100, 100, 100, 100]))

function _merge_ae_filter!(df)
    # filter out simple vae as its properties are shown in the separate table focused on ae
    filter!(x -> (x.modelname != "vae_simple"), df)

    # filter out ocsvm_nu
    filter!(x -> (x.modelname != "ocsvm_nu"), df)

    # merge autoencoders
    apply_aliases!(df, col="modelname", d=AE_MERGE)

    # filter only reconstruction_samples score
    filter!(x -> ~(x.modelname in (["aae", "vae", "wae"])) || occursin("_score=reconstruction-sampled", x.parameters), df)

    # filter only default parameters of "ocsvm" - rbf + nu = 0.5, tune gamma
    filter!(x -> (x.modelname != "ocsvm_rbf") || occursin("nu=0.5", x.parameters), df)
    df
end

function _filter_ensembles!(df)
    filter!(x -> occursin("ignore_nan=true_method=mean", x.parameters), df)
end

df_tabular = load(datadir("evaluation/tabular_eval.bson"))[:df];
df_tabular_ens = load(datadir("evaluation_ensembles/tabular_eval.bson"))[:df];
@info "Loaded results from tabular evaluation."

function basic_summary_table(df; suffix="", prefix="", downsample=Dict{String, Int}())
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)

    for metric in [:auc, :tpr_5]
        agg_names = ["maxmean", "meanmax"]
        agg_funct = [aggregate_stats_max_mean, aggregate_stats_mean_max]
        for (name, agg) in zip(agg_names, agg_funct)
            val_metric = _prefix_symbol("val", metric)
            tst_metric = _prefix_symbol("tst", metric)

            df_agg = agg(df, val_metric, downsample=downsample)

            df_agg["model_type"] = copy(df_agg["modelname"])
            apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
            apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
            apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

            sort!(df_agg, (:dataset, :model_type, :modelname))
            rt = rank_table(df_agg, tst_metric)
            rt[end-2, 1] = "\$\\sigma_1\$"
            rt[end-1, 1] = "\$\\sigma_{10}\$"
            rt[end, 1] = "rnk"

            filename = "$(projectdir())/paper/tables/$(prefix)_$(metric)_$(metric)_$(name)$(suffix).tex"
            open(filename, "w") do io
                print_rank_table(io, rt; backend=:tex)
            end
        end
    end
end

function basic_summary_table_tabular(df; suffix="", downsample=Dict{String, Int}())
    basic_summary_table(df; prefix="tabular", suffix=suffix, downsample=downsample)
end

# basic_summary_table_tabular(copy(df_tabular))
# basic_summary_table_tabular(_merge_ae_filter!(copy(df_tabular)), suffix="_merge_ae")
basic_summary_table_tabular(_merge_ae_filter!(copy(df_tabular)), suffix="_merge_ae_down", downsample=DOWNSAMPLE)
# basic_summary_table_tabular(copy(df_tabular_ens), suffix="_ensembles")
basic_summary_table_tabular(_filter_ensembles!(copy(df_tabular_ens)), suffix="_ensembles_merge_ae_down")
@info "basic_summary_table_tabular"

# generic for both images and tabular data just change prefix and filter input
function plot_knowledge_ranks(df; prefix="tabular", suffix="", format="pdf", downsample=Dict{String, Int}())
    for (mn, metric) in zip(["AUC", "TPR@5"],[:auc, :tpr_5])
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        for (name, agg) in zip(
                    ["maxmean", "meanmax"],
                    [aggregate_stats_max_mean, aggregate_stats_mean_max])

            n = (prefix == "tabular") ? 4 : 1
            for (ctype, cnames, criterions) in zip(
                                ["pat", "pac", "patn"],
                                [PAT_METRICS_NAMES[n:end], PAC_METRICS_NAMES, PATN_METRICS_NAMES],
                                [_prefix_symbol.("val", PAT_METRICS[n:end]), _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("val", PATN_METRICS)])
                ranks = []
                metric_means = []
                extended_criterions = vcat(criterions, [val_metric, tst_metric])
                extended_cnames = vcat(cnames, ["\$$(mn)_{val}\$", "\$$(mn)_{tst}\$"])

                for criterion in extended_criterions
                    df_agg = agg(df, criterion; downsample=downsample)
                    sort!(df_agg, (order(:dataset), order(:modelname)))
                    rt = rank_table(df_agg, tst_metric)
                    push!(ranks, rt[end:end, :])
                    push!(metric_means, mean(Matrix(rt[1:end-3, 2:end]), dims=1))
                end

                df_ranks = vcat(ranks...)
                metric_mean = vcat(metric_means...)

                models = names(df_ranks)[2:end]

                a = PGFPlots.Axis([PGFPlots.Plots.Linear(
                                1:length(extended_criterions),
                                df_ranks[:, i + 1]) for (i, m) in enumerate(models)],
                        ylabel="avg. rnk",
                        style="xtick=$(_pgf_array(1:length(extended_criterions))),
                            xticklabels=$(_pgf_array(extended_cnames)),
                            width=6cm, height=4cm, scale only axis=true,
                            x tick label style={rotate=50,anchor=east}")
                b = PGFPlots.Axis([PGFPlots.Plots.Linear(
                                1:length(extended_criterions),
                                metric_mean[:,i],
                                legendentry=m) for (i, m) in enumerate(models)],
                        ylabel="avg. $mn",
                        legendStyle = "at={(0.5,1.02)}, anchor=south",
                        style="width=6cm, height=3cm, scale only axis=true,
                        xtick=$(_pgf_array(1:length(criterions))),
                        xticklabels={}, legend columns = -1")
                g = PGFPlots.GroupPlot(1, 2, groupStyle = "vertical sep = 0.5cm")
                push!(g, b); push!(g, a);
                filename = "$(projectdir())/paper/figures/$(prefix)_knowledge_rank_$(ctype)_$(metric)_$(name)$(suffix).$(format)"
                PGFPlots.save(filename, g; include_preamble=false)
            end
        end
    end
end

function plot_knowledge_tabular_repre(df, models; suffix="", format="pdf", downsample=Dict{String, Int}())
    filter!(x -> (x.modelname in models), df)
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    apply_aliases!(df, col="modelname", d=MODEL_ALIAS)
    plot_knowledge_ranks(df; prefix="tabular", suffix=suffix, format=format, downsample=downsample)
end

function plot_knowledge_tabular_type(df; suffix="", format="pdf")
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    apply_aliases!(df, col="modelname", d=MODEL_TYPE)
    plot_knowledge_ranks(df; prefix="tabular", suffix=suffix, format=format)
end

representatives=["ocsvm", "wae", "MAF", "fmgan", "vae_ocsvm", "vae+ocsvm"]
# plot_knowledge_tabular_repre(copy(df_tabular), representatives; format="tex", suffix="_representatives")
# plot_knowledge_tabular_repre(copy(df_tabular_ens), representatives; suffix="_representatives_ensembles", format="tex")
plot_knowledge_tabular_repre(
    _merge_ae_filter!(copy(df_tabular)), representatives; suffix="_representatives_merge_ae_down", format="tex", downsample=DOWNSAMPLE)
@info "plot_knowledge_tabular_repre"

# these ones are really costly as there are 12 aggregations over all models
# plot_knowledge_tabular_type(copy(df_tabular); format="tex")
# plot_knowledge_tabular_type(copy(df_tabular_ens), suffix="_ensembles", format="tex")
@info "plot_knowledge_tabular_type"

function rank_comparison_agg(df, tm=("AUC", :auc); prefix="tabular", suffix="")
    mn, metric = tm
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)

    ranks = []
    agg_names = ["maxmean", "meanmax", "maxmean10", "meanmax10"]
    agg_funct = [aggregate_stats_max_mean, aggregate_stats_mean_max, aggregate_stats_max_mean_top_10, aggregate_stats_mean_max_top_10]
    for (name, agg) in zip(agg_names, agg_funct)
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        df_agg = agg(df, val_metric)

        df_agg["model_type"] = copy(df_agg["modelname"])
        apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
        apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
        apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

        sort!(df_agg, (:dataset, :model_type, :modelname))

        rt = rank_table(df_agg, tst_metric)
        models = names(rt)[2:end]
        rt["agg"] = name

        select!(rt, vcat(["agg"], models))
        push!(ranks, rt[end:end,:])
    end
    df_ranks = reduce(vcat, ranks)

    hl_best_rank = LatexHighlighter(
                    (data, i, j) -> (data[i,j] == minimum(df_ranks[i, 2:end])),
                    ["color{red}","textbf"])

    filename = "$(projectdir())/paper/tables/$(prefix)_aggcomp_$(metric)$(suffix).tex"
    open(filename, "w") do io
        pretty_table(
            io, df_ranks,
            backend=:latex,
            formatters=ft_printf("%.1f"),
            highlighters=(hl_best_rank),
            nosubheader=true,
            tf=latex_booktabs
        )
    end
end

# superseeded by per_seed table
# rank_comparison_agg(copy(df_tabular), ("AUC", :auc))
# rank_comparison_agg(copy(df_tabular), ("TPR@5", :tpr_5))
# @info "rank_comparison_agg"

function rank_comparison_metric(df, ta=("maxmean", aggregate_stats_max_mean); suffix="", downsample=Dict{String, Int}())
    name, agg = ta
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)

    ranks = []
    for (mn, metric) in zip(["AUC", "TPR@5"],[:auc, :tpr_5])
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        df_agg = agg(df, val_metric, downsample=downsample)

        df_agg["model_type"] = copy(df_agg["modelname"])
        apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
        apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
        apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

        sort!(df_agg, (:dataset, :model_type, :modelname))

        rt = rank_table(df_agg, tst_metric)
        models = names(rt)[2:end]
        rt["crit"] = mn

        select!(rt, vcat(["crit"], models))
        push!(ranks, rt[end:end,:])
    end
    df_ranks = reduce(vcat, ranks)

    hl_best_rank = LatexHighlighter(
                    (data, i, j) -> (data[i,j] == minimum(df_ranks[i, 2:end])),
                    ["color{red}","textbf"])

    filename = "$(projectdir())/paper/tables/tabular_metriccomp_$(name)$(suffix).tex"
    open(filename, "w") do io
        pretty_table(
            io, df_ranks,
            backend=:latex,
            formatters=ft_printf("%.1f"),
            highlighters=(hl_best_rank),
            nosubheader=true,
            tf=latex_booktabs
        )
    end
end

# superseeded by per_seed table
# rank_comparison_metric(copy(df_tabular), ("maxmean", aggregate_stats_max_mean))
# rank_comparison_metric(copy(df_tabular), ("meanmax", aggregate_stats_mean_max))
# rank_comparison_metric(_merge_ae_filter!(copy(df_tabular)), ("maxmean", aggregate_stats_max_mean), suffix="_merge_ae_down", downsample=DOWNSAMPLE)
# rank_comparison_metric(_merge_ae_filter!(copy(df_tabular)), ("meanmax", aggregate_stats_mean_max), suffix="_merge_ae_down", downsample=DOWNSAMPLE)
# @info "rank_comparison_metric"

function comparison_ensemble(df, df_ensemble, tm=("AUC", :auc); downsample=Dict{String,Int}(), prefix="", suffix="")
    mn, metric = tm

    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    apply_aliases!(df_ensemble, col="modelname", d=MODEL_MERGE)

    models = unique(df_ensemble.modelname)
    filter!(x -> x.modelname in models, df)

    val_metric = _prefix_symbol("val", metric)
    tst_metric = _prefix_symbol("tst", metric)

    function _rank(d)
        df_agg = aggregate_stats_max_mean(d, val_metric; downsample=downsample)
        df_agg["model_type"] = copy(df_agg["modelname"])
        apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
        apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
        apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)
        sort!(df_agg, (:dataset, :model_type, :modelname))
        rank_table(df_agg, tst_metric)
    end

    rt = _rank(df)
    rt[end, 1] = "baseline"
    rt_ensemble = _rank(df_ensemble)
    rt_ensemble[end, 1] = "ensembles"

    df_ranks = vcat(rt[end:end, :], rt_ensemble[end:end, :])
    dif = Matrix(rt_ensemble[1:end-3, 2:end]) - Matrix(rt[1:end-3, 2:end])
    mean_dif = round.(mean(dif, dims=1), digits=2)
    push!(df_ranks, ["avg. change", mean_dif...])

    hl_best_rank = LatexHighlighter(
                    (data, i, j) -> (i < 3) && (data[i,j] == minimum(df_ranks[i, 2:end])),
                    ["color{red}","textbf"])

    hl_best_dif = LatexHighlighter(
                    (data, i, j) -> (i == 3) && (data[i,j] == maximum(df_ranks[i, 2:end])),
                    ["color{blue}","textbf"])

    f_float = (v, i, j) -> (j > 1) && (i == 3) ? ft_printf("%.2f")(v,i,j) : ft_printf("%.1f")(v,i,j)

    filename = "$(projectdir())/paper/tables/$(prefix)_ensemblecomp_$(metric)$(suffix).tex"
    open(filename, "w") do io
        pretty_table(
            io, df_ranks,
            backend=:latex,
            formatters=f_float,
            highlighters=(hl_best_rank, hl_best_dif),
            nosubheader=true,
            tf=latex_booktabs
        )
    end
end

function comparison_tabular_ensemble(df, df_ensemble, tm; downsample=DOWNSAMPLE, suffix="")
    comparison_ensemble(
        df, df_ensemble, tm;
        downsample=downsample, prefix="tabular", suffix=suffix)
end

# comparison_tabular_ensemble(copy(df_tabular), copy(df_tabular_ens), ("AUC", :auc))
# comparison_tabular_ensemble(copy(df_tabular), copy(df_tabular_ens), ("TPR@5", :tpr_5))
comparison_tabular_ensemble(
    _merge_ae_filter!(copy(df_tabular)),
    _filter_ensembles!(copy(df_tabular_ens)),
    ("AUC", :auc); downsample=DOWNSAMPLE, suffix="_merge_ae_down_mean")
comparison_tabular_ensemble(
    _merge_ae_filter!(copy(df_tabular)),
    _filter_ensembles!(copy(df_tabular_ens)),
    ("TPR@5", :tpr_5); downsample=DOWNSAMPLE, suffix="_merge_ae_down_mean")
@info "comparison_tabular_ensemble"

# join baseline and ensembles, allows to dig into where they help
# compare only those models for which ensembles are computed
# baseline = copy(df_tabular);
# models = unique(df_tabular_ens.modelname);
# ensembles = copy(vcat(filter(x -> (x.modelname in models), baseline), df_tabular_ens));
# comparison_tabular_ensemble(
#     baseline,
#     ensembles, ("AUC", :auc);
#     suffix="_only_improve")
# comparison_tabular_ensemble(
#     baseline,
#     ensembles, ("TPR@5", :tpr_5);
#     suffix="_only_improve")
# @info "comparison_tabular_ensemble_only_improve"

# how the choice of creation criterion and size of an ensembles affect results
function ensemble_sensitivity(df, df_ensemble; prefix="tabular", suffix="")
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    apply_aliases!(df_ensemble, col="modelname", d=MODEL_MERGE)

    for metric in [:auc, :tpr_5]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        function _rank(df_agg)
            df_agg["model_type"] = copy(df_agg["modelname"])
            apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
            apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
            apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

            sort!(df_agg, (:dataset, :model_type, :modelname))

            _rt = rank_table(df_agg, tst_metric)
            models = names(_rt)[2:end]
            _rt, models
        end

        ranks = []

        # computing baseline
        df_agg = aggregate_stats_max_mean(df, val_metric)
        rt, models = _rank(df_agg)
        rt["criterion-size"] = "baseline"
        select!(rt, vcat(["criterion-size"], models))
        push!(ranks, rt[end:end, :])

        # computing different ensembles
        for (cn, c) in zip(["AUC", "TPR@5"], [:val_auc, :val_tpr_5])
            for s in [5, 10]
                dff = filter(x -> startswith(x.parameters, "criterion=$(c)") && endswith(x.parameters, "_size=$(s)"), df_ensemble)
                df_agg = aggregate_stats_max_mean(dff, val_metric)
                rt_ensemble, models = _rank(df_agg)
                rt_ensemble["criterion-size"] = "$(cn)-$(s)"
                select!(rt_ensemble, vcat(["criterion-size"], models))

                dif = Matrix(rt_ensemble[1:end-3, 2:end]) - Matrix(rt[1:end-3, 2:end])
                mean_dif = round.(mean(dif, dims=1), digits=2)

                push!(rt_ensemble, ["avg. change", mean_dif...])
                push!(ranks, rt_ensemble[end-1:end, :])
            end
        end

        df_ranks = reduce(vcat, ranks)
        hl_best_rank = LatexHighlighter(
                        (data, i, j) -> ((i == 1) || (i%2 == 0)) && (data[i,j] == minimum(df_ranks[i, 2:end])),
                        ["color{red}","textbf"])

        hl_best_dif = LatexHighlighter(
                        (data, i, j) -> (i != 1) && (i%2 == 1) && (data[i,j] == maximum(df_ranks[i, 2:end])),
                        ["color{blue}","textbf"])

        f_float = (v, i, j) -> (i != 1) && (i%2 == 1) ? ft_printf("%.2f")(v,i,j) : ft_printf("%.1f")(v,i,j)

        filename = "$(projectdir())/paper/tables/$(prefix)_ensemble_size_$(metric)$(suffix).tex"
        open(filename, "w") do io
            pretty_table(
                io, df_ranks,
                backend=:latex,
                formatters=f_float,
                highlighters=(hl_best_rank, hl_best_dif),
                nosubheader=true,
                tf=latex_booktabs)
        end
    end
end

ensemble_sensitivity(
    _merge_ae_filter!(copy(df_tabular)),
    _filter_ensembles!(copy(df_tabular_ens)); prefix="tabular")
@info "ensemble_sensitivity_images"

# training time rank vs avg rank
function plot_tabular_fit_time(df; time_col=:fit_t, suffix="", format="pdf", downsample=Dict{String,Int}())
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    df["model_type"] = copy(df["modelname"])
    apply_aliases!(df, col="model_type", d=MODEL_TYPE)
    apply_aliases!(df, col="modelname", d=MODEL_ALIAS)

    # add first stage time for encoding and training of encoding
    df["fit_t"] .+= df["fs_fit_t"] .+ df["fs_fit_t"]
    # add all eval times together
    df["total_eval_t"] = df["tr_eval_t"] .+ df["tst_eval_t"] .+ df["val_eval_t"]

    for (mn, metric) in zip(["AUC", "TPR@5"],[:auc, :tpr_5])
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        # only meanmax to be consistent
        for (name, agg) in zip(
                    ["meanmax"],
                    [aggregate_stats_mean_max])

            df_agg = agg(df, val_metric, add_col=:total_eval_t, downsample=downsample)
            rt = rank_table(df_agg, tst_metric)

            # time is computed only for the best models
            # maybe should be averaged over all samples
            df_agg[time_col] .= -df_agg[time_col]
            rtt = rank_table(df_agg, time_col)

            models = names(rt)[2:end]
            x = Vector(rt[end, 2:end])
            y = Vector(rtt[end, 2:end])
            # labels cannot be shifted uniformly
            a = PGFPlots.Axis([
                    PGFPlots.Plots.Scatter(x, y),
                    [PGFPlots.Plots.Node(m, xx + 0.5, yy + 0.7) for (m, xx, yy) in zip(models, x, y)]...],
                    xlabel="avg. rnk",
                    ylabel="avg. time rnk")
            filename = "$(projectdir())/paper/figures/tabular_$(time_col)_vs_$(metric)_$(name)$(suffix).$(format)"
            PGFPlots.save(filename, a; include_preamble=false)
        end
    end
end

# plot_tabular_fit_time(copy(df_tabular); time_col=:fit_t, format="tex")
# plot_tabular_fit_time(copy(df_tabular); time_col=:tr_eval_t, format="tex")
# plot_tabular_fit_time(copy(df_tabular); time_col=:val_eval_t, format="tex")
# plot_tabular_fit_time(copy(df_tabular); time_col=:tst_eval_t, format="tex")
# plot_tabular_fit_time(copy(df_tabular); time_col=:total_eval_t, format="tex")

plot_tabular_fit_time(_merge_ae_filter!(copy(df_tabular)); time_col=:fit_t, format="tex", suffix="_merge_ae_down", downsample=DOWNSAMPLE)
plot_tabular_fit_time(_merge_ae_filter!(copy(df_tabular)); time_col=:tr_eval_t, format="tex", suffix="_merge_ae_down", downsample=DOWNSAMPLE)
plot_tabular_fit_time(_merge_ae_filter!(copy(df_tabular)); time_col=:val_eval_t, format="tex", suffix="_merge_ae_down", downsample=DOWNSAMPLE)
plot_tabular_fit_time(_merge_ae_filter!(copy(df_tabular)); time_col=:tst_eval_t, format="tex", suffix="_merge_ae_down", downsample=DOWNSAMPLE)
plot_tabular_fit_time(_merge_ae_filter!(copy(df_tabular)); time_col=:total_eval_t, format="tex", suffix="_merge_ae_down", downsample=DOWNSAMPLE)
@info "plot_tabular_fit_time"

# show basic table only for autoencoders
function tabular_autoencoders_investigation(df; split_vamp=false, suffix="")
    df["model_type"] = copy(df["modelname"])
    apply_aliases!(df, col="model_type", d=MODEL_TYPE)
    filter!(x -> (x.model_type == "autoencoders"), df)
    apply_aliases!(df, col="modelname", d=MODEL_ALIAS)
    apply_aliases!(df, col="dataset", d=DATASET_ALIAS)

    # split aae_full and wae_full based on used prior
    if split_vamp
        aaefv_mask = (df.modelname .== "aaef") .& occursin.("prior=vamp", df.parameters)
        df[aaefv_mask, :modelname] .=  "aaefv"
        waefv_mask = (df.modelname .== "waef") .& occursin.("prior=vamp", df.parameters)
        df[waefv_mask, :modelname] .=  "waefv"
    else
        apply_aliases!(df, col="modelname", d=Dict("aaev" => "aae", "waev" => "wae"))
    end

    # define further splitting of models based on parameters
    jc_mask = occursin.("jacodeco", df.parameters)
    df[jc_mask, :modelname] .=  df[jc_mask, :modelname] .*"-jc"

    lm_mask = occursin.("latent-mean", df.parameters)
    df[lm_mask, :modelname] .=  "disregard"

    l_mask = occursin.("latent_", df.parameters)
    df[l_mask, :modelname] .=  "disregard"

    ls_mask = occursin.("latent-sampled", df.parameters)
    df[ls_mask, :modelname] .=  "disregard"

    rm_mask = occursin.("reconstruction-mean", df.parameters)
    df[rm_mask, :modelname] .=  df[rm_mask, :modelname] .*"-rm"

    r_mask = occursin.("reconstruction_", df.parameters)
    df[r_mask, :modelname] .=  "disregard"

    rs_mask = occursin.("reconstruction-sampled", df.parameters)
    df[rs_mask, :modelname] .=  df[rs_mask, :modelname] .*"-rs"

    d_mask = occursin.("disc_", df.parameters)
    df[d_mask, :modelname] .=  df[d_mask, :modelname] .*"-d"

    el_mask = occursin.("elbo_", df.parameters)
    df[el_mask, :modelname] .=  df[el_mask, :modelname] .*"-el"

    mse_mask = occursin.("score=mse_", df.parameters)
    df[mse_mask, :modelname] .=  df[mse_mask, :modelname] .*"-rm"

    filter!(x -> (x.modelname != "disregard"), df)

    for metric in [:auc, :tpr_5]
        for (name, agg) in zip(
                    ["maxmean", "meanmax"],
                    [aggregate_stats_max_mean, aggregate_stats_mean_max])
            val_metric = _prefix_symbol("val", metric)
            tst_metric = _prefix_symbol("tst", metric)

            df_agg = agg(df, val_metric)
            sort!(df_agg, (:dataset, :modelname))
            rt = rank_table(df_agg, tst_metric)

            rt[end-2, 1] = "\$\\sigma_1\$"
            rt[end-1, 1] = "\$\\sigma_{10}\$"
            rt[end, 1] = "rnk"

            filename = "$(projectdir())/paper/tables/tabular_ae_only_$(metric)_$(metric)_$(name)$(suffix).tex"
            open(filename, "w") do io
                print_rank_table(io, rt; backend=:tex)
            end

            # due to its width we opted for box plot
            # recompute ranks for each dataset
            ranks = compute_ranks(rt[1:end-3, 2:end])
            models = names(rt)[2:end]

            #### grouping by method
            mgroups = ones(Int, length(models))
            mgroups[startswith.(models, "avae")] .= 1
            mgroups[startswith.(models, "gano")] .= 2
            mgroups[startswith.(models, "aae-")] .= 3
            mgroups[startswith.(models, "aaev")] .= 4
            mgroups[startswith.(models, "aaef-")] .= 5
            mgroups[startswith.(models, "aaefv")] .= 6
            mgroups[startswith.(models, "vae-")] .= 7
            mgroups[startswith.(models, "vaef")] .= 8
            mgroups[startswith.(models, "vaes")] .= 9
            mgroups[startswith.(models, "wae-")] .= 10
            mgroups[startswith.(models, "waev")] .= 11
            mgroups[startswith.(models, "waef-")] .= 12
            mgroups[startswith.(models, "waefv")] .= 13

            sm = sortperm(mgroups)
            mgroups = mgroups[sm]

            models = models[sm]
            ranks = ranks[sm, :]

            #### grouping by score, 5 types + 1 for other scores
            groups = 6*ones(Int, length(models))
            groups[endswith.(models, "-rm")] .= 1
            groups[endswith.(models, "-rs")] .= 2
            groups[endswith.(models, "-jc")] .= 3
            groups[endswith.(models, "-d")] .= 4
            groups[endswith.(models, "-el")] .= 5
            ####

            # compute statistics for boxplot
            rmin, rmax = maximum(ranks, dims=2), maximum(ranks, dims=2)
            rmean = mean(ranks, dims=2)
            rmedian = median(ranks, dims=2)
            rlowq = quantile.(eachrow(ranks), 0.25)
            rhighq = quantile.(eachrow(ranks), 0.75)

            a = pgf_boxplot_grouped_colorpos(rlowq, rhighq, rmedian, rmin, rmax, models, groups, mgroups; h="12cm", w="8cm")
            filename = "$(projectdir())/paper/figures/tabular_ae_only_box_$(metric)_$(name)$(suffix).tex"
            open(filename, "w") do f
                write(f, a)
            end
        end
    end
end

# tabular_autoencoders_investigation(copy(df_tabular))
tabular_autoencoders_investigation(copy(df_tabular); split_vamp=true, suffix="_split_vamp")
@info "tabular_autoencoders_investigation"

# does crossvalidation matters?
function per_seed_ranks_tabular(df; suffix="", downsample=Dict{String, Int}())
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)

    for metric in [:auc, :tpr_5]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        function _rank(df_agg)
            df_agg["model_type"] = copy(df_agg["modelname"])
            apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
            apply_aliases!(df_agg, col="dataset", d=DATASET_ALIAS)
            apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

            sort!(df_agg, (:dataset, :model_type, :modelname))

            rt = rank_table(df_agg, tst_metric)
            models = names(rt)[2:end]
            rt, models
        end

        ranks = []
        for (name, agg) in zip(
                ["maxmean", "meanmax"],
                [aggregate_stats_max_mean, aggregate_stats_mean_max])
            df_agg = agg(df, val_metric; downsample=downsample)
            rt, models = _rank(df_agg)
            rt["exp"] = name
            select!(rt, vcat(["exp"], models))
            push!(ranks, rt[end:end,:])
        end

        for seed in 1:5
            dff = filter(x -> (x.seed == seed), df)
            # silence the warnings for insufficient number of seeds
            df_agg = aggregate_stats_max_mean(dff, val_metric; verbose=false, downsample=downsample)
            rt, models = _rank(df_agg)
            rt["exp"] = "seed=$seed"
            select!(rt, vcat(["exp"], models))
            push!(ranks, rt[end:end,:])
        end

        df_ranks = reduce(vcat, ranks)
        hl_best_rank = LatexHighlighter(
                        (data, i, j) -> (data[i,j] == minimum(df_ranks[i, 2:end])),
                        ["color{red}","textbf"])

        filename = "$(projectdir())/paper/tables/tabular_per_seed_ranks_$(metric)$(suffix).tex"
        open(filename, "w") do io
            pretty_table(
                io, df_ranks,
                backend=:latex,
                formatters=ft_printf("%.1f"),
                highlighters=(hl_best_rank),
                nosubheader=true,
                tf=latex_booktabs)
        end
    end
end

# per_seed_ranks_tabular(copy(df_tabular))
# per_seed_ranks_tabular(_merge_ae_filter!(copy(df_tabular)), suffix="_merge_ae")
per_seed_ranks_tabular(_merge_ae_filter!(copy(df_tabular)), suffix="_merge_ae_down", downsample=DOWNSAMPLE)
# per_seed_ranks_tabular(copy(df_tabular_ens); suffix="_ensembles")
@info "per_seed_ranks_tabular"

######################################################################################
#######################               IMAGES                ##########################
######################################################################################
df_images = load(datadir("evaluation/images_eval.bson"))[:df];
df_images_loi = load(datadir("evaluation/images_leave-one-in_eval.bson"))[:df];
df_images_ens = load(datadir("evaluation_ensembles/images_eval.bson"))[:df];
@info "Loaded results from images"

function filter_img_models!(df)
    filter!(x -> (x.modelname != "aae_ocsvm") && (x.modelname != "fAnoGAN-GP"), df)
end

function basic_summary_table_images(df; suffix="")
    filter!(x -> x.seed == 1, df)
    basic_summary_table(df; prefix="images", suffix=suffix)
end

# basic_summary_table_images(copy(df_images))
# basic_summary_table_images(copy(df_images_ens), suffix="_ensembles")

basic_summary_table_images(filter_img_models!(copy(df_images)), suffix="_filter")
basic_summary_table_images(copy(df_images_loi), suffix="_loi")
basic_summary_table_images(filter_img_models!(copy(df_images_ens)), suffix="_ensembles_filter")
@info "basic_summary_table_images"

function rank_comparison_agg_images(df, tm=("AUC", :auc); suffix="")
    filter!(x -> (x.seed == 1), df)
    rank_comparison_agg(df, tm; prefix="images", suffix=suffix)
end

# rank_comparison_agg_images(copy(df_images), ("AUC", :auc))
# rank_comparison_agg_images(copy(df_images), ("TPR@5", :tpr_5))
@info "rank_comparison_agg_images"

function comparison_images_ensemble(df, df_ensemble, tm=("AUC", :auc); suffix="")
    filter!(x -> (x.seed == 1), df)
    filter!(x -> (x.seed == 1), df_ensemble)
    comparison_ensemble(df, df_ensemble, tm; prefix="images", suffix=suffix)
end

# comparison_images_ensemble(copy(df_images), copy(df_images_ens), ("AUC", :auc))
# comparison_images_ensemble(copy(df_images), copy(df_images_ens), ("TPR@5", :tpr_5))
comparison_images_ensemble(
    filter_img_models!(copy(df_images)),
    _filter_ensembles!(filter_img_models!(copy(df_images_ens))),
    ("AUC", :auc); suffix="_filter")
comparison_images_ensemble(
    filter_img_models!(copy(df_images)),
    _filter_ensembles!(filter_img_models!(copy(df_images_ens))),
    ("TPR@5", :tpr_5); suffix="_filter")

@info "comparison_images_ensemble"

function ensemble_sensitivity_images(df, df_ensemble)
    filter!(x -> (x.seed == 1), df)
    filter!(x -> (x.seed == 1), df_ensemble)
    ensemble_sensitivity(df, df_ensemble; prefix="images")
end

ensemble_sensitivity_images(
        filter_img_models!(copy(df_images)),
        _filter_ensembles!(filter_img_models!(copy(df_images_ens))))

@info "ensemble_sensitivity_images"

# join baseline and ensembles, allows to dig into where they help
# compare only those models for which ensembles are computed
# baseline_img = copy(df_images);
# models_img = unique(df_images_ens.modelname);
# ensembles_img = copy(vcat(filter(x -> (x.modelname in models_img), baseline_img), df_images_ens));
# comparison_images_ensemble(
#     baseline_img,
#     ensembles_img, ("AUC", :auc);
#     suffix="_only_improve")
# comparison_images_ensemble(
#     baseline_img,
#     ensembles_img, ("TPR@5", :tpr_5);
#     suffix="_only_improve")
# @info "comparison_images_ensemble_only_improve"

# basic tables when anomaly_class is treated as separate dataset
function basic_tables_images_per_ac(df; suffix="")
    dff = filter(x -> (x.seed == 1), df)
    apply_aliases!(dff, col="dataset", d=DATASET_ALIAS)
    for d in unique(dff.dataset)
        mask = (dff.dataset .== d)
        dff[mask, :dataset] .= dff[mask, :dataset] .* ":" .* convert_anomaly_class.(dff[mask, :anomaly_class], d)
    end
    select!(dff, Not(:anomaly_class))

    apply_aliases!(dff, col="modelname", d=MODEL_MERGE)
    for metric in [:auc, :tpr_5]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        # does not play well with the aggregation when anomaly_class is present
        df_agg = aggregate_stats_mean_max(dff; min_samples=1)

        df_agg["model_type"] = copy(df_agg["modelname"])
        apply_aliases!(df_agg, col="modelname", d=MODEL_ALIAS)
        apply_aliases!(df_agg, col="model_type", d=MODEL_TYPE)

        sort!(df_agg, (:dataset, :model_type, :modelname))
        rt = rank_table(df_agg, tst_metric)
        rt[end-2, 1] = "\$\\sigma_1\$"
        rt[end-1, 1] = "\$\\sigma_{10}\$"
        rt[end, 1] = "rnk"

        filename = "$(projectdir())/paper/tables/images_per_ac_$(metric)_$(metric)$(suffix).tex"
        open(filename, "w") do io
            print_rank_table(io, rt; backend=:tex)
        end
    end
end

# these results are part of supplementary
# df_images_vae = vcat(load(datadir("evaluation/vae_seeds_eval.bson"))[:df], filter(x -> (x.modelname == "vae"), df_images));
# df_images_ganomaly = filter(x -> (x.modelname == "Conv-GANomaly"), df_images)
# function split_seeds!(df)
#     df["modelname"] = string.(df["seed"]) .* df["modelname"]
#     df
# end
# basic_summary_table(split_seeds!(copy(df_images_vae)), prefix="images", suffix="_vae_seed")
# basic_summary_table(split_seeds!(copy(df_images_ganomaly)), prefix="images", suffix="_gano_seed")
# basic_tables_images_per_ac(split_seeds!(copy(df_images_vae)), suffix="_vae_seed")
# basic_tables_images_per_ac(split_seeds!(copy(df_images_ganomaly)), suffix="_gano_seed")

# basic_tables_images_per_ac(copy(df_images))
# basic_tables_images_per_ac(copy(df_images_ens), suffix="_ensembles")
basic_tables_images_per_ac(filter_img_models!(copy(df_images)), suffix="_filter")
basic_tables_images_per_ac(copy(df_images_loi), suffix="_loi")
basic_tables_images_per_ac(filter_img_models!(copy(df_images_ens)), suffix="_ensembles_filter")
@info "basic_tables_images_per_ac"

function plot_knowledge_images_repre(df, models; suffix="", format="pdf")
    filter!(x -> (x.seed == 1), df)
    filter!(x -> (x.modelname in models), df)
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    apply_aliases!(df, col="modelname", d=MODEL_ALIAS)
    plot_knowledge_ranks(df; prefix="images", suffix=suffix, format=format)
end

function plot_knowledge_images_type(df; suffix="", format="pdf")
    filter!(x -> (x.seed == 1), df)
    apply_aliases!(df, col="modelname", d=MODEL_MERGE)
    apply_aliases!(df, col="modelname", d=MODEL_TYPE)
    plot_knowledge_ranks(df; prefix="images", suffix=suffix, format=format)
end

representatives=["ocsvm", "aae", "fmgan", "vae_ocsvm"]
plot_knowledge_images_repre(copy(df_images), representatives; format="tex", suffix="_representatives")
# plot_knowledge_images_repre(copy(df_images_ens), representatives; suffix="_representatives_ensembles", format="tex")
@info "plot_knowledge_images_repre"

# these ones are really costly as there are 12 aggregations over all models
# plot_knowledge_images_type(copy(df_images); format="tex")
# plot_knowledge_images_type(copy(df_images_ens), suffix="_ensembles", format="tex")
# @info "plot_knowledge_images_type"

# this is just specific figure for the main text
function plot_knowledge_combined(df_tab, df_img; format="pdf")
    filter!(x -> (x.seed == 1), df_img);
    apply_aliases!(df_img, col="modelname", d=MODEL_MERGE);
    apply_aliases!(df_img, col="modelname", d=MODEL_TYPE);

    apply_aliases!(df_tab, col="modelname", d=MODEL_MERGE);
    apply_aliases!(df_tab, col="modelname", d=MODEL_TYPE);

    for (mn, metric) in [collect(zip(["AUC", "TPR@5"],[:auc, :tpr_5]))[1]]
        val_metric = _prefix_symbol("val", metric)
        tst_metric = _prefix_symbol("tst", metric)

        for (ctype, cnames, criterions) in [collect(zip(
                            ["pat", "pac", "patn"],
                            [PAT_METRICS_NAMES, PAC_METRICS_NAMES, PATN_METRICS_NAMES],
                            [_prefix_symbol.("val", PAT_METRICS), _prefix_symbol.("val", PAC_METRICS), _prefix_symbol.("val", PATN_METRICS)]))[end]]

            function _rank(df, criterions, agg)
                ranks, metric_means = [], []
                for criterion in criterions
                    df_agg = agg(df, criterion)
                    sort!(df_agg, (order(:dataset), order(:modelname)))
                    rt = rank_table(df_agg, tst_metric)
                    push!(ranks, rt[end:end, 2:end])
                    push!(metric_means, mean(Matrix(rt[1:end-3, 2:end]), dims=1))
                end
                vcat(ranks...), vcat(metric_means...)
            end

            function _plot(df_ranks, metric_mean, criterions, cnames, models)
                a = PGFPlots.Axis([PGFPlots.Plots.Linear(
                            1:length(criterions),
                            df_ranks[:, i]) for (i, m) in enumerate(models)],
                    ylabel="avg. rnk",
                    style="xtick=$(_pgf_array(1:length(criterions))),
                        xticklabels=$(_pgf_array(cnames)),
                        width=5cm, height=3cm, scale only axis=true,
                        x tick label style={rotate=50,anchor=east}")
                b = PGFPlots.Axis([PGFPlots.Plots.Linear(
                            1:length(criterions),
                            metric_mean[:, i],
                                legendentry=m) for (i, m) in enumerate(models)],
                    ylabel="avg. $mn",
                    legendStyle = "at={(0.1,1.15)}, anchor=west",
                    style="width=5cm, height=3cm, scale only axis=true,
                    xtick=$(_pgf_array(1:length(criterions))),
                    xticklabels={}, legend columns = -1")
                a, b
            end

            extended_criterions = (ctype == "pat") ? criterions[4:end] : criterions
            extended_criterions = vcat(extended_criterions, [val_metric, tst_metric])
            extended_cnames = (ctype == "pat") ? cnames[4:end] : cnames
            extended_cnames = vcat(extended_cnames, ["\$$(mn)_{val}\$", "\$$(mn)_{tst}\$"])

            ranks_tab, metric_means_tab = _rank(df_tab, extended_criterions, aggregate_stats_mean_max)
            p = [1,2,4,5,3] # make sure that flows are last
            models_tab = names(ranks_tab)[p]
            select!(ranks_tab, p)
            metric_means_tab = metric_means_tab[:, p]
            a_tab, b_tab = _plot(ranks_tab, metric_means_tab, extended_criterions, extended_cnames, models_tab)

            extended_criterions = vcat(criterions, [val_metric, tst_metric])
            extended_cnames = vcat(cnames, ["\$$(mn)_{val}\$", "\$$(mn)_{tst}\$"])

            ranks_img, metric_means_img = _rank(df_img, extended_criterions, aggregate_stats_max_mean)
            models_img = names(ranks_img)
            a_img, b_img = _plot(ranks_img, metric_means_img, extended_criterions, extended_cnames, models_img)
            a_img.ylabel = ""
            b_img.ylabel = ""
            for b in b_img.plots
                b.legendentry=nothing
            end

            g = PGFPlots.GroupPlot(2, 2, groupStyle = "vertical sep = 0.5cm, horizontal sep = 1.0cm")
            push!(g, b_tab); push!(g, b_img);
            push!(g, a_tab); push!(g, a_img);

            filename = "$(projectdir())/paper/figures/combined_knowledge_rank_$(ctype)_$(metric).$(format)"
            PGFPlots.save(filename, g; include_preamble=false)
        end
    end
end

# no filter applied here as we aggregate all the models into groups
# shows the best possible performance
plot_knowledge_combined(copy(df_tabular), copy(df_images); format="tex")
@info "plot_combine_knowledge_type"

@info "----------------- DONE ---------------------"


text:
function multiply_slice_matrix_left!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    N=size(bond_ch,2) #number of bonds in checkerboard
    num_ch=mc.p.num_ch
    time_slice = cld(slice,num_ch) # imagenary time
    cb = mod1(slice,num_ch) # number of checkerboard in time slice

	if cb==4
		diag_terms = mc.diag_terms
		M .= M.*diag_terms
	end

    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
            bond_index = bond_ch[1,n,cb]
            if mc.conf[bond_index,time_slice]==1
                hopping_mat = mc.hopping_term_up
            else
                hopping_mat=mc.hopping_term_down
            end

                B =hopping_mat*[M[i,:]' ; M[j,:]']

	        M[i,:] .= B[1,:]
	        M[j,:] .= B[2,:]
	    end
	end

    nothing
end
function multiply_slice_matrix_right!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    N=size(bond_ch,2) #number of bonds in checkerboard
    num_ch=mc.p.num_ch
    time_slice = cld(slice,num_ch) # imagenary time
    cb = mod1(slice,num_ch) # number of checkerboard in time slice
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
            bond_index = bond_ch[1,n,cb]
            if mc.conf[bond_index,time_slice]==1
                hopping_mat = mc.hopping_term_up
            else
                hopping_mat = mc.hopping_term_down
            end

                B = [M[i,:]  M[j,:]]*hopping_mat

	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
	if cb==4
		diag_terms = mc.diag_terms
		M .= M.*diag_terms'
	end
    nothing
end

function multiply_slice_matrix_inv_left!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    N=size(bond_ch,2) #number of bonds in checkerboard
    num_ch=mc.p.num_ch
    time_slice = cld(slice,num_ch) # imagenary time
    cb = mod1(slice,num_ch) # number of checkerboard in time slice
    h_inv = zeros(Float64,2,2)
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
            bond_index = bond_ch[1,n,cb]
            if mc.conf[bond_index,time_slice]==1
                hopping_mat = mc.hopping_term_up
            else
                hopping_mat=mc.hopping_term_down
            end
	        h_inv[1,1]=hopping_mat[1,1]
	        h_inv[2,2]=hopping_mat[2,2]
	        h_inv[1,2]=-hopping_mat[1,2]
	        h_inv[2,1]=-hopping_mat[2,1]
	        B = h_inv*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end
	if cb==4
		diag_terms_inv = mc.diag_terms_inv
		M .= M.*diag_terms_inv
	end
    nothing
end
function multiply_slice_matrix_inv_right!(mc::DQMC_bond, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    N=size(bond_ch,2) #number of bonds in checkerboard
    num_ch=mc.p.num_ch
    time_slice = cld(slice,num_ch) # imagenary time
    cb = mod1(slice,num_ch) # number of checkerboard in time slice
    if cb==4
		diag_terms_inv = mc.diag_terms_inv
		M .= M.*diag_terms_inv'
	end
	h_inv = zeros(Float64,2,2)
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
            bond_index = bond_ch[1,n,cb]
            if mc.conf[bond_index,time_slice]==1
                hopping_mat = mc.hopping_term_up
            else
                hopping_mat=mc.hopping_term_down
            end
            h_inv[1,1]=hopping_mat[1,1]
           h_inv[2,2]=hopping_mat[2,2]
           h_inv[1,2]=-hopping_mat[1,2]
           h_inv[2,1]=-hopping_mat[2,1]
	        B = [M[i,:]  M[j,:]]*h_inv
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
    nothing
end
function multiply_daggered_slice_matrix_left!(mc::DQMC_bond, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    N=size(bond_ch,2) #number of bonds in checkerboard
    num_ch=mc.p.num_ch
    time_slice = cld(slice,num_ch) # imagenary time
    cb = mod1(slice,num_ch) # number of checkerboard in time slice
	if cb==4
		diag_terms = mc.diag_terms
		M .= M.*diag_terms
	end
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
            bond_index = bond_ch[1,n,cb]
            if mc.conf[bond_index,time_slice]==1
                hopping_mat = mc.hopping_term_up
            else
                hopping_mat=mc.hopping_term_down
            end
	        B = adjoint(hopping_mat)*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end
    nothing
end
function multiply_daggered_slice_matrix_right!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    N=size(bond_ch,2) #number of bonds in checkerboard
    num_ch=mc.p.num_ch
    time_slice = cld(slice,num_ch) # imagenary time
    cb = mod1(slice,num_ch) # number of checkerboard in time slice
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
            bond_index = bond_ch[1,n,cb]
            if mc.conf[bond_index,time_slice]==1
                hopping_mat = mc.hopping_term_up
            else
                hopping_mat=mc.hopping_term_down
            end
	        B = [M[i,:]  M[j,:]]*adjoint(hopping_mat)
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
	if cb==4
		diag_terms = mc.diag_terms
		M .= M.*diag_terms'
	end
    nothing
end
#=
function multiply_slice_matrix_left!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
	cb = mod1(slice,4) # !!4 is the number of checkerboards - replace later with a variable!!
	if cb==4
		diag_terms = mc.diag_terms
		M = M.*diag_terms
	end

    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        B = hopping_mat[:,:,n,slice,cb]*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end

    nothing
end
function multiply_slice_matrix_right!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
	cb = mod1(slice,4) # !!4 is the number of checkerboards - replace later in a variable!!
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        B = [M[i,:]  M[j,:]]*hopping_mat[:,:,n,slice,cb]
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
	if cb==4
		diag_terms = mc.diag_terms
		M = M.*diag_terms'
	end
    nothing
end

function multiply_slice_matrix_inv_left!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
	cb = mod1(slice,4) # !!4 is the number of checkerboards - replace later in a variable!!
    h_inv = zeros(Float64,2,2)
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]

	        h_inv[1,1]=hopping_mat[1,1,n,slice,cb]
	        h_inv[2,2]=hopping_mat[2,2,n,slice,cb]
	        h_inv[1,2]=-hopping_mat[1,2,n,slice,cb]
	        h_inv[2,1]=-hopping_mat[2,1,n,slice,cb]
	        B = h_inv*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end
	if cb==4
		diag_terms_inv = mc.diag_terms_inv
		M = M.*diag_terms_inv
	end
    nothing
end
function multiply_slice_matrix_inv_right!(mc::DQMC_bond, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
	cb = mod1(slice,4) # !!4 is the number of checkerboards - replace later in a variable!!
	if cb==4
		diag_terms_inv = mc.diag_terms_inv
		M = M.*diag_terms_inv'
	end
	h_inv = zeros(Float64,2,2)
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        h_inv[1,1]=hopping_mat[1,1,n,slice,cb]
	        h_inv[2,2]=hopping_mat[2,2,n,slice,cb]
	        h_inv[1,2]=-hopping_mat[1,2,n,slice,cb]
	        h_inv[2,1]=-hopping_mat[2,1,n,slice,cb]
	        B = [M[i,:]  M[j,:]]*h_inv
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
    nothing
end
function multiply_daggered_slice_matrix_left!(mc::DQMC_bond, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
	cb = mod1(slice,4) # !!4 is the number of checkerboards - replace later in a variable!!
	if cb==4
		diag_terms = mc.diag_terms
		M = M.*diag_terms
	end
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        B = adjoint(hopping_mat[:,:,n,slice,cb])*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end
    nothing
end
function multiply_daggered_slice_matrix_right!(mc::DQMC_bond,  slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
	cb = mod1(slice,4) # !!4 is the number of checkerboards - replace later in a variable!!
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        B = [M[i,:]  M[j,:]]*adjoint(hopping_mat[:,:,n,slice,cb])
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
	if cb==4
		diag_terms = mc.diag_terms
		M = M.*diag_terms'
	end
    nothing
end
=#
#=
"""
	slice_matrix(mc::DQMC_CBFalse, m::Model, slice::Int, power::Float64=1.)

Direct calculation of effective slice matrix, i.e. no checkerboard.
Calculates `Beff(slice) = exp(−1/2∆tauT)exp(−1/2∆tauT)exp(−∆tauV(slice))`.
"""
function slice_matrix(mc::DQMC_CBFalse, m::Model, slice::Int,
					power::Float64=1.)
	eT = mc.s.hopping_matrix_exp
	eTinv = mc.s.hopping_matrix_exp_inv
	eV = mc.s.eV

	interaction_matrix_exp!(mc, m, eV, mc.conf, slice, power)

	if power > 0
		return eT * eT * eV
	else
		return eV * eTinv * eTinv
	end
end
=#

# CheckerboardTrue
#=
const DQMC_CBTrue = DQMC{M, CheckerboardTrue} where M

function slice_matrix(mc::DQMC_CBTrue, m::Model, slice::Int, power::Float64=1.)
  M = eye(heltype(mc), m.flv*m.l.sites)
  if power > 0
    multiply_slice_matrix_left!(mc, m, slice, M)
  else
    multiply_slice_matrix_inv_left!(mc, m, slice, M)
  end
  return M
end
=#
#=
function multiply_slice_matrix_inv_left!(mc::DQMC_bond, cb::Int, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
    h_inv = zeros(Float64,2,2)
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]

	        h_inv[1,1]=hopping_mat[1,1,n,slice,cb]
	        h_inv[2,2]=hopping_mat[2,2,n,slice,cb]
	        h_inv[1,2]=-hopping_mat[1,2,n,slice,cb]
	        h_inv[2,1]=-hopping_mat[2,1,n,slice,cb]
	        B = h_inv*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end
    nothing
end
function multiply_slice_matrix_inv_right!(mc::DQMC_bond, cb::Int, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
    h_inv = zeros(Float64,2,2)
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        h_inv[1,1]=hopping_mat[1,1,n,slice,cb]
	        h_inv[2,2]=hopping_mat[2,2,n,slice,cb]
	        h_inv[1,2]=-hopping_mat[1,2,n,slice,cb]
	        h_inv[2,1]=-hopping_mat[2,1,n,slice,cb]
	        B = [M[i,:]  M[j,:]]*h_inv
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
    nothing
end
function multiply_daggered_slice_matrix_left!(mc::DQMC_bond, cb::Int, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        B = adjoint(hopping_mat[:,:,n,slice,cb])*[M[i,:]' ; M[j,:]']
	        M[i,:] = B[1,:]
	        M[j,:] = B[2,:]
	    end
	end
    nothing
end
function multiply_daggered_slice_matrix_right!(mc::DQMC_bond,  cb::Int, slice::Int, M::AbstractMatrix{T}) where T<:Number
    bond_ch = mc.model.bond_checkerboard
    hopping_mat = mc.hopping_mat#[:,:,n,slice,cb]
    N=size(hopping_mat,3) #number of bonds in checkerboard
    @inbounds @views begin
	    for n in 1:N
	        i = bond_ch[2,n,cb]
	        j = bond_ch[3,n,cb]
	        B = [M[i,:]  M[j,:]]*adjoint(hopping_mat[:,:,n,slice,cb])
	        M[:,i] = B[:,1]
	        M[:,j] = B[:,2]
	    end
	end
    nothing
end
=#

ext: jl

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
using TableFrames
using NamedTuples
using Base.Test

@testset "DataFrame" begin
    df = DataFrame(:id => collect(1:5), :val => collect(10:10:50))
    @test df[:id] == collect(1:5)
    @test df[:id, 2] == 2
    @test df[:val, 2] == 20
    r = df[2]
    @test r.id == 2
    @test r.val == 20

    df[:val, 2] = 100
    @test df[:val, 2] == 100
    df[2] = record(:id => 10, :val => 75)
    r = df[2]
    @test r.id == 10
    @test r.val == 75

    rs = collect(rows(df))
    @test length(rs) == 5
    @test isa(rs, Vector)
    @test isa(rs[1], NamedTuple)

    df_str = sprint(show, df)
end


text:
function inverse_characteristic_matrix_ellipsoid(q0::Float64, q1::Float64, q2::Float64, q3::Float64, r1::Float64, r2::Float64, r3::Float64)

	a11::Float64 = (4.0 * r1^2 * r2^2 * q0^2 * q2^2 - 8.0 * r1^2 * r2^2 * q0 * q1 * q2 * q3 + 4.0 * r1^2 * r2^2 * q1^2 * q3^2 + 4.0 * r1^2 * r3^2 * q0^2 * q3^2 + 8.0 * r1^2 * r3^2 * q0 * q1 * q2 * q3 + 4.0 * r1^2 * r3^2 * q1^2 * q2^2 + r2^2 * r3^2 * q0^4 + 2.0 * r2^2 * r3^2 * q0^2 * q1^2 - 2.0 * r2^2 * r3^2 * q0^2 * q2^2 - 2.0 * r2^2 * r3^2 * q0^2 * q3^2 + r2^2 * r3^2 * q1^4 - 2.0 * r2^2 * r3^2 * q1^2 * q2^2 - 2.0 * r2^2 * r3^2 * q1^2 * q3^2 + r2^2 * r3^2 * q2^4 + 2.0 * r2^2 * r3^2 * q2^2 * q3^2 + r2^2 * r3^2 * q3^4)/(r1^2 * r2^2 * r3^2)
	a12::Float64 = -(2.0 * (2.0 * r1^2 * r2^2 * q0^2 * q1 * q2 - 2.0 * r1^2 * r2^2 * q0 * q1^2 * q3 + 2.0 * r1^2 * r2^2 * q0 * q2^2 * q3 - 2.0 * r1^2 * r2^2 * q1 * q2 * q3^2 - r1^2 * r3^2 * q0^3 * q3 - r1^2 * r3^2 * q0^2 * q1 * q2 + r1^2 * r3^2 * q0 * q1^2 * q3 - r1^2 * r3^2 * q0 * q2^2 * q3 + r1^2 * r3^2 * q0 * q3^3 + r1^2 * r3^2 * q1^3 * q2 - r1^2 * r3^2 * q1 * q2^3 + r1^2 * r3^2 * q1 * q2 * q3^2 + r2^2 * r3^2 * q0^3 * q3 - r2^2 * r3^2 * q0^2 * q1 * q2 + r2^2 * r3^2 * q0 * q1^2 * q3 - r2^2 * r3^2 * q0 * q2^2 * q3 - r2^2 * r3^2 * q0 * q3^3 - r2^2 * r3^2 * q1^3 * q2 + r2^2 * r3^2 * q1 * q2^3 + r2^2 * r3^2 * q1 * q2 * q3^2))/(r1^2 * r2^2 * r3^2)
	a13::Float64 = (2.0 * (- r1^2 * r2^2 * q0^3 * q2 + r1^2 * r2^2 * q0^2 * q1 * q3 + r1^2 * r2^2 * q0 * q1^2 * q2 + r1^2 * r2^2 * q0 * q2^3 - r1^2 * r2^2 * q0 * q2 * q3^2 - r1^2 * r2^2 * q1^3 * q3 - r1^2 * r2^2 * q1 * q2^2 * q3 + r1^2 * r2^2 * q1 * q3^3 - 2.0 * r1^2 * r3^2 * q0^2 * q1 * q3 - 2.0 * r1^2 * r3^2 * q0 * q1^2 * q2 + 2.0 * r1^2 * r3^2 * q0 * q2 * q3^2 + 2.0 * r1^2 * r3^2 * q1 * q2^2 * q3 + r2^2 * r3^2 * q0^3 * q2 + r2^2 * r3^2 * q0^2 * q1 * q3 + r2^2 * r3^2 * q0 * q1^2 * q2 - r2^2 * r3^2 * q0 * q2^3 - r2^2 * r3^2 * q0 * q2 * q3^2 + r2^2 * r3^2 * q1^3 * q3 - r2^2 * r3^2 * q1 * q2^2 * q3 - r2^2 * r3^2 * q1 * q3^3))/(r1^2 * r2^2 * r3^2)
    a21::Float64 = -(2.0 * (2.0 * r1^2 * r2^2 * q0^2 * q1 * q2 - 2.0 * r1^2 * r2^2 * q0 * q1^2 * q3 + 2.0 * r1^2 * r2^2 * q0 * q2^2 * q3 - 2.0 * r1^2 * r2^2 * q1 * q2 * q3^2 - r1^2 * r3^2 * q0^3 * q3 - r1^2 * r3^2 * q0^2 * q1 * q2 + r1^2 * r3^2 * q0 * q1^2 * q3 - r1^2 * r3^2 * q0 * q2^2 * q3 + r1^2 * r3^2 * q0 * q3^3 + r1^2 * r3^2 * q1^3 * q2 - r1^2 * r3^2 * q1 * q2^3 + r1^2 * r3^2 * q1 * q2 * q3^2 + r2^2 * r3^2 * q0^3 * q3 - r2^2 * r3^2 * q0^2 * q1 * q2 + r2^2 * r3^2 * q0 * q1^2 * q3 - r2^2 * r3^2 * q0 * q2^2 * q3 - r2^2 * r3^2 * q0 * q3^3 - r2^2 * r3^2 * q1^3 * q2 + r2^2 * r3^2 * q1 * q2^3 + r2^2 * r3^2 * q1 * q2 * q3^2))/(r1^2 * r2^2 * r3^2)
	a22::Float64 = (4.0 * r1^2 * r2^2 * q0^2 * q1^2 + 8.0 * r1^2 * r2^2 * q0 * q1 * q2 * q3 + 4.0 * r1^2 * r2^2 * q2^2 * q3^2 + r1^2 * r3^2 * q0^4 - 2.0 * r1^2 * r3^2 * q0^2 * q1^2 + 2.0 * r1^2 * r3^2 * q0^2 * q2^2 - 2.0 * r1^2 * r3^2 * q0^2 * q3^2 + r1^2 * r3^2 * q1^4 - 2.0 * r1^2 * r3^2 * q1^2 * q2^2 + 2.0 * r1^2 * r3^2 * q1^2 * q3^2 + r1^2 * r3^2 * q2^4 - 2.0 * r1^2 * r3^2 * q2^2 * q3^2 + r1^2 * r3^2 * q3^4 + 4.0 * r2^2 * r3^2 * q0^2 * q3^2 - 8.0 * r2^2 * r3^2 * q0 * q1 * q2 * q3 + 4.0 * r2^2 * r3^2 * q1^2 * q2^2)/(r1^2 * r2^2 * r3^2)
	a23::Float64 = -(2.0 * (- r1^2 * r2^2 * q0^3 * q1 - r1^2 * r2^2 * q0^2 * q2 * q3 + r1^2 * r2^2 * q0 * q1^3 + r1^2 * r2^2 * q0 * q1 * q2^2 - r1^2 * r2^2 * q0 * q1 * q3^2 + r1^2 * r2^2 * q1^2 * q2 * q3 + r1^2 * r2^2 * q2^3 * q3 - r1^2 * r2^2 * q2 * q3^3 + r1^2 * r3^2 * q0^3 * q1 - r1^2 * r3^2 * q0^2 * q2 * q3 - r1^2 * r3^2 * q0 * q1^3 + r1^2 * r3^2 * q0 * q1 * q2^2 - r1^2 * r3^2 * q0 * q1 * q3^2 + r1^2 * r3^2 * q1^2 * q2 * q3 - r1^2 * r3^2 * q2^3 * q3 + r1^2 * r3^2 * q2 * q3^3 + 2.0 * r2^2 * r3^2 * q0^2 * q2 * q3 - 2.0 * r2^2 * r3^2 * q0 * q1 * q2^2 + 2.0 * r2^2 * r3^2 * q0 * q1 * q3^2 - 2.0 * r2^2 * r3^2 * q1^2 * q2 * q3))/(r1^2 * r2^2 * r3^2)
    a31::Float64 = (2.0 * (- r1^2 * r2^2 * q0^3 * q2 + r1^2 * r2^2 * q0^2 * q1 * q3 + r1^2 * r2^2 * q0 * q1^2 * q2 + r1^2 * r2^2 * q0 * q2^3 - r1^2 * r2^2 * q0 * q2 * q3^2 - r1^2 * r2^2 * q1^3 * q3 - r1^2 * r2^2 * q1 * q2^2 * q3 + r1^2 * r2^2 * q1 * q3^3 - 2.0 * r1^2 * r3^2 * q0^2 * q1 * q3 - 2.0 * r1^2 * r3^2 * q0 * q1^2 * q2 + 2.0 * r1^2 * r3^2 * q0 * q2 * q3^2 + 2.0 * r1^2 * r3^2 * q1 * q2^2 * q3 + r2^2 * r3^2 * q0^3 * q2 + r2^2 * r3^2 * q0^2 * q1 * q3 + r2^2 * r3^2 * q0 * q1^2 * q2 - r2^2 * r3^2 * q0 * q2^3 - r2^2 * r3^2 * q0 * q2 * q3^2 + r2^2 * r3^2 * q1^3 * q3 - r2^2 * r3^2 * q1 * q2^2 * q3 - r2^2 * r3^2 * q1 * q3^3))/(r1^2 * r2^2 * r3^2)
	a32::Float64 = -(2.0 * (- r1^2 * r2^2 * q0^3 * q1 - r1^2 * r2^2 * q0^2 * q2 * q3 + r1^2 * r2^2 * q0 * q1^3 + r1^2 * r2^2 * q0 * q1 * q2^2 - r1^2 * r2^2 * q0 * q1 * q3^2 + r1^2 * r2^2 * q1^2 * q2 * q3 + r1^2 * r2^2 * q2^3 * q3 - r1^2 * r2^2 * q2 * q3^3 + r1^2 * r3^2 * q0^3 * q1 - r1^2 * r3^2 * q0^2 * q2 * q3 - r1^2 * r3^2 * q0 * q1^3 + r1^2 * r3^2 * q0 * q1 * q2^2 - r1^2 * r3^2 * q0 * q1 * q3^2 + r1^2 * r3^2 * q1^2 * q2 * q3 - r1^2 * r3^2 * q2^3 * q3 + r1^2 * r3^2 * q2 * q3^3 + 2.0 * r2^2 * r3^2 * q0^2 * q2 * q3 - 2.0 * r2^2 * r3^2 * q0 * q1 * q2^2 + 2.0 * r2^2 * r3^2 * q0 * q1 * q3^2 - 2.0 * r2^2 * r3^2 * q1^2 * q2 * q3))/(r1^2 * r2^2 * r3^2)
	a33::Float64 = (r1^2 * r2^2 * q0^4 - 2.0 * r1^2 * r2^2 * q0^2 * q1^2 - 2.0 * r1^2 * r2^2 * q0^2 * q2^2 + 2.0 * r1^2 * r2^2 * q0^2 * q3^2 + r1^2 * r2^2 * q1^4 + 2.0 * r1^2 * r2^2 * q1^2 * q2^2 - 2.0 * r1^2 * r2^2 * q1^2 * q3^2 + r1^2 * r2^2 * q2^4 - 2.0 * r1^2 * r2^2 * q2^2 * q3^2 + r1^2 * r2^2 * q3^4 + 4.0 * r1^2 * r3^2 * q0^2 * q1^2 - 8.0 * r1^2 * r3^2 * q0 * q1 * q2 * q3 + 4.0 * r1^2 * r3^2 * q2^2 * q3^2 + 4.0 * r2^2 * r3^2 * q0^2 * q2^2 + 8.0 * r2^2 * r3^2 * q0 * q1 * q2 * q3 + 4.0 * r2^2 * r3^2 * q1^2 * q3^2)/(r1^2 * r2^2 * r3^2)

	return (a11, a12, a13, a21, a22, a23, a31, a32, a33)
end


text:
"""

    plot_posterior_check(y, yrep, point_est, prob, plot_legend, args...; kwargs...)

Draw a density plot of a random sample of draws from the response variable posterior distribution against the density estimation of the actual data to visualise model fit.

Usage:
```julia-repl
plot_posterior_check(y, yrep, point_est, prob, plot_legend)
```

Details:

The `yrep` matrix should have individual draws from the posterior distribution in rows and unique values (matching the length of the `y` vector in length) in columns.

Arguments:

- `y` : The vector of response variable values.
- `yrep` : The Draws x Values matrix of posterior predictions.
- `point_est` : The type of point estimate to use.
- `prob` : The probability of the credible interval to calculate.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_posterior_check(y::Array, yrep, point_est::String = "median", prob::Float64 = 0.95, plot_legend::Bool = true, args...; kwargs...)

    # Check point estimate argument

    (point_est == "mean" || point_est == "median") || error("`point_est` should be a String of either 'mean' or 'median'.")

    # Check prob argument

    prob > 0 || error("`prob` should be a single Float64 value of 0 < prob < 1.")
    prob < 1 || error("`prob` should be a single Float64 value of 0 < prob < 1.")
    quantileRange = generatequantile(prob)

    # Check object sizes

    size(y, 2) == 1 || error("`y` should have a dimension of N x 1.")
    length1 = size(y, 1)

    if isa(yrep, Array)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    elseif isa(yrep, DataFrame)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    else
        error("`yrep` should be an object of class `Array` or `DataFrame` containing Draws x Values matrix information.")
    end

    length1 == length2 || error("Number of columns in `yrep` should match the length of vector `y`.")

    #-------- Draw plot --------

    gr() # gr backend for graphics
    actualcolour = cgrad(:blues)[1]
    drawscolour = cgrad(:blues)[2]
    Random.seed!(123) # Fix seed for reproducibility
    myPlot = plot()

    # Check if yrep matrix is all integers or numerics
    # to infer discrete/continuous model for plot

    if all(isa.(yrep, Integer)) && all(isa.(y, Integer))
        plotType = "discrete"
    else
        plotType = "continuous"
    end

    # Wrangle from wide to long

    if isa(yrep, Array)
        tmp = DataFrame(yrep)
    else
        tmp = yrep
    end

    #----------
    # Draw plot
    #----------

    if plotType == "discrete"

        #-----------------------
        # Convert to proportions
        #-----------------------

        # y vector

        a = countmap(y)
        b = DataFrame(hcat([[key, val] for (key, val) in a]...)')
        b = rename(b, :x1 => :value)
        b = rename(b, :x2 => :tally)
        b.props = b.tally / sum(b.tally)

        # yrep array - compute proportions for each simulated draw, median and intervals

        tmp.iteration = rownumber.(eachrow(tmp)) # Add iteration IDs
        tmp = stack(tmp, 1:size(tmp,2)-1)
        tmp = combine(groupby(tmp, [:iteration, :value]), nrow => :count)

        tmp = combine(groupby(tmp, :iteration), :count => (x -> x ./ sum(x)) => :props, :value => :value)

        if point_est == "median"
            tmp = combine(groupby(tmp, [:value]), :props => median => :centre, :props => (x -> quantile(x, quantileRange[1])) => :lower, :props => (x -> quantile(x, quantileRange[2])) => :upper)
        else
            tmp = combine(groupby(tmp, [:value]), :props => mean => :centre, :props => (x -> quantile(x, quantileRange[1])) => :lower, :props => (x -> quantile(x, quantileRange[2])) => :upper)
        end

        # Get margins to extend vertically from point estimate for yerror bars on plot

        tmp.lower_margin = tmp.centre - tmp.lower
        tmp.upper_margin = tmp.upper - tmp.centre

        # Draw plot

        myPlot = plot(b[:,1], b[:,3], seriestype = :bar, fillalpha = 0.8,
                     xlabel = "Value", ylabel = "Proportion",
                     label = "y", fill = actualcolour,
                     title = "Posterior Predictive Check",
                     size = (600, 600), legend = plot_legend)

        # Add draw median and 95% credible intervals

        plot!(tmp[:,1], tmp[:,2], seriestype = :scatter, color = drawscolour,
              yerror = (tmp[:,5], tmp[:,6]),
              label = "yrep", legend = plot_legend,
              markerstrokecolor = drawscolour, markersize = 5)

    else
        # Plot posterior draws

        Random.seed!(123)

        for n in 1:nrows

            tmp2 = DataFrame(tmp[n, :])
            tmp2 = stack(tmp2, 1:length2)

            # Add iteration to the plot

            if n == nrows
                plot!(tmp2[!, :value], linealpha = 0.4, xlabel = "", ylabel = "", label = "yrep",
                seriestype = :density, color = drawscolour, size = (400, 400))
            else
                plot!(tmp2[!, :value], linealpha = 0.4, xlabel = "", ylabel = "", label = "",
                seriestype = :density, color = drawscolour, size = (400, 400))
            end
        end

        # Plot actual data

        plot!(y, linealpha = 1, xlabel = "Value", ylabel = "Density", label = "y",
            seriestype = :density, color = actualcolour, legend = plot_legend,
            title = "Posterior Predictive Check", linewidth = 2)
    end

    return myPlot
end

"""

    plot_hist_check(y, yrep, point_est, plot_legend, args...; kwargs...)

Draw a plot with a binned histogram of posterior-predicted response variable values against a measure of centrality of the actual data to visualise model fit.

Usage:
```julia-repl
plot_hist_check(y, yrep, point_est, plot_legend)
```

Arguments:

- `y` : The vector of response variable values.
- `yrep` : The Draws x Values matrix of posterior predictions.
- `point_est` : The type of point estimate to use.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_hist_check(y::Array, yrep, point_est::String = "median", plot_legend::Bool = true, args...; kwargs...)

    # Check point estimate argument

    (point_est == "mean" || point_est == "median") || error("`point_est` should be a String of either 'mean' or 'median'.")

    # Check object sizes

    size(y, 2) == 1 || error("`y` should have a dimension of N x 1.")
    length1 = size(y, 1)

    if isa(yrep, Array)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    elseif isa(yrep, DataFrame)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    else
        error("`yrep` should be an object of class `Array` or `DataFrame` containing Draws x Values matrix information.")
    end

    length1 == length2 || error("Number of columns in `yrep` should match the length of vector `y`.")

    # Compute median for real data vector

    if point_est == "median"
        m = median(y)
    else
        m = mean(y)
    end

    # Wrangle posterior draws from wide to long and compute point estimate for each value

    if isa(yrep, Array)
        tmp = DataFrame(yrep)
    else
        tmp = yrep
    end

    tmp = stack(tmp, 1:length2)

    if point_est == "median"
        tmp = combine(groupby(tmp, :variable), :value => median)
    else
        tmp = combine(groupby(tmp, :variable), :value => mean)
    end

    # Set up graphics helpers

    drawscolour = cgrad(:blues)[1]
    actualcolour = cgrad(:blues)[2]

    # Draw plot

    myPlot = plot(tmp[!, :value_median], seriestype = :histogram, fillalpha = 0.5,
                  xlabel = "Value", ylabel = "", label = "yrep",
                  color = drawscolour, title = "Posterior Predictive Check", size = (400, 400),
                  legend = plot_legend)

    plot!([m], seriestype = "vline", color = actualcolour, label = "y", linewidth = 2.5)

    return myPlot
end

"""

    plot_ecdf_check(y, yrep, plot_legend, args...; kwargs...)

Draw an empirical cumulative density function plot of a random sample of draws from the response variable posterior distribution against the density estimation of the actual data to visualise model fit.

Usage:
```julia-repl
plot_ecdf_check(y, yrep, plot_legend)
```

Details:

The `yrep` matrix should have individual draws from the posterior distribution in rows and unique values (matching the length of the `y` vector in length) in columns.

Arguments:

- `y` : The vector of response variable values.
- `yrep` : The Draws x Values matrix of posterior predictions.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_ecdf_check(y::Array, yrep, plot_legend::Bool = true, args...; kwargs...)

    # Check object sizes

    size(y, 2) == 1 || error("`y` should have a dimension of N x 1.")
    length1 = size(y, 1)

    if isa(yrep, Array)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    elseif isa(yrep, DataFrame)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    else
        error("`yrep` should be an object of class `Array` or `DataFrame` containing Draws x Values matrix information.")
    end

    length1 == length2 || error("Number of columns in `yrep` should match the length of vector `y`.")

    #-------- Draw plot --------

    gr() # gr backend for graphics
    drawscolour = cgrad(:blues)[1]
    actualcolour = cgrad(:blues)[2]
    Random.seed!(123) # Fix seed for reproducibility
    myPlot = plot()

    # Wrangle from wide to long

    if isa(yrep, Array)
        tmp = DataFrame(yrep)
    else
        tmp = yrep
    end

    # Compute ECDF and plot posterior draws

    for n in 1:nrows

        tmp2 = DataFrame(tmp[n, :])
        tmp2 = stack(tmp2, 1:length2)

        # ECDF calculation

        myECDF = StatsBase.ecdf(tmp2[!, :value])

        # Add iteration to the plot

        if n == nrows
            plot!(myECDF, linealpha = 0.4, xlabel = "", ylabel = "", label = "yrep",
            color = drawscolour, size = (400, 400))
        else
            plot!(myECDF, linealpha = 0.4, xlabel = "", ylabel = "", label = "",
            color = drawscolour, size = (400, 400))
        end
    end

    # Compute ECDF of actual data and add to plot

    yECDF = StatsBase.ecdf(y)

    if plot_legend == true
        plot!(yECDF, linealpha = 1, xlabel = "Value", ylabel = "Density", label = "y",
        color = actualcolour, legend = :bottomright, title = "Posterior Predictive Check",
        linewidth = 1.5)
    else
        plot!(yECDF, linealpha = 1, xlabel = "Value", ylabel = "Density", label = "y",
        color = actualcolour, legend = false, title = "Posterior Predictive Check", linewidth = 1.5)
    end

    return myPlot
end


text:
abstract type InitialConditionChronology end
abstract type InitialConditionType end

"""
    InterProblemChronology()

    Type struct to select an information sharing model between stages that uses results from the most recent stage executed to calculate the initial conditions. This model takes into account solutions from stages defined finer resolutions
"""

struct InterProblemChronology <: InitialConditionChronology end

"""
    InterProblemChronology()

    Type struct to select an information sharing model between stages that uses results from the same recent stage to calculate the initial conditions. This model ignores solutions from stages defined finer resolutions.
"""
struct IntraProblemChronology <: InitialConditionChronology end

######################### Initial Conditions Definitions#####################################
struct DevicePower <: InitialConditionType end
struct DeviceStatus <: InitialConditionType end
struct InitialTimeDurationOn <: InitialConditionType end
struct InitialTimeDurationOff <: InitialConditionType end
struct InitialEnergyLevel <: InitialConditionType end
struct InitialEnergyLevelUp <: InitialConditionType end
struct InitialEnergyLevelDown <: InitialConditionType end
struct AreaControlError <: InitialConditionType end


text:
# Load packages and helper functions
using MPI, Statistics
include("../parallelhelpers.jl")

MPI.Init()
comm = MPI.COMM_WORLD
my_rank = MPI.Comm_rank(comm)
n_procs = MPI.Comm_size(comm)
root = 0

# Load the file names
dataLoc = "../../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Distribute filenames by rank
myfnames = fnames[my_rank+1:n_procs:length(fnames)]

# Calculate the word counts
mycounts =  getcounts.(myfnames)

# Gather all counts onto Process 0

# First we have to serialize the counts so we can send it
msg = MPI.serialize(mycounts)

# Since each msg size is different we have to get the length of each one
counts = Int32.(MPI.Allgather(length(msg),comm))

# Then we can do an Allgather so everyone has all the counts
res = MPI.Allgatherv(msg, counts, comm)

# Deserialize the counts
idx = [0; cumsum(counts)]
res = [res[idx[i]+1:idx[i+1]] for i in 1:length(idx)-1]
allcounts = vcat(MPI.deserialize.(res)...)

# Merge them to find the overall counts
overallcounts = merge(+,allcounts...)

# Calculate the normalized counts
normcount = getnormcount.(mycounts,Ref(overallcounts))

# Print out the top 5 words for each file and their counts
for i=1:length(normcount)
    top5 = sort(collect(normcount[i]), by=x->x[2], rev=true)[1:min(5,length(normcount[i])),:]
    println(top5)
end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
module diffprog

end # module


text:
"""
    curve_length(f::Function,a::Real,b::Real;rtol=1e-13) -> Float64

Calculates the length of the curve given by `f`, where `f` is function of one parameter, on interval <`a`, `b`> of said parameter.

# Arguments
- `f::Function`: function of the curve to be measured,
- `a::Real`: left-bound of the interval to be integrated on,
- `b::Real`: right-bound of the interval to be integrated on,
- `rtol::Real`: default `rtol = 1e-13`, relative tolerance,

# Keywords

# Returns
- `Float64`: length of the curve defined by the `f` on the interval <`a`, `b`>.

# Throws
"""
function curve_length(f::Function,a::Real,b::Real;rtol=1e-13)
    df = x -> ForwardDiff.derivative(f,x)
    dr = x -> sqrt(1+df(x)^2)
    l,_ = quadgk(dr,a,b;rtol=rtol)
    return l
end

"""
    curve_fraction(f::Function,a::Real,b::Real,r::Real) -> Float64

Retunrs the position of the fraction of the length from the interval <`a`, `b`>, where the length of the curve is given by one parameter function `f`.

# Arguments
- `f::Function`: function of the curve,
- `a::Real`: left-bound of the interval to be integrated on,
- `b::Real`: right-bound of the interval to be integrated on,
- `r::Real`: fraction of the length of the curve given by `f` from the left,
- `rtol:Real`: default `rtol = 1e-13`, relative tolerance,

# Keywords

# Returns
- `Float64`: position of the fraction of the length from the interval <`a`, `b`>.

# Throws
"""
function curve_fraction(f::Function,a::Real,b::Real,r::Real;rtol=1e-13)
    df = x -> ForwardDiff.derivative(f,x)
    dr = x -> sqrt(1.0 + df(x)^2)

    l,_ = quadgk(dr,a,b;rtol=rtol)
    lf = l*r

    xo = (a+b)/2
    xn = 0.0
    eps_n = 1.0
    while eps_n > xo*rtol
        l,_ = quadgk(dr,a,xo;rtol=rtol)
        xn = xo - (l - lf)/dr(xo)
        eps_n = abs(xo-xn)
        xo = xn
    end

    return xo
end


text:
# MonoidAlternative
# =================

# combine Exceptions
@test neutral(Exception) == MultipleExceptions()
@test ErrorException("hi") ⊕ ErrorException("ho") == MultipleExceptions(ErrorException("hi"), ErrorException("ho"))

# orelse Try
@test Try(3) ⊘ Try(4) == Try(3)  # take the first non-ErrorException("error")
@test Try(ErrorException("error")) ⊘ Try(4) == Try(4)  # take the first non-ErrorException("error")
@test Try(ErrorException("error")) ⊘ Try(4) == Try(4)  # take the first non-ErrorException("error")
@test (@Try error("error")) ⊘ Try(4) == Try(4)  # take the first non-ErrorException("error")

# combine Try
@test Try("hi") ⊕ Try("ho") == Try("hiho")
@test Try(ErrorException("error")) ⊕ Try(4) == Try(4)
@test Try(ErrorException("error")) ⊕ Try(ErrorException("exception")) == Try(MultipleExceptions(ErrorException("error"), ErrorException("exception")))
@test flip_left_right(Try(ErrorException("error"))) ⊕ flip_left_right(Try(ErrorException("exception"))) == flip_left_right(Try(MultipleExceptions(ErrorException("error"), ErrorException("exception"))))

# FunctorApplicativeMonad
# =======================

@test map(Try(3)) do x
  x*x
end == Try(9)

@test map(Try(ErrorException("error"))) do x
  x*x
end == Try(ErrorException("error"))

@test mapn(Try(3), Try("hi")) do x, y
  "$x, $y"
end == Try("3, hi")

@test mapn(Try(3), Try(ErrorException("error"))) do x, y
  "$x, $y"
end isa Const{<:Exception}

h = @syntax_flatmap begin
  a = Try(3)
  b = Try("hi")
  @pure "$a, $b"
end
@test h == Try("3, hi")

h = @syntax_flatmap begin
  a = Try(3)
  b = Try(ErrorException("error"))
  @pure "$a, $b"
end
@test h == Try(ErrorException("error"))

@test pure(Try, 4) == Try(4)


text:
module QuantumSDPs

import Convex

using LinearAlgebra
import LinearAlgebra.I

if VERSION < v"1.2pre"
    (I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))
end

include("types/QuantumVariable.jl")
include("choi2so.jl")

end # module


text:
function activepowervariables(m::JuMP.AbstractModel, devices::Array{R,1}, time_periods::Int64) where {A <: JumpExpressionMatrix, R <: PowerSystems.RenewableGen}

    on_set = [d.name for d in devices]

    t = 1:time_periods

    p_re = @variable(m, p_re[on_set,t] >= 0.0, start = 0.0)

    return p_re

end

function reactivepowervariables(m::JuMP.AbstractModel, devices::Array{R,1}, time_periods::Int64) where {A <: JumpExpressionMatrix, R <: PowerSystems.RenewableGen}

    on_set = [d.name for d in devices]

    t = 1:time_periods

    q_re = @variable(m, q_re[on_set,t], start = 0.0) # Power output of generators

    return q_re
end


text:
# run: ~/.julia/bin/mpiexecjl -n 4 julia --project diffusion_1D_mpi.jl
using Plots, Printf, MAT
import MPI

# enable plotting by default
if !@isdefined do_save; do_save = true end

# MPI functions
@views function update_halo(A, neighbors_x, comm)
    if neighbors_x[1] >= 0 # MPI_PROC_NULL?
        sendbuf = A[2]
        recvbuf = zeros(1)
        MPI.Send(sendbuf,  neighbors_x[1], 0, comm)
        MPI.Recv!(recvbuf, neighbors_x[1], 1, comm)
        A[1] = recvbuf[1]
    end
    if neighbors_x[2] >= 0 # MPI_PROC_NULL?
        sendbuf = A[end-1]
        recvbuf = zeros(1)
        MPI.Send(sendbuf,  neighbors_x[2], 1, comm)
        MPI.Recv!(recvbuf, neighbors_x[2], 0, comm)
        A[end] = recvbuf[1]
    end
    return
end

@views function diffusion_1D_mpi(; do_save=false)
    # MPI
    MPI.Init()
    dims        = [0]
    comm        = MPI.COMM_WORLD
    nprocs      = MPI.Comm_size(comm)
    MPI.Dims_create!(nprocs, dims)
    comm_cart   = MPI.Cart_create(comm, dims, [0], 1)
    me          = MPI.Comm_rank(comm_cart)
    coords      = MPI.Cart_coords(comm_cart)
    neighbors_x = MPI.Cart_shift(comm_cart, 0, 1)
    if (me==0) println("nprocs=$(nprocs), dims[1]=$(dims[1])") end
    # Physics
    lx    = 10.0
    λ     = 1.0
    nt    = 100
    # Numerics
    nx    = 32                 # local
    nx_g  = dims[1]*(nx-2) + 2 # global
    # Derived numerics
    dx    = lx/nx_g            # global
    dt    = dx^2/λ/2.1
    # Array allocation
    qHx   = zeros(nx-1)
    # Initial condition
    x0    = coords[1]*(nx-2)*dx
    xc    = [x0 + ix*dx - dx/2 - 0.5*lx  for ix=1:nx]
    H     = exp.(.-xc.^2)
    t_tic = 0.0
    # Time loop
    for it = 1:nt
        if (it==11) t_tic = Base.time() end
        qHx        .= -λ*diff(H)/dx
        H[2:end-1] .= H[2:end-1] - dt*diff(qHx)/dx
        update_halo(H, neighbors_x, comm_cart)
    end
    t_toc = Base.time()-t_tic
    if (me==0) @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", t_toc, round((2/1e9*nx*sizeof(lx))/(t_toc/(nt-10)), sigdigits=2)) end
    if do_save file = matopen("$(@__DIR__)/H_$(me).mat", "w"); write(file, "H", Array(H)); close(file) end
    MPI.Finalize()
    return
end

diffusion_1D_mpi(; do_save=do_save)


text:
module WoodburyMatrices

using LinearAlgebra
import LinearAlgebra: det, logdet, logabsdet, ldiv!, mul!, adjoint, transpose, diag
import Base: +, *, \, inv, convert, copy, show, similar, axes, size
using SparseArrays

export AbstractWoodbury, Woodbury, SymWoodbury

abstract type AbstractWoodbury{T} <: Factorization{T} end

safeinv(A) = inv(A)
safeinv(A::SparseMatrixCSC) = safeinv(Matrix(A))

myldiv!(A, B)       = ldiv!(A, B)
myldiv!(dest, A, B) = ldiv!(dest, A, B)
if VERSION <= v"1.4.0-DEV.635"
    myldiv!(A::Diagonal, B)       = (B .= A.diag .\ B)
    myldiv!(dest, A::Diagonal, B) = (dest .= A.diag .\ B)
end

include("woodbury.jl")
include("symwoodbury.jl")
include("sparsefactors.jl")
include("deprecated.jl")

# Traits and algorithms expressible in terms of AbstractWoodbury

size(W::AbstractWoodbury) = size(W.A)
size(W::AbstractWoodbury, d) = size(W.A, d)
axes(W::AbstractWoodbury) = axes(W.A)
axes(W::AbstractWoodbury, d) = _axes(W.A, d)
_axes(F::Factorization, dim::Int) = hasmethod(axes, Tuple{typeof(F), Int}) ? axes(F, dim) : Base.OneTo(size(F, dim))
_axes(A::AbstractArray, dim::Int) = axes(A, dim)

Base.Matrix(W::AbstractWoodbury{T}) where {T} = Matrix{T}(W)
Base.Matrix{T}(W::AbstractWoodbury) where {T} = convert(Matrix{T}, W)
Base.Array(W::AbstractWoodbury) = Matrix(W)
Base.Array{T}(W::AbstractWoodbury) where T = Matrix{T}(W)

convert(::Type{Matrix{T}}, W::AbstractWoodbury) where {T} = convert(Matrix{T}, Matrix(W.A) + W.U*W.C*W.V)

# This is a slow hack, but generally these matrices aren't sparse.
SparseArrays.sparse(W::AbstractWoodbury) = sparse(Matrix(W))

# Multiplication

*(W::AbstractWoodbury, B::AbstractMatrix)=W.A*B + W.U*(W.C*(W.V*B))

function *(W::AbstractWoodbury, x::AbstractVector)
    # A reduced-allocation optimization (using temp storage for the multiplications)
    mul!(W.tmpN1, W.A, x)
    mul!(W.tmpk1, W.V, x)
    mul!(W.tmpk2, W.C, W.tmpk1)
    mul!(W.tmpN2, W.U, W.tmpk2)
    return W.tmpN1 + W.tmpN2
end

# Division

function _ldiv(W::AbstractWoodbury, R::AbstractMatrix)
    AinvR = W.A\R
    return AinvR - W.A\(W.U*(W.Cp*(W.V*AinvR)))
end

\(W::AbstractWoodbury, R::AbstractMatrix) = _ldiv(W, R)
\(W::AbstractWoodbury, D::Diagonal) = _ldiv(W, D)

ldiv!(W::AbstractWoodbury, B::AbstractVector) = ldiv!(B, W, B)

function ldiv!(dest::AbstractVector, W::AbstractWoodbury, B::AbstractVector)
    @noinline throwdmm(W, B) = throw(DimensionMismatch("Vector length $(length(B)) must match matrix size $(size(W,1))"))

    length(B) == size(W, 1) || throwdmm(W, B)
    _ldiv!(dest, W, W.A, B)
    return dest
end

function _ldiv!(dest, W, A::Union{Factorization,Diagonal}, B)
    myldiv!(W.tmpN1, A, B)
    mul!(W.tmpk1, W.V, W.tmpN1)
    mul!(W.tmpk2, W.Cp, W.tmpk1)
    mul!(W.tmpN2, W.U, W.tmpk2)
    myldiv!(A, W.tmpN2)
    for i = 1:length(W.tmpN2)
        @inbounds dest[i] = W.tmpN1[i] - W.tmpN2[i]
    end
    return dest
end
_ldiv!(dest, W, A, B) = _ldiv!(dest, W, lu(A), B)

det(W::AbstractWoodbury) = det(W.A)*det(W.C)/det(W.Cp)
logdet(W::AbstractWoodbury) = logdet(W.A) + logdet(W.C) - logdet(W.Cp)
function logabsdet(W::AbstractWoodbury)
    lad_A = logabsdet(W.A)
    lad_C = logabsdet(W.C)
    lad_Cp = logabsdet(W.Cp)
    lad = lad_A[1] + lad_C[1] - lad_Cp[1]
    s = lad_A[2] * lad_C[2] / lad_Cp[2]
    return lad, s
end

*(W::AbstractWoodbury, α::Real) = α*W
+(M::AbstractMatrix, W::AbstractWoodbury) = W + M

end


text:
using NumberIntervals, Turmeric, StaticArrays

rts = roots(sin, NumberInterval(-5, 5))

rts = roots(sin, NumberInterval(-5, 6), Bisection())
rts = roots(sin, NumberInterval(-5, 6), Newton())

# 2D:
f(x, y) = SVector(x^2 + y^2 - 1, y - x)
f(X) = f(X...)

rts = roots(f, SVector(NumberInterval(-5, 5), NumberInterval(-5, 5)))

# complex:
x = NumberInterval(-5, 6)
Xc = Complex(x, x)
f(z) = z^3 - 1

rts = roots(f, Xc, Bisection())
rts = roots(f, Xc, Newton())
rts = roots(f, Xc)

# From R docs:

# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden

function g(x)
    (x1, x2, x3) = x
    SVector(    x1^2 + x2^2 + x3^2 - 1,
                x1^2 + x3^2 - 0.25,
                x1^2 + x2^2 - 4x3
            )
end

X = NumberInterval(-5, 5)
@time rts = roots(g, SVector(X, X, X))

h(xv) = ((x,y) = xv; SVector(2*x - y - exp(-x), -x + 2*y - exp(-y)))

rts = roots(h, SVector(X, X), Bisection())

# Dennis-Schnabel:
h(xv) = ((x, y) = xv; SVector(x^2 + y^2 - 2, exp(x - 1) + y^3 - 2))

rts = roots(h, SVector(X, X), Bisection())

# Test suites:

# http://folk.uib.no/ssu029/Pdf_file/Testproblems/testprobRheinboldt03.pdf

# http://www.mat.univie.ac.at/~neum/glopt/test.html

# http://titan.princeton.edu/TestProblems/

# http://www-sop.inria.fr/saga/POL/

##  MINPACK benchmarks: https://github.com/JuliaNLSolvers/NLsolve.jl/blob/master/test/minpack.jl

rosenbrock(x, y) = SVector( 1 - x, 10 * (y - x^2) )
rosenbrock(X) = rosenbrock(X...)

rosenbrock(xx) = ( (x, y) = xx; SVector( 1 - x, 1000 * (y - x^2) ) )
X = NumberInterval(-1e5, 1e5)

rts = roots(rosenbrock, SVector(X, X))

# and other files in NLsolve test suite

# Testing unconstrained optimization software, Moré, Garbow, Hillstrom ACM Trans. Math. Soft. 7 (1), 17-41 (1981)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AccountApi <: SwaggerApi
    client::Swagger.Client
end

"""

Checks whether the specified account name is available or taken.
Param: location::String (required)
Param: parameters::CheckNameAvailabilityParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NameAvailabilityInformation
"""
function accountCheckNameAvailability(_api::AccountApi, location::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", NameAvailabilityInformation, "/subscriptions/{subscriptionId}/providers/Microsoft.DataLakeStore/locations/{location}/checkNameAvailability", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Creates the specified Data Lake Store account.
Param: resourceGroupName::String (required)
Param: name::String (required)
Param: parameters::DataLakeStoreAccount (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DataLakeStoreAccount
"""
function accountCreate(_api::AccountApi, resourceGroupName::String, name::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", DataLakeStoreAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{name}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "name", name)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified Data Lake Store account.
Param: resourceGroupName::String (required)
Param: name::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function accountDelete(_api::AccountApi, resourceGroupName::String, name::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{name}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "name", name)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Attempts to enable a user managed Key Vault for encryption of the specified Data Lake Store account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function accountEnableKeyVault(_api::AccountApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/enableKeyVault", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified Data Lake Store account.
Param: resourceGroupName::String (required)
Param: name::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DataLakeStoreAccount
"""
function accountGet(_api::AccountApi, resourceGroupName::String, name::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{name}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "name", name)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists the Data Lake Store accounts within the subscription. The response includes a link to the next page of results, if any.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: top::Int32
Param: skip::Int32
Param: select::String
Param: orderby::String
Param: count::Bool
Return: DataLakeStoreAccountListResult
"""
function accountList(_api::AccountApi, api_version::String, subscriptionId::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    Swagger.validate_param("top", "accountList", :minimum, top, 1, false)

    Swagger.validate_param("skip", "accountList", :minimum, skip, 1, false)

    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreAccountListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.DataLakeStore/accounts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$skip", skip)  # type Int32
    Swagger.set_param(_ctx.query, "\$select", select)  # type String
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "\$count", count)  # type Bool
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists the Data Lake Store accounts within a specific resource group. The response includes a link to the next page of results, if any.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: top::Int32
Param: skip::Int32
Param: select::String
Param: orderby::String
Param: count::Bool
Return: DataLakeStoreAccountListResult
"""
function accountListByResourceGroup(_api::AccountApi, resourceGroupName::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    Swagger.validate_param("top", "accountListByResourceGroup", :minimum, top, 1, false)

    Swagger.validate_param("skip", "accountListByResourceGroup", :minimum, skip, 1, false)

    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreAccountListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$skip", skip)  # type Int32
    Swagger.set_param(_ctx.query, "\$select", select)  # type String
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "\$count", count)  # type Bool
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Updates the specified Data Lake Store account information.
Param: resourceGroupName::String (required)
Param: name::String (required)
Param: parameters::DataLakeStoreAccountUpdateParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DataLakeStoreAccount
"""
function accountUpdate(_api::AccountApi, resourceGroupName::String, name::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", DataLakeStoreAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{name}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "name", name)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export accountCheckNameAvailability, accountCreate, accountDelete, accountEnableKeyVault, accountGet, accountList, accountListByResourceGroup, accountUpdate