text:
"""
    read_file(file::String)
    read_file(path::Array{String,1}, file<:File; kwargs...)
    read_file(path::Array{String,1}, file::XLSXInput)
    read_file(path::String, x::T) where T<:File
    read_file(editor::T) where T<:Edit
This method to reads .csv mapping files required for editing. These files must be stored in
the `data/coremaps` directory. It returns a .csv file.

# Arguments
- `path::Array{String,1}` or `path::String`: Path to file *location*;
    does not include file name.
- `file::String`: Full path to file, including file name.
- `file<:File`: A SLiDE DataType used to store information about a file. Options include:
    - [`SLiDE.CSVInput`](@ref)
    - [`SLiDE.XLSXInput`](@ref)
    - [`SLiDE.GAMSInput`](@ref)
    - [`SLiDE.SetInput`](@ref)
    - [`SLiDE.DataInput`](@ref)
- `editor<:Edit`: A SLiDE DataType used to store information about an edit to make in a
    DataFrame. Specifically, this function might be called for edit types that include the
    field `file` in reference to
    - [`SLiDE.Group`](@ref)
    - [`SLiDE.Map`](@ref)

# Returns
- `df::DataFrame`: If the input is a csv or xlsx file, this method will return a DataFrame.
- `yml::Dict{Any,Any}`: If the input file is a yaml file, this method will return a
    dictionary. All keys that correspond with SLiDE DataStream DataTypes will be converted
    to (lists of) those types.
"""
function read_file(path::Array{String,1}, file::GAMSInput; remove_notes::Bool=false)
    filepath = joinpath(path..., file.name)
    xf = readlines(filepath)
    df = split_gams(xf; colnames=file.col)
    return df
end

function read_file(path::Array{String,1}, file::CSVInput; remove_notes::Bool=false)
    filepath = joinpath(path..., file.name)
    df = SLiDE._read_csv(filepath)

    if remove_notes && size(df, 1) > 1
        df = SLiDE._remove_header(df, filepath)
        df = SLiDE._remove_footer(df)
        df = SLiDE._remove_empty(df)
    end
    return df
end

function read_file(path::Array{String,1}, file::XLSXInput; remove_notes::Bool=false)
    filepath = joinpath(path..., file.name)
    xf = XLSX.readdata(filepath, file.sheet, file.range)

    # Delete rows containing only missing values.
    xf = xf[[!all(row) for row in eachrow(ismissing.(xf))],:]
    return convert_type(DataFrame, xf)
end

function read_file(path::Array{String,1}, file::DataInput; remove_notes::Bool=false)
    df = read_file(path, convert_type(CSVInput, file))
    df = edit_with(df, Rename.(propertynames(df), file.col))

    (:value in propertynames(df)) && (df[!,:value] .= convert_type.(Float64, df[:,:value]))
    return df
end

function read_file(path::Array{String,1}, file::SetInput; remove_notes::Bool=false)
    filepath = joinpath(path..., file.name)
    df = _read_csv(filepath)
    return sort(df[:,1])
end

function read_file(path::String, file::T; remove_notes::Bool=false) where T <: File
    return read_file([path], file; remove_notes=remove_notes)
end

function read_file(editor::T) where T <: Edit
    DIR = joinpath("data", "coremaps")
    df = read_file(joinpath(DIR, editor.file))
    return df
end

function read_file(file::String; colnames=[])
    file = joinpath(file)

    if occursin(".map", file) | occursin(".set", file)
        # return gams_to_dataframe(readlines(file); colnames=colnames)
        return split_gams(readlines(file); colnames=colnames)
    end

    if occursin(".yml", file) | occursin(".yaml", file)
        y = YAML.load(open(file))
        # Here, we first list all sub-subtypes of DataStream (DataTypes that are used in
        # editing datasource files). Then, we find where they overlap with keys in the
        # dictionary read from the YAML file.
        TYPES = string.([IU.subtypes.(IU.subtypes(DataStream))...; IU.subtypes(CGE)])
        KEYS = intersect(TYPES, collect(keys(y)))
        [y[k] = load_from(datatype(k), y[k]) for k in KEYS]
        return y

    elseif occursin(".csv", file)
        df = _read_csv(file)
        return df
    end
end

"""
"""
function _read_csv(filepath::String; header::Int=1)
    df = CSV.read(filepath, DataFrame,
        silencewarnings=true,
        ignoreemptylines=true,
        comment="#",
        missingstrings=["","\xc9","..."];
        header=header)

    (header > 0 && isempty(df)) && (df = _read_csv(filepath; header=0))

    return df
end

"""
This function removes footer rows. These are identified by rows at the bottom of the sheet
that are empty with the exception of the first column.
"""
function _remove_footer(df::DataFrame)
    N = size(df, 1)
    while all(ismissing.(values(df[N,2:end])))
        N -= 1
    end
    return df[1:N,:]
end

"""
"""
function _remove_header(df::DataFrame, filepath::String)
    cols = propertynames(df)
    missing_cols = .|(occursin.(:Column, cols), occursin.(:missing, cols))

    # If only one column was read, maybe this is because there was a problem?
    if size(df, 2) == 1
        # @warn("Removing header WITH reading ($filepath)")
        xf = DelimitedFiles.readdlm(filepath, ',', Any, '\n')
        HEAD = _find_header(xf)
        df = _read_csv(filepath; header=HEAD)

    elseif sum(missing_cols) > 0
        # @warn("Removing header without reading ($filepath)")
        HEAD = _find_header(df)
        df = rename(df[HEAD + 1:end,:], Pair.(cols, Symbol.(ensurearray(df[HEAD,:]))))
    end
    return df
end

"""
"""
function _find_header(df::DataFrame)
    HEAD = 1
    while all(ismissing.(values(df[HEAD,2:end])))
        HEAD += 1
    end
    return HEAD
end

function _find_header(xf::Array{Any,2})
    HEAD = 1
    while all(length.(xf[HEAD,2:end]) .== 0)
        HEAD += 1
    end
    return HEAD
end

"""
Ensure there are no empty columns (those containing all missing values)
"""
function _remove_empty(df::DataFrame)
    LAST = size(df, 2)
    while eltype(df[:,LAST]) === Missing
        LAST -= 1
    end
    return df[:,1:LAST]
end
# _remove_empty(df::DataFrame) = df[:,eltype.(eachcol(df)) .!== Missing]

"""
    split_gams(xf::Array{String,1}; colnames = false)
This function converts a GAMS map or set to a DataFrame, expanding sets into multiple rows.

# Arguments
- `xf::Array{String,1}`: A list of rows of text from a .map or a .set input file.

# Keywords
- `colnames = []`: A user has the option to specify the column propertynames of the output
    DataFrame. If none are specified, a default of `[missing, missing_1, ...]` will be used,
    consistent with the default column headers for `CSV.read()` if column propertynames are missing.

# Returns
- `df::DataFrame`: A dataframe representation of the GAMS map or set.
"""
function split_gams(x::String)
    matches = collect(eachmatch(r"((?<=\").*(?=\")|\((?>[^()]|(?R))*\)|[\w\d]+)", x))
    return [_expand_set(matches[jj][1]) for jj in 1:length(matches)]
end

function split_gams(xf::Array{String,1}; colnames=[])
    xf = xf[xf .!== ""]                             # blank lines
    xf = xf[match.(r"^\s*SET", xf) .=== nothing]    # set definitions
    xf = xf[match.(r"^\s*\*", xf) .=== nothing]     # commented lines
    length(xf) == 0 && (return nothing)

    data = split_gams.(xf)

    if length(unique(length.(data))) == 1
        COLS = length(data[1])
        cols = isempty(colnames) ? _generate_id(COLS) : colnames
        data = if all(typeof.(data) .== Array{String,1})
            DataFrame(permutedims(hcat(data...)), cols)
        else
            vcat([DataFrame(Pair.(cols, row)) for row in data]...)
        end
    end

    return data
end

"""
"""
function _expand_set(x)
    m = match(r"^\((.*)\)$", x)
    x = (m !== nothing) ? string.(strip.(split(m[1], ","))) : string(x)
    return x
end

"""
    read_from(path::String)
This function reads information as specified by the path argument.

# Arguments
- `path::String` to a directory containing files to read *or* to a yaml file with
    information on what to read and how.

# Returns
- `d::Dict` of file contents.
"""
function read_from(path::String; ext=".csv", run_bash::Bool=false)
    d = if any(occursin.([".yml",".yaml"], path))
        _read_from_yaml(path; run_bash=run_bash)
    elseif isdir(path)
        _read_from_dir(path; ext=ext, run_bash=run_bash)
    else
        @error("Cannot read from $path. Function input must point to an existing directory or yaml file.")
    end
    return d
