text:
## Mimi UI
using Dates
using CSVFiles

function dataframe_or_scalar(m::Model, comp_name::Symbol, item_name::Symbol)
    dims = dim_names(m, comp_name, item_name)
    return length(dims) > 0 ? getdataframe(m, comp_name, item_name) : m[comp_name, item_name]
end

##
## 1. Generate the VegaLite spec for a variable or parameter
##

# Get spec
function _spec_for_item(m::Model, comp_name::Symbol, item_name::Symbol; interactive::Bool=true)
    dims = dim_names(m, comp_name, item_name)
    if length(dims) > 2
        # Drop references to singleton dimensions
        dims = tuple([dim for dim in dims if dim_count(m, dim) != 1]...)
    end

    # Control flow logic selects the correct plot type based on dimensions
    # and dataframe fields
    if length(dims) == 0
        value = m[comp_name, item_name]
        name = "$comp_name : $item_name = $value"
        spec = createspec_singlevalue(name)
    elseif length(dims) > 2
        @warn("$comp_name.$item_name has > 2 indexed dimensions, not yet implemented in explorer")
        return nothing
    else
        name = "$comp_name : $item_name"
        df = getdataframe(m, comp_name, item_name)
        dffields = map(string, names(df))         # convert to string once before creating specs

        # a 'time' field necessitates a line plot
        if "time" in dffields

            # need to reorder the df to have 'time' as the first dimension
            ti = findfirst(isequal("time"), dffields)
            if ti != 1
                fields1, fields2 = dffields[1:ti-1], dffields[ti+1:end]
                dffields = ["time", fields1..., fields2...]
                df = df[:, [Symbol(name) for name in dffields]]
            end

            if length(dffields) > 2
                spec = createspec_multilineplot(name, df, dffields, dims, interactive=interactive)
            else
                spec = createspec_lineplot(name, df, dffields, interactive=interactive)
            end

        #otherwise we are dealing with a barplot
        else
            spec = createspec_barplot(name, df, dffields)
        end
    end

    return spec

end

function _spec_for_sim_item(sim_inst::SimulationInstance, comp_name::Symbol, item_name::Symbol, results::DataFrame; model_index::Int = 1, interactive::Bool=true)

    # Control flow logic selects the correct plot type based on dimensions
    # and dataframe fields
    m = sim_inst.models[model_index]
    dims = dim_names(m, comp_name, item_name)
    if length(dims) > 2
        # Drop references to singleton dimensions
        dims = tuple([dim for dim in dims if dim_count(m, dim) != 1]...)
    end

    dffields = map(string, names(results))         # convert to string once before creating specs

    name = "$comp_name : $item_name"

    if length(dims) == 0 # histogram
        spec = createspec_histogram(name, results, dffields; interactive = interactive)
    elseif length(dims) > 2
        @warn("$name has > 2 indexed dimensions, not yet implemented in explorer")
        return nothing
    else

        # check if there are too many dimensions to map and if so, error
        if length(dffields) > 4
            error()

        # a 'time' field necessitates a trumpet plot
        elseif dffields[1] == "time"
            if length(dffields) > 3
                spec =createspec_multitrumpet(name, results, dffields; interactive = interactive)
            else
                spec = createspec_trumpet(name, results, dffields; interactive = interactive)
            end
        #otherwise we are dealing with layered histograms
        else
            spec = createspec_multihistogram(name,results, dffields; interactive = interactive)
        end
    end

    return spec

end

# Create the list of variables and parameters
function menu_item_list(model::Model)
    all_menuitems = []

    for comp_name in map(nameof, compdefs(model))
        items = vcat(variable_names(model, comp_name), parameter_names(model, comp_name))

        for item_name in items
            menu_item = _menu_item(model, comp_name, item_name)
            if menu_item !== nothing
                push!(all_menuitems, menu_item)
            end
        end
    end

    # Return sorted list so that the UI list of items will be in alphabetical order
    return sort(all_menuitems, by = x -> lowercase(x["name"]))
end

function menu_item_list(sim_inst::SimulationInstance)
    all_menuitems = []
    for datum_key in sim_inst.sim_def.savelist

        menu_item = _menu_item(sim_inst, datum_key)
        if menu_item !== nothing
            push!(all_menuitems, menu_item)
        end
    end

    # Return sorted list so that the UI list of items will be in alphabetical order
    return sort(all_menuitems, by = x -> lowercase(x["name"]))
end

function _menu_item(m::Model, comp_name::Symbol, item_name::Symbol)
    dims = dim_names(m, comp_name, item_name)
    if length(dims) > 2
        # Drop references to singleton dimensions
        dims = tuple([dim for dim in dims if dim_count(m, dim) != 1]...)
    end

    if length(dims) == 0
        value = m[comp_name, item_name]
        name = "$comp_name : $item_name = $value"
    elseif length(dims) > 2
        @warn("$comp_name.$item_name has > 2 indexed dimensions, not yet implemented in explorer")
        return nothing
    else
        name = "$comp_name : $item_name"          # the name is needed for the list label
    end

    menu_item = Dict("name" => name, "comp_name" => comp_name, "item_name" => item_name)
    return menu_item
end

function _menu_item(sim_inst::SimulationInstance, datum_key::Tuple{Symbol, Symbol})
    (comp_name, item_name) = datum_key
    dims = dim_names(sim_inst.models[1], comp_name, item_name)
    if length(dims) > 2
        # Drop references to singleton dimensions
        dims = tuple([dim for dim in dims if dim_count(m, dim) != 1]...)
    end

    if length(dims) > 2
        @warn("$comp_name.$item_name has >2 graphing dims, not yet implemented in explorer")
        return nothing
    else
        name = "$comp_name : $item_name"          # the name is needed for the list label
    end

    menu_item = Dict("name" => "$item_name", "comp_name" => comp_name, "item_name" => item_name)
    return menu_item
end

##
## 2. Create individual specs
##

# Specs for explore(m::Model)
function createspec_lineplot(name, df, dffields; interactive::Bool=true)
    interactive ? createspec_lineplot_interactive(name, df, dffields) : createspec_lineplot_static(name, df, dffields)
end

function createspec_lineplot_interactive(name, df, dffields)
    datapart = getdatapart(df, dffields, :line) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "type" => "line",
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => "$name (use bottom plot for interactive selection)",
            "data"=> Dict("values" => datapart),
            "vconcat" => [
                Dict(
                    # "transform" => [Dict("filter" => Dict("selection" => "brush"))],
                    "width" => _plot_width,
                    "height" => _plot_height,
                    "mark" => Dict("type" => "line", "point" => true),
                    "encoding" => Dict(
                        "x" => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear",
                            "axis" => Dict("title"=> ""),
                            "scale" => Dict("domain" => Dict("selection" => "brush"))
                        ),
                        "y" => Dict(
                            "field" => dffields[2],
                            "type" => "quantitative"
                        )
                    )
                ), Dict(
                    "width" => _plot_width,
                    "height" => _slider_height,
                    "mark" => Dict("type" => "line", "point" => true),
                    "selection" => Dict("brush" => Dict("type" => "interval", "encodings" => ["x"])),
                    "encoding" => Dict(
                        "x" => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear"
                        ),
                        "y" => Dict(
                            "field" => dffields[2],
                            "type" => "quantitative",
                            "axis" => Dict("tickCount" => 3, "grid" => false
                            )
                        )
                    )
                )
            ]
        )
    )
    return spec
end

function createspec_lineplot_static(name, df, dffields)
    datapart = getdatapart(df, dffields, :line) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "type" => "line",
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data"=> Dict("values" => datapart),
            "mark" => Dict("type" => "line"),
            "encoding" => Dict(
                "x" => Dict(
                    "field" => dffields[1],
                    "type" => "temporal",
                    "timeUnit" => "utcyear"
                ),
                "y" => Dict(
                    "field" => dffields[2],
                    "type" => "quantitative"
                )
            ),
            "width" => _plot_width,
            "height" => _plot_height
        )
    )
    return spec
end

function createspec_multilineplot(name, df, dffields, multidims; interactive::Bool=true)
    strmultidims = [String(dim) for dim in multidims]
    interactive ? createspec_multilineplot_interactive(name, df, dffields, strmultidims) : createspec_multilineplot_static(name, df, dffields, strmultidims)
end

function createspec_multilineplot_interactive(name, df, dffields, strmultidims)
    datapart = getdatapart(df, dffields, :multiline) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => "$name (use bottom plot for interactive selection)",
            "data"  => Dict("values" => datapart),
            "vconcat" => [
                Dict(
                    # "transform" => [Dict("filter" => Dict("selection" => "brush"))],
                    "mark" => Dict("type" => "line", "point" => true),
                    "encoding" => Dict(
                        "x"     => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear",
                            "axis" => Dict("title"=> ""),
                            "scale" => Dict("domain" => Dict("selection" => "brush"))
                            ),
                        "y"     => Dict(
                            "field" => dffields[3],
                            "type" => "quantitative"
                            ),
                        "color" => Dict("field" => strmultidims[findfirst(strmultidims .!= "time")],
                                        "type" => "nominal",
                            "scale" => Dict("scheme" => "category20"))
                    ),
                    "width"  => _plot_width,
                    "height" => _plot_height
                ), Dict(
                    "width" => _plot_width,
                    "height" => _slider_height,
                    "mark" => Dict("type" => "line", "point" => true),
                    "selection" => Dict("brush" => Dict("type" => "interval", "encodings" => ["x"])),
                    "encoding" => Dict(
                        "x" => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear"
                        ),
                        "y" => Dict(
                            "field" => dffields[3],
                            "type" => "quantitative",
                            "axis" => Dict("tickCount" => 3, "grid" => false)
                        ),
                        "color" => Dict("field" => strmultidims[findfirst(strmultidims .!= "time")],
                                        "type" => "nominal",
                            "scale" => Dict("scheme" => "category20")
                        )
                    )
                )
            ]
        ),
    )
    return spec
end