end

"""
    _read_from_dir(dir::String; kwargs...)
This function reads all of the files from the input directory and returns the contents of
those of the specified extension.

# Arguments
- `dir::String`: Relative path to directory to read

# Keywords
- `ext::String = ".csv"`: File extension to read and return
- `run_bash::Bool = false`: If there's a shell script in `dir`, run it to generate/update
    directory contents before reading the files.

# Returns
- `d::Dict{Symbol,Any}` of file contents where the key references the source file name.
"""
function _read_from_dir(dir::String; ext::String=".csv", run_bash::Bool=false)
    files = readdir(dir)

    # If the path contains both a .gdx file and a bash shell script, assume that the script
    # is there to execute "gdxdump" on the shell files.
    run_bash && (files = _run_bash(dir, files))

    @info("Reading $ext files from $dir.")
    files = Dict(_inp_key(f) => f for f in files if occursin(ext, f))
    d = Dict(k => read_file(joinpath(dir, f)) for (k, f) in files)

    _delete_empty!(d)

    return d
end

"""
    _delete_empty!(d::Dict)
This function removes dictionary entries with empty values.
"""
function _delete_empty!(d::Dict)
    for k in keys(d)
        if isempty(d[k])
            @warn("Removing empty entry with key $k from the dictionary.")
            delete!(d, k)
        end
    end
    return d
end

_delete_empty!(d::Any) = d

"""
    _read_from_yaml(path::String)
"""
function _read_from_yaml(path::String; run_bash::Bool=false)
    y = read_file(path)

    # If the yaml file includes the key "Path", this indicates that the yaml file
    if "Path" in keys(y)
        # Look for shell scripts in this path, and run them if they are there.
        # If none are found, nothing will happen.
        run_bash && _run_bash(joinpath(SLIDE_DIR, ensurearray(y["Path"])...))

        files = ensurearray(values(find_oftype(y, File)))
        inp = "Input" in keys(y) ? y["Input"] : files
        d = _read_from_yaml(y["Path"], inp)
        d = _edit_from_yaml(d, y, inp)
    else
        inp = ensurearray(values(find_oftype(y, CGE)))
        d = _read_from_yaml(inp)
    end
    return d
end

function _read_from_yaml(path::String, files::Dict)
    path = joinpath(SLIDE_DIR, path)
    d = Dict(_inp_key(k) => read_file(joinpath(path, f)) for (k, f) in files)
    return _delete_empty!(d)
end

function _read_from_yaml(path::String, files::Array{T,1}) where {T <: File}
    path = joinpath(SLIDE_DIR, path)
    d = Dict(_inp_key(f) => read_file(path, f) for f in files)
    return _delete_empty!(d)
end

_read_from_yaml(path::Array{String,1}, file::Any) = _read_from_yaml(joinpath(path...), file)
_read_from_yaml(lst::Array{Parameter,1}) = _delete_empty!(Dict(_inp_key(x) => x for x in lst))

"""
"""
_edit_from_yaml(d::Dict, editor::Dict, files::Array) = d

_edit_from_yaml(d::Dict{Symbol,DataFrame}, editor::Dict, files::Dict) = d
# function _edit_from_yaml(d::Dict{Symbol,DataFrame}, editor::Dict, files::Dict)
#     # For when the yaml file only points to files to read and doesn't contain any edits.
#     return Dict(k => edit_with(y) for (k,y) in d)
# end

function _edit_from_yaml(d::Dict{Symbol,Dict{Any,Any}}, editor, files)
    # For when the yaml file read contains a list of yaml files, each containing edits.
    # See: src/readfiles/build/shareinp.yml.
    return Dict(k => edit_with(y) for (k, y) in d)
    # !!!! Maybe use _read_from_yaml again here hahah.
end

function _edit_from_yaml(d::Dict{Symbol,DataFrame}, editor::Dict, files::Array{T,1}) where {T <: File}
    # For when the yaml file defines files AND edits that should be made to ALL data inputs.
    [d[_inp_key(f)] = edit_with(d[_inp_key(f)], editor, f) for f in files]
    return d
end

function _edit_from_yaml(d::Dict{Symbol,DataFrame}, editor::Dict, files::Array{DataInput,1})
    # For when the yaml file defines files AND edits that should be made to ALL data inputs.
    # This is the same as when editing an array of files EXCEPT that we order the columns
    # as specified. See: dev/readfiles/1_data_out.yml.
    [d[_inp_key(f)] = select(edit_with(d[_inp_key(f)], editor, f), f.col) for f in files]
    return d
end

"""
    _run_bash(dir::String, files)
    _run_bash(path::String)
This function runs a shell script if one is found.

# Arguments
- `path::String`: Relative path to a specific shell script to run or to a directory that
    might contain a shell script.
- `file::String` or `files::Array{String,1}`: A list of files in the specified directory.

# Returns
- `files::Array{String,1}`: An updated list of files in the specified directory after
    running the shell script if a list of files is given as an argument.
"""
function _run_bash(path::String, files::Array{String,1})
    scripts = files[occursin.(".sh", files)]
    isempty(scripts) && return

    # Save the current directory so we can return to it later. Enter the directory
    # containing the bash file(s), run it/them, and return to the original directory.
    # We default to iterating over a loop of an array of files, even if that array contains
    # only one file, to minimize changing directories.
    curr_dir = pwd()
    cd(path)
    for s in scripts
        @info("Running bash script $s in $path.")
        run(`bash $s`)
    end
    cd(curr_dir)
    files = readdir(path)
    return readdir(path)
end

function _run_bash(path::String)
    if isdir(path)
        _run_bash(path, readdir(path))
    elseif isfile(path) && (path[end - 2:end] == ".sh")
        dir = joinpath(splitpath(path)[1:end - 1]...)
        file = splitpath(path)[end]
        _run_bash(dir, file)
    end
end

_run_bash(path::String, file::String) = _run_bash(path, ensurearray(file))

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
function parsetext!(wdg::Widgets.AbstractWidget, name = "function"; text = observe(wdg), on = text, parse = Meta.parse, default = (args...) -> nothing)
    f = default
    name = Symbol(name)
    wdg[name] = Observable{Any}(f)
    InteractBase.on(on) do s
        update_function!(wdg[name], text[]; parse = parse)
    end
    wdg
end

function update_function!(func::AbstractObservable, s; parse = Meta.parse)
    try
        sym = gensym()
        @eval Main $sym = $(parse(s))
        func[] = getproperty(Main, sym)
    catch
    end
end

"""
`Undo(obs::Observable{T}, stack = T[obs[]]; stacksize = 10) where {T}`

Return a `Undo` object that stores up to `stacksize` past occurences of `Observable` `obs` in `Vector` `stack`.
(::Undo)() sets the observable `obs` back to its previous state as far as memory permits.
"""
struct Undo{T, F}
    obs::Observable{T}
    stack::Vector{T}
    stacksize::Int
    exclude::F

    function Undo(obs::Observable{T}, stack = T[obs[]]; stacksize = 10) where {T}
        exclude = on(obs) do val
            push!(stack, val)
            length(stack) > stacksize && popfirst!(stack)
        end
        new{T, typeof(exclude)}(obs, stack, stacksize, exclude)
    end
end

Widgets.observe(u::Undo) = u.obs

(u::Undo)() = undo!(u.obs, u.stack, u.exclude)

function undo!(obs::Observable, stack, exclude)
    pop!(stack)
    isempty(stack) && error("Stack is finished, cannot undo any more")
    Observables.setexcludinghandlers(obs, last(stack), t -> t != exclude)
    obs
end

skipnan(v) = Iterators.filter(!isnan,v)

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Augmented _Gradients_ Transect Data
#
# Here we (1) retrieve [SCOPE-Gradients](http://scope.soest.hawaii.edu/data/gradients/data/) cruise data from the [Simons' CMAP](https://cmap.readthedocs.io/en/latest/) data base and (2) collocate with mean dymanic topography / SSH.
#
# <img src="cbiomes-01.png" alt="Drawing" style="height: 150px;"/>

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Set up tools
#
# [PyCmap](https://github.com/simonscmap/pycmap) is the Python API that we will use in Julia, via [PyCall.jl](https://github.com/JuliaPy/PyCall.jl), to query the CMAP data base. [Plots.jl](http://docs.juliaplots.org/latest/) is a common `Julia` plotting package and `helper_functions.jl` adds a few convenience functions.

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# _Pre-requisites:_
#
# - 1. _install the [PyCmap](https://github.com/simonscmap/pycmap) python package and its dependencies using `pip`_
# - 2. _compile [PyCall.jl](https://github.com/simonscmap/pycmap) using external python distribution that installed `PyCmap`_
# - 3. _obtain your own API key from [Simons' CMAP](https://simonscmap.com) (free; takes `<30s`)_
# - 4. _import `pycmap` (via `pycall`), `Plots`, and `helper functions`_
#
# _You may need to replace `your-own-API-key` (as outline below) with your own API key from [Simons' CMAP](https://simonscmap.com) and uncomment the command below._

# + {"cell_style": "split", "slideshow": {"slide_type": "skip"}}
if false
    run(`pip install pycmap`) #pycmap is used via PyCall later
    run(pipeline(`which python`,"whichpython.txt")) #external python path
    ENV["PYTHON"]=readline("whichpython.txt")
    import Pkg; Pkg.build("PyCall")
end

# + {"cell_style": "split", "slideshow": {"slide_type": "skip"}}
using PyCall
PyCmap = pyimport("pycmap")
#cmap = PyCmap.API(token="your-own-API-key")
cmap = PyCmap.API()

using Plots
include("helper_functions.jl")

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Get Data Catalog
#
# _Downloading the catalog is a simple way to verify that cmap is all set-up. The commented `df.to_csv` command writes the content of `df` to a new `catalog.csv` file. Alternatively, `Pandas.jl` can be used as also shown._

# + {"slideshow": {"slide_type": "skip"}}
df = cmap.get_catalog();
#df.to_csv("catalog.csv")

#df=Pandas.DataFrame(cmap.get_catalog());
#to_csv(df,"catalog.csv")

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Download Data (`.csv` Files)
#
# A simple method is to download data from `CMAP` and store it to a `CSV` file which any software can then reload. Below, `cmap_helpers.tables` provide CMAP `table` lists for [SCOPE-Gradients](http://scope.soest.hawaii.edu/data/gradients/data/) cruise data, which `cmap.get_dataset` downloads one at a time. _Alternatively one can use the computer's memory (later slides)._

# + {"slideshow": {"slide_type": "skip"}}
pth="../samples/gradients/"
!isdir("$pth") ? mkdir("$pth") : nothing

Γ=1:3
γ=filter(x -> occursin("tblKM1906",x), readdir(pth))
!isempty(γ) ? Γ = [] : nothing

for g in Γ
    list0=cmap_helpers.tables("G$g")
    for i in list0
        df=cmap.get_dataset(i)
        df.to_csv("$pth$i.csv")
    end
end

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Read Data (& Meta-Data)
#
# As an example below we read, and then plot, the `LISST` data collected during the `Gradients 3` cruise.

# + {"slideshow": {"slide_type": "skip"}}
s=cmap_helpers.get("tblKM1906_Gradients3_uway_optics","LISST_small")
m=cmap_helpers.get("tblKM1906_Gradients3_uway_optics","LISST_medium")
l=cmap_helpers.get("tblKM1906_Gradients3_uway_optics","LISST_large")

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Collocate Ancillary Data
#
# Here we interpolate a `sea surface height` climatology estimate ([Forget et al 2015](http://doi.org/10.5194/gmd-8-3071-2015)) along the ship track. _Any number of commonly available methods can readily be used to interpolate gridded estimates to observed locations. For `MITgcm` output in our example, we use `MeshArrays.jl` method._

# + {"slideshow": {"slide_type": "skip"}}
ssh=cbiomes_helpers.myinterp(pth,"SSH",s["lon"],s["lat"])
ssh=merge(ssh,Dict("Unit" => "m", "Variable" => "SSH", "Long_Name" => "Sea Surface Height (Mean Dynamic Topography)"))

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Plot v Latitude
#
#  LISST data collected in Gradients 3 and the sea surface height climatology.

# + {"slideshow": {"slide_type": "skip"}}
t=1:10:length(s["lat"])
scatter(s["lat"][t],s["val"][t],marker = 1.5,label=s["Long_Name"],
    xlabel="°N",ylabel=s["Unit"], title="Gradients 3 (subset)")
scatter!(m["lat"][t],m["val"][t],marker = 1.5,label=m["Long_Name"])
scatter!(l["lat"][t],l["val"][t],marker = 1.5,label=l["Long_Name"])
plot!(ssh["lat"][t],(1.0.-ssh["val"][t]).*5.0,linecolor=:black,label="(1.- ssh in m) * 5")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Gradients As Function Of Station ID
#
# `Gradients1`, `Gradients2`, and `Gradients3` are plotted one after the other. _Note the strong covariations, in each case, with the ship moving back and forth across the gyre in each case._

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
s1=cmap_helpers.get("tblKOK1606_Gradients1_uway_optics","LISST_small")
m1=cmap_helpers.get("tblKOK1606_Gradients1_uway_optics","LISST_medium")
l1=cmap_helpers.get("tblKOK1606_Gradients1_uway_optics","LISST_large")
ssh1=cbiomes_helpers.myinterp(pth,"SSH",s1["lon"],s1["lat"])

t=1:5:length(s1["lat"])
scatter(s1["val"][t],marker = 2,label=s["Long_Name"],
    ylabel=s["Unit"], title="Gradients 1 (subset)")
scatter!(m1["val"][t],marker = 2,label=m["Long_Name"])
scatter!(l1["val"][t],marker = 2,label=l["Long_Name"])
plot!((1.0.-ssh1["val"][t]).*20.0,linecolor=:black,label="(1.- ssh in m) * 20")
plot!(s1["lat"][t]./10.,linecolor=:black,linestyle = :dash,label="°N / 10")

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
s2=cmap_helpers.get("tblMGL1704_Gradients2_uway_optics","LISST_small")
m2=cmap_helpers.get("tblMGL1704_Gradients2_uway_optics","LISST_medium")
l2=cmap_helpers.get("tblMGL1704_Gradients2_uway_optics","LISST_large")
ssh2=cbiomes_helpers.myinterp(pth,"SSH",s2["lon"],s2["lat"])

t=1:5:length(s2["lat"])
scatter(s2["val"][t],marker = 2,label=s["Long_Name"],
    ylabel=s["Unit"], title="Gradients 2 (subset)")
scatter!(m2["val"][t],marker = 2,label=m["Long_Name"])
scatter!(l2["val"][t],marker = 2,label=l["Long_Name"])
plot!((1.0.-ssh2["val"][t]).*5.0,linecolor=:black,label="(1.- ssh in m) * 5")
plot!(s2["lat"][t]./10.0,linecolor=:black,linestyle = :dash,label="°N / 10")

# + {"slideshow": {"slide_type": "subslide"}}
t=1:10:length(s["lat"])
plot(s["val"][t],marker = 2,label=s["Long_Name"],
    ylabel=s["Unit"], title="Gradients 3 (subset)")
plot!(m["val"][t],marker = 2,label=m["Long_Name"])
plot!(l["val"][t],marker = 2,label=l["Long_Name"])
plot!((1.0.-ssh["val"][t]).*5.0,linecolor=:black,label="(1.- ssh in m) * 5")
plot!(s["lat"][t]./10.0,linecolor=:black,linestyle = :dash,label="(lat in °N) / 5")
# -

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


text:
using Plots
include("analysis_functions.jl")

# plotting shortcut to make MHC and MHC+DOS curves
function plot_comparison(dos_func, eη_max, E_min, E_max, average_dos; kT=.026, λ=.26, length=500, plot_title="")
    eη_range = range(-eη_max, eη_max, length=length)
    MHC_k = [compute_k_MHC(E_min, E_max, λ, eη, average_dos; kT=kT) for eη in eη_range]
    MHC_DOS_k = [compute_k_MHC_DOS(E_min, E_max, λ, eη, dos_func; kT=kT) for eη in eη_range]
    if any(MHC_k.<0) | any(MHC_DOS_k.<0)
        println("negative k...uh-oh...")
        println(sum(MHC_k.<0))
        println(sum(MHC_DOS_k.<0))
    end
    plot(eη_range, hcat(log10.(abs.(MHC_k)), log10.(abs.(MHC_DOS_k))), label=["MHC" "MHC+DOS"], legend=:bottomright, xlabel="eη", ylabel="log(k)", title=plot_title, color=[:purple :green])
end