function createspec_multilineplot_static(name, df, dffields, strmultidims)
    datapart = getdatapart(df, dffields, :multiline) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data"  => Dict("values" => datapart),

            "mark" => Dict("type" => "line"),
            "encoding" => Dict(
                "x"     => Dict(
                    "field" => dffields[1],
                    "type" => "temporal",
                    "timeUnit" => "utcyear"
                    ),
                "y"     => Dict(
                    "field" => dffields[3],
                    "type" => "quantitative"
                    ),
                "color" => Dict("field" => strmultidims[findfirst(strmultidims .!= "time")], "type" => "nominal",
                    "scale" => Dict("scheme" => "category20"))
            ),
            "width"  => _plot_width,
            "height" => _plot_height
        ),
    )
    return spec
end

function createspec_barplot(name, df, dffields)
    datapart = getdatapart(df, dffields, :bar) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "type" => "bar",
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data"=> Dict("values" => datapart),
            "mark" => "bar",
            "encoding" => Dict(
                "x" => Dict("field" => dffields[1], "type" => "ordinal"),
                "y" => Dict("field" => dffields[2], "type" => "quantitative" )
                ),
            "width"  => _plot_width,
            "height" => _plot_height
        )
    )
    return spec
end

function createspec_singlevalue(name)

    datapart = [];
    spec = Dict(
        "name" => name,
        "type" => "singlevalue",
        "VLspec" => Dict()
    )
    return spec
end

# Specs for explore(sim_inst::SimulationInstance)

function createspec_trumpet(name, df, dffields; interactive::Bool=true)
    df_reduced = trumpet_df_reduce(df, :trumpet) #reduce the dataframe down to only the data needed for max, min, and mean lines
    interactive ? createspec_trumpet_interactive(name, df_reduced, dffields) : createspec_trumpet_static(name, df_reduced, dffields)
end

# https://vega.github.io/vega-lite/examples/layer_line_errorband_ci.html
function createspec_trumpet_static(name, df, dffields)

    datapart = getdatapart(df, dffields, :trumpet) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data"=> Dict("values" => datapart),
            "encoding" => Dict(
                "x"     => Dict(
                    "field" => dffields[1],
                    "type" => "temporal",
                    "timeUnit" => "utcyear",
                )
            ),
            "layer" => [
                Dict(
                    "mark" => "line",
                    "encoding" => Dict(
                        "y" => Dict(
                            "aggregate" => "mean",
                            "field" => dffields[2],
                            "type" => "quantitative"
                        )
                    )
                ),
                Dict(
                    "mark" => "area",
                    "encoding" => Dict(
                        "y" => Dict(
                            "aggregate" => "max",
                            "field" => dffields[2],
                            "type" => "quantitative",
                            "title" => "$(dffields[2])"
                        ),
                        "y2" => Dict(
                            "aggregate" => "min",
                            "field" => dffields[2]
                        ),
                        "opacity" => Dict(
                            "value" => 0.5
                        )
                    )
                )

            ]
        ),
        "width"  => _plot_width,
        "height" => _plot_height
    )
    return spec
end

function createspec_trumpet_interactive(name, df, dffields)

    datapart = getdatapart(df, dffields, :trumpet) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => "$name (use bottom plot for interactive selection)",
            "data"=> Dict("values" => datapart),
            "vconcat" => [
                Dict(
                    "width" => _plot_width,
                    "height" => _plot_height,
                    "encoding" => Dict(
                        "x"     => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear",
                            "scale" => Dict("domain" => Dict("selection" => "brush"))
                        )
                    ),
                    "layer" => [
                        Dict(
                            "mark" => "line",
                            "encoding" => Dict(
                                "y" => Dict(
                                    "aggregate" => "mean",
                                    "field" => dffields[2],
                                    "type" => "quantitative"
                                )
                            )
                        ),
                        Dict(
                            "mark" => "area",
                            "encoding" => Dict(
                                "y" => Dict(
                                    "aggregate" => "max",
                                    "field" => dffields[2],
                                    "type" => "quantitative",
                                    "title" => "$(dffields[2])"
                                ),
                                "y2" => Dict(
                                    "aggregate" => "min",
                                    "field" => dffields[2]
                                ),
                                "opacity" => Dict(
                                    "value" => 0.5
                                )
                            )
                        )

                    ]
                ),
                Dict(
                    "width" => _plot_width,
                    "height" => _slider_height,
                    "encoding" => Dict(
                        "x"     => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear"
                        )
                    ),

                    "layer" => [
                        Dict(
                            "mark" => "line",
                            "encoding" => Dict(
                                "y" => Dict(
                                    "aggregate" => "mean",
                                    "field" => dffields[2],
                                    "type" => "quantitative"
                                )
                            )
                        ),
                        Dict(
                            "mark" => "area",
                            "selection" => Dict("brush" => Dict("type" => "interval", "encodings" => ["x"])),
                            "encoding" => Dict(
                                "y" => Dict(
                                    "aggregate" => "max",
                                    "field" => dffields[2],
                                    "type" => "quantitative",
                                    "title" => "$(dffields[2])"
                                ),
                                "y2" => Dict(
                                    "aggregate" => "min",
                                    "field" => dffields[2]
                                ),
                                "opacity" => Dict(
                                    "value" => 0.5
                                )
                            )
                        )

                    ]
                )
            ]
        )
    )
    return spec
end

function createspec_multitrumpet(name, df, dffields; interactive::Bool = true)
    df_reduced = trumpet_df_reduce(df, :multitrumpet) #reduce the dataframe down to only the data needed for max, min, and mean lines
    interactive ? createspec_multitrumpet_interactive(name, df, dffields) : createspec_multitrumpet_static(name, df, dffields)
end

function createspec_multitrumpet_interactive(name, df, dffields)
    datapart = getdatapart(df, dffields, :multitrumpet) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => "$name (use top plot for interactive selection)",
            "data"=> Dict("values" => datapart),
            "vconcat" => [

            # summary graphic of all mean lines
            Dict(
                "mark" => "line",
                "width"  => _plot_width,
                "height" => _slider_height,
                "selection" => Dict("brush" => Dict("type" => "interval", "encodings" => ["x"])),
                "encoding" => Dict(
                    "x" => Dict(
                        "field" => dffields[1],
                        "type" => "temporal",
                        "timeUnit" => "utcyear"
                    ),
                    "y" => Dict(
                        "aggregate" => "mean",
                        "field" => dffields[3],
                        "type" => "quantitative",
                        "title" => "Mean $(dffields[3])"
                    ),
                    "color" => Dict("field" => dffields[2], "type" => "nominal",
                    "scale" => Dict("scheme" => "category20"))
                )
            ),

            # faceted rows
            Dict(
                "facet" => Dict("row" => Dict("field" => dffields[2], "type" => "nominal")),
                "spec" => Dict(
                        Dict(
                            "width" => _plot_width,
                            "height" => _plot_height / 2,
                            "encoding" => Dict(
                                "x"     => Dict(
                                    "field" => dffields[1],
                                    "type" => "temporal",
                                    "timeUnit" => "utcyear",
                                    "scale" => Dict("domain" => Dict("selection" => "brush"))
                                ),
                                "color" => Dict("field" => dffields[2], "type" => "nominal",
                                "scale" => Dict("scheme" => "category20"))
                            ),
                            "layer" => [
                                Dict(
                                    "mark" => "line",
                                    "encoding" => Dict(
                                        "y" => Dict(
                                            "aggregate" => "mean",
                                            "field" => dffields[3],
                                            "type" => "quantitative"
                                        ),
                                        "color" => Dict("field" => dffields[2], "type" => "nominal",
                                        "scale" => Dict("scheme" => "category20"))
                                    )
                                ),
                                Dict(
                                    "mark" => "area",
                                    "encoding" => Dict(
                                        "y" => Dict(
                                            "aggregate" => "max",
                                            "field" => dffields[3],
                                            "type" => "quantitative",
                                            "title" => "Mean $(dffields[3]) with Max/Min"
                                        ),
                                        "y2" => Dict(
                                            "aggregate" => "min",
                                            "field" => dffields[3]
                                        ),
                                        "opacity" => Dict(
                                            "value" => 0.5
                                        ),
                                        "color" => Dict("field" => dffields[2], "type" => "nominal",
                                        "scale" => Dict("scheme" => "category20"))
                                    )
                                )
                            ]
                        )
                )
            )
            ]
        )
    )
    return spec
end

function createspec_multitrumpet_static(name, df, dffields)
    datapart = getdatapart(df, dffields, :multitrumpet) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data" => Dict("values" => datapart),
            "vconcat" => [

            # summary graphic of all mean lines
            Dict(
                "width"  => _plot_width,
                "height" => _slider_height,
                "mark" => "line",
                "encoding" => Dict(
                    "x" => Dict(
                        "field" => dffields[1],
                        "type" => "temporal",
                        "timeUnit" => "utcyear"
                    ),
                    "y" => Dict(
                        "aggregate" => "mean",
                        "field" => dffields[3],
                        "type" => "quantitative",
                        "title" => "Mean $(dffields[3])"
                    ),
                    "color" => Dict("field" => dffields[2], "type" => "nominal",
                    "scale" => Dict("scheme" => "category20"))
                )
            ),
            # faceted rows
            Dict(
                "facet" => Dict("row" => Dict("field" => dffields[2], "type" => "nominal")),
                "spec" => Dict(
                    "width"  => _plot_width,
                    "height" => _plot_height / 2,
                    "encoding" => Dict(
                        "x"     => Dict(
                            "field" => dffields[1],
                            "type" => "temporal",
                            "timeUnit" => "utcyear"
                        ),
                        "color" => Dict("field" => dffields[2], "type" => "nominal",
                        "scale" => Dict("scheme" => "category20"))
                    ),
                    "layer" => [
                        Dict(
                            "mark" => "line",
                            "encoding" => Dict(
                                "y" => Dict(
                                    "aggregate" => "mean",
                                    "field" => dffields[3],
                                    "type" => "quantitative"
                                ),
                                "color" => Dict("field" => dffields[2], "type" => "nominal",
                                "scale" => Dict("scheme" => "category20"))
                            )
                        ),
                        Dict(
                            "mark" => "area",
                            "encoding" => Dict(
                                "y" => Dict(
                                    "aggregate" => "max",
                                    "field" => dffields[3],
                                    "type" => "quantitative",
                                    "title" => "Mean $(dffields[3]) with Min/Max"
                                ),
                                "y2" => Dict(
                                    "aggregate" => "min",
                                    "field" => dffields[3]
                                ),
                                "opacity" => Dict(
                                    "value" => 0.5
                                ),
                                "color" => Dict("field" => dffields[2], "type" => "nominal",
                                "scale" => Dict("scheme" => "category20"))
                            )
                        )
                    ]
                )
            )
            ]
        )
    )
    return spec