# if you want to visualize each term in integrand separately
function plot_integrand_components(dos_func, E_min, E_max, average_dos; eη=0, kT=.026, λ=.26, length=500, plot_title="")
    E_range = range(E_min, E_max, length=length)
    dos_vals = dos_func.(E_range)
    p1 = plot(E_range, dos_vals, color=:black, ylabel="electrode states", label="electrons/holes", title=plot_title)

    fd_vals_e = fermi_dirac.(E_range; kT=kT)
    fd_vals_h = 1 .- fd_vals_e
    p2 = plot(E_range, hcat(fd_vals_e, fd_vals_h), label=["electrons" "holes"], ylabel="occupation", legend=:right)

    ox_dos_vals = marcus_integrand.(E_range, λ, eη, true; kT=kT)
    red_dos_vals = marcus_integrand.(E_range, λ, eη, false; kT=kT)
    p3 = plot(E_range, hcat(red_dos_vals, ox_dos_vals), label=["reduction" "oxidation"], ylabel="electrolyte states")

    plot(p1, p2, p3, layout=(3,1))
end

function plot_fits(exp_data, dos_f, avg_dos, E_min, E_max, MHC_λ, MHC_A, MHC_DOS_λ, MHC_DOS_A; plot_title="")
    # compute k's
    V = exp_data[:,1]
    V_mag = 1.1 * maximum(abs.(V))
    V_range = range(-V_mag, V_mag, length=200)
    MHC_k = [MHC_A*compute_k_MHC(E_min, E_max, MHC_λ, V, avg_dos) for V in V_range]
    MHC_DOS_k = [MHC_DOS_A*compute_k_MHC_DOS(E_min, E_max, MHC_DOS_λ, V, dos_f) for V in V_range]

    # scatter plot of experimental data, lines for fits
    xs = Vector[V, V_range, V_range]
    ys = Vector[exp_data[:,2], MHC_k, MHC_DOS_k]
    plot(xs, ys, seriestype=[:scatter :line :line ], label=["experiment" string("MHC: λ=",round(MHC_λ, digits=3), "; A=", round(MHC_A, digits=2)) string("MHC+DOS: λ=", round(MHC_DOS_λ, digits=3), "; MHC_DOS_A=", round(MHC_DOS_A, digits=2))], xlabel="V", ylabel="log(k or I)", yscale=:log10, leg=:bottomright, title=plot_title)
end


text:
module KMeans

using Query, DataFrames, Clustering

export call

function score(centers, X)
  dmat = Distances.pairwise(Distances.SqEuclidean(), centers, X)
  sum(findmin(dmat; dims=1)[1])
end

function kmeans_internal(values, X, cluster_num)
  result = kmeans(X, cluster_num; maxiter=300, display=:none)
  df = DataFrame(index=1:length(values), data=values)
  centers = mapslices(x->[x], result.centers, dims=1)[:]
  clusters = df |>
               @groupby(result.assignments[_.index]) |>
               @map(x -> x |> @map(typeof(_.data) <: Array ? _.data : _.data.value) |> collect) |>
               collect
  (centers, clusters)
end

function call(values, cluster_num)
  X = hcat(values...)
  centers, clusters = (nothing, nothing)
  score = nothing
  for i in 1:100
    centers_, clusters_ = kmeans_internal(values, X, cluster_num)
    if isnothing(score) || score > score(centers_, clusters_)
      centers, clusters = (centers_, clusters_)
    end
  end

  Dict(:clusters => clusters, :centers => centers)
end

end


text:
module REPLHistory

using InteractiveUtils
using REPL

export history

function history(n = 5)
    out = ""
    n = n - 1
    list = [];
    for l in reverse(readlines(REPL.find_hist_file()))[2:end]
        l = strip(l)
        if n != 0
            (l[1] != '#') && (n = n - 1)
            push!(list, l);
        else
            break
        end
    end
    next = ""
    for l in reverse(list)
        l = next * l
        if l == "# mode: shell"
            next = "   #;"
        elseif l == "# mode: pkg"
            next = "   #]"
        elseif l == "# mode: help"
            next = "   #?"
        else
            (l[1] != '#') && (out = out * l * "\n")
            next = ""
        end
    end
    return out
end

end


text:
# TODO:
# * add transform and transfom! functions
# * add NT (or better name) to column selector passing NamedTuple
#   (also in other places: filter, combine)
# * add select/select!/transform/transform! for GroupedDataFrame

# normalize_selection function makes sure that whatever input format of idx is it
# will end up in one of four canonical forms
# 1) Int
# 2) AbstractVector{Int}
# 3) Pair{Int, Pair{ColRename, Symbol}}
# 4) Pair{Int, <:Pair{<:Base.Callable, Symbol}}
# 5) Pair{AbstractVector{Int}, <:Pair{<:Base.Callable, Symbol}}
# 6) Pair{Int, Pair{ByRow, Symbol}}
# 7) Pair{AbstractVector{Int}, Pair{ByRow, Symbol}}

"""
    ColRename

A singleton type indicating that column renaming operation was requested in `select`.
"""
struct ColRename end

"""
    ByRow

A type used for selection operations to signal that the wrapped function should
be applied to each element (row) of the selection.
"""
struct ByRow{T}
    fun::T
end

(f::ByRow)(cols::AbstractVector...) = f.fun.(cols...)

# add a method to funname defined in other/utils.jl
funname(row::ByRow) = funname(row.fun)

normalize_selection(idx::AbstractIndex, sel) =
    try
        idx[sel]
    catch e
        if e isa MethodError && e.f === getindex && e.args === (idx, sel)
            throw(ArgumentError("Unrecognized column selector: $sel"))
        else
            rethrow(e)
        end
    end

function normalize_selection(idx::AbstractIndex, sel::ColumnIndex)
    c = idx[sel]
    return c => ColRename() => _names(idx)[c]
end

function normalize_selection(idx::AbstractIndex, sel::Pair{<:ColumnIndex, Symbol})
    c = idx[first(sel)]
    return c => ColRename() => last(sel)
end

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:Any,<:Pair{<:Union{Base.Callable, ByRow}, Symbol}})
    rawc = first(sel)
    if rawc isa AbstractVector{Int}
        c = rawc
    elseif rawc isa AbstractVector{Symbol}
        c = [idx[n] for n in rawc]
    else
        c = try
                idx[rawc]
            catch e
                if e isa MethodError && e.f === getindex && e.args === (idx, rawc)
                    throw(ArgumentError("Unrecognized column selector: $rawc"))
                else
                    rethrow(e)
                end
            end
    end
    if length(c) == 0 && first(last(sel)) isa ByRow
        throw(ArgumentError("at least one column must be passed to a " *
                            "`ByRow` transformation function"))
    end
    return c => last(sel)
end

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:ColumnIndex,<:Union{Base.Callable, ByRow}})
    c = idx[first(sel)]
    fun = last(sel)
    newcol = Symbol(_names(idx)[c], "_", funname(fun))
    return c => fun => newcol
end

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:Any, <:Union{Base.Callable,ByRow}})
    rawc = first(sel)
    if rawc isa AbstractVector{Int}
        c = rawc
    elseif rawc isa AbstractVector{Symbol}
        c = [idx[n] for n in rawc]
    else
        c = try
                idx[rawc]
            catch e
                if e isa MethodError && e.f === getindex && e.args === (idx, rawc)
                    throw(ArgumentError("Unrecognized column selector: $rawc"))
                else
                    rethrow(e)
                end
            end
    end
    if length(c) == 0 && last(sel) isa ByRow
        throw(ArgumentError("at least one column must be passed to a " *
                            "`ByRow` transformation function"))
    end
    fun = last(sel)
    if length(c) > 3
        newcol = Symbol(join(@views(_names(idx)[c[1:2]]), '_'), "_etc_", funname(fun))
    elseif isempty(c)
        newcol = Symbol(funname(fun))
    else
        newcol = Symbol(join(view(_names(idx), c), '_'), '_', funname(fun))
    end
    return c => fun => newcol
end

function select_transform!(nc::Pair{Int, Pair{ColRename, Symbol}},
                           df::AbstractDataFrame, newdf::DataFrame,
                           transformed_cols::Dict{Symbol, Any}, copycols::Bool)
    col_idx, (_, newname) = nc
    # it is allowed to request column tranformation only once
    @assert !hasproperty(newdf, newname)
    newdf[!, newname] = copycols ? df[:, col_idx] : df[!, col_idx]
    # mark that column transformation was applied
    # nothing is not possible otherwise as a value in this dict
    transformed_cols[newname] = nothing
end

function select_transform!(nc::Pair{<:Union{Int, AbstractVector{Int}},
                                    <:Pair{<:Union{Base.Callable, ByRow}, Symbol}},
                           df::AbstractDataFrame, newdf::DataFrame,
                           transformed_cols::Dict{Symbol, Any}, copycols::Bool)
    col_idx, (fun, newname) = nc
    @assert !hasproperty(newdf, newname)
    cdf = eachcol(df)
    if col_idx isa Int
        res = fun(df[!, col_idx])
    else
        res = fun((cdf[i] for i in col_idx)...)
    end
    if res isa Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}
        throw(ArgumentError("return value from function $fun " *
                            "of type $(typeof(res)) is currently not allowed."))
    end
    if res isa AbstractVector
        respar = parent(res)
        if copycols && !(fun isa ByRow) &&
            (res isa SubArray || any(i -> respar === parent(cdf[i]), col_idx))
            newdf[!, newname] = copy(res)
        else
            newdf[!, newname] = res
        end
    else
        newdf[!, newname] = [res]
    end
    transformed_cols[newname] = nothing
end

"""
    select!(df::DataFrame, inds...)

Mutate `df` in place to retain only columns specified by `inds...` and return it.

Arguments passed as `inds...` can be any index that is allowed for column indexing.
In particular, regular expressions, `All`, `Between`, and `Not` selectors are supported.

Columns can be renamed using the `old_column => new_column_name` syntax,
and transformed using the `old_column => fun => new_column_name` syntax.
`new_column_name` must be a `Symbol`, and `fun` a function or a type. If `old_column`
is a `Symbol` or an integer then `fun` is applied to the corresponding column vector.
Otherwise `old_column` can be any column indexing syntax, in which case `fun`
will be passed the column vectors specified by `old_column` as separate arguments.
If `fun` returns a value of type other than `AbstractVector` then it will be wrapped
into a 1-element vector, unless its type is one of `AbstractDataFrame`, `NamedTuple`,
`DataFrameRow`, `AbstractMatrix`, in which case an error is thrown as currently these
return types are not allowed.

To apply `fun` to each row instead of whole columns, it can be wrapped in a `ByRow`
struct. In this case if `old_column` is a `Symbol` or an integer then `fun` is applied
to each element (row) of `old_column` using broadcasting. Otherwise `old_column` can be
any column indexing syntax, in which case `fun` will be passed one argument for each of
the columns specified by `old_column`. If `ByRow` is used it is not allowed for
`old_column` to select an empty set of columns.

Column transformation can also be specified using the short `old_column => fun` form.
In this case, `new_column_name` is automatically generated as `\$(old_column)_\$(fun)`.
Up to three column names are used for multiple input columns and they are joined
using `_`; if more than three columns are passed then the name consists of the
first two names and `etc` suffix then, e.g. `[:a,:b,:c,:d] => fun` produces
the new column name `:a_b_etc_fun`.

If a collection of column names is passed to `select!` then requesting duplicate column
names in target data frame are accepted (e.g. `select!(df, [:a], :, r"a")` is allowed)
and only the first occurrence is used. In particular a syntax to move column `:col`
to the first position in the data frame is `select!(df, :col, :)`.
On the contrary, output column names of renaming, transformation and single column
selection operations must be unique, so e.g. `select!(df, :a, :a => :a)` or
`select!(df, :a, :a => ByRow(sin) => :a)` are not allowed.

Note that including the same column several times in the data frame via renaming
or transformations that return the same object without copying will create column aliases.
An example of such a situation is `select!(df, :a, :a => :b, :a => identity => :c)`.

# Examples
```jldoctest
julia> df = DataFrame(a=1:3, b=4:6)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> select!(df, 2)
3×1 DataFrame
│ Row │ b     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 4     │
│ 2   │ 5     │
│ 3   │ 6     │

julia> df = DataFrame(a=1:3, b=4:6)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> select!(df, :a => ByRow(sin) => :c, :b)
3×2 DataFrame
│ Row │ c        │ b     │
│     │ Float64  │ Int64 │
├─────┼──────────┼───────┤
│ 1   │ 0.841471 │ 4     │
│ 2   │ 0.909297 │ 5     │
│ 3   │ 0.14112  │ 6     │

julia> select!(df, :, [:c, :b] => (c,b) -> c .+ b .- sum(b)/length(b))
3×3 DataFrame
│ Row │ c        │ b     │ c_b_function │
│     │ Float64  │ Int64 │ Float64      │
├─────┼──────────┼───────┼──────────────┤
│ 1   │ 0.841471 │ 4     │ -0.158529    │
│ 2   │ 0.909297 │ 5     │ 0.909297     │
│ 3   │ 0.14112  │ 6     │ 1.14112      │
```

"""
function select!(df::DataFrame, inds::AbstractVector{Int})
    if isempty(inds)
        empty!(_columns(df))
        empty!(index(df))
        return df
    end
    indmin, indmax = extrema(inds)
    if indmin < 1
        throw(ArgumentError("indices must be positive"))
    end
    if indmax > ncol(df)
        throw(ArgumentError("indices must not be greater than number of columns"))
    end
    if !allunique(inds)
        throw(ArgumentError("indices must not contain duplicates"))
    end
    copy!(_columns(df), _columns(df)[inds])
    x = index(df)
    copy!(_names(x), _names(df)[inds])
    empty!(x.lookup)
    for (i, n) in enumerate(x.names)
        x.lookup[n] = i
    end
    return df
end

select!(df::DataFrame, c::Int) = select!(df, [c])
select!(df::DataFrame, c::Union{AbstractVector{<:Integer}, AbstractVector{Symbol},
                                Colon, All, Not, Between, Regex}) =
    select!(df, index(df)[c])

function select!(df::DataFrame, cs...)
    newdf = select(df, cs..., copycols=false)
    copy!(_columns(df), _columns(newdf))
    x = index(df)
    copy!(_names(x), _names(newdf))
    empty!(x.lookup)
    for (i, n) in enumerate(x.names)
        x.lookup[n] = i
    end
    return df
end