end

function createspec_histogram(name, df, dffields; interactive::Bool = true)
    interactive ? createspec_histogram_interactive(name, df, dffields) : createspec_histogram_static(name, df, dffields)
end

function createspec_histogram_static(name, df, dffields)
    datapart = getdatapart(df, dffields, :histogram) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "type" => "histogram",
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data"=> Dict("values" => datapart),
            "mark" => Dict("type" => "bar"),
            "encoding" => Dict(
                "x" => Dict(
                    "field" => dffields[1],
                    "type" => "quantitative",
                    "bin" => Dict("maxbins" => 15)
                ),
                "y" => Dict(
                    "aggregate" => "count",
                    "type" => "quantitative",
                    "title" => "count"
                )
            ),
            "width" => _plot_width,
            "height" => _plot_height
        )
    )
    return spec
end

function createspec_histogram_interactive(name, df, dffields) # for now the same as static version
    createspec_histogram_static(name, df, dffields)
end

function createspec_multihistogram(name, df, dffields; interactive::Bool = true)
    interactive ? createspec_multihistogram_interactive(name, df, dffields) : createspec_multihistogram_static(name, df, dffields)
end

function createspec_multihistogram_static(name, df, dffields)
    datapart = getdatapart(df, dffields, :multihistogram) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => name,
            "data"  => Dict("values" => datapart),
            "vconcat" => [
                # layered histogram of all
                Dict(
                    "mark" => Dict("type" => "bar"),
                    "width"  => _plot_width,
                    "height" => _slider_height,
                    "encoding" => Dict(
                        "x"     => Dict(
                            "field" => dffields[2],
                            "type" => "quantitative",
                            "bin" => Dict("maxbins" => 15)
                            ),
                        "y"     => Dict(
                            "field" => dffields[3],
                            "aggregate" => "count",
                            "type" => "quantitative",
                            "stack" => nothing
                            ),
                        "color" => Dict("field" => dffields[1], "type" => "nominal",
                            "scale" => Dict("scheme" => "category20")),
                        "opacity" => Dict("value" => 0.7)
                    )
                ),

                # faceted rows
                Dict(
                    "facet" => Dict("row" => Dict("field" => dffields[1], "type" => "nominal")),
                    "spec" => Dict(
                        Dict(
                            "mark" => Dict("type" => "bar"),
                            "width"  => _plot_width,
                            "height" => _plot_height / 2,
                            "encoding" => Dict(
                                "x"     => Dict(
                                    "field" => dffields[2],
                                    "type" => "quantitative",
                                    "bin" => Dict("maxbins" => 15)
                                    ),
                                "y"     => Dict(
                                    "field" => dffields[3],
                                    "aggregate" => "count",
                                    "type" => "quantitative",
                                    "stack" => true
                                    ),
                                "color" => Dict("field" => dffields[1], "type" => "nominal",
                                "scale" => Dict("scheme" => "category20"))
                            )
                        )
                    )
                )
            ]
        )
    )
    return spec
end

function createspec_multihistogram_interactive(name, df, dffields)
    datapart = getdatapart(df, dffields, :multihistogram) #returns JSONtext type
    spec = Dict(
        "name"  => name,
        "VLspec" => Dict(
            "\$schema" => "https://vega.github.io/schema/vega-lite/v3.json",
            "description" => "plot for a specific component variable pair",
            "title" => "$name (use top plot for interactive selection)",
            "data"  => Dict("values" => datapart),
            "vconcat" => [

                # layered histogram of all
                Dict(
                    "mark" => Dict("type" => "bar"),
                    "width"  => _plot_width,
                    "height" => _slider_height,
                    "selection" => Dict("brush" => Dict("type" => "interval", "encodings" => ["x"])),
                    "encoding" => Dict(
                        "x"     => Dict(
                            "field" => dffields[2],
                            "type" => "quantitative",
                            "bin" => Dict("maxbins" => 15)
                        ),
                        "y"     => Dict(
                            "field" => dffields[3],
                            "aggregate" => "count",
                            "type" => "quantitative",
                            "stack" => nothing
                        ),
                        "color" => Dict("field" => dffields[1], "type" => "nominal",
                            "scale" => Dict("scheme" => "category20")),
                        "opacity" => Dict("value" => 0.7)
                    )
                ),

                # faceted rows
                Dict(
                    "facet" => Dict("row" => Dict("field" => dffields[1], "type" => "nominal")),
                    "spec" => Dict(
                        Dict(
                            "mark" => Dict("type" => "bar"),
                            "width"  => _plot_width,
                            "height" => _plot_height / 2,
                            "encoding" => Dict(
                                "x"     => Dict(
                                    "field" => dffields[2],
                                    "type" => "quantitative",
                                    "bin" => Dict("maxbins" => 15),
                                    "scale" => Dict("domain" => Dict("selection" => "brush"))
                                    ),
                                "y"     => Dict(
                                    "field" => dffields[3],
                                    "aggregate" => "count",
                                    "type" => "quantitative"
                                    ),
                                "color" => Dict("field" => dffields[1], "type" => "nominal",
                                "scale" => Dict("scheme" => "category20"))
                            )
                        )
                    )
                )
            ]
        )
    )
    return spec
end

# Plot constants
global const _plot_width  = 450
global const _plot_height = 410
global const _slider_height = 90

##
## Helper functions
##

# Various functions to get the JSONtext of the data
function getdatapart(df, dffields, plottype::Symbol)

    # start the main string
    sb = StringBuilder()
    append!(sb, "[");

    # get the specific string for this type of data
    datasb = StringBuilder()
    numrows = length(df[!, 1]);

    # loop over rows and create a dictionary for each row
    if plottype == :multitrumpet #4D with 3 indices
        cols = (df[!, 1], df[!, 2], df[!, 3], df[!, 4])
        datastring = getdatapart_4d(cols, dffields, numrows, datasb)
    elseif plottype == :multiline || plottype == :trumpet #3D with 2 indices, one of which is time
        cols = (df[!, 1], df[!, 2], df[!, 3])
        datastring = getdatapart_3d_time(cols, dffields, numrows, datasb)
    elseif plottype == :multihistogram #3D with 2 indices, none of which is time
        cols = (df[!, 1], df[!, 2], df[!, 3])
        datastring = getdatapart_3d(cols, dffields, numrows, datasb)
    elseif plottype == :line  #2D with 1 index, one of which is time
        cols = (df[!, 1], df[!, 2])
        datastring = getdatapart_2d_time(cols, dffields, numrows, datasb)
    else # :bar and :histogram
        cols = (df[!, 1], df[!, 2])
        datastring = getdatapart_2d(cols, dffields, numrows, datasb)
    end

    append!(sb, datastring * "]");
    datapart = String(sb)

    return JSON.JSONText(datapart)
end