"""
    select(df::AbstractDataFrame, inds...; copycols::Bool=true)

Create a new data frame that contains columns from `df` specified by `inds` and return it.

Arguments passed as `inds...` can be any index that is allowed for column indexing.
In particular, regular expressions, `All`, `Between`, and `Not` selectors are supported.

Also if `df` is a `DataFrame` or `copycols=true` then column renaming and transformations
are supported.

Columns can be renamed using the `old_column => new_column_name` syntax,
and transformed using the `old_column => fun => new_column_name` syntax.
`new_column_name` must be a `Symbol`, and `fun` a function or a type. If `old_column`
is a `Symbol` or an integer then `fun` is applied to the corresponding column vector.
Otherwise `old_column` can be any column indexing syntax, in which case `fun`
will be passed the column vectors specified by `old_column` as separate arguments.
If `fun` returns a value of type other than `AbstractVector` then it will be wrapped
into a 1-element vector, unless its type is one of `AbstractDataFrame`, `NamedTuple`,
`DataFrameRow`, `AbstractMatrix`, in which case an error is thrown as currently these
return types are not allowed.

To apply `fun` to each row instead of whole columns, it can be wrapped in a `ByRow`
struct. In this case if `old_column` is a `Symbol` or an integer then `fun` is applied
to each element (row) of `old_column` using broadcasting. Otherwise `old_column` can be
any column indexing syntax, in which case `fun` will be passed one argument for each of
the columns specified by `old_column`. If `ByRow` is used it is not allowed for
`old_column` to select an empty set of columns.

Column transformation can also be specified using the short `old_column => fun` form.
In this case, `new_column_name` is automatically generated as `\$(old_column)_\$(fun)`.
Up to three column names are used for multiple input columns and they are joined
using `_`; if more than three columns are passed then the name consists of the
first two names and `etc` suffix then, e.g. `[:a,:b,:c,:d] => fun` produces
the new column name `:a_b_etc_fun`.

If a collection of column names is passed to `select!` then requesting duplicate column
names in target data frame are accepted (e.g. `select!(df, [:a], :, r"a")` is allowed)
and only the first occurrence is used. In particular a syntax to move column `:col`
to the first position in the data frame is `select!(df, :col, :)`.
On the contrary, output column names of renaming, transformation and single column
selection operations must be unique, so e.g. `select!(df, :a, :a => :a)` or
`select!(df, :a, :a => ByRow(sin) => :a)` are not allowed.

If `df` is a `DataFrame` a new `DataFrame` is returned.
If `copycols=false`, then the returned `DataFrame` shares column vectors with `df` where possible.
If `copycols=true` (the default), then the returned `DataFrame` will not share columns with `df`.
The only exception for this rule is the `old_column => fun => new_column` transformation
when `fun` returns a vector that is not allocated by `fun` but is neither a `SubArray` nor one
of the input vectors.
In such a case a new `DataFrame` might contain aliases. Such a situation can only happen
with transformations which returns vectors other than their inputs, e.g. with
`select(df, :a => (x -> c) => :c1, :b => (x -> c) => :c2)`  when `c` is a vector object
or with `select(df, :a => (x -> df.c) => :c2)`.

If `df` is a `SubDataFrame` and `copycols=true` then a `DataFrame` is returned and
the same copying rules apply as for a `DataFrame` input:
this means in particular that selected columns will be copied.
If `copycols=false`, a `SubDataFrame` is returned without copying columns.

Note that including the same column several times in the data frame via renaming or
transformations that return the same object when `copycols=false` will create column
aliases. An example of such a situation is
`select(df, :a, :a => :b, :a => identity => :c, copycols=false)`.

# Examples
```jldoctest
julia> df = DataFrame(a=1:3, b=4:6)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> select(df, :b)
3×1 DataFrame
│ Row │ b     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 4     │
│ 2   │ 5     │
│ 3   │ 6     │

julia> select(df, Not(:b)) # drop column :b from df
3×1 DataFrame
│ Row │ a     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 1     │
│ 2   │ 2     │
│ 3   │ 3     │

julia> select(df, :a => :c, :b)
3×2 DataFrame
│ Row │ c     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> select(df, :a => ByRow(sin) => :c, :b)
3×2 DataFrame
│ Row │ c        │ b     │
│     │ Float64  │ Int64 │
├─────┼──────────┼───────┤
│ 1   │ 0.841471 │ 4     │
│ 2   │ 0.909297 │ 5     │
│ 3   │ 0.14112  │ 6     │

julia> select(df, :, [:a, :b] => (a,b) -> a .+ b .- sum(b)/length(b))
3×3 DataFrame
│ Row │ a     │ b     │ a_b_function │
│     │ Int64 │ Int64 │ Float64      │
├─────┼───────┼───────┼──────────────┤
│ 1   │ 1     │ 4     │ 0.0          │
│ 2   │ 2     │ 5     │ 2.0          │
│ 3   │ 3     │ 6     │ 4.0          │
```

"""
select(df::DataFrame, inds::AbstractVector{Int}; copycols::Bool=true) =
    DataFrame(_columns(df)[inds], Index(_names(df)[inds]),
              copycols=copycols)
select(df::DataFrame, c::Union{AbstractVector{<:Integer}, AbstractVector{Symbol},
                               Colon, All, Not, Between, Regex}; copycols::Bool=true) =
    select(df, index(df)[c], copycols=copycols)
select(df::DataFrame, c::ColumnIndex; copycols::Bool=true) =
    select(df, [c], copycols=copycols)

select(df::DataFrame, cs...; copycols::Bool=true) =
    _select(df, [normalize_selection(index(df), c) for c in cs], copycols)

function _select(df::AbstractDataFrame, normalized_cs, copycols::Bool)
    @assert !(df isa SubDataFrame && copycols==false)
    newdf = DataFrame()
    # the role of transformed_cols is the following
    # * make sure that we do not use the same target column name twice in transformations;
    #   note though that it can appear in no-transformation selection like
    #   `select(df, :, :a => ByRow(sin) => :a), where :a is produced both by `:`
    #   and by `:a => ByRow(sin) => :a`
    # * make sure that if some column is produced by transformation like
    #   `:a => ByRow(sin) => :a` and it appears earlier or later in non-transforming
    #   selection like `:` or `:a` then the transformation is computed and inserted
    #   in to the target data frame once and only once the first time the target column
    #   is requested to be produced.
    #
    # For example in:
    #
    # julia> df = DataFrame(a=1:2, b=3:4)
    # 2×2 DataFrame
    # │ Row │ a     │ b     │
    # │     │ Int64 │ Int64 │
    # ├─────┼───────┼───────┤
    # │ 1   │ 1     │ 3     │
    # │ 2   │ 2     │ 4     │
    #
    # julia> select(df, :, :a => ByRow(sin) => :a, :a, 1)
    # 2×2 DataFrame
    # │ Row │ a        │ b     │
    # │     │ Float64  │ Int64 │
    # ├─────┼──────────┼───────┤
    # │ 1   │ 0.841471 │ 3     │
    # │ 2   │ 0.909297 │ 4     │
    #
    # we compute column :a immediately when we process `:` although it is specified
    # later by `:a=>sin=>:a` because we know from `transformed_cols` variable that
    # it will be computed later via a transformation
    transformed_cols = Dict{Symbol, Any}()
    for nc in normalized_cs
        if nc isa Pair
            newname = last(last(nc))
            @assert newname isa Symbol
            if haskey(transformed_cols, newname)
                throw(ArgumentError("duplicate target column name $newname passed"))
            end
            transformed_cols[newname] = nc
        end
    end
    for nc in normalized_cs
        if nc isa AbstractVector{Int}
            allunique(nc) || throw(ArgumentError("duplicate column names selected"))
            for i in nc
                newname = _names(df)[i]
                # as nc is a multiple column selection without transformations
                # we allow duplicate column names with selections applied earlier
                # and ignore them for convinience, to allow for e.g. select(df, :x1, :)
                if !hasproperty(newdf, newname)
                    if haskey(transformed_cols, newname)
                        # if newdf does not have a column newname
                        # but a column transformation was requested for this column
                        # then apply the transformation immediately
                        # in such a case nct may not be nothing, as if it were
                        # nothing then newname should be preasent in newdf already
                        nct = transformed_cols[newname]
                        @assert nct !== nothing
                        select_transform!(nct, df, newdf, transformed_cols, copycols)
                    else
                        newdf[!, newname] = copycols ? df[:, i] : df[!, i]
                    end
                end
            end
        else
            # nc is normalized so it has a form src_cols => fun => Symbol
            newname = last(last(nc))
            if hasproperty(newdf, newname)
                # it is possible that the transformation has already been applied
                # via multiple column selection, like in select(df, :, :x1 => :y1)
                # but then transformed_cols[newname] must be nothing
                @assert transformed_cols[newname] === nothing
            else
                select_transform!(nc, df, newdf, transformed_cols, copycols)
            end
        end
    end
    return newdf
end

select(dfv::SubDataFrame, ind::ColumnIndex; copycols::Bool=true) =
    select(dfv, [ind], copycols=copycols)
select(dfv::SubDataFrame, inds::Union{AbstractVector{<:Integer}, AbstractVector{Symbol},
                                      Colon, All, Not, Between, Regex}; copycols::Bool=true) =
    copycols ? dfv[:, inds] : view(dfv, :, inds)

function select(dfv::SubDataFrame, inds...; copycols::Bool=true)
    if copycols
        return _select(dfv, [normalize_selection(index(dfv), c) for c in inds], true)
    else
        # we do not support transformations here
        # newinds should not be large so making it Vector{Any} should be OK
        newinds = []
        seen_single_column = Set{Int}()
        for ind in inds
            if ind isa ColumnIndex
                ind_idx = index(dfv)[ind]
                if ind_idx in seen_single_column
                    throw(ArgumentError("selecting the same column multiple times using" *
                                        " Symbol or integer is not allowed ($ind was " *
                                        "passed more than once"))
                else
                    push!(seen_single_column, ind_idx)
                end
            else
                newind = normalize_selection(index(dfv), ind)
                if newind isa Pair
                    throw(ArgumentError("transforming and renaming columns of a " *
                                        "SubDataFrame is not allowed when `copycols=false`"))
                end
                push!(newinds, newind)
            end
        end
        return view(dfv, :, All(newinds...))
    end
end


text:
module PrivateExample4
export hello, domath

"""
    hello(who::String)

Return "Hello, `who`".
"""
hello4(who::String) = "Hello, $who"

"""
    domath(x::Number)

Return `x + 5`.
"""
domath4(x::Number) = x + 5

end


text:
# simple XGBoost example
# with tuning along the lines of this post:
# https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/

# Warning: This is a demonstration only. The small size of the data
# set and the large number of hyperparameters being tuned implies that
# overfitting is likely.

# Declare variables const for a speed boost

# using Revise
using StatsBase
using MLJ
using Random
using Plots
plotly()

# load a task (data plus learning objective):
task = load_crabs()
shuffle!(MersenneTwister(134), task)

# list models matching the task:
models(task)

# assuming XGBoost package is in load path, we can
# load the model code for the applicable XGBoost model:
@load XGBoostClassifier

## EXTRACT A HOLDOUT TEST SET

train, test = partition(1:nrows(task), 0.7)

testtask = task[test]
task = task[train]

## BALANCED TARGET VALUES?

countmap(task.y)
# B:O is 67:73, so yes.

## DEFINE THE XGBOOST MODEL

# a time-dependent random seed is obtained with seed=-1:
boost = XGBoostClassifier(num_round=10, eta=0.05)

## DETERMINE NUMBER OF ITERATIONS FOR DEFAULT LEARNING RATE

# wrap model in task:
mach = machine(boost, task)

# range for forest.n:
r = range(boost, :num_round, lower=10, upper=1000)

# generate a learning curve, cross entropy loss vs num_round:
curve = learning_curve!(mach,
                        resampling=CV(),
                        nested_range=(num_round=r,),
                        measure=cross_entropy,
                        resolution=100)
plot(curve.parameter_values, curve.measurements)

boost.num_round = 300

## TUNE max_depth AND min_child_weight

# to check order of paramaters:
params(boost)

r1 = range(boost, :max_depth, lower=3, upper=10)
r2 = range(boost, :min_child_weight, lower=0, upper=5)

# define nested ranges, matching pattern of params(forest):
nested_ranges = (max_depth=r1, min_child_weight=r2)

tuned1 = TunedModel(model=boost,
                   tuning=Grid(resolution=8),
                   resampling=CV(),
                   nested_ranges=nested_ranges,
                   measure=cross_entropy)

# tune and plot the results:
mach = machine(tuned1, task)
fit!(mach)
plot(mach)

# fine tune:
tuned1.nested_ranges = (max_depth=range(boost, :max_depth, lower=3, upper=7),
                        min_child_weight = range(boost, :min_child_weight, lower=0, upper=2))
fit!(mach)
plot(mach)

# extract optimal model:
boost = fitted_params(mach).best_model
boost.min_child_weight # 0.57
boost.max_depth        # 3

## TUNE gamma

# We'll use a learning curve to examine the effect of gamma:
mach = machine(boost, task)
curve = learning_curve!(mach,
                        resampling=CV(),
                        nested_range=(gamma=range(boost, :gamma, lower=0, upper=10),),
                        measure=cross_entropy,
                        resolution=30)
plot(curve.parameter_values, curve.measurements)
# gamma appears to have no effect on performance
boost.gamma = 0

## RECALIBRATE NUMBER OF ITERATIONS

# generate a learning curve, cross entropy loss vs num_round:
curve = learning_curve!(mach,
                        resampling=CV(),
                        nested_range=(num_round=r,), # as before
                        measure=cross_entropy,
                        resolution=100)
plot(curve.parameter_values, curve.measurements)

boost.num_round = curve.parameter_values[argmin(curve.measurements)] # 560

## TUNE subsample AND colsample_bytree

r1 = range(boost, :subsample, lower=0.6, upper=1.0)
r2 = range(boost, :colsample_bytree, lower=0.6, upper=1.0)

# define nested ranges, matching pattern of params(forest):
nested_ranges = (subsample=r1, colsample_bytree=r2)

tuned1 = TunedModel(model=boost,
                   tuning=Grid(resolution=8),
                   resampling=CV(),
                   nested_ranges=nested_ranges,
                   measure=cross_entropy)

# tune and plot the results:
mach = machine(tuned1, task)
fit!(mach)
plot(mach)

# fine tune:
tuned1.nested_ranges = (subsample=range(boost, :subsample, lower=0.5, upper=0.85),
                        colsample_bytree=range(boost, :colsample_bytree, lower=0.8, upper=1.0))
fit!(mach)
plot(mach)

# extract optimal model:
boost = fitted_params(mach).best_model
boost.subsample          # 0.7
boost.colsample_bytree   # 1

## LOWER THE LEARNING RATE AND ADD MORE TREES

boost.eta = 0.005

# generate a learning curve, cross entropy loss vs num_round:
mach = machine(boost, task)

r = range(boost, :num_round, lower=10, upper=7000)

curve = learning_curve!(mach,
                        resampling=CV(),
                        nested_range=(num_round=r,), # as before
                        measure=cross_entropy,
                        resolution=30, verbosity=1)
plot(curve.parameter_values, curve.measurements)

boost.num_round = curve.parameter_values[argmin(curve.measurements)] # 6277

## TRAIN ON ALL THE DATA AND EVALUATE ON THE HOLDOUT SET

fit!(mach)
yhat = predict(mach, testtask.X)
cross_entropy(yhat, testtask.y)
misclassification_rate(yhat, testtask.y)


text:
using LAP_julia: load_H, make_flow_from_H, gen_rand_points, get_valid_landmarks
using DataFrames, CSV

## create a birl like structure for with the oxford affine dataset
base_path = "/Users/MrTrololord/Downloads/mri_head_homography/"
target_images_path = joinpath(base_path, "images", "target")
target_images_fnames = filter(x -> !isfile(x), readdir(target_images_path))
target_images_paths = map(fname -> joinpath(target_images_path, fname), target_images_fnames)

COUNT = 500
max_disp = 40
# flow_gen_func(img_size) = gen_quad_flow(img_size, max_disp)
flow_gen_func(img_size) = gen_homo_flow(img_size, max_disp)

for target_img_path in target_images_paths
    k = parse(Int, target_img_path[findfirst(r"[0-9]+", target_img_path)])
    target_img = load(target_img_path)
    # create save paths
    target_land_path = joinpath(base_path, "landmarks", "target", "landmarks$(k).csv")
    source_land_path = joinpath(base_path, "landmarks", "source", "landmarks$(k).csv")
    source_img_path = joinpath(base_path, "images", "source", "img$(k).png")
    # generate rand inds in target
    target_inds = gen_rand_points(target_img, COUNT, "Gridded")
    cur_flow = flow_gen_func(size(target_img))
    # edit target inds and save
    target_inds = LAP_julia.get_valid_landmarks(cur_flow, target_inds)
    LAP_julia.save_landmarks(target_inds, target_land_path)
    # warp target inds and save
    source_points = LAP_julia.move_landmarks(transpose(LAP_julia.inds_to_points(target_inds)), cur_flow)
    # rng = extrema.(indices_spatial(cur_flow))
    # @info source_points[:,1] size(cur_flow) rng
    # @assert all(LAP_julia.is_in_bounds.(source_points[:,1], rng[1]...) .& LAP_julia.is_in_bounds.(source_points[:,2], rng[2]...))
    LAP_julia.save_landmarks(source_points, source_land_path)
    # LAP_julia.save_shift_landmarks(target_land_path, cur_flow.*(-1), source_land_path)
    # create source img
    source = warp_img(target_img, -real.(cur_flow), -imag.(cur_flow))
    save(source_img_path, source)
    @info k target_img_path size(target_inds)
end

# create loc table
tmp_df = CSV.read("/Volumes/davidkunz/Documents/bakalarka/birl_small/dataset_small_repaired.csv") |> DataFrame

df = deepcopy(tmp_df)
delete!(df, 1:size(tmp_df,1))

# rename!(df, Symbol.(replace.(string.(names(df)), Ref(" "=> "_"))))
# rename!(df, Symbol.(names(df)))

names(df)
eltypes(df)

abc = [:a , :b, :c, :d, :e, :f, :g, :h, :i, :j, :k]
new_df = DataFrame(a = eltypes(df)[1][],
                   b = eltypes(df)[2][],
                   c = eltypes(df)[3][],
                   d = eltypes(df)[4][],
                   e = eltypes(df)[5][],
                   f = eltypes(df)[6][],
                   g = eltypes(df)[7][],
                   h = eltypes(df)[8][],
                   i = eltypes(df)[9][],
                   j = eltypes(df)[10][],
                   k = eltypes(df)[11][])

function rename_dict_keys(from_name, to_name, dict)
    new_dict = Dict()
    rename_dict = Dict(from_name .=> to_name)
    for (key, value) in rename_dict
        new_dict[value] = dict[key]
    end
    return new_dict
end

experiment_count = count(x -> !isfile(x), readdir(target_images_path))