function getdatapart_4d(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{
            \"" * dffields[1]  * "\":\"" * string(Date(cols[1][i])) * "\",
            \"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\",
            \"" * dffields[3] * "\":\"" * string(cols[3][i]) * "\",
            \"" * dffields[4] * "\":\"" * string(cols[4][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function getdatapart_3d_time(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{\"" * dffields[1]  * "\":\"" * string(Date(cols[1][i]))
            * "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\",\""
            * dffields[3] * "\":\"" * string(cols[3][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function getdatapart_3d(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{\"" * dffields[1]  * "\":\"" * string(cols[1][i])
            * "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\",\""
            * dffields[3] * "\":\"" * string(cols[3][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function getdatapart_2d_time(cols, dffields, numrows, datasb)
    for i = 1:numrows
        append!(datasb, "{\"" * dffields[1]  * "\":\"" * string(Date(cols[1][i]))
            * "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\"}")

        if i != numrows
            append!(datasb, ",")
        end
    end

    return String(datasb)
end

function getdatapart_2d(cols, dffields, numrows, datasb)
    for i = 1:numrows

        append!(datasb, "{\"" * dffields[1] * "\":\"" * string(cols[1][i]) *
            "\",\"" * dffields[2] * "\":\"" * string(cols[2][i]) * "\"}") #end of dictionary

        if i != numrows
            append!(datasb, ",")
        end
    end
    return String(datasb)
end

function trumpet_df_reduce(df, plottype::Symbol)

    if plottype == :trumpet
        col_index = 2
    else
        col_index = 3
    end

    col = names(df)[col_index]
    groupby_keys = []
    for i = 1:length(names(df))
        i != col_index && push!(groupby_keys,  names(df)[i])
    end

    df_new = by(df, groupby_keys, value = col => minimum)
    append!(df_new, by(df, groupby_keys, value = col => maximum))
    append!(df_new, by(df, groupby_keys, value = col => mean))
    rename!(df_new, :value => col)

    if plottype == :trumpet
        reorder_cols = [groupby_keys[1], col, groupby_keys[2:end]...]
    else
        reorder_cols = [groupby_keys[1:2]..., col, groupby_keys[3:end]...]
    end

    return df_new[:, reorder_cols]
end


text:
const DATA_FILE = joinpath(@__DIR__, "dummy.csv")
const NUM_DIMENSIONS, NUM_OBSERVATIONS = size(load_data(DATA_FILE)[1])
const SOLVER = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

experiment = Dict{Symbol, Any}(
    :hash => 1,
    :data_file => DATA_FILE,
    :data_set_name => "Dummy",
    :output_file => joinpath(@__DIR__, "example_pool_qs.json"),
    :log_dir => joinpath(@__DIR__, "example_pool_qs.log"),
    :split_strategy_name => "Sf",
    :initial_pool_strategy_name => "Pu",
    :model => Dict(:type => :SVDDneg,
                   :param => Dict{Symbol, Any}(),
                   :init_strategy => SimpleCombinedStrategy(RuleOfThumbScott(), BoundedTaxErrorEstimate(0.05, 0.02, 0.98))),
    :query_strategy => Dict(:type => :RandomOutlierPQs,
                            :param => Dict{Symbol, Any}()),
    :split_strategy => OneClassActiveLearning.DataSplits(trues(NUM_OBSERVATIONS)),
    :oracle => Dict{Symbol, Any}(
        :type => :PoolOracle,
        :param => Dict{Symbol, Any}()
    ),
    :param => Dict(:num_al_iterations => 10,
                   :solver => SOLVER,
                   :initial_pools => fill(:U, NUM_OBSERVATIONS),
                   :adjust_K => true,
                   :initial_pool_resample_version => 1))

ext: jl


text:
module ADI

using Parameters
using Reexport
using Statistics
using UnPack

@reexport using HCIToolbox

export reconstruct,
       decompose,
       Median,
       PCA,
       GreeDS,
       TPCA,
       NMF,
       SingleSDI,
       DoubleSDI,
       SliceSDI

"""
    ADI.ADIAlgorithm

This abstract type is used for defining ADI algorithms. See the extended help (`??ADIAlgorithm`) for interface details.

# Extended help
## Interface
To extend `ADIAlgorithm` you may implement the following

| function | default | description |
|----------|---------|-------------|
| [`reconstruct`](@ref) | | Subroutine for creating the full reconstructed cube with the PSF |
| `(::ADIAlgorithm)` | subtracts output of `reconstruct`, then derotates and collapses | Subroutine for returning the reduced residual cube |
"""
abstract type ADIAlgorithm end

"""
    reconstruct(::ADIAlgorithm, cube, angles, [cube_ref]; kwargs...)

Reconstruct the PSF approximation for the given algorithm, using `cube_ref` as the reference cube if given.
"""
function reconstruct end

"""
    (::ADIAlgorithm)(cube, angles, [cube_ref]; kwargs...)

Fully process an ADI data cube using [`reconstruct`](@ref) and collapsing the residuals. Keyword arguments will be passed to `HCIToolbox.collapse!`.
"""
function (alg::ADIAlgorithm)(cube::AbstractArray{T,3}, angles, args...; kwargs...) where T
    S = reconstruct(alg, cube, angles, args...; kwargs...)
    residual_cube = cube .- S
    return collapse!(residual_cube, angles; kwargs...)
end

include("median.jl")

"""
    ADI.LinearAlgorithm <: ADI.ADIAlgorithm

This abstract type is used for defining linear ADI algorithms. See the extended help (`??LinearAlgorithm`) for interface details.

# Extended help
## Interface
To extend `LinearAlgorithm` you may implement the following

| function | default | description |
|----------|---------|-------------|
| [`decompose`](@ref) | | Subroutine for fitting the linear basis and coefficients as unrolled matrices |
| [`reconstruct`](@ref) | Computes the inner product of the design matrix and weights from `decompose` | Subroutine for creating the full reconstructed cube with the PSF |
| `(::LinearAlgorithm)` | subtracts output of `reconstruct`, then derotates and collapses | Subroutine for returning the reduced residual cube |
"""
abstract type LinearAlgorithm <: ADIAlgorithm end

"""
    ADI.decompose(::LinearAlgorithm, cube, angles, [cube_ref]; kwargs...)
"""
function decompose end

function reconstruct(alg::LinearAlgorithm, cube::AbstractArray{T,3}, angles, args...; kwargs...) where T
    # assumed sizes are (n, Npx) (n, M)
    basis, weights = decompose(alg, cube, angles, args...; kwargs...)
    return reconstruct(alg, basis, weights; kwargs...)
end

function reconstruct(alg::LinearAlgorithm, basis::AbstractMatrix, weights::AbstractMatrix; kwargs...)
    # assumed sizes are (n, Npx) (n, M)
    return weights * basis |> expand
end
# The core decomposition routines
include("pca.jl")
include("greeds.jl")
include("nmf.jl")
include("sdi.jl")

include("metrics/Metrics.jl")
@reexport using .Metrics

end


text:
export K2KrylovParams

"""
Type to use the K2 formulation with a Krylov method, using the package
[`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
The outer constructor

    K2KrylovParams(; uplo = :L, kmethod = :minres, preconditioner = :Identity,
                   rhs_scale = true, form_mat = false,
                   atol0 = 1.0e-4, rtol0 = 1.0e-4,
                   atol_min = 1.0e-10, rtol_min = 1.0e-10,
                   ρ0 = sqrt(eps()) * 1e5, δ0 = sqrt(eps()) * 1e5,
                   ρ_min = 1e2 * sqrt(eps()), δ_min = 1e2 * sqrt(eps()),
                   memory = 20)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
The available methods are:
- `:minres`
- `:minres_qlp`
- `:symmlq`

The list of available preconditioners for this solver is displayed here: [`RipQP.PreconditionerDataK2`](@ref).
"""
mutable struct K2KrylovParams <: AugmentedParams
  uplo::Symbol
  kmethod::Symbol
  preconditioner::Symbol
  rhs_scale::Bool
  form_mat::Bool
  atol0::Float64
  rtol0::Float64
  atol_min::Float64
  rtol_min::Float64
  ρ0::Float64
  δ0::Float64
  ρ_min::Float64
  δ_min::Float64
  mem::Int
end

function K2KrylovParams(;
  uplo::Symbol = :L,
  kmethod::Symbol = :minres,
  preconditioner::Symbol = :Identity,
  rhs_scale::Bool = true,
  form_mat::Bool = false,
  atol0::T = 1.0e-4,
  rtol0::T = 1.0e-4,
  atol_min::T = 1.0e-10,
  rtol_min::T = 1.0e-10,
  ρ0::T = sqrt(eps()) * 1e5,
  δ0::T = sqrt(eps()) * 1e5,
  ρ_min::T = 1e2 * sqrt(eps()),
  δ_min::T = 1e2 * sqrt(eps()),
  mem::Int = 20,
) where {T <: Real}
  if uplo == :L && form_mat
    error("matrix can only be created if uplo == :U")
  end
  return K2KrylovParams(
    uplo,
    kmethod,
    preconditioner,
    rhs_scale,
    form_mat,
    atol0,
    rtol0,
    atol_min,
    rtol_min,
    ρ0,
    δ0,
    ρ_min,
    δ_min,
    mem,
  )
end

mutable struct MatrixTools{T}
  diag_Q::SparseVector{T, Int} # Q diag
  diagind_K::Vector{Int}
end

mutable struct PreallocatedDataK2Krylov{
  T <: Real,
  S,
  M <: Union{LinearOperator{T}, AbstractMatrix{T}},
  MT <: Union{MatrixTools{T}, Int},
  Pr <: PreconditionerDataK2,
  Ksol <: KrylovSolver,
} <: PreallocatedDataAugmentedKrylov{T, S}
  pdat::Pr
  D::S                                  # temporary top-left diagonal
  rhs::S
  rhs_scale::Bool
  regu::Regularization{T}
  δv::Vector{T}
  K::M # augmented matrix
  mt::MT
  KS::Ksol
  atol::T
  rtol::T
  atol_min::T
  rtol_min::T
end

function opK2prod!(
  res::AbstractVector{T},
  nvar::Int,
  Q::Union{AbstractMatrix{T}, AbstractLinearOperator{T}},
  D::AbstractVector{T},
  A::Union{AbstractMatrix{T}, AbstractLinearOperator{T}},
  δv::AbstractVector{T},
  v::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  @views mul!(res[1:nvar], Q, v[1:nvar], -α, β)
  res[1:nvar] .+= α .* D .* v[1:nvar]
  if uplo == :U
    @views mul!(res[1:nvar], A, v[(nvar + 1):end], α, one(T))
    @views mul!(res[(nvar + 1):end], A', v[1:nvar], α, β)
  else
    @views mul!(res[1:nvar], A', v[(nvar + 1):end], α, one(T))
    @views mul!(res[(nvar + 1):end], A, v[1:nvar], α, β)
  end
  res[(nvar + 1):end] .+= @views (α * δv[1]) .* v[(nvar + 1):end]
end

function PreallocatedData(
  sp::K2KrylovParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}

  # init Regularization values
  D = similar(fd.c, id.nvar)
  if iconf.mode == :mono
    regu = Regularization(T(sp.ρ0), T(sp.δ0), T(sp.ρ_min), T(sp.δ_min), :classic)
    D .= -T(1.0e0) / 2
  else
    regu =
      Regularization(T(sp.ρ0), T(sp.δ0), T(sqrt(eps(T)) * 1e0), T(sqrt(eps(T)) * 1e0), :classic)
    D .= -T(1.0e-2)
  end
  δv = [regu.δ] # put it in a Vector so that we can modify it without modifying opK2prod!
  if sp.form_mat
    diag_Q = get_diag_Q(fd.Q.data.colptr, fd.Q.data.rowval, fd.Q.data.nzval, id.nvar)
    K = Symmetric(create_K2(id, D, fd.Q.data, fd.A, diag_Q, regu), fd.uplo)
    diagind_K = get_diag_sparseCSC(K.data.colptr, id.ncon + id.nvar)
    mt = MatrixTools(diag_Q, diagind_K)
  else
    K = LinearOperator(
      T,
      id.nvar + id.ncon,
      id.nvar + id.ncon,
      true,
      true,
      (res, v, α, β) -> opK2prod!(res, id.nvar, fd.Q, D, fd.A, δv, v, α, β, fd.uplo),
    )
    mt = 0
  end

  rhs = similar(fd.c, id.nvar + id.ncon)
  KS = init_Ksolver(K, rhs, sp)
  pdat = eval(sp.preconditioner)(id, fd, regu, D, K)

  return PreallocatedDataK2Krylov(
    pdat,
    D,
    rhs,
    sp.rhs_scale,
    regu,
    δv,
    K, #K
    mt,
    KS,
    T(sp.atol0),
    T(sp.rtol0),
    T(sp.atol_min),
    T(sp.rtol_min),
  )
end

function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK2Krylov{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}
  pad.rhs .= dd
  if pad.rhs_scale
    rhsNorm = kscale!(pad.rhs)
  end
  ksolve!(pad.KS, pad.K, pad.rhs, pad.pdat.P, verbose = 0, atol = pad.atol, rtol = pad.rtol)
  update_kresiduals_history!(res, pad.K, pad.KS.x, pad.rhs)
  if pad.rhs_scale
    kunscale!(pad.KS.x, rhsNorm)
  end

  dd .= pad.KS.x

  return 0
end

function update_pad!(
  pad::PreallocatedDataK2Krylov{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end

  if pad.atol > pad.atol_min
    pad.atol /= 10
  end
  if pad.rtol > pad.rtol_min
    pad.rtol /= 10
  end

  pad.D .= -pad.regu.ρ
  pad.D[id.ilow] .-= pt.s_l ./ itd.x_m_lvar
  pad.D[id.iupp] .-= pt.s_u ./ itd.uvar_m_x
  pad.δv[1] = pad.regu.δ
  if typeof(pad.K) <: Symmetric{T, SparseMatrixCSC{T, Int}}
    pad.D[pad.mt.diag_Q.nzind] .-= pad.mt.diag_Q.nzval
    pad.K.data.nzval[view(pad.mt.diagind_K, 1:(id.nvar))] = pad.D
    pad.K.data.nzval[view(pad.mt.diagind_K, (id.nvar + 1):(id.ncon + id.nvar))] .= pad.regu.δ
  end

  update_preconditioner!(pad.pdat, pad, itd, pt, id, fd, cnts)

  return 0
end


text:
function getSMulKernel{T <: FloatingPoint}(::Type{T})
    nType = T == Float64 ? "double" : "float"

    return "
        #if defined(cl_khr_fp64)  // Khronos extension available?
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
        #elif defined(cl_amd_fp64)  // AMD extension available?
        #pragma OPENCL EXTENSION cl_amd_fp64 : enable
        #endif

        #define number $nType

        __kernel void smul(
                      const number s,
                      __global const number *a,
                      __global number *out) {

        int i = get_global_id(0);

        out[i] = s * a[i];
    }
"
end

function smulCL!{T <: FloatingPoint}(
    s :: Real, a_buff,
    out_buff,
    d1 :: Int64, d2 :: Int64,
    ctx, queue, program, :: Type{T})

    k = Kernel(program, "smul")

    call(queue, k, d1 * d2, nothing, convert(T, s), a_buff, out_buff)
end


text:
using Bukdu

struct WelcomeController <: ApplicationController
    conn::Conn
end

function index(c::WelcomeController)
    render(JSON, "Hello World")
end

routes() do
    get("/", WelcomeController, index)
end

Bukdu.start(8080)

Base.JLOptions().isinteractive==0 && wait()


text:
## fixme ** should work
# or add_dirs in project
VW_DIR="/home/graham_s/VirtualWorlds/projects/"
push!(LOAD_PATH,joinpath(VW_DIR,"TBComponents.jl","src"))
push!(LOAD_PATH,joinpath(VW_DIR,"DDIMeta.jl","src"))


text:
export  lex_plist,
        NUMBER, STRING, IDENT,
        UNKNOWN, ERROR, EOF,
        BEGIN_TAG, END_TAG, CLOSE_TAG, END_AND_CLOSE_TAG,
        HEXBINARY, LPAREN, RPAREN, LBRACE, RBRACE, COMMA, EQUAL, SEMICOLON

@enum(TokenType,
      NUMBER, STRING, IDENT,                                            # Generic
      UNKNOWN, ERROR, EOF,                                              # Control

      # XML
      BEGIN_TAG,         # <tag
      END_TAG,           # >
      CLOSE_TAG,         # </tag>
      END_AND_CLOSE_TAG, # />
      TEXT,              # Text such as "bar" inside tags: <foo>bar</foo>

      # PList
      HEXBINARY,
      LPAREN = Int('('),
      RPAREN = Int(')'),
      LBRACE = Int('{'),
      RBRACE = Int('}'),
      COMMA  = Int(','),
      EQUAL  = Int('='),
      SEMICOLON = Int(';'))

include("token.jl")
include("lexer.jl")
include("basiclexer.jl")

function lex_plist(input::String)
    lex(input, lex_basic)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Materials, Tensors

dummyel = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
dummyip = first(get_integration_points(dummyel))
mutable struct Simulator
    stresses :: Vector{Vector{Float64}}
    strains :: Vector{Vector{Float64}}
    times :: Vector{Float64}
    material :: M where {M<:AbstractMaterial}
end

function Simulator(material)
    return Simulator([],[],[],material)
end

function initialize!(simulator, strains, times)
    simulator.strains = strains
    simulator.times = times
    return nothing
end

function run!(simulator)
    material = simulator.material
    times = simulator.times
    strains = simulator.strains
    t_n = times[1]
    strain_n = strains[1]
    push!(simulator.stresses, copy(tovoigt(material.variables.stress)))
    for i in 2:length(times)
        strain = strains[i]
        t = times[i]
        dstrain = strain - strain_n
        dt = t - t_n
        material.ddrivers.strain = fromvoigt(SymmetricTensor{2,3}, dstrain)
        material.ddrivers.time = dt
        integrate_material!(material)
        material_postprocess_increment!(material, dummyel, dummyip, t)
        push!(simulator.stresses, copy(tovoigt(material.variables.stress)))
        strain_n = strain
        t_n = t
    end
    return nothing
end


text:
function get_partition_code_dict(dict_forward::Bool=true)
    if dict_forward == true
        return partition_code_dict = Dict("intergenic"=>1,"periexonic"=>2,"exon"=>3)
    else
        return code_partition_dict = Dict(1=>"intergenic", 2=>"periexonic", 3=>"exon")
    end
end

function make_padded_df(position_fasta::String, gff3_path::String, genome_path::String, genome_index_path::String, pad::Integer)
    position_reader = FASTA.Reader(open((position_fasta),"r"))
    genome_reader = open(FASTA.Reader, genome_path, index=genome_index_path)
    scaffold_df = build_scaffold_df(gff3_path)
    position_df = DataFrame(SeqID = String[], Start=Int[], End=Int[], PadSeq = LongSequence[], PadStart=Int[], RelStart=Int[], SeqOffset=Int[])
    scaffold_seq_dict = build_scaffold_seq_dict(genome_path, genome_index_path)

    for entry in position_reader
        scaffold = FASTA.identifier(entry)

        if scaffold != "MT"
            desc_array = split(FASTA.description(entry))
            pos_start = parse(Int, desc_array[2])
            pos_end = parse(Int, desc_array[4])
            scaffold_end = scaffold_df.End[findfirst(isequal(scaffold), scaffold_df.SeqID)]

            pad_start=max(1,pos_start-pad)
            pad_length= pos_start - pad_start
            seq_offset = pad - pad_length
            padded_seq = fetch_sequence(scaffold, scaffold_seq_dict, pad_start, pos_end, '+')

            if !hasambiguity(padded_seq)
                push!(position_df, [scaffold, pos_start, pos_end, padded_seq, pad_start, pad_length, seq_offset])
            end
        end
    end

    close(position_reader)
    close(genome_reader)
    return position_df
end

function add_partition_masks!(position_df::DataFrame, gff3_path::String, perigenic_pad::Integer=500, columns::Tuple{Symbol,Symbol,Symbol}=(:SeqID, :PadSeq, :PadStart))
    partitions=["exon", "periexonic", "intergenic"]
    partition_coords_dict = partition_genome_coordinates(gff3_path, perigenic_pad)
    partitioned_scaffolds = divide_partitions_by_scaffold(partition_coords_dict)
    maskcol = [zeros(Int64,0,0) for i in 1:size(position_df,1)]
    position_df.MaskMatrix=maskcol

    @Threads.threads for entry in eachrow(position_df)
        scaffold = entry[columns[1]]
        maskLength = length(entry[columns[2]])
        seqStart = entry[columns[3]]

        scaffold_coords_dict = Dict{String,DataFrame}()

        for ((partition, part_scaffold), df) in partitioned_scaffolds
            if scaffold == part_scaffold
                scaffold_coords_dict[partition] = df
            end
        end

        entry.MaskMatrix=mask_sequence_by_partition(maskLength, seqStart, scaffold_coords_dict)
    end

end
                #add_partition_masks!() SUBFUNCTIONS
                function divide_partitions_by_scaffold(partition_coords_dict::Dict{String, DataFrame})
                    scaffold_coords_dict = Dict{Tuple{String, String}, DataFrame}()
                    for (partition_id, partition_df) in partition_coords_dict
                        for scaffold_subframe in groupby(partition_df, :SeqID)
                            scaffold_id = scaffold_subframe.SeqID[1]
                            scaffold_df = copy(scaffold_subframe)
                            scaffold_coords_dict[(partition_id, scaffold_id)] = scaffold_df
                        end
                    end
                    return scaffold_coords_dict
                end

                function mask_sequence_by_partition(maskLength::Integer, seqStart::Integer, scaffold_coords_dict::Dict{String, DataFrame})
                    partition_code_dict = get_partition_code_dict()
                    seqMask = zeros(Integer, (maskLength, 2))
                    position = seqStart
                    while position <= seqStart+maskLength
                        position_partition, partition_extent, position_strand = find_position_partition(position, scaffold_coords_dict)

                        partition_code = partition_code_dict[position_partition]
                        mask_position = position - seqStart + 1
                        seqMask[mask_position:min(maskLength,mask_position + partition_extent),1] .= partition_code
                        if position_strand == '+'
                            seqMask[mask_position:min(maskLength,mask_position + partition_extent),2] .= 1
                        elseif position_strand == '-'
                            seqMask[mask_position:min(maskLength,mask_position + partition_extent),2] .= -1
                        else
                            seqMask[mask_position:min(maskLength,mask_position + partition_extent),2] .= 0
                        end

                        position += partition_extent + 1
                    end

                    return seqMask
                end

                function find_position_partition(position::Integer, partition_dict::Dict{String, DataFrame})
                    foundPos = false
                    position_partition_id = ""
                    three_prime_extent = 0
                    sample_strand = 0
                    for (partition_id, partition) in partition_dict
                        hitindex = findfirst(x->x>=position, partition.End)
                        if hitindex != nothing
                            if position >= partition.Start[hitindex]
                                foundPos=true
                                position_partition_id = partition_id
                                three_prime_extent = partition.End[hitindex] - position
                                if partition_id == "exon" || partition_id == "periexonic"
                                    sample_strand = partition.Strand[hitindex]
                                end
                            end
                        end
                    end

                    if foundPos == false
                        throw(DomainError("Position $position not found among partition coordinates!"))
                    else
                        return position_partition_id, three_prime_extent, sample_strand
                    end
                end

#function to partition a genome into coordinate sets of:
#merged exons
#"periexonic" sequences (genes with 5' and 3' boundaries projected -/+perigenic_pad bp, minus exons) - includes promoter elements, introns, 3' elements
#intergenic sequences (everything else)
#given a valid gff3
function partition_genome_coordinates(gff3_path::String, perigenic_pad::Integer=500)
    # construct dataframes of scaffolds and metacoordinates
    scaffold_df = build_scaffold_df(gff3_path)

    #partition genome into intragenic, periexonic, and exonic coordinate sets
    #assemble exonic featureset
    exon_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[], Strand=Char[])
    build_feature_df!(gff3_path, "CDS", "MT", exon_df)

    #project exon coordinates onto the scaffold bitwise, merging overlapping features and returning the merged dataframe
    merged_exon_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[], Strand=Char[])
    @showprogress 1 "Partitioning exons..." for scaffold_subframe in DataFrames.groupby(exon_df, :SeqID) # for each scaffold subframe that has exon features
        scaffold_id  = scaffold_subframe.SeqID[1] #get the scaffold id
        scaffold_bitarray = init_scaffold_bitarray(scaffold_id, scaffold_df, false, stranded=true) #init a stranded bitarray of scaffold length
        project_features_to_bitarray!(scaffold_subframe,scaffold_bitarray) #project features as Trues on bitarray of falses
        merged_subframe = get_feature_df_from_bitarray(scaffold_id, scaffold_bitarray) #get a feature df from the projected bitarray
        append!(merged_exon_df,merged_subframe) #append the merged scaffold df to the overall merged df
    end

    #assemble gene featureset
    gene_df =  DataFrame(SeqID = String[], Start = Integer[], End = Integer[], Strand = Char[])
    build_feature_df!(gff3_path, "gene", "MT", gene_df)

    perigenic_pad > 0 && add_pad_to_coordinates!(gene_df, scaffold_df, perigenic_pad) #if a perigenic pad is specified (to capture promoter/downstream elements etc in the periexonic set), apply it to the gene coords

    #build intergenic coordinate set by subtracting gene features from scaffold bitarrays
    intergenic_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[])
    @showprogress 1 "Partitioning intergenic regions..." for scaffold_subframe in DataFrames.groupby(scaffold_df, :SeqID) # for each scaffold subframe that has exon features
        scaffold_id  = scaffold_subframe.SeqID[1] #get the scaffold id
        scaffold_bitarray = init_scaffold_bitarray(scaffold_id, scaffold_df, true) #init a bitarray of scaffold length
        if any(isequal(scaffold_id),gene_df.SeqID) #if any genes on the scafold
            scaffold_genes=gene_df[findall(isequal(scaffold_id), gene_df.SeqID), :] #get the gene rows by finding the scaffold_id
            subtract_features_from_bitarray!(scaffold_genes,scaffold_bitarray) #subtract the gene positions from the scaffold bitarray
        end
        intragenic_subframe = get_feature_df_from_bitarray(scaffold_id, scaffold_bitarray) #get a feature df from the projected bitarray
        append!(intergenic_df,intragenic_subframe) #append the merged scaffold df to the overall merged df
    end

    #build periexonic set by projecting gene coordinates onto the scaffold bitwise, then subtracting the exons
    periexonic_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[], Strand=Char[])
    @showprogress 1 "Partitioning periexonic regions..." for gene_subframe in DataFrames.groupby(gene_df, :SeqID) # for each scaffold subframe that has exon features
        scaffold_id  = gene_subframe.SeqID[1] #get the scaffold id
        scaffold_bitarray = init_scaffold_bitarray(scaffold_id, scaffold_df, false, stranded=true) #init a bitarray of scaffold length
        project_features_to_bitarray!(gene_subframe,scaffold_bitarray) #project features as Trues on bitarray of falses
        if any(isequal(scaffold_id),merged_exon_df.SeqID)
            scaffold_exons=merged_exon_df[findall(isequal(scaffold_id), merged_exon_df.SeqID), :]
            subtract_features_from_bitarray!(scaffold_exons,scaffold_bitarray)
        end
        periexonic_subframe = get_feature_df_from_bitarray(scaffold_id, scaffold_bitarray) #get a feature df from the projected bitarray
        append!(periexonic_df,periexonic_subframe) #append the merged scaffold df to the overall merged df
    end

    return Dict("exon"=>merged_exon_df, "periexonic"=>periexonic_df, "intergenic"=>intergenic_df)
end

            #partition_genome_coordinates() SUBFUNCTIONS
            #BITARRAY SCAFFOLD REPRESENTATION SUBFUNCTIONS
            function init_scaffold_bitarray(scaffold_id::String, scaffold_df, value::Bool; stranded::Bool=false)
                scaffold_length = scaffold_df.End[findfirst(isequal(scaffold_id), scaffold_df.SeqID)]
                if stranded
                    value ? (return rscaffold_bitarray = trues(scaffold_length,2)) : (return rscaffold_bitarray = falses(scaffold_length,2))
                else
                    value ? (return rscaffold_bitarray = trues(scaffold_length,1)) : (return rscaffold_bitarray = falses(scaffold_length,1))
                end
            end

            function project_features_to_bitarray!(scaffold_feature_sf::SubDataFrame, scaffold_bitarray::BitArray)
                @inbounds for item in eachrow(scaffold_feature_sf)
                    scaffold_bitarray[item.Start:item.End,1] = [true for base in item.Start:item.End]
                    if size(scaffold_bitarray)[2] == 2 #if the bitarray is stranded
                        if item.Strand == '+'
                            scaffold_bitarray[item.Start:item.End,2] = [true for base in item.Start:item.End]
                        end
                    end
                end
            end

            function subtract_features_from_bitarray!(scaffold_feature_sf::DataFrame, scaffold_bitarray::BitArray)
                @inbounds for item in eachrow(scaffold_feature_sf)
                    scaffold_bitarray[item.Start:item.End,1] = [false for base in item.Start:item.End]
                end
            end

            function get_feature_df_from_bitarray(scaffold_id::String, scaffold_bitarray::BitArray)
                size(scaffold_bitarray)[2] == 2 ? scaffold_feature_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[], Strand=Char[]) : scaffold_feature_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[])

                new_feature_start = findnext(view(scaffold_bitarray,:,1), 1)
                while new_feature_start != nothing # while new features are still found on the bitarray
                        if size(scaffold_bitarray)[2] == 2 #if stranded, get strand info
                            scaffold_bitarray[new_feature_start,2] == 1 ? new_feature_strand = '+' : new_feature_strand = '-'
                        end
                        if findnext(!eval,view(scaffold_bitarray,:,1),new_feature_start) != nothing
                            new_feature_end = findnext(!eval,view(scaffold_bitarray,:,1),new_feature_start)-1 #find next false after feature start and subtract 1 for feature end
                        else
                            new_feature_end = size(scaffold_bitarray)[1] #if none is found, the end of the feature is the end of hte scaffold
                        end
                        size(scaffold_bitarray)[2] == 2 ? push!(scaffold_feature_df,[scaffold_id, new_feature_start, new_feature_end, new_feature_strand]) : push!(scaffold_feature_df,[scaffold_id, new_feature_start, new_feature_end]) #push stranded feature info as appropriate
                        new_feature_start = findnext(view(scaffold_bitarray,:,1), new_feature_end+1)
                end
                return scaffold_feature_df
            end


text:
const L_system_instructions = @SVector ['F', # Go forward
                                        '-', # Turn right 90°
                                        '+'  # Turn left 90°
                                       ]
struct L_System
    axiom::Char
    rules::Dict{Char, String}
end

function apply_rules(L::L_System, N::Int64; production_sequence::String = "")
    if production_sequence == ""
        production_sequence = string(L.axiom)
    end
    if N == 0
        return production_sequence
    else
        new_production_sequence = String[]
        for i ∈  production_sequence
            if i ∈  L_system_instructions
                push!(new_production_sequence, string(i))
            else
                push!(new_production_sequence, L.rules[i])
            end
        end
        return apply_rules(L, N-1, production_sequence = join(new_production_sequence))
    end
end

# For rectangular domains, with angles constrained to n(π/2) where n = 0,1,2,...
# θ₀ = 0 ⟹     0°
# θ₀ = 1 ⟹    90°
# θ₀ = 2 ⟹   180°
# θ₀ = 3 ⟹   270°
function generate_points(production_sequence::String, p₀::Point_2D,
                         θ₀::Int64, Δx::Float64, Δy::Float64)
    points = [p₀]
    p = p₀
    θ = θ₀
    for i ∈  production_sequence
        if i == 'F'
            p = p + Point_2D(cos(θ*π/2)*Δx, sin(θ*π/2)*Δy)
            push!(points, p)
        elseif i == '+'
            θ = mod(θ + 1, 4)
        elseif i == '-'
            θ = mod(θ - 1, 4)
        end
    end
    return points
end

function hilbert_curve()
    # Note Hilbert curve returns 2^(n^2-1) points
    return L_System('A', Dict('A' => "+BF-AFA-FB+", 'B' => "-AF+BFB+FA-"))
end

# Generate the vector of points making up the Hilbert curve that fill a bounding box bb
# npoints is the target number of points, where the number of points returned will always
# be greater than or equal to npoints
function hilbert_curve(bb::Rectangle_2D, npoints::Int64)
    L = hilbert_curve()
    # Note Hilbert curve returns 2^(2n) points
    n = ceil(Int64, log(npoints)/(2log(2)))
    npoints_actual = 2^(2n)
    Δx = width(bb)/(sqrt(npoints_actual) - 1)
    Δy = height(bb)/(sqrt(npoints_actual) - 1)
    production_sequence = apply_rules(L, n)
    return generate_points(production_sequence, bb.bl, 0, Δx, Δy)
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
struct Ionosphere <: Tabular end
url(::Ionosphere) = "https://archive.ics.uci.edu/ml/machine-learning-databases/ionosphere/ionosphere.data"
checksum(::Ionosphere) = "46d52186b84e20be52918adb93e8fb9926b34795ff7504c24350ae0616a04bbd"
prep(::Ionosphere) = path -> preprocess(path, Ionosphere())
target(::Ionosphere) = 35
size(::Ionosphere) = (351, 0, 0)
problem(::Ionosphere) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
module ReverseDiffX

export NotTracked

using LinearAlgebra
using ..ReverseDiff
using StaticArrays
using Distributions
using PDMats
using ..ForwardDiff

using Base.Broadcast: BroadcastStyle, ArrayStyle, Broadcasted, broadcasted
using ..ForwardDiff: Dual
using ..ReverseDiff: SpecialInstruction, value, value!, deriv, track, record!,
                     tape, unseed!, @grad, TrackedReal, TrackedVector,
                     TrackedMatrix, TrackedArray
using ..DistributionsAD: DistributionsAD

import SpecialFunctions, NaNMath
import ..DistributionsAD: turing_chol, symm_turing_chol, adapt_randn,
    simplex_logpdf
import Base.Broadcast: materialize
import StatsFuns: logsumexp

const TrackedVecOrMat{V,D} = Union{TrackedVector{V,D},TrackedMatrix{V,D}}
const RDBroadcasted{F, T} = Broadcasted{<:Any, <:Any, F, T}

import Random

import Distributions: logpdf,
                      _logpdf,
                      loglikelihood,
                      Gamma,
                      MvNormal,
                      MvLogNormal,
                      Dirichlet,
                      Wishart,
                      InverseWishart,
                      PoissonBinomial,
                      isprobvec

using ..DistributionsAD: TuringPoissonBinomial,
                         TuringMvNormal,
                         TuringMvLogNormal,
                         TuringWishart,
                         TuringInverseWishart,
                         TuringDirichlet,
                         TuringScalMvNormal,
                         TuringDiagMvNormal,
                         TuringDenseMvNormal,
                         VectorOfMultivariate,
                         FillVectorOfMultivariate

###############
## logsumexp ##
###############

logsumexp(x::TrackedArray; dims=:) = track(logsumexp, x, dims = dims)
@grad function logsumexp(x::AbstractArray; dims)
    x_value = value(x)
    lse = logsumexp(x_value; dims=dims)
    return lse, Δ -> (Δ .* exp.(x_value .- lse),)
end

############
## linalg ##
############

function LinearAlgebra.cholesky(A::Symmetric{<:Any, <:TrackedMatrix}; check=true)
    uplo = A.uplo == 'U' ? (:U) : (:L)
    factors, info = symm_turing_chol(parent(A), check, uplo)
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end
function LinearAlgebra.cholesky(A::TrackedMatrix; check=true)
    factors, info = turing_chol(A, check)
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end

function symm_turing_chol(x::TrackedArray{V,D}, check, uplo) where {V,D}
    tp = tape(x)
    x_value = value(x)
    (factors,info), back = DistributionsAD.symm_turing_chol_back(x_value, check, uplo)
    C = Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
    out = track(C.factors, D, tp)
    record!(tp, SpecialInstruction, symm_turing_chol, (x, check, uplo), out, (back, issuccess(C)))
    return out, C.info
end
function turing_chol(x::TrackedArray{V,D}, check) where {V,D}
    tp = tape(x)
    x_value = value(x)
    (factors,info), back = DistributionsAD.turing_chol_back(x_value, check)
    C = Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
    out = track(C.factors, D, tp)
    record!(tp, SpecialInstruction, turing_chol, (x, check), out, (back, issuccess(C)))
    return out, C.info
end

for f in (:turing_chol, :symm_turing_chol)
    @eval begin
        @noinline function ReverseDiff.special_reverse_exec!(
            instruction::SpecialInstruction{typeof($f)},
        )
            output = instruction.output
            instruction.cache[2] || throw(PosDefException(C.info))
            input = instruction.input
            input_deriv = deriv(input[1])
            P = instruction.cache[1]
            input_deriv .+= P((factors = deriv(output),))[1]
            unseed!(output)
            return nothing
        end
    end
end

@noinline function ReverseDiff.special_forward_exec!(
    instruction::SpecialInstruction{typeof(turing_chol)},
)
    output, input = instruction.output, instruction.input
    factors = turing_chol(value.(input)...)[1]
    value!(output, factors)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(
    instruction::SpecialInstruction{typeof(symm_turing_chol)},
)
    output, input = instruction.output, instruction.input
    factors = symm_turing_chol(value.(input)...)[1]
    value!(output, factors)
    return nothing
end

adapt_randn(rng::Random.AbstractRNG, x::TrackedArray, dims...) = adapt_randn(rng, value(x), dims...)

function PoissonBinomial(p::TrackedArray{<:Real}; check_args=true)
    return TuringPoissonBinomial(p; check_args = check_args)
end

Gamma(α::TrackedReal, θ::Real; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::Real, θ::TrackedReal; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::TrackedReal, θ::TrackedReal; check_args=true) = pgamma(α, θ, check_args = check_args)
pgamma(α, θ; check_args=true) = Gamma(promote(α, θ)..., check_args = check_args)
Gamma(α::T; check_args=true) where {T <: TrackedReal} = Gamma(α, one(T), check_args = check_args)
function Gamma(α::T, θ::T; check_args=true) where {T <: TrackedReal}
    check_args && Distributions.@check_args(Gamma, α > zero(α) && θ > zero(θ))
    return Gamma{T}(α, θ)
end

# Work around to stop TrackedReal of Inf and -Inf from producing NaN in the derivative
function Base.minimum(d::LocationScale{T}) where {T <: TrackedReal}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * minimum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(minimum)(d.ρ))
    end
end
function Base.maximum(d::LocationScale{T}) where {T <: TrackedReal}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * maximum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(maximum)(d.ρ))
    end
end

## MvNormal

for (f, T) in (
    (:_logpdf, :TrackedVector),
    (:logpdf, :TrackedMatrix),
    (:loglikelihood, :TrackedMatrix),
)
    @eval begin
        function $f(d::MvNormal{<:Real,<:PDMats.ScalMat}, x::$T{<:Real})
            return $f(TuringScalMvNormal(d.μ, sqrt(d.Σ.value)), x)
        end
        function $f(d::MvNormal{<:Real,<:PDMats.PDiagMat}, x::$T{<:Real})
            return $f(TuringDiagMvNormal(d.μ, sqrt.(d.Σ.diag)), x)
        end
        function $f(d::MvNormal{<:Real,<:PDMats.PDMat}, x::$T{<:Real})
            return $f(TuringDenseMvNormal(d.μ, d.Σ.chol), x)
        end

        function $f(d::MvLogNormal{<:Real,<:PDMats.ScalMat}, x::$T{<:Real})
            return $f(
                TuringMvLogNormal(TuringScalMvNormal(d.normal.μ, sqrt(d.normal.Σ.value))),
                x,
            )
        end
        function $f(d::MvLogNormal{<:Real,<:PDMats.PDiagMat}, x::$T{<:Real})
            return $f(
                TuringMvLogNormal(TuringDiagMvNormal(d.normal.μ, sqrt.(d.normal.Σ.diag))),
                x,
            )
        end
        function $f(d::MvLogNormal{<:Real,<:PDMats.PDMat}, x::$T{<:Real})
            return $f(
                TuringMvLogNormal(TuringDenseMvNormal(d.normal.μ, d.normal.Σ.chol)),
                x,
            )
        end
    end
end

# zero mean, dense covariance
MvNormal(A::TrackedMatrix) = TuringMvNormal(A)

# zero mean, diagonal covariance
MvNormal(σ::TrackedVector) = TuringMvNormal(σ)

# dense mean, dense covariance
MvNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)

# dense mean, diagonal covariance
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvNormal(m, D)
end

# dense mean, diagonal covariance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvNormal(m, σ)

# dense mean, constant variance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedReal) = TuringMvNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Real) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedReal) = TuringMvNormal(m, σ)

# dense mean, constant variance
function MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvNormal(m, A)
end
function MvNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvNormal(m, A)
end
function MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvNormal(m, A)
end

# zero mean,, constant variance
MvNormal(d::Int, σ::TrackedReal) = TuringMvNormal(d, σ)

# zero mean, dense covariance
MvLogNormal(A::TrackedMatrix) = TuringMvLogNormal(TuringMvNormal(A))

# zero mean, diagonal covariance
MvLogNormal(σ::TrackedVector) = TuringMvLogNormal(TuringMvNormal(σ))

# dense mean, dense covariance
MvLogNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))

# dense mean, diagonal covariance
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:TrackedReal, <:TrackedVector{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end

# dense mean, diagonal covariance
MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))

# dense mean, constant variance
function MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedReal)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::TrackedVector{<:Real}, σ::Real)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedReal)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end

# dense mean, constant variance
function MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end

# zero mean,, constant variance
MvLogNormal(d::Int, σ::TrackedReal) = TuringMvLogNormal(TuringMvNormal(d, σ))

# Dirichlet

Dirichlet(alpha::AbstractVector{<:TrackedReal}) = TuringDirichlet(alpha)
Dirichlet(d::Integer, alpha::TrackedReal) = TuringDirichlet(d, alpha)

function _logpdf(d::Dirichlet, x::AbstractVector{<:TrackedReal})
    return _logpdf(TuringDirichlet(d), x)
end
function logpdf(d::Dirichlet, x::AbstractMatrix{<:TrackedReal})
    return logpdf(TuringDirichlet(d), x)
end
function loglikelihood(d::Dirichlet, x::AbstractMatrix{<:TrackedReal})
    return loglikelihood(TuringDirichlet(d), x)
end

for func_header in [
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::Real, x::AbstractVector{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractVector)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractVector{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractVector{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractVector{<:TrackedReal})),

    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::Real, x::AbstractMatrix{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractMatrix)),
    :(simplex_logpdf(alpha::AbstractVector, lmnB::TrackedReal, x::AbstractMatrix{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::Real, x::AbstractMatrix{<:TrackedReal})),
    :(simplex_logpdf(alpha::AbstractVector{<:TrackedReal}, lmnB::TrackedReal, x::AbstractMatrix{<:TrackedReal})),
]
    @eval $func_header = track(simplex_logpdf, alpha, lmnB, x)
end
@grad function simplex_logpdf(alpha, lmnB, x::AbstractVector)
    _alpha = value(alpha)
    _lmnB = value(lmnB)
    _x = value(x)
    simplex_logpdf(_alpha, _lmnB, _x), Δ -> begin
        (Δ .* log.(_x), -Δ, Δ .* (_alpha .- 1) ./ _x)
    end
end
@grad function simplex_logpdf(alpha, lmnB, x::AbstractMatrix)
    _alpha = value(alpha)
    _lmnB = value(lmnB)
    _x = value(x)
    simplex_logpdf(_alpha, _lmnB, _x), Δ -> begin
        (log.(_x) * Δ, -sum(Δ), ((_alpha .- 1) ./ _x) * Diagonal(Δ))
    end
end

Distributions.Wishart(df::TrackedReal, S::Matrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::AbstractMatrix{<:TrackedReal}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractMatrix{<:TrackedReal}) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::TrackedMatrix) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::TrackedMatrix) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::AbstractPDMat{<:TrackedReal}) = TuringWishart(df, S)
Distributions.Wishart(df::TrackedReal, S::AbstractPDMat{<:TrackedReal}) = TuringWishart(df, S)

Distributions.InverseWishart(df::TrackedReal, S::Matrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractMatrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::AbstractMatrix{<:TrackedReal}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractMatrix{<:TrackedReal}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::TrackedMatrix) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::TrackedMatrix) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::AbstractPDMat{<:TrackedReal}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::TrackedReal, S::AbstractPDMat{<:TrackedReal}) = TuringInverseWishart(df, S)

function _logpdf(d::Wishart, X::TrackedMatrix)
    return _logpdf(TuringWishart(d), X)
end
function logpdf(d::Wishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringWishart(d), X)
end
function loglikelihood(d::Wishart, X::AbstractArray{<:TrackedMatrix})
    return loglikelihood(TuringWishart(d), X)
end

function _logpdf(d::InverseWishart, X::TrackedMatrix)
    return _logpdf(TuringInverseWishart(d), X)
end
function logpdf(d::InverseWishart, X::AbstractArray{<:TrackedMatrix})
    return logpdf(TuringInverseWishart(d), X)
end
function loglikelihood(d::InverseWishart, X::AbstractArray{<:TrackedMatrix})
    return loglikelihood(TuringInverseWishart(d), X)
end

# isprobvec

function isprobvec(p::TrackedArray{<:Real})
    pdata = value(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end
function isprobvec(p::SubArray{<:TrackedReal, 1, <:TrackedArray{<:Real}})
    pdata = value(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end
end


text:
#
# changes from 2021/2 files to remove the "20 uplift" from October 2021
# See:
# Department for Work and Pensions. 2021. ‘EXPLANATORY MEMORANDUM TO THE UNIVERSAL CREDIT (EXTENSION OF CORONAVIRUS MEASURES) REGULATIONS 2021 2021 No. 313’. https://www.legislation.gov.uk/uksi/2021/313/pdfs/uksiem_20210313_en.pdf.
#

sys.uc.taper = 55.0

## FIXME JSA/ESA/IS/HB no uplift ??? Think so.

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


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

name = "iperf"
version = v"3.7.0"

# Collection of sources required to complete build
sources = [
    "https://github.com/esnet/iperf/archive/3.7.tar.gz" =>
    "c349924a777e8f0a70612b765e26b8b94cc4a97cc21a80ed260f65e9823c8fc5",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/iperf-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-profiling
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(p -> !isa(p, Windows), supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libiperf", :libiperf),
    ExecutableProduct("iperf3", :iperf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

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


text:
module VITA2D

export Model, Config, Gradient, Thermo, StochasticMatrix, SGDParams
export half, third
export warmup!, mc, pmc
export SGDParams
export sgd

using Parameters
using DelimitedFiles
using Dates
using Random
using LinearAlgebra
using Statistics
using Distributed

set_zero_subnormals(true)

# @everywhere using Pkg
# @everywhere Pkg.activate("/home/mbeach/VITA2D.jl")
# @everywhere Pkg.instantiate()
# @everywhere Pkg.update()
# @everywhere Pkg.resolve()

include("types.jl")
include("utils.jl")
include("sr.jl")
include("metro.jl")
include("update.jl")
include("mc.jl")
include("parallel.jl")
include("parallel_optimise.jl")

end # module


text:
include("Frame.jl")

"""
ADのバイナリ－物理量の変換公式は
物理量 = ボードのAD変換レンジ*(サンプリングバイナリデータ)/2^(分解能) - ボードのAD変換レンジの半分。
また、PIDの信号を高さ情報として取り込んでいるので、

バイナリデータが大きい　＝　PIDの出力電圧が高い　＝　試料に対して押し込んでいる　＝　高さが低い

という関係になる。
したがって、高さの最大最小は

高さ最値小 = バイナリ最大値から計算した高さ
高さ最大値 = バイナリ最小値から計算した高さ
"""
function binaryToPhysicalQuantity(data, header, chanelType)
    # 電圧データに変換
    cc = header.adRange / header.AdResolution
    adUiniRange = header.adRange / 2

    # 電圧データを高さor位相データに変換するための乗数。ErrorとPhaseの場合はTopo像とは逆符号になる点に注意。
    multiplier = Float64(0)
    if chanelType == "topography"
        multiplier = header.piezoGainZ * header.piezoExtensionZ
    elseif chanelType == "error"
        multiplier = -1.0 * header.sensorSensitivity
    elseif chanelType == "phase"
        multiplier = if header.phaseSensitivity != 0 phaseSensitivity else -1.0 end
    else
        # ここには来ない
        @assert false
    end

    # nm -> angstrom(要議論)
    unitConversion = 10

    for y in 1:header.pixelY, x in 1:header.pixelX
        data[y, x] = (adUiniRange - data[y, x] * cc) * multiplier * unitConversion
    end
end

function readFrameHeader(io::IOStream)
    number       = Int64(read(io, Int32))
    maxData      = Int64(read(io, Int16))
    minData      = Int64(read(io, Int16))
    offsetX      = Int64(read(io, Int16))
    offsetY      = Int64(read(io, Int16))
    tiltX        = Float64(read(io, Float32))
    tiltY        = Float64(read(io, Float32))
    isStimulated = read(io, Bool)
    booked1      = Int64(read(io, Int8))
    booked2      = Int64(read(io, Int16))
    booked3      = Int64(read(io, Int32))
    booked4      = Int64(read(io, Int32))

    return FrameHeader(number, maxData, minData, offsetX, offsetY, tiltX, tiltY, isStimulated, booked1, booked2, booked3, booked4)
end

function readImage(io::IOStream, header, chanelType)
    data = zeros(header.pixelY, header.pixelX)
    # TODO: 平均化回数が1じゃないことがあるらしい(今は不必要？)
    for y in 1:header.pixelY, x in 1:header.pixelX
        data[y, x] = Int64(read(io, Int16))
    end
    binaryToPhysicalQuantity(data, header, chanelType)
    return data
end

function readFrame(io::IOStream, header)
    frameHeader = readFrameHeader(io)

    data = readImage(io, header, header.dataType1ch)

    if header.dataType2ch == "none"
        return Frame(frameHeader, data, nothing)
    end

    subData = readImage(io, header, header.dataType2ch)

    return Frame(frameHeader, data, subData)
end

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
lines = readlines("input.txt")

function syntax_check(lines)
    B = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>') # brackets
    ES = Dict(')' => 3, ']' => 57, '}' => 1197, '>' => 25137) # error scores
    CS = Dict('(' => 1, '[' => 2, '{' => 3, '<' => 4) # completion scores

    lefts = Char[]
    error_score = 0
    completion_score = 0
    completion_scores = Int[]

    for line in lines
        for c in line
            if c in keys(B)
                push!(lefts, c)
            else
                if c == B[pop!(lefts)]
                    continue
                else
                    error_score += ES[c]
                    empty!(lefts)
                    break
                end
            end
        end
        if !isempty(lefts)
            while !isempty(lefts)
                completion_score = 5*completion_score + CS[pop!(lefts)]
            end
            push!(completion_scores, completion_score)
            completion_score -= completion_score
        end
    end
    return error_score, completion_scores
end

error_score, completion_scores = syntax_check(lines)
@info "Syntax error score = $error_score"
@info "Median completion score = $(sort!(completion_scores)[end÷2+1])"