for k in 1:experiment_count
    row = Dict{String}{Any}()
    row["Column1"] = k
    # images
    row["Target image"] = joinpath("images", "target", "img$k.png")
    row["Source image"] = joinpath("images", "source", "img$k.png")
    # lnds
    row["Target landmarks"] = joinpath("landmarks", "target", "landmarks$k.csv")
    row["Source landmarks"] = joinpath("landmarks", "source", "landmarks$k.csv")
    # sizes
    source_img = load(joinpath(base_path, "images", "source", "img$k.png"))
    row["Image size [pixels]"] = string(size(source_img))
    row["Image diagonal [pixels]"] = round(sqrt(sum(size(source_img).^2)), digits=1)
    row["status"] = "training"
    # empty values
    row["Warped target landmarks"] = missing
    row["Warped source landmarks"] = missing
    row["Execution time [minutes]"] = missing

    @info rename_dict_keys(names(tmp_df), abc, row)
    push!(new_df, rename_dict_keys(names(tmp_df), abc, row))
end

new_df

# rename back
renamed = rename!(new_df, names(tmp_df))

CSV.write(joinpath(base_path, "loc_table.csv"), renamed)


text:
const seps = "cfhistuCFHISTU"
const sep_div = 3.5
const guard_div = 12
const minalphalen = 16

function consistent_shuffle(astr::String, salt::String)
    len_salt = length(salt)
    len_salt == 0 && return astr

    str = copy(convert(Vector{UInt8}, astr))
    i, idx, i_sum = length(str) - 1, 0, 0
    while i > 0
        idx %= len_salt
        asc = Int(salt[idx+1])
        i_sum += asc
        j = (asc + idx + i_sum) % i
        temp = str[j+1]
        str[j+1] = str[i+1]
        str[i+1] = temp

        i = i - 1
        idx = idx + 1
    end
    return ascii(String(str))
end

function hash_num(input::Int, alphabet::String)
    alphalen = length(alphabet)
    hash = UInt8[]
    push!(hash,alphabet[(input%alphalen)+1])
    input = div(input,alphalen)
    while input > 0
        push!(hash,alphabet[(input%alphalen)+1])
        input = div(input,alphalen)
    end
    convert(String,reverse(hash))
end

function unhash_num(input::SubString{String}, alphabet::String)
    num = 0
    alphalen = length(alphabet)
    inputlen = length(input)
    for i in 1:inputlen
        pos = search(alphabet, input[i])-1
        num += pos * alphalen^(inputlen-i)
    end
    num
end

immutable Hashid
    salt::String
    min_length::Int
    alphabet::Any
    separators::Any
    guards::Any

    function Hashid(salt::String = "", min_length::Int = 0, alphabet::String = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890")
        alphabet=replace(alphabet,r"\s+","")    # remove space characters
        alpha_separators = unique(intersect(alphabet, seps))
        alpha = ascii(String(setdiff(alphabet,alpha_separators)))
        if length(alpha)+length(alpha_separators) < minalphalen
            error("This algorithm requires an alphabet whose length >= $minalphalen")
        end
        if min_length < 0
            error("Minimum hash length must be >= 0")
        end
        alpha_len = length(alpha)
        sep_len = length(alpha_separators)

        alpha_separators = consistent_shuffle(ascii(String(alpha_separators)), salt)
        if alpha_separators == "" || alpha_len/sep_len > sep_div
            new_len = Int(ceil(alpha_len/sep_div))
            (new_len) == 1 && (new_len = 2)

            if new_len > sep_len
                diff = new_len - sep_len
                alpha_separators *= alpha[1:diff]
                alpha = alpha[diff+1:end]
                alpha_len = length(alpha)
            else
                alpha_separators = alpha_separators[1:new_len]
            end
        end

        alpha = consistent_shuffle(alpha, salt)
        guard_count = Int(ceil(length(alpha)/guard_div))

        if alpha_len < 3
            guards = alpha_separators[1:guard_count]
            alpha_separators = alpha_separators[guard_count+1:end]
        else
            guards = alpha[1:guard_count]
            alpha = alpha[guard_count+1:end]
        end

        new(salt, min_length, alpha, alpha_separators, guards)
    end
end

function _encode(h::Hashid, nums::Int...)
    numberHashInt = 0
    for i in 1:length(nums)
        if nums[i] < 0
            throw(DomainError())
        end
        numberHashInt += (nums[i] % (i+99))
    end
    alphabet = h.alphabet
    lottery = ret = string(alphabet[(numberHashInt % length(h.alphabet))+1])
    buffer = ""; ret_str = string(ret)
    for i in 1:length(nums)
        num = nums[i]
        buffer = lottery * h.salt * alphabet
        alphabet = consistent_shuffle(alphabet, buffer[1:length(alphabet)])
        last = hash_num(num, alphabet)
        ret_str *= last

        if i < length(nums)
            num %= Int(last[1])+i
            seps_idx = num % length(h.separators)+1
            ret_str = string(ret_str,h.separators[seps_idx])
        end
    end
    if length(ret_str) < h.min_length
        guard_idx = ((numberHashInt + Int(ret_str[1])) % length(h.guards))+1
        guard = h.guards[guard_idx]
        ret_str = string(guard,ret_str)
        if length(ret_str) < h.min_length
            guard_idx = ((numberHashInt + Int(ret_str[2])) % length(h.guards))+1
            guard = h.guards[guard_idx]
            ret_str = string(ret_str, guard)
        end
    end
    halflen = div(length(alphabet),2)
    while length(ret_str) < h.min_length
        alphabet = consistent_shuffle(alphabet, alphabet)
        ret_str = alphabet[halflen+1:end] * ret_str * alphabet[1:halflen]
        excess = length(ret_str) - h.min_length
        if excess > 0
            start = div(excess,2)
            ret_str = ret_str[start:start + h.min_length]
        end
    end
    ret_str
end

function _decode(h::Hashid, hash::String)
    ret = Int[]
    guards = Array{Char,1}(h.guards)
    separators = Array{Char,1}(h.separators)
    hash_array = split(hash, guards)
    op = ""
    alphabet = h.alphabet
    for tmp in hash_array
        op *= tmp * ", "
    end
    hash_breakdown = hash_array[2 <= length(hash_array) <= 3 ? 2 : 1]
    lottery = hash_breakdown[1]
    hash_array = split(hash_breakdown[2:end], separators)
    for sub_hash in hash_array
        buffer = string(lottery, h.salt, alphabet)
        alphabet = consistent_shuffle(alphabet, buffer[1:length(alphabet)])
        push!(ret, unhash_num(sub_hash, alphabet))
    end
    ret
end

function encode(h::Hashid, nums::Int...)
    length(nums) == 0 ? "" : _encode(h, nums...)
end

function encode(h::Hashid, nums::Vector{Int})
    _encode(h, nums...)
end

function decode(h::Hashid, e::String)
    e == "" ? Int[] : _decode(h, e)
end


text:
using LinearAlgebra
using Statistics

"""
Apply the Local Ensemble Transform Kalman Filter algorithm.

- `A`: state variable ensemble matrix
- `d`: observation vector
- `HA`: model-predicted observation ensemble matrix
- `E`: observation error ensemble matrix
- `xs`: indices of state vector that map to observation local patches
- `ys`: indices of observation local patches
- `ρ`: covariance inflation parameter (optional, default value of 1.05)
- `diagR`: force observation error covariance to be diagonal, i.e. no spatial correlation (optional, default is false)

"""
function letkf(A::Matrix, d::Vector, HA::Matrix, E::Matrix, xs::Vector, ys::Vector; ρ=1.05, diagR=false)
    Aa = zeros(size(A))
    ndim, nens = size(A)
    nobs = length(d)
    if size(E) == (nobs, nobs)
        R = E
    else
        if diagR
            R = diagm(0 => diag((1 / (nens - 1)) * E * E'))
        else
            R = (1 / (nens - 1)) * E * E'
        end
    end
    Y = HA .- mean(HA, dims=2)
    X = A .- mean(A, dims=2)
    for (lx, ly) in zip(xs, ys)
        Xₗ = X[lx, :]
        Yₗ = Y[ly, :]
        Rₗ = R[ly, ly]
        C = Yₗ' * pinv(Rₗ)
        P = pinv((nens - 1) / ρ * I + C * Yₗ)
        W = real(sqrt((nens - 1) * P))
        w = P * C * (d[ly] .- mean(HA[ly, :], dims=2))
        W = W .+ w
        Xₗᵃ = Xₗ * W .+ mean(A[lx, :], dims=2)
        Aa[lx, :] = Xₗᵃ
    end
    Aa
end