text:
"""
```
plot_history_and_forecast(m, var, class, input_type, cond_type;
    title = "", plot_handle = plot(), kwargs...)

plot_history_and_forecast(m, vars, class, input_type, cond_type;
    forecast_string = "", use_bdd = :unbdd,
    plotroot = figurespath(m, \"forecast\"), titles = [],
    plot_handles = fill(plot(), length(vars)), verbose = :low,
    add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
    kwargs...)

plot_history_and_forecast(ms::Vector, vars::Vector{Symbol}, class::Symbol,
    input_type::Symbol, cond_type::Symbol;
    forecast_string::String = "",
    use_bdd::Symbol = :bdd,
    modal_line::Bool = false, untrans::Bool = false,
    fourquarter::Bool = false,
    plotroot::String = figurespath(ms[1], "forecast"),
    titles::Vector{String} = String[],
    plot_handles::Vector{Plots.Plot} = Plots.Plot[plot() for i = 1:length(vars)],
    verbose::Symbol = :low, names = Dict{Symbol, Vector{String}}(),
    start_date::Date = Date(2019,3,31), end_date::Date = iterate_quarters(date_forecast_start(m), 12),
    add_new_model::Bool = false, new_model::AbstractDSGEModel = Model1002(),
    new_cond_type::Symbol = :full, outfile_end::String =  "",
    new_forecast_string::String = forecast_string,
    add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
    forecast_strings::Vector{String} = repeat([forecast_string], length(ms)),
    kwargs...)
```

Plot history and forecast for `var` or `vars`. If these correspond to a
full-distribution forecast, you can specify the `bands_style` and `bands_pcts`.

### Inputs

- `m::AbstractDSGEModel`
- `var::Symbol` or `vars::Vector{Symbol}`: variable(s) to be plotted,
  e.g. `:obs_gdp` or `[:obs_gdp, :obs_nominalrate]`
- `class::Symbol`
- `input_type::Symbol`
- `cond_type::Symbol`

### Keyword Arguments

- `forecast_string::String = ""`
- `use_bdd::Symbol = :unbdd`: specifies which combination of means and bands to use
    a. `:bdd` -> bounded bands and bounded means (from `:bddforecastobs`, etc.)
    b. `:bdd_and_unbdd` -> bounded bands (from `:bddforecastobs`, etc.) and unbounded means (from `:forecastobs`, etc.)
    c. `:unbdd` -> unbounded bands and unbounded means (from `:forecastobs`, etc.)
- `modal_line::Bool = false`: if true, the modal line is plotted instead of the mean.
- `untrans::Bool = false`: whether to plot untransformed (model units) history and forecast
- `fourquarter::Bool = false`: whether to plot four-quarter history and forecast
- `plotroot::String = figurespath(m, "forecast")`: if nonempty, plots will be saved in that directory
- `title::String = ""` or `titles::Vector{String} = []`
- `plot_handle::Plot` or `plot_handles::Vector{Plot}`: existing plot(s) on which
  to overlay new forecast plot(s)
- `verbose::Symbol = :low`
- `names::Dict{Symbol, Vector{String}}`: Legend names for lines
- `start_date::Date` = Date(2019,3,31)`: Date plot should start at
- `end_date::Date` = iterate_quarters(date_forecast_start(m), 12)`: Date plot should end at
- `add_new_model::Bool` = false: Whether history should be plotted using a new model
- `new_model::AbstractDSGEModel` = Model1002(): The new model to use to plot history
- `new_cond_type::Symbol` = :full: Cond_type of new model
- `outfile_end::String = ""`: String to append to the end of file name of plot
- `new_forecast_string::String = forecast_string`: Forecast string for new model
- `add_trendline::Bool`: Whether add a trendline or not
- `trend_vals::Vector{Float64}`: The values to multiply last historical value by
    for the trendline (should start with 0.0 to keep last historical value the same)
- `forecast_strings::Vector{String}`: Forecast string to use for each give model.

See `?histforecast` or `?histforecast_vector' for additional keyword arguments, all of which can be passed
into `plot_history_and_forecast`.

### Output
- `p::Plot` or `plots::OrderedDict{Symbol, Plot}`
"""
function plot_history_and_forecast(m::AbstractDSGEModel, var::Symbol, class::Symbol,
                                   input_type::Symbol, cond_type::Symbol;
                                   title::String = "", plot_handle::Plots.Plot = plot(),
                                   kwargs...)

    plots = plot_history_and_forecast(m, [var], class, input_type, cond_type;
                                      titles = isempty(title) ? String[] : [title],
                                      plot_handles = Plots.Plot[plot_handle],
                                      kwargs...)
    return plots[var]
end

function plot_history_and_forecast(m::AbstractDSGEModel, vars::Vector{Symbol}, class::Symbol,
                                   input_type::Symbol, cond_type::Symbol;
                                   forecast_string::String = "",
                                   use_bdd::Symbol = :bdd,
                                   modal_line::Bool = false, untrans::Bool = false,
                                   fourquarter::Bool = false,
                                   plotroot::String = figurespath(m, "forecast"),
                                   titles::Vector{String} = String[],
                                   plot_handles::Vector{Plots.Plot} = Plots.Plot[plot() for i = 1:length(vars)],
                                   verbose::Symbol = :low, outfile_end::String = "",
                                   add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
                                   trend_start_date::Date = date_forecast_start(m),
                                   kwargs...)
    # Determine output_vars
    if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    elseif untrans
        hist_prod  = :histut
        fcast_prod = :forecastut
    elseif fourquarter
        hist_prod  = :hist4q
        fcast_prod = :forecast4q
    else
        hist_prod  = :hist
        fcast_prod = :forecast
    end

    # Read in MeansBands
    hist  = read_mb(m, input_type, cond_type, Symbol(hist_prod, class), forecast_string = forecast_string)
    fcast = read_mb(m, input_type, cond_type, Symbol(fcast_prod, class), forecast_string = forecast_string,
                    use_bdd = use_bdd, modal_line = modal_line)

    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(m, var, class, detexify = detexify_title), vars)
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title, plot_handle) in zip(vars, titles, plot_handles)

        # Call recipe
        plots[var] = plot(plot_handle)
        histforecast!(var, hist, fcast;
                      ylabel = series_ylabel(m, var, class, untrans = untrans,
                                             fourquarter = fourquarter),
                      title = title, add_trendline = add_trendline, trend_vals = trend_vals,
                      trend_start_date = trend_start_date, kwargs...)
        # Save plot
        if !isempty(plotroot)
            output_file = get_forecast_filename(plotroot, filestring_base(m), input_type, cond_type,
                                                Symbol(fcast_prod, "_", detexify(var)),
                                                forecast_string = forecast_string * outfile_end,
                                                fileformat = plot_extension())
            @show output_file
            DSGE.save_plot(plots[var], output_file, verbose = verbose)
        end
    end
    return plots
end

function plot_history_and_forecast(ms::Vector, vars::Vector{Symbol}, class::Symbol,
                                   input_type::Symbol, cond_type::Symbol;
                                   forecast_string::String = "",
                                   use_bdd::Symbol = :bdd,
                                   modal_line::Bool = false, untrans::Bool = false,
                                   fourquarter::Bool = false,
                                   plotroot::String = figurespath(ms[1], "forecast"),
                                   titles::Vector{String} = String[],
                                   plot_handles::Vector{Plots.Plot} = Plots.Plot[plot() for i = 1:length(vars)],
                                   verbose::Symbol = :low, names = Dict{Symbol, Vector{String}}(),
                                   start_date::Date = date_forecast_start(m), end_date::Date = iterate_quarters(date_forecast_start(m), 12),
                                   add_new_model::Bool = false, new_model::AbstractDSGEModel = Model1002(),
                                   new_cond_type::Symbol = :full, outfile_end::String = "",
                                   new_forecast_string::String = forecast_string,
                                   add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
                                   forecast_strings::Vector{String} = repeat([forecast_string], length(ms)),
                                   trend_start_date::Date = start_date,
                                   kwargs...)
    # Determine output_vars
    if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    elseif untrans
        hist_prod  = :histut
        fcast_prod = :forecastut
    elseif fourquarter
        hist_prod  = :hist4q
        fcast_prod = :forecast4q
    else
        hist_prod  = :hist
        fcast_prod = :forecast
    end

    # Read in MeansBands
    hist = Vector{MeansBands}(undef, length(ms))
    fcast = Vector{MeansBands}(undef, length(ms))
    for i in 1:length(ms)
        hist[i]  = read_mb(ms[i], input_type, cond_type, Symbol(hist_prod, class), forecast_string = forecast_strings[i])
        fcast[i] = read_mb(ms[i], input_type, cond_type, Symbol(fcast_prod, class), forecast_string = forecast_strings[i],
                        use_bdd = use_bdd, modal_line = modal_line)
    end

    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(ms[1], var, class, detexify = detexify_title), vars)
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title, plot_handle) in zip(vars, titles, plot_handles)
        new_data = []
        if add_new_model
            new_data, tmp = read_forecast_output(new_model, :full, new_cond_type, Symbol(:hist, class), var,
                                                 forecast_string = new_forecast_string)
            new_data = vec(mean(new_data,dims=1)) .* 4
        end

        # Call recipe
        plots[var] = plot(plot_handle)
        histforecast_vector!(var, hist, fcast;
                      ylabel = series_ylabel(ms[1], var, class, untrans = untrans,
                                             fourquarter = fourquarter),
                      title = title, names = names, start_date = start_date, end_date = end_date,
                             add_new_model = add_new_model, new_data = new_data,
                             add_trendline = add_trendline, trend_vals = trend_vals,
                             trend_start_date = trend_start_date, kwargs...)
        # Save plot
        if !isempty(plotroot)
            output_file = get_forecast_filename(plotroot, filestring_base(ms[1]), input_type, cond_type,
                                                Symbol(fcast_prod, "_", detexify(var)),
                                                forecast_string = forecast_strings[1] * outfile_end,
                                                fileformat = plot_extension())

            save_plot(plots[var], output_file, verbose = verbose)
        end
    end
    return plots
end

@userplot HistForecast

"""
```
histforecast(var, hist, forecast;
    start_date = hist.means[1, :date], end_date = forecast.means[end, :date],
    names = Dict{Symbol, String}(), colors = Dict{Symbol, Any}(),
    alphas = Dict{Symbol, Float64}(), styles = Dict{Symbol, Symbol}(),
    bands_pcts = union(which_density_bands(hist, uniquify = true),
                       which_density_bands(forecast, uniquify = true)),
    bands_style = :fan, label_bands = false, transparent_bands = true,
    tick_size = 2)
```

User recipe called by `plot_history_and_forecast`.

### Inputs

- `var::Symbol`: e.g. `obs_gdp`
- `hist::MeansBands`
- `forecast::MeansBands`

### Keyword Arguments

- `start_date::Date`
- `end_date::Date`
- `names::Dict{Symbol, String}`: maps keys `[:hist, :forecast, :bands]` to
  labels. If a key is missing from `names`, a default value will be used
- `colors::Dict{Symbol, Any}`: maps keys `[:hist, :forecast, :bands]` to
  colors
- `alphas::Dict{Symbol, Float64}`: maps keys `[:hist, :forecast, :bands]` to
  transparency values (between 0.0 and 1.0)
- `styles::Dict{Symbol, Symbol}`: maps keys `[:hist, :forecast, :bands]` to
  linestyles
- `bands_pcts::Vector{String}`: which bands percentiles to plot
- `bands_style::Symbol`: either `:fan` or `:line`
- `label_bands::Bool`
- `transparent_bands::Bool`
- `tick_size::Int`: x-axis (time) tick size in units of years
- `add_new_model::Bool`: Adding history from another model?
- `new_data::Array{Flaot64,1}`: The new history to plot

Additionally, all Plots attributes (see docs.juliaplots.org/latest/attributes)
are supported as keyword arguments.
"""
histforecast

@recipe function f(hf::HistForecast;
                   start_date::Date = hf.args[2].means[1, :date],
                   end_date::Date = hf.args[3].means[end, :date],
                   names = Dict{Symbol, String}(),
                   colors = Dict{Symbol, Any}(),
                   alphas = Dict{Symbol, Float64}(),
                   styles = Dict{Symbol, Symbol}(),
                   bands_pcts = union(which_density_bands(hf.args[2], uniquify = true),
                                      which_density_bands(hf.args[3], uniquify = true)),
                   bands_style = :fan,
                   label_bands = false,
                   transparent_bands = true,
                   add_new_model::Bool = false, new_data::Array{Float64,1} = [],
                   add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
                   trend_start_date::Date = start_date, tick_size = 2)

    # Error checking
    if length(hf.args) != 3 || typeof(hf.args[1]) != Symbol ||
        typeof(hf.args[2]) != MeansBands || typeof(hf.args[3]) != MeansBands

        error("histforecast must be given a Symbol and two MeansBands. Got $(typeof(hf.args))")
    end

    for dict in [names, colors, styles, alphas]
        bad_keys = setdiff(keys(dict), [:hist, :forecast, :bands])
        if !isempty(bad_keys)
            error("Invalid key(s) in $dict: $bad_keys")
        end
    end

    # Concatenate MeansBands
    var, hist, forecast = hf.args
    combined = cat(hist, forecast)
    dates = combined.means[!, :date]
    start_trend_ind = findfirst(x -> x == trend_start_date, dates)
    if isnothing(start_trend_ind)
        @warn "trend_start_date ($(trend_start_date))is not located in the history or forecast, " *
            "so we are plotting the trend from the start_date ($(start_date))"
        trend_start_date = start_date
        start_trend_ind = findfirst(x -> x == trend_start_date, dates)
    end
    start_trend_val = combined.means[start_trend_ind, var]
    end_ind = findfirst(dates .== end_date)

    # Assign date ticks
    date_ticks = Base.filter(x -> start_date <= x <= end_date,    dates)
    date_ticks = Base.filter(x -> Dates.month(x) == 3,            date_ticks)
    date_ticks = Base.filter(x -> Dates.year(x) % tick_size == 0, date_ticks)
    xticks --> (map(Dates.value, date_ticks), map(Dates.year, date_ticks))

    # Bands
    sort!(bands_pcts, rev = true) # s.t. non-transparent bands will be plotted correctly
    inds = findall(start_date .<= combined.bands[var][!, :date] .<= end_date)

    for (i, pct) in enumerate(bands_pcts)
        seriestype := :line

        x = combined.bands[var][inds, :date]
        lb = combined.bands[var][inds, Symbol(pct, " LB")]
        ub = combined.bands[var][inds, Symbol(pct, " UB")]

        bands_color = haskey(colors, :bands) ? colors[:bands] : :blue
        bands_alpha = haskey(alphas, :bands) ? alphas[:bands] : 0.1
        bands_linestyle = haskey(styles, :bands) ? styles[:bands] : :solid

        if bands_style == :fan
            @series begin
                if transparent_bands
                    fillcolor := bands_color
                    fillalpha := bands_alpha
                else
                    if typeof(bands_color) in [Symbol, String]
                        bands_color = parse(Colorant, bands_color)
                    end
                    fillcolor := weighted_color_mean(bands_alpha*i, bands_color, colorant"white")
                    fillalpha := 1
                end
                linealpha  := 0
                fillrange  := ub
                label      := label_bands ? "$pct Bands" : ""
                x, lb
            end
        elseif bands_style == :line
            # Lower bound
            @series begin
                linewidth --> 2
                linecolor :=  bands_color
                linestyle :=  bands_linestyle
                label     :=  label_bands ? "$pct LB" : ""
                x, lb
            end

            # Upper bound
            @series begin
                linewidth --> 2
                linecolor :=  bands_color
                linestyle :=  bands_linestyle
                label     :=  label_bands ? "$pct UB" : ""
                x, ub
            end
        else
            error("bands_style must be either :fan or :line. Got $bands_style")
        end
    end

    # Mean history
    @series begin
        seriestype :=  :line
        linewidth  --> 2
        linecolor  :=  haskey(colors, :hist) ? colors[:hist] : :black
        if VERSION >= v"1.3"
            seriesalpha := haskey(alphas, :hist) ? alphas[:hist] : 1.0
        else
            alpha       := haskey(alphas, :hist) ? alphas[:hist] : 1.0
        end
        linestyle  :=  haskey(styles, :hist) ? styles[:hist] : :solid
        label      :=  haskey(names, :hist) ? names[:hist] : "History"

        if add_new_model
            inds = intersect(findall(start_date .<= dates .<= end_date), 1:length(new_data))
            combined.means[inds, :date], new_data[inds]
        else
            inds = intersect(findall(start_date .<= dates .<= end_date),
                             findall(hist.means[1, :date] .<= dates .<= hist.means[end, :date]))
            combined.means[inds, :date], combined.means[inds, var]
        end
    end

    # Mean forecast
    @series begin
        seriestype :=  :line
        linewidth  --> 2
        linecolor  :=  haskey(colors, :forecast) ? colors[:forecast] : :red
        if VERSION >= v"1.3"
            seriesalpha := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
        else
            alpha       := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
        end
        linestyle  :=  haskey(styles, :forecast) ? styles[:forecast] : :solid
        label      :=  haskey(names, :forecast) ? names[:forecast] : "Forecast"

        inds = intersect(findall(start_date .<= dates .<= end_date),
                         findall(hist.means[end, :date] .<= dates .<= forecast.means[end, :date]))
        combined.means[inds, :date], combined.means[inds, var]
    end

    # Add trendline
    if add_trendline
        @series begin
            seriestype :=  :line
            linewidth  --> 2
            # linecolor  :=  haskey(colors, :forecast) ? colors[:forecast] : :red
            if VERSION >= v"1.3"
                seriesalpha := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
            else
                alpha       := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
            end
            linestyle  :=  :dash # haskey(styles, :forecast) ? styles[:forecast] : :solid
            label      :=  "Trend"

            dates[start_trend_ind:end_ind], start_trend_val .+ trend_vals[1:(end_ind-start_trend_ind+1)]
        end
    end
end

# HistForecast but for vector without bands
@userplot HistForecast_vector
"""
```
histforecast_vector(var, hist, forecast;
    start_date::Date = hf.args[2][1].means[1, :date],
    end_date::Date = hf.args[3][1].means[end, :date],
    names = Dict{Symbol, Vector{String}}(),
    colors = Dict{Symbol, Any}(),
    alphas = Dict{Symbol, Float64}(),
    styles = Dict{Symbol, Symbol}(),
    add_new_model::Bool = false, new_data::Array{Float64,1} = [],
    bands_pcts = union(which_density_bands(hf.args[2][1], uniquify = true),
        which_density_bands(hf.args[3][1], uniquify = true)),
    bands_style = :fan,
    label_bands = false,
    transparent_bands = true,
    add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
    tick_size = 2)
```

User recipe called by `plot_history_and_forecast`.

### Inputs

- `var::Symbol`: e.g. `obs_gdp`
- `hist::Vector{MeansBands}`
- `forecast::Vector{MeansBands}`

### Keyword Arguments

- `start_date::Date`
- `end_date::Date`
- `names::Dict{Symbol, String}`: maps keys `[:hist, :forecast, :bands]` to
  labels. If a key is missing from `names`, a default value will be used
- `colors::Dict{Symbol, Any}`: maps keys `[:hist, :forecast, :bands]` to
  colors
- `alphas::Dict{Symbol, Float64}`: maps keys `[:hist, :forecast, :bands]` to
  transparency values (between 0.0 and 1.0)
- `styles::Dict{Symbol, Symbol}`: maps keys `[:hist, :forecast, :bands]` to
  linestyles
- `bands_pcts::Vector{String}`: which bands percentiles to plot
- `bands_style::Symbol`: either `:fan` or `:line`
- `label_bands::Bool`
- `transparent_bands::Bool`
- `tick_size::Int`: x-axis (time) tick size in units of years
- `add_new_model::Bool`: Adding history from another model?
- `new_data::Array{Float64,1}`: The new history to plot
- `add_trendline::Bool`: Whether add a trendline or not
- `trend_vals::Vector{Float64}`: The values to multiply last historical value by
    for the trendline (should start with 0.0 to keep last historical value the same)
- `plot_all_histories::Bool`: Whether to plot histories for each model (generally used
    for pseudo observables like Natural Rate.

Additionally, all Plots attributes (see docs.juliaplots.org/latest/attributes)
are supported as keyword arguments.
"""
histforecast_vector

@recipe function f(hf::HistForecast_vector;
                   start_date::Date = hf.args[2][1].means[1, :date],
                   end_date::Date = hf.args[3][1].means[end, :date],
                   names = Dict{Symbol, Vector{String}}(),
                   colors = Dict{Symbol, Any}(),
                   alphas = Dict{Symbol, Float64}(),
                   styles = Dict{Symbol, Symbol}(),
                   add_new_model::Bool = false, new_data::Array{Float64,1} = [],
                   bands_pcts = union(which_density_bands(hf.args[2][1], uniquify = true),
                                      which_density_bands(hf.args[3][1], uniquify = true)),
                   bands_style = :fan,
                   label_bands = false,
                   transparent_bands = true,
                   add_trendline::Bool = false, trend_vals::Vector{Float64} = [1.0],
                   trend_start_date::Date = start_date, plot_all_histories::Bool = false,
                   tick_size = 2)

    # Error checking
    if length(hf.args) != 3 || typeof(hf.args[1]) != Symbol ||
        typeof(hf.args[2][1]) != MeansBands || typeof(hf.args[3][1]) != MeansBands

        error("histforecast must be given a Symbol and two MeansBands. Got $(typeof(hf.args))")
    end

    for dict in [names, colors, styles, alphas]
        bad_keys = setdiff(keys(dict), [:hist, :forecast, :bands])
        if !isempty(bad_keys)
            error("Invalid key(s) in $dict: $bad_keys")
        end
    end

    # Concatenate MeansBands
    var, hist, forecast = hf.args
    combined = cat(hist[1], forecast[1])
    dates = combined.means[!, :date]
    start_trend_ind = findfirst(dates .== trend_start_date)
    if isnothing(start_trend_ind)
        @warn "trend_start_date ($(trend_start_date))is not located in the history or forecast, " *
            "so we are plotting the trend from the start_date ($(start_date))"
        trend_start_date = start_date
        start_trend_ind = findfirst(x -> x == trend_start_date, dates)
    end
    start_trend_val = combined.means[start_trend_ind, var]
    end_ind = findfirst(dates .== end_date)

    # Assign date ticks
    date_ticks = Base.filter(x -> start_date <= x <= end_date,    dates)
    date_ticks = Base.filter(x -> Dates.month(x) == 3,            date_ticks)
    date_ticks = Base.filter(x -> Dates.year(x) % tick_size == 0, date_ticks)
    xticks --> (map(Dates.value, date_ticks), map(Dates.year, date_ticks))

    # Bands
    sort!(bands_pcts, rev = true) # s.t. non-transparent bands will be plotted correctly
    inds = findall(start_date .<= combined.bands[var][!, :date] .<= end_date)
#=
    for hfs_ind in 1:length(hist)
        combined = cat(hist[hfs_ind], forecast[hfs_ind])
        inds = findall(start_date .<= combined.bands[var][!, :date] .<= end_date)

        for (i, pct) in enumerate(bands_pcts)
            seriestype := :line

            x = combined.bands[var][inds, :date]
            lb = combined.bands[var][inds, Symbol(pct, " LB")]
            ub = combined.bands[var][inds, Symbol(pct, " UB")]

            #bands_color = haskey(colors, :bands) ? colors[:bands] : :blue
            bands_alpha = haskey(alphas, :bands) ? alphas[:bands] : 0.1
            bands_linestyle = haskey(styles, :bands) ? styles[:bands] : :solid

            if bands_style == :fan
                @series begin
                    if transparent_bands
                        #fillcolor := bands_color
                        fillalpha := bands_alpha
                    else
                        #if typeof(bands_color) in [Symbol, String]
                        #    bands_color = parse(Colorant, bands_color)
                        #end
                        #fillcolor := weighted_color_mean(bands_alpha*i, bands_color, colorant"white")
                        fillalpha := 1
                    end
                    linealpha  := 0
                    fillrange  := ub
                    label      := label_bands ? "$pct Bands" : ""
                    x, lb
                end
            elseif bands_style == :line
                # Lower bound
                @series begin
                    linewidth --> 2
                    #linecolor :=  bands_color
                    linestyle :=  bands_linestyle
                    label     :=  label_bands ? "$pct LB" : ""
                    x, lb
                end

                # Upper bound
                @series begin
                    linewidth --> 2
                    #linecolor :=  bands_color
                    linestyle :=  bands_linestyle
                    label     :=  label_bands ? "$pct UB" : ""
                    x, ub
                end
            else
                error("bands_style must be either :fan or :line. Got $bands_style")
            end
        end
    end
=#
    # Mean history
    end_ind = plot_all_histories && !add_new_model ? length(hist) : 1
    for i in 1:end_ind
        combined = cat(hist[i], forecast[i])
        dates = combined.means[!, :date]
        inds = findall(start_date .<= combined.bands[var][!, :date] .<= end_date)

        @series begin
            seriestype :=  :line
            linewidth  --> 2
            linecolor  :=  haskey(colors, :hist) ? colors[:hist][i] : palette(:default)[i]
            if VERSION >= v"1.3"
                seriesalpha := haskey(alphas, :hist) ? alphas[:hist] : 1.0
            else
                alpha       := haskey(alphas, :hist) ? alphas[:hist] : 1.0
            end
            linestyle  :=  haskey(styles, :hist) ? styles[:hist] : :dash
            label      :=  haskey(names, :hist) ? names[:hist][i] : "History"

            if add_new_model
            inds = intersect(findall(start_date .<= dates .<= end_date), 1:length(new_data))
                combined.means[inds, :date], new_data[inds]
            else
                inds = intersect(findall(start_date .<= dates .<= end_date),
                                 findall(hist[i].means[1, :date] .<= dates .<= hist[i].means[end, :date]))
                combined.means[inds, :date], combined.means[inds, var]
            end
        end
    end

    # Add other mean forecasts too
    for hfs_ind in 1:length(hist)
        combined = cat(hist[hfs_ind], forecast[hfs_ind])
        dates = combined.means[!, :date]
        inds = findall(start_date .<= combined.bands[var][!, :date] .<= end_date)

        for (i, pct) in enumerate(bands_pcts)
            seriestype := :line

            x = combined.bands[var][inds, :date]
            lb = combined.bands[var][inds, Symbol(pct, " LB")]
            ub = combined.bands[var][inds, Symbol(pct, " UB")]

            bands_color = haskey(colors, :bands) ? colors[:bands] : palette(:default)[hfs_ind]
            # get_color_palette(:auto, plot_color(:white), 17)[hfs_ind] # :blue
            bands_alpha = haskey(alphas, :bands) ? alphas[:bands] : 0.1
            bands_linestyle = haskey(styles, :bands) ? styles[:bands] : :solid

            if bands_style == :fan
                @series begin
                    if transparent_bands
                        fillcolor := bands_color
                        fillalpha := bands_alpha
                    else
                        if typeof(bands_color) in [Symbol, String]
                            bands_color = parse(Colorant, bands_color)
                        end
                        fillcolor := weighted_color_mean(bands_alpha*i, bands_color, colorant"white")
                        fillalpha := 1
                    end
                    linealpha  := 0
                    fillrange  := ub
                    label      := label_bands ? "$pct Bands" : ""
                    x, lb
                end
            elseif bands_style == :line
                # Lower bound
                @series begin
                    linewidth --> 2
                    linecolor :=  bands_color
                    linestyle :=  bands_linestyle
                    label     :=  label_bands ? "$pct LB" : ""
                    x, lb
                end

                # Upper bound
                @series begin
                    linewidth --> 2
                    linecolor :=  bands_color
                    linestyle :=  bands_linestyle
                    label     :=  label_bands ? "$pct UB" : ""
                    x, ub
                end
            else
                error("bands_style must be either :fan or :line. Got $bands_style")
            end
        end

        @series begin
        seriestype :=  :line
        linewidth  --> 2
        linecolor  :=  haskey(colors, :forecast) ? colors[:forecast] : palette(:default)[hfs_ind]
            # get_color_palette(:auto, plot_color(:white), 17)[hfs_ind] # :red
        if VERSION >= v"1.3"
            seriesalpha := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
        else
            alpha       := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
        end
        linestyle  :=  haskey(styles, :forecast) ? styles[:forecast] : :solid
        label      :=  haskey(names, :forecast) ? names[:forecast][hfs_ind] : "Forecast"

        inds = intersect(findall(start_date .<= dates .<= end_date),
                         findall(hist[hfs_ind].means[end, :date] .<= dates .<= forecast[hfs_ind].means[end, :date]))
        combined.means[inds, :date], combined.means[inds, var]
        end
    end

    # Add trendline
    if add_trendline
        @series begin
            seriestype :=  :line
            linewidth  --> 2
            # linecolor  :=  haskey(colors, :forecast) ? colors[:forecast] : :red
            if VERSION >= v"1.3"
                seriesalpha := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
            else
                alpha       := haskey(alphas, :forecast) ? alphas[:forecast] : 1.0
            end
            linestyle  :=  :dash # haskey(styles, :forecast) ? styles[:forecast] : :solid
            label      :=  "Trend"

            dates[start_trend_ind:end_ind], start_trend_val .+ trend_vals[1:(end_ind-start_trend_ind+1)]
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ b438ab92-fefd-11ea-1054-632292d0965c
begin
	using Flux
	using Flux: @epochs
	using DataFrames
	using Gadfly
	using LinearAlgebra
	using Random
end

# ╔═╡ d6e28796-fefb-11ea-1b08-bbe1e9f98682
md"# 线性回归

## 基本

- 模型
- 训练
- 数据
- 损失函数
- 优化
***
"

# ╔═╡ a3a08c06-fefc-11ea-38d3-2526f9a1ccd2
md"$y = Xw + b$"

# ╔═╡ b3f67ca0-fefc-11ea-17cb-373174529e5c
md"
根据误差，利用梯度对参数$W$和$b$进行更新。

$W = W - \frac{lr}{batchsize} \frac{\partial L(W, b)}{\partial W}$

其中$batchsize$是每批数据训练大小，$lr$是学习率。"

# ╔═╡ 5a0fac6a-fefd-11ea-1daa-135ba73382ac
md"**例子**"

# ╔═╡ 924bbece-fefe-11ea-3cdf-75c5fe8eb392
begin
	#prepare data
	Random.seed!(123)
	num_inputs = 2
	num_inputs = 2
	num_examples = 1000
	True_w = [2.5, 3.0]
	True_b = 4.2
	features = randn(Float64, (num_examples, num_inputs))
	labels = features * True_w .+ True_b
	labels += randn(size(labels)) * 0.05
	regData = []
	for i in 1:num_examples
  		push!(regData, (features[i,:], labels[i]))
	end

end

# ╔═╡ 6293416a-feff-11ea-08f3-ab1625da01ab
x1 = reshape(features[:,1], (1, 1000))

# ╔═╡ 67ef56b2-feff-11ea-0d3a-5fcfab8cd1d8
x2 = reshape(features[:, 2], (1, 1000))

# ╔═╡ 00a23774-ff00-11ea-2e13-e96bf9905d60
#define Model
m = Chain(Dense(2,1))

# ╔═╡ 0632bd28-ff00-11ea-3f4f-214aaf6ecf73
#define Loss function
loss(x,y) = Flux.mse(m(x), y)

# ╔═╡ 1706aba0-ff00-11ea-12b7-159a1cc6ca0a
#opt
opt = Descent(0.3)

# ╔═╡ 32f1e2a8-ff00-11ea-2c18-95d2536c5068
@epochs 100 Flux.train!(loss, Flux.params(m), regData, opt)

# ╔═╡ 537447da-ff00-11ea-1c38-a30b89cc8c10
m([x1; x2])

# ╔═╡ 63707fde-ff00-11ea-2485-2b716b32b584
labels

# ╔═╡ c351b254-ff00-11ea-2b8a-677da5b21bd2
plot(x=m([x1;x2]), y=labels, Guide.xlabel("Predict Value"), Guide.ylabel("reality value"), Guide.title("Pre vs Real"))

# ╔═╡ 951be72a-ff01-11ea-1199-3dd0e9831890
md"
详细步骤：👇

```julia
batch_size = 10
#init w, b
w = rand(2)
b = rand(1)

#learning rate
lr = 0.3

#epoch
epoch = 3

linreg(X, w, b) = X*w .+ b
#define Loss
function Loss(ŷ, y)
	sum((ŷ - y) .^ 2)
end

function train_model(batch_size, features, labels, w, b, lr, epoch)
	dt_row, dt_col = size(features)
	indices = 1:dt_row
	Rond_indices = Random.shuffle(indices)
	for j in 1:epoch
		for i in range(1, dt_row, step = batch_size)
			j = Rond_indices[i:min(i+batch_size, dt_row)]
			X = features[j,:]
			y = labels[j,:]
			L = Loss(linreg(X, w, b), y)
			gs = gradient(() -> Loss(linreg(X, w, b), y), Flux.params(w, b))
			w -= lr/batch_size * gs[w]
			b -= lr/batch_size * gs[b]
		end
	end
	return w,b

end

ŵ, b̂ = train_model(batch_size, features, labels, w, b, lr, epoch)

pre_y = features * ŵ .+ b̂

print(pre_y)

```"

# ╔═╡ Cell order:
# ╟─d6e28796-fefb-11ea-1b08-bbe1e9f98682
# ╟─a3a08c06-fefc-11ea-38d3-2526f9a1ccd2
# ╟─b3f67ca0-fefc-11ea-17cb-373174529e5c
# ╟─5a0fac6a-fefd-11ea-1daa-135ba73382ac
# ╠═b438ab92-fefd-11ea-1054-632292d0965c
# ╠═924bbece-fefe-11ea-3cdf-75c5fe8eb392
# ╠═6293416a-feff-11ea-08f3-ab1625da01ab
# ╠═67ef56b2-feff-11ea-0d3a-5fcfab8cd1d8
# ╠═00a23774-ff00-11ea-2e13-e96bf9905d60
# ╠═0632bd28-ff00-11ea-3f4f-214aaf6ecf73
# ╠═1706aba0-ff00-11ea-12b7-159a1cc6ca0a
# ╠═32f1e2a8-ff00-11ea-2c18-95d2536c5068
# ╠═537447da-ff00-11ea-1c38-a30b89cc8c10
# ╠═63707fde-ff00-11ea-2485-2b716b32b584
# ╠═c351b254-ff00-11ea-2b8a-677da5b21bd2
# ╟─951be72a-ff01-11ea-1199-3dd0e9831890

ext: jl

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
using Optim

# Wavelength range [nm]
λ = 400:1000
# Angle of incidence [degrees]
θ = [5.]
# Polarisation (p->1.0, s->0.0, weighted->between 0.0 and 1.0)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Components of EMA
n1 = RIdb.air(beam.λ)
n2 = RIdb.silicon(beam.λ)
n3 = RIdb.glass(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(n1), # 1
    ModelFit(:looyenga; N=(ninc=n1, nhost=n2)), # 2
    ModelFit(:looyenga; N=(ninc=n1, nhost=n2)), # 3
    LayerTMMO(n3), # 4
]

# Set the order of the layers (build the system) with the layers
order = [1, # incident medium
         2, 3, 2, 3, 2, 3, 2, 3, # top DBR of the MC
         2, 2, # defect of the MC
         3, 2, 3, 2, 3, 2, 3, 2, # bottom DBR of the MC
         4, # substrate
]

# Absolute reflectance spectrum to fit stored in Utils
Rexp = SpectraDB.fp_spectrum(beam.λ)

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=false,
);

# Seeds for each ModelFit layer defined above (No alpha)
seed = [
    [119.0, 0.8], # layers[2]
    [76.0, 0.5], # layers[3]
]

solOptim = fit_tmm_optics(
    Reflectance(Rexp), seed, beam, layers;
    order=order, options=options, alg=SAMIN(), lb=0.5.*seed, ub=1.5.*seed,
)

# Seeds for each ModelFit layer defined above plus the alpha (with alpha)
seed2 = [
    [120, 0.8],
    [83.0, 0.5],
    [0.995], # alpha
]

solOptim2 = fit_tmm_optics(
    Reflectance(Rexp), seed2, beam, layers;
    order=order, options=options, alg=SAMIN(), lb=0.5.*seed2, ub=1.5.*seed2, alpha=true,
)


text:
#!/bin/julia
function bind_pipes()
    # UNIX named pipe for echo
    game_pipe = ARGS[1]
    while true
        s = read(`cat $game_pipe`, String)
        # println("Msg from $game_pipe : $s")
    end
end

function bind_write()
    pipe = ARGS[2]
    while true
        run(`./teste/sc.sh w $pipe`; wait=true)
        run(`echo \"w\" $pipe`)
        sleep(0.2)
    end
end

@async bind_pipes()
@async bind_write()

wait()

max_stars_repo_name: GabrielTrettel/tetris

max_issues_repo_name: GabrielTrettel/tetris

max_forks_repo_name: GabrielTrettel/tetris


text:
# Routines to calculate the azimuth, backazimuth and great circle distance
# from headers in traces.

"Earth elliposoid semi-axes in WGS84"
const earth_r_major_WGS84 = 6378137.0000
const earth_r_minor_WGS84 = 6356752.3142
"Flattening of the Earth in WGS84"
const f_WGS84 = (earth_r_major_WGS84 - earth_r_minor_WGS84)/earth_r_major_WGS84

"""
    _great_circle(lon0, lat0, lon1, lat1, f=f_WGS84) -> gcarc, az, baz

Return the great-circle distance, `gcarc`, forward azimuth `az` and backazimuth `baz`
between two points, all specified in degrees.
"""
function _great_circle(lon0, lat0, lon1, lat1, f=f_WGS84)
    lon0, lat0, lon1, lat1 = Float64(lon0), Float64(lat0), Float64(lon1), Float64(lat1)
    gcarc, az, baz = vincentydist(f, 1.0, deg2rad(lat0), deg2rad(lon0),
        deg2rad(lat1), deg2rad(lon1))
    rad2deg(gcarc), rad2deg(az), rad2deg(baz)
end

"""
    update_great_circle!(s::SACtr)

If all headers `evlo`, `evla`, `stlo` and `stla` are set, update the values of
`az`, `baz` and `gcarc`.
"""
function update_great_circle!(s::SACtr)
    any([s.evlo, s.evla, s.stlo, s.stla] .== sac_rnull) && return
    s.gcarc, s.az, s.baz = _great_circle(s.evlo, s.evla, s.stlo, s.stla)
end

# The following functions have been taken from
#   https://github.com/acrosby/GreatCircle.jl
# commit
#   c4860071103c1895a86561c0a0a772bf89705fba
#
# See the file LICENSE.md for details of the licence.
#
# -----------------------------------------------------------------------
# | Algorithms from Geocentric Datum of Australia Technical Manual |
# | |
# | http://www.anzlic.org.au/icsm/gdatum/chapter4.html |
# | |
# | This page last updated 11 May 1999 |
# | |
# | Computations on the Ellipsoid |
# | |
# | There are a number of formulae that are available |
# | to calculate accurate geodetic positions, |
# | azimuths and distances on the ellipsoid. |
# | |
# | Vincenty's formulae (Vincenty, 1975) may be used |
# | for lines ranging from a few cm to nearly 20,000 km, |
# | with millimetre accuracy. |
# | The formulae have been extensively tested |
# | for the Australian region, by comparison with results |
# | from other formulae (Rainsford, 1955 & Sodano, 1965). |
# | |
# | * Inverse problem: azimuth and distance from known |
# | latitudes and longitudes |
# | * Direct problem: Latitude and longitude from known |
# | position, azimuth and distance. |
# | * Sample data |
# | * Excel spreadsheet |
# | |
# | Vincenty's Inverse formulae |
# | Given: latitude and longitude of two points |
# | (phi1, lembda1 and phi2, lembda2), |
# | Calculate: the ellipsoidal distance (s) and |
# | forward and reverse azimuths between the points (alpha12, alpha21). |
# | |
# -----------------------------------------------------------------------
"""
    vincentypt(flattening, semimajor_radius, lat, lon, azimuth, distance) -> lat, lon, backazimuth

Return lat and long of projected point and backazimuth,
given a reference point and a distance and azimuth to project.

Coordinates and azimuth are in radians.
"""
function vincentypt(f::T, a::T, phi1::T, lembda1::T, alpha12::T, s::T) where {T<:AbstractFloat}
    two_pi = 2π
    if alpha12 < 0.0
        alpha12 = alpha12 + two_pi
    end
    if alpha12 > two_pi
        alpha12 = alpha12 - two_pi
    end

    b = a .* (1.0 - f)

    TanU1 = (1 - f) .* tan(phi1)
    U1 = atan( TanU1 )
    sigma1 = atan( TanU1, cos(alpha12) )
    Sinalpha = cos(U1) .* sin(alpha12)
    cosalpha_sq = 1.0 - Sinalpha .* Sinalpha

    u2 = cosalpha_sq .* (a .* a - b .* b ) ./ (b .* b)
    A = 1.0 + (u2 ./ 16384) .* (4096 + u2 .* (-768 + u2 .* (320 - 175 .* u2) ) )
    B = (u2 ./ 1024) .* (256 + u2 .* (-128 + u2 .* (74 - 47 .* u2) ) )

    # Starting with the approximation
    sigma = (s ./ (b .* A))

    # Not moving anywhere. We can return the location that was passed in.
    if sigma == 0
        return phi1, lembda1, alpha12
    end

    last_sigma = 2.0 .* sigma + 2.0 # something impossible

    # Iterate the following three equations
    # until there is no significant change in sigma
    # two_sigma_m , delta_sigma
    while ( abs( (last_sigma - sigma) ./ sigma) > 1.0e-9 )
        global two_sigma_m = 2 .* sigma1 + sigma
        delta_sigma = B .* sin(sigma) .* ( cos(two_sigma_m) + (B./4) .* (cos(sigma) .* (-1 + 2 .* cos(two_sigma_m).^2 - (B./6) .* cos(two_sigma_m) .* (-3 + 4 .* sin(sigma).^2) .* (-3 + 4 .* cos(two_sigma_m).^2 ))))
        last_sigma = sigma
        sigma = (s ./ (b .* A)) + delta_sigma
    end

    phi2 = atan( (sin(U1) .* cos(sigma) + cos(U1) .* sin(sigma) .* cos(alpha12) ), ((1-f) .* sqrt( Sinalpha.^2 + (sin(U1) .* sin(sigma) - cos(U1) .* cos(sigma) .* cos(alpha12)).^2)))

    lembda = atan( (sin(sigma) .* sin(alpha12 )), (cos(U1) .* cos(sigma) - sin(U1) .* sin(sigma) .* cos(alpha12)))

    C = (f./16) .* cosalpha_sq .* (4 + f .* (4 - 3 .* cosalpha_sq ))

    omega = lembda - (1-C) .* f .* Sinalpha .* (sigma + C .* sin(sigma) .* (cos(two_sigma_m) + C .* cos(sigma) .* (-1 + 2 .* cos(two_sigma_m).^2 )))

    lembda2 = lembda1 + omega

    alpha21 = atan( Sinalpha, (-sin(U1) .* sin(sigma) + cos(U1) .* cos(sigma) .* cos(alpha12)))
    alpha21 = alpha21 + pi

    if alpha21 < 0.0
        alpha21 = alpha21 + two_pi
    end
    if alpha21 > two_pi
        alpha21 = alpha21 - two_pi
    end
    return phi2, lembda2, alpha21

end

function vincentypt(f::T, a::T, phi1::Array{T,1}, lembda1::Array{T,1}, alpha12::Array{T,1}, s::Array{T,1}) where {T<:AbstractFloat}
    lenphi = length(phi1)
    @assert lenphi == length(lembda1)
    @assert lenphi == length(alpha12)
    @assert lenphi == length(s)
    phi2 = zeros(phi1)
    lembda2 = zeros(phi1)
    alpha21 = zeros(phi1)
    for i=1:lenphi
        @inbounds phi2[i], lembda2[i], alpha21[i] = vincentypt(f, a, phi1[i], lembda1[i], alpha12[i], s[i])
    end
    return phi2, lembda2, alpha21
end

vincentypt(f::T, a::T, phi1::Array{T,1}, lembda1::Array{T,1}, alpha12::Array{T,1}, s::T) where {T<:AbstractFloat} = vincentypt(f, a, phi1, lembda1, alpha12, ones(length(phi1)) * s)
vincentypt(f::T, a::T, phi1::Array{T,1}, lembda1::Array{T,1}, alpha12::T, s::Array{T,1}) where {T<:AbstractFloat} = vincentypt(f, a, phi1, lembda1, ones(length(phi1)) * alpha12, s)
vincentypt(f::T, a::T, phi1::Array{T,1}, lembda1::Array{T,1}, alpha12::T, s::T) where {T<:AbstractFloat} = vincentypt(f, a, phi1, lembda1, ones(length(phi1)) * alpha12, ones(length(phi1)) * s)

"""
    vincentydist(flattening, semimajor_radius, lat1, lon1, lat2, lon2) -> dist, azimuth, backazimuth

Returns the distance between two geographic points on the ellipsoid,
and the forward and reverse azimuths between these points.

Latitude, longitude and (back)azimuth are in radians; distances are in the units
provided for the semimajor axis length.
"""
function vincentydist(f::T, a::T, phi1::T, lembda1::T, phi2::T, lembda2::T) where {T<:AbstractFloat}
    two_pi = 2π
    if (abs( phi2 - phi1 ) < 1e-8) & ( abs( lembda2 - lembda1) < 1e-8 )
        return 0.0, 0.0, 0.0
    end

    b = a .* (1.0 - f)

    TanU1 = (1 - f) .* tan( phi1 )
    TanU2 = (1 - f) .* tan( phi2 )

    U1 = atan(TanU1)
    U2 = atan(TanU2)

    lembda = lembda2 - lembda1
    last_lembda = -4000000.0 # an impossibe value
    omega = lembda

    # Iterate the following equations until there is no significant change in lembda
    alpha, sigma, Sin_sigma, Cos2sigma_m, Cos_sigma, sqr_sin_sigma = -999999., -999999., -999999., -999999., -999999., -999999.
    while ( (last_lembda < -3000000.0) | (lembda != 0) ) & ( abs( (last_lembda - lembda)./lembda) > 1.0e-9 )
        sqr_sin_sigma = (cos(U2) .* sin(lembda)).^2 + ( (cos(U1) .* sin(U2) - sin(U1) .* cos(U2) .* cos(lembda) )).^2
        Sin_sigma = sqrt( sqr_sin_sigma )
        Cos_sigma = sin(U1) .* sin(U2) + cos(U1) .* cos(U2) .* cos(lembda)
        sigma = atan( Sin_sigma, Cos_sigma )

        Sin_alpha = cos(U1) .* cos(U2) .* sin(lembda) ./ sin(sigma)

        if (Sin_alpha >= 1)# & (Sin_alpha .== 1.0)
            Sin_alpha = 1.0
        elseif (Sin_alpha <= -1)# & (Sin_alpha .== -1.0)
            Sin_alpha = -1.0
        end

        alpha = asin( Sin_alpha )
        Cos2sigma_m = cos(sigma) - (2 .* sin(U1) .* sin(U2) ./ cos(alpha).^2 )
        C = (f./16) .* cos(alpha).^2 .* (4 + f .* (4 - 3 .* cos(alpha).^2))
        last_lembda = lembda
        lembda = omega + (1-C) .* f .* sin(alpha) .* (sigma + C .* sin(sigma) .* (Cos2sigma_m + C .* cos(sigma) .* (-1 + 2 .* Cos2sigma_m.^2 )))
    end

    u2 = cos(alpha).^2 .* (a.*a-b.*b) ./ (b.*b)

    A = 1 + (u2./16384) .* (4096 + u2 .* (-768 + u2 .* (320 - 175 .* u2)))

    B = (u2./1024) .* (256 + u2 .* (-128+ u2 .* (74 - 47 .* u2)))

    delta_sigma = B .* Sin_sigma .* (Cos2sigma_m + (B./4) .* (Cos_sigma .* (-1 + 2 .* Cos2sigma_m.^2 ) - (B./6) .* Cos2sigma_m .* (-3 + 4 .* sqr_sin_sigma) .* (-3 + 4 .* Cos2sigma_m.^2 )))

    s = b .* A .* (sigma - delta_sigma)

    alpha12 = atan( (cos(U2) .* sin(lembda)), (cos(U1) .* sin(U2) - sin(U1) .* cos(U2) .* cos(lembda)))
    alpha21 = atan( (cos(U1) .* sin(lembda)), (-sin(U1) .* cos(U2) + cos(U1) .* sin(U2) .* cos(lembda)))

    if alpha12 < 0.0
        alpha12 = alpha12 + two_pi
    end
    if alpha12 > two_pi
        alpha12 = alpha12 - two_pi
    end
    alpha21 = alpha21 + two_pi ./ 2.0

    if alpha21 < 0.0
        alpha21 = alpha21 + two_pi
    end
    if alpha21 > two_pi
        alpha21 = alpha21 - two_pi
    end
    return s, alpha12, alpha21
end

function vincentydist(f::T, a::T, phi1::Array{T,1}, lembda1::Array{T,1}, phi2::Array{T,1}, lembda2::Array{T,1}) where {T<:AbstractFloat}
    lenphi = length(phi1)
    @assert lenphi == length(lembda1)
    @assert lenphi == length(phi2)
    @assert lenphi == length(lembda2)
    s = zeros(phi1)
    alpha12 = zeros(phi1)
    alpha21 = zeros(phi1)
    for i=1:lenphi
        @inbounds s[i], alpha12[i], alpha21[i] = vincentydist(f, a, phi1[i], lembda1[i], phi2[i], lembda2[i])
    end
    return s, alpha12, alpha21
end

vincentydist(f::T, a::T, phi1::T, lembda1::Array{T,1}, phi2::T, lembda2::Array{T,1}) where {T<:AbstractFloat} = vincentydist(f, a, phi1 * ones(lembda1), lembda1, phi2 * ones(lembda1), lembda2)
vincentydist(f::T, a::T, phi1::Array{T,1}, lembda1::T, phi2::Array{T,1}, lembda2::T) where {T<:AbstractFloat} = vincentydist(f, a, phi1, lembda1 * ones(phi1), phi2, lembda2 * ones(phi1))


text:
#L1DistLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

l1fn = "L(r)=r"
l1dist = latexify(l1fn)
Plots.plot(L1DistLoss(), c=:olive, linewidth=2, annotation=(1.5, .4, L"L\left( r \right) = |r|", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
include("drop_clock.jl")
include("lock_clock.jl")
include("hold_clock.jl")
include("push_clock.jl")
include("drag_clock.jl")


text:
# SymbolicAlgebra.jl

# [[file:~/Documents/Julia/scrap.org::*SymbolicAlgebra.jl][SymbolicAlgebra.jl:1]]
#_____________________________________________
# Addition
Base.:(+)(x::T, y::T) where {T<:Symbolic} = promote(T)(:+, stripiden.([x, y])) |> simplify
Base.:+(x::Symbolic) = x #|> simplify

#_____________________________________________
# Subtraction
Base.:(-)(x::T, y::T) where {T<:Symbolic} = promote(T)(:+, stripiden.([x, -y]))  |> simplify
Base.:(-)(x::T) where {T<:Symbolic} = promote(T)(:*, [-1, x])  |> simplify

#_____________________________________________
# Multiplication
Base.:(*)(x::T, y::T) where {T<:Symbolic} = ((x == y) ? promote(T)(:^, stripiden.([x, 2])) :
                                                        promote(T)(:*, stripiden.([x, y])))  |> simplify
LinearAlgebra.dot(a::T, b::T) where {T<:Symbolic} = a*b
Base.zero(a::Symbolic) = 0

#_____________________________________________
# Division
Base.:(/)(x::T, y::T) where {T<:Symbolic} = promote(T)(:*, stripiden.([x, y^-1]))  |> simplify

#_____________________________________________
# Exponents
Base.:(^)(x::T, y::T) where {T<:Symbolic}   = promote(T)(:^, stripiden.([x, y]))  |> simplify
Base.:(^)(x::T, y::Int) where {T<:Symbolic} = promote(T)(:^, stripiden.([x, y]))  |> simplify

#_____________________________________________
# other
Base.inv(x::T) where {T<:Symbolic} = promote(T)(:^, stripiden.([x, -1]))  |> simplify

Base.:( \ )(x::T,y::T) where {T<:Symbolic} = inv(x)*y

Base.abs(x::T) where {T<:Symbolic} = sqrt(x^2)

Base.conj(x::Union{AbstractSymExpr,AbstractSym}) = x

SymNum = Union{Symbolic,Number}

function Base.atan(x::T, y::V) where {T<:SymNum,V<:SymNum}
   promote_SymForm(x,y)(:atan, stripiden.([x,y]))
end

function Base.hypot(x::T, y::V) where {T<:Symbolic,V<:Symbolic}
   promote_SymForm(x,y)(:hypot, stripiden.([x,y]))
end
function Base.hypot(x::T, y::V) where {T<:Number,V<:Symbolic}
   promote_SymForm(x,y)(:hypot, stripiden.([x,y]))
end
function Base.hypot(x::T, y::V) where {T<:Symbolic,V<:Number}
   promote_SymForm(x,y)(:hypot, stripiden.([x,y]))
end

function Base.max(x::T, y::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(x,y)(:max, stripiden.([x,y]))
end

function Base.min(x::T, y::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(x,y)(:min, stripiden.([x,y]))
end

function Base.:<(x::T, y::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(x,y)(:<, stripiden.([x,y]))
end

function SpecialFunctions.besselj(ν::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(ν,x)(:besselj, stripiden.([ν, x]))
end

function SpecialFunctions.besseli(ν::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(ν,x)(:besseli, stripiden.([ν, x]))
end

function SpecialFunctions.bessely(ν::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(ν,x)(:bessely, stripiden.([ν, x]))
end

function SpecialFunctions.besselk(ν::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(ν,x)(:besselk, stripiden.([ν, x]))
end

function SpecialFunctions.hankelh1(ν::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(ν,x)(:hankelh1, stripiden.([ν, x]))
end

function SpecialFunctions.hankelh2(ν::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(ν,x)(:hankelh2, stripiden.([ν, x]))
end

function SpecialFunctions.polygamma(m::T, x::V) where {T<:SymNum,V<:SymNum}
    promote_SymForm(m,x)(:polygamma, stripiden.([m, x]))
end
function SpecialFunctions.polygamma(m::Int, x::V) where {V<:SymNum}
    promote_SymForm(m,x)(:polygamma, stripiden.([m, x]))
end

function SpecialFunctions.beta(a::T, b::V) where {T<:Number,V<:Symbolic}
    promote_SymForm(a,b)(:beta, stripiden.([a, b]))
end
function SpecialFunctions.beta(a::T, b::V) where {T<:Symbolic,V<:Number}
    promote_SymForm(a,b)(:beta, stripiden.([a, b]))
end
function SpecialFunctions.beta(a::T, b::V) where {T<:Symbolic,V<:Symbolic}
    promote_SymForm(a,b)(:beta, stripiden.([a, b]))
end

function SpecialFunctions.lbeta(a::T, b::V) where {T<:Number,V<:Symbolic}
    promote_SymForm(a,b)(:lbeta, stripiden.([a, b]))
end
function SpecialFunctions.lbeta(a::T, b::V) where {T<:Symbolic,V<:Number}
    promote_SymForm(a,b)(:lbeta, stripiden.([a, b]))
end
function SpecialFunctions.lbeta(a::T, b::V) where {T<:Symbolic,V<:Symbolic}
    promote_SymForm(a,b)(:lbeta, stripiden.([a, b]))
end

#_____________________________________________
# More math functions
for (M, f, arity) in DiffRules.diffrules()
    if arity == 1 && (M == :Base || M == :SpecialFunctions) && f ∉ [:inv, :+, :-, :abs]
        @eval begin
            $M.$f(x::T) where {T<:Symbolic} = promote(T)(Symbol($f), stripiden.([x]))  |> simplify
        end
    end
end

# This is different from a possible `isnegative` in that `-x` could be
# positive, if `x` is negative. This function returns `true` if there
# is a factor `-1` present.
isnegated(::Number) = false
isnegated(r::Real) = r < 0
isnegated(s::SymExpr) = s.op.name == :(*) && findfirst(isequal(-1), s.args) !== nothing

isdenominator(a::Number) = false
isdenominator(a::SymExpr) = (a.op.name == :(^) && isnegated(a.args[2]))

Base.numerator(s::Sym) = s
Base.denominator(s::Sym) = 1

function Base.numerator(expr::SymExpr)
    isdenominator(expr) && return 1
    if expr.op.name == :(*)
        args = filter(a -> !isdenominator(a), expr.args)
        length(args) > 1 ? prod(args) : args[1]
    else
        expr
    end
end

function Base.denominator(expr::SymExpr)
    # # This would be the simplest implemenation, if simplification
    # # worked (e.g. inv(xy/sin(x)) becomes (x * y * sin(x) ^ -1) ^ -1
    # # at the moment):
    # return numerator(inv(expr))
    if expr.op.name == :(*)
        den = filter(isdenominator, expr.args)
        isempty(den) && return 1
        # Negate exponents of all factors that should be in the
        # denominator.
        args = map(inv, den)#  do a
        #     a.args[1] ^ (-a.args[2])
        # end
        length(args) > 1 ? prod(args) : args[1]
    elseif expr.op.name == :(^)
        numerator(inv(expr))
    else
        1
    end
end

# SymbolicAlgebra.jl:1 ends here


text:
using KernelAbstractions, CUDA, SimpleDiffEq
using StaticArrays

function loop(u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du1 = σ*(u[2]-u[1])
    du2 = u[1]*(ρ-u[3]) - u[2]
    du3 = u[1]*u[2] - β*u[3]
    return SVector{3}(du1, du2, du3)
end
const func = ODEFunction(loop)
u0 = 10ones(Float32,3)
const su0= SVector{3}(u0)
const dt = 1f-1
const tspan = (0.0f0, 10.0f0)

const odeoop = ODEProblem{false}(loop, SVector{3}(u0), (0.0f0, 10.0f0),  Float32[10, 28, 8/3])
sol2 = solve(odeoop,GPUSimpleTsit5(),dt=dt)
cps = Array([@SVector [10f0,28f0,8/3f0] for i in 1:32])
ps = CuArray([@SVector [10f0,28f0,8/3f0] for i in 1:32])
CUDA.allowscalar(false)

function f(p)
    solve(odeoop,GPUSimpleTsit5(),dt=dt)[end]
end

map(f,cps)

@kernel function f_kernel(ps)
    solve(odeoop,GPUSimpleTsit5(),dt=dt)[end]
end

version = CUDADevice()
wait(version, f_kernel(version)(ps;ndrange=length(ps)))

function DiffEqBase.solve(prob::ODEProblem,
                          alg::GPUSimpleTsit5;
                          dt = 0.1f0)
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p
  t = tspan[1]
  tf = prob.tspan[2]
  #ts = tspan[1]:dt:tspan[2]
  us = MVector{101,typeof(u0)}(undef)
  us[1] = u0
  u = u0
  k7 = f(u, p, t)

  cs, as, btildes, rs = SimpleDiffEq._build_atsit5_caches(eltype(u0))
  c1, c2, c3, c4, c5, c6 = cs
  a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
  a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
  btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = btildes

  # FSAL
  for i in 2:101
      uprev = u; k1 = k7; t = tspan[1] + dt * i
      tmp = uprev+dt*a21*k1
      k2 = f(tmp, p, t+c1*dt)
      tmp = uprev+dt*(a31*k1+a32*k2)
      k3 = f(tmp, p, t+c2*dt)
      tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
      k4 = f(tmp, p, t+c3*dt)
      tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
      k5 = f(tmp, p, t+c4*dt)
      tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
      k6 = f(tmp, p, t+dt)
      u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
      k7 = f(u, p, t+dt)
      us[i] = u
  end

  SArray(us)
end

# Works now!
function f(p)
    solve(odeoop,GPUSimpleTsit5(),dt=dt)[end]
end

map(f,cps)

@kernel function f_kernel(ps)
    solve(odeoop,GPUSimpleTsit5(),dt=dt)[end]
end

version = CUDADevice()
wait(version, f_kernel(version)(ps;ndrange=length(ps)))

# Fails to change parameters
function f(p)
    remake(odeoop,p=p)
    solve(odeoop,GPUSimpleTsit5(),dt=dt)[end]
end

map(f,cps)

@kernel function f_kernel(ps)
    i = @index(Global)
    remake(odeoop,p=ps[i])
    solve(odeoop,GPUSimpleTsit5(),dt=dt)[end]
end

version = CUDADevice()
wait(version, f_kernel(version)(ps;ndrange=length(ps)))


text:
@testset "Test lagrange_to_state and state_to_lagrange" begin

    atol = 1000*eps()

    sensors = -11.0:0.5:10.0

    Nv = 10
    Nx = 3*Nv
    U = randn(ComplexF64)
    # Define the state
    x = rand(Nx)
    x0 = deepcopy(x)

    config = let Nv = Nv,
                 U = U,
                 ss = sensors, Δt = 5e-3, δ = 1e-1,
                 ϵX = 1e-3, ϵΓ = 1e-3,
                 β = 1.0,
                 ϵY = 1e-2
        VortexConfig(Nv, U, ss, Δt, δ, ϵX, ϵΓ, β, ϵY)
    end

    blobs = state_to_lagrange(x0, config)
    blobs0 = deepcopy(blobs)
    for i=1:Nv
        bi₊ = blobs[1][i]
        bi₋ = blobs[2][i]

        @test isapprox(bi₋.z, conj(bi₊.z), atol = atol)
        @test isapprox(bi₋.S, -bi₊.S, atol = atol)
    end

    xbis = lagrange_to_state(blobs0, config)

    @test isapprox(xbis, x, atol = atol)
end


text:
# Base32 Decoder
# ==============

struct Base32Decoder <: Codec
    table::CodeTable32
    state::State
    buffer::Buffer
end

Base32Decoder(table::CodeTable32) = Base32Decoder(table, State(), Buffer(7))

"""
    Base32Decoder(;hex::Bool=false, ignore::String=$(escape_string(whitespace)))

Create a base32 decoding codec.

Arguments
---------
- `hex`: use extended hex alphabet (Table 4 of RFC4648)
- `ignore`: ASCII characters that will be ignored while decoding
"""
function Base32Decoder(;hex::Bool=false, ignore::String=whitespace)
    if hex
        table = BASE32_HEX
    else
        table = BASE32_STD
    end
    table = copy(table)
    ignorechars!(table, ignore)
    return Base32Decoder(table)
end

const Base32DecoderStream{S} = TranscodingStream{Base32Decoder,S} where S<:IO

"""
    Base32DecoderStream(stream::IO; kwargs...)

Create a base32 decoding stream (see `Base32Decoder` for `kwargs`).
"""
function Base32DecoderStream(stream::IO; kwargs...)
    return TranscodingStream(Base32Decoder(;kwargs...), stream)
end

function TranscodingStreams.startproc(
        codec :: Base32Decoder,
        state :: Symbol,
        error :: Error)
    start!(codec.state)
    return :ok
end

function TranscodingStreams.process(
        codec  :: Base32Decoder,
        input  :: Memory,
        output :: Memory,
        error  :: Error)
    table = codec.table
    state = codec.state
    buffer = codec.buffer

    # Check if we can encode data.
    if !is_running(state)
        error[] = ArgumentError("decoding is already finished")
        return 0, 0, :error
    elseif output.size < 3
        # Need more output space.
        return 0, 0, :ok
    end

    # Load the first bytes.
    i = j = 0
    while buffer.size < 7 && i + 1 ≤ input.size
        buffer[buffer.size+=1] = input[i+=1]
    end
    c1 = c2 = c3 = c4 = c5 = c6 = c7 = c8 = BASE32_CODEIGN
    if buffer.size ≥ 1
        c1 = decode(table, buffer[1])
    end
    if buffer.size ≥ 2
        c2 = decode(table, buffer[2])
    end
    if buffer.size ≥ 3
        c3 = decode(table, buffer[3])
    end
    if buffer.size ≥ 4
        c4 = decode(table, buffer[4])
    end
    if buffer.size ≥ 5
        c5 = decode(table, buffer[5])
    end
    if buffer.size ≥ 6
        c6 = decode(table, buffer[6])
    end
    if buffer.size ≥ 7
        c7 = decode(table, buffer[7])
    end
    empty!(buffer)

    # Start decoding loop.
    #  01234567 89012345 67890123 45678901 23456789
    # +--------+--------+--------+--------+--------+
    # |< 1 >< 2| >< 3 ><| 4 >< 5 |>< 6 >< |7 >< 8 >|
    # +--------+--------+--------+--------+--------+
    status = :ok
    @inbounds while true
        if c1 > 0x1f || c2 > 0x1f || c3 > 0x1f || c4 > 0x1f ||
           c5 > 0x1f || c6 > 0x1f || c7 > 0x1f || c8 > 0x1f
           i, j, status = decode32_irregular(
               table,
               c1, c2, c3, c4, c5, c6, c7, c8,
               input,  i,
               output, j,
               error)
        else
            output[j+1] = c1 << 3 | c2 >> 2
            output[j+2] = c2 << 6 | c3 << 1 | c4 >> 4
            output[j+3] = c4 << 4 | c5 >> 1
            output[j+4] = c5 << 7 | c6 << 2 | c7 >> 3
            output[j+5] = c7 << 5 | c8
            j += 5
        end
        if i + 8 ≤ input.size && j + 5 ≤ output.size && status == :ok
            c1 = decode(table, input[i+1])
            c2 = decode(table, input[i+2])
            c3 = decode(table, input[i+3])
            c4 = decode(table, input[i+4])
            c5 = decode(table, input[i+5])
            c6 = decode(table, input[i+6])
            c7 = decode(table, input[i+7])
            c8 = decode(table, input[i+8])
            i += 8
        else
            break
        end
    end

    # Epilogue.
    if status == :end || status == :error
        finish!(state)
    end
    return i, j, status
end

function decode32_irregular(table, c1, c2, c3, c4, c5, c6, c7, c8, input, i, output, j, error)
    # Skip ignored chars.
    while true
        if c1 == BASE32_CODEIGN
            c1, c2, c3, c4, c5, c6, c7 = c2, c3, c4, c5, c6, c7, c8
        elseif c2 == BASE32_CODEIGN
            c2, c3, c4, c5, c6, c7 = c3, c4, c5, c6, c7, c8
        elseif c3 == BASE32_CODEIGN
            c3, c4, c5, c6, c7 = c4, c5, c6, c7, c8
        elseif c4 == BASE32_CODEIGN
            c4, c5, c6, c7 = c5, c6, c7, c8
        elseif c5 == BASE32_CODEIGN
            c5, c6, c7 = c6, c7, c8
        elseif c6 == BASE32_CODEIGN
            c6, c7 = c7, c8
        elseif c7 == BASE32_CODEIGN
            c7 = c8
        elseif c8 == BASE32_CODEIGN
            # pass
        else
            break
        end
        if i + 1 ≤ input.size
            c8 = decode(table, input[i+=1])
        else
            c8 = BASE32_CODEEND
            break
        end
    end

    # Write output.
    if max(c1, c2, c3, c4, c5, c6, c7, c8) ≤ 0x1f
        output[j+=1] = c1 << 3 | c2 >> 2
        output[j+=1] = c2 << 6 | c3 << 1 | c4 >> 4
        output[j+=1] = c4 << 4 | c5 >> 1
        output[j+=1] = c5 << 7 | c6 << 2 | c7 >> 3
        output[j+=1] = c7 << 5 | c8
        status = :ok
    elseif max(c1, c2, c3, c4, c5, c6, c7) ≤ 0x1f && c8 == BASE32_CODEPAD
        output[j+=1] = c1 << 3 | c2 >> 2
        output[j+=1] = c2 << 6 | c3 << 1 | c4 >> 4
        output[j+=1] = c4 << 4 | c5 >> 1
        output[j+=1] = c5 << 7 | c6 << 2 | c7 >> 3
        status = :end
    elseif max(c1, c2, c3, c4, c5) ≤ 0x1f && c6 == c7 == c8 == BASE32_CODEPAD
        output[j+=1] = c1 << 3 | c2 >> 2
        output[j+=1] = c2 << 6 | c3 << 1 | c4 >> 4
        output[j+=1] = c4 << 4 | c5 >> 1
        status = :end
    elseif max(c1, c2, c3, c4) ≤ 0x1f && c5 == c6 == c7 == c8 == BASE32_CODEPAD
        output[j+=1] = c1 << 3 | c2 >> 2
        output[j+=1] = c2 << 6 | c3 << 1 | c4 >> 4
        status = :end
    elseif max(c1, c2) ≤ 0x1f && c3 == c4 == c5 == c6 == c7 == c8 == BASE32_CODEPAD
        output[j+=1] = c1 << 3 | c2 >> 2
        status = :end
    elseif c1 == c2 == c3 == c4 == c5 == c6 == c7 == BASE32_CODEIGN && c8 == BASE32_CODEEND
        status = :end
    else
        error[] = DecodeError("invalid base32 data")
        status = :error
    end

    return i, j, status
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const VisionFeaturePrint_Scene_SceneVersion = (;[
    Symbol("SCENE_VERSION_INVALID") => Int32(0),
    Symbol("SCENE_VERSION_1") => Int32(1),
]...)

mutable struct VisionFeaturePrint_Scene <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function VisionFeaturePrint_Scene(; kwargs...)
        obj = new(meta(VisionFeaturePrint_Scene), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct VisionFeaturePrint_Scene
const __meta_VisionFeaturePrint_Scene = Ref{ProtoMeta}()
function meta(::Type{VisionFeaturePrint_Scene})
    ProtoBuf.metalock() do
        if !isassigned(__meta_VisionFeaturePrint_Scene)
            __meta_VisionFeaturePrint_Scene[] = target = ProtoMeta(VisionFeaturePrint_Scene)
            allflds = Pair{Symbol,Union{Type,String}}[:version => Int32]
            meta(target, VisionFeaturePrint_Scene, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_VisionFeaturePrint_Scene[]
    end
end
function Base.getproperty(obj::VisionFeaturePrint_Scene, name::Symbol)
    if name === :version
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

const VisionFeaturePrint_Objects_ObjectsVersion = (;[
    Symbol("OBJECTS_VERSION_INVALID") => Int32(0),
    Symbol("OBJECTS_VERSION_1") => Int32(1),
]...)

mutable struct VisionFeaturePrint_Objects <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function VisionFeaturePrint_Objects(; kwargs...)
        obj = new(meta(VisionFeaturePrint_Objects), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct VisionFeaturePrint_Objects
const __meta_VisionFeaturePrint_Objects = Ref{ProtoMeta}()
function meta(::Type{VisionFeaturePrint_Objects})
    ProtoBuf.metalock() do
        if !isassigned(__meta_VisionFeaturePrint_Objects)
            __meta_VisionFeaturePrint_Objects[] = target = ProtoMeta(VisionFeaturePrint_Objects)
            fnum = Int[1,100]
            allflds = Pair{Symbol,Union{Type,String}}[:version => Int32, :output => Base.Vector{AbstractString}]
            meta(target, VisionFeaturePrint_Objects, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_VisionFeaturePrint_Objects[]
    end
end
function Base.getproperty(obj::VisionFeaturePrint_Objects, name::Symbol)
    if name === :version
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :output
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    else
        getfield(obj, name)
    end
end

mutable struct VisionFeaturePrint <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function VisionFeaturePrint(; kwargs...)
        obj = new(meta(VisionFeaturePrint), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct VisionFeaturePrint
const __meta_VisionFeaturePrint = Ref{ProtoMeta}()
function meta(::Type{VisionFeaturePrint})
    ProtoBuf.metalock() do
        if !isassigned(__meta_VisionFeaturePrint)
            __meta_VisionFeaturePrint[] = target = ProtoMeta(VisionFeaturePrint)
            fnum = Int[20,21]
            allflds = Pair{Symbol,Union{Type,String}}[:scene => VisionFeaturePrint_Scene, :objects => VisionFeaturePrint_Objects]
            oneofs = Int[1,1]
            oneof_names = Symbol[Symbol("VisionFeaturePrintType")]
            meta(target, VisionFeaturePrint, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_VisionFeaturePrint[]
    end
end
function Base.getproperty(obj::VisionFeaturePrint, name::Symbol)
    if name === :scene
        return (obj.__protobuf_jl_internal_values[name])::VisionFeaturePrint_Scene
    elseif name === :objects
        return (obj.__protobuf_jl_internal_values[name])::VisionFeaturePrint_Objects
    else
        getfield(obj, name)
    end
end

export VisionFeaturePrint_Scene_SceneVersion, VisionFeaturePrint_Scene, VisionFeaturePrint_Objects_ObjectsVersion, VisionFeaturePrint_Objects, VisionFeaturePrint


text:
using Test
using RigidBodyDynamics
using StrandbeestRobot
using NBInclude

@testset "Initial configuration" begin
    @testset "floating base" begin
        mechanism = StrandbeestRobot.mechanism(floating=true)
        state = MechanismState(mechanism)
        cost = StrandbeestRobot.loop_joint_error(mechanism)
        StrandbeestRobot.solve_initial_configuration!(state)
        for joint in tree_joints(mechanism)
            @test RigidBodyDynamics.is_configuration_normalized(joint, configuration(state, joint); atol=1e-8)
        end
        @test cost(configuration(state)) <= 1e-9
    end

    @testset "fixed base" begin
        mechanism = StrandbeestRobot.mechanism(floating=false)
        state = MechanismState(mechanism)
        cost = StrandbeestRobot.loop_joint_error(mechanism)
        StrandbeestRobot.solve_initial_configuration!(state)
        @test cost(configuration(state)) <= 1e-9
    end
end

@testset "notebooks" begin
    notebook_dir = joinpath(@__DIR__, "..", "notebooks")
    @nbinclude(joinpath(notebook_dir, "Basic simulation.ipynb"); regex = r"^((?!\#NBSKIP).)*$"s)
    @nbinclude(joinpath(notebook_dir, "Passive walking.ipynb"); regex = r"^((?!\#NBSKIP).)*$"s)
    @nbinclude(joinpath(notebook_dir, "Powered walking.ipynb"); regex = r"^((?!\#NBSKIP).)*$"s)
end


text:
"""
bretonnet_silbert.jl

Bretonnet-Silbert pseudopotential
"""

immutable BretonnetSilbert <: PseudoPotential
  z::Vector{Float64} # valence of core ion
  rc::Vector{Float64} # pseudopotential core radius
  a::Vector{Float64} # softness parameters
end

function BretonnetSilbert(z::Real, rc::Real, a::Real)
  BretonnetSilbert([z], [rc], [a])
end

function formfactor(nfe::NFE{BretonnetSilbert})::Vector{Function}
  ρ::Float64 = nfe.ρ

  N::Int = length(nfe.pseudo.z)
  ret = Vector{Function}(N)

  for i in 1:N
    z::Float64 = nfe.pseudo.z[i]
    rc::Float64 = nfe.pseudo.rc[i]
    a::Float64 = nfe.pseudo.a[i]

    B₁::Float64 = (z/rc) * (1 - 2a/rc) * exp(rc/a)
    B₂::Float64 = (2z/rc) * (a/rc - 1) * exp(0.5rc/a)

    function ω(q)::Float64
      J₁(q) = 2 - exp(-rc/a) * (
      (rc * (1 + a^2 * q^2) / a + 1 - a^2 * q^2) * sin(q*rc) / (a*q) +
      (2 + rc * (1 + a^2 * q^2) / a) * cos(q*rc) )

      J₂(q) = 2 - exp(-rc/2a) * (
      (rc * (1 + 4 * a^2 * q^2) / 2a + 1 - 4 * a^2 * q^2) * sin(q*rc) / (2a*q) +
      (2 + rc * (1 + 4 * a^2 * q^2) / 2a) * cos(q*rc) )

      4π * ρ * a^3 * (
      B₁ * J₁(q) / (1 + a^2 * q^2)^2 + 8B₂ * J₂(q) / (1 + 4 * a^2 * q^2)^2) -
      (4π * z * ρ / q^2) * cos(q*rc)
    end

    ret[i] = ω
  end

  return ret
end


text:
function apply_kernel!(img::ImgArr, kernel::Array{Int, 2})
	c, w, h = size(img)
	kh, kw = size(kernel)
	x_off = fld(kw - 1, 2)
	y_off = fld(kh - 1, 2)
	cpyImg::ImgArr = ImgArr(undef, (3, w, h))
	for i in CartesianIndices((w, h))
		x, y = Tuple(i)
		acc = Int[0, 0, 0]
		kersum::Int = 0
		for j in CartesianIndices((max(1, x - x_off) : min(w, x + x_off), max(1, y - y_off) : min(h, y + y_off)))
			kx, ky = Tuple(j)
			curker = kernel[ky + y_off - y + 1, kx + x_off - x + 1]
			kersum += abs(Int(curker))
			acc .+= img[1:3, kx, ky] * curker
		end
		if kersum > 0
			map!(acc -> fld.(acc, kersum), acc, acc)
		end
		cpyImg[:, x, y] .= abs.(acc)
	end
	img[1:3, :, :] .= cpyImg
end

function apply_kernel!(img::ImgFlt, kernel::Array{Float32, 2})
	c, w, h = size(img)
	kh, kw = size(kernel)
	x_off = fld(kw - 1, 2)
	y_off = fld(kh - 1, 2)
	cpyImg::ImgFlt = ImgFlt(undef, (3, w, h))
	acc::Array{Float32, 1} = Float32[0, 0, 0]
	kerInd = (- x_off : x_off, - y_off : y_off)
	for x in 1 : w
		if x % 100 == 0
			println(x)
		end
		for y in 1 : h
			acc .= 0
			for kx in kerInd[1]
				for ky in kerInd[2]
					acc .+= img[:, clamp(kx + x, 1, w), clamp(ky + y, 1, h)] * kernel[ky + y_off + 1, kx + x_off + 1]
				end
			end
			cpyImg[:, x, y] .= acc
		end
	end
	img .= cpyImg
end

function apply_kernel!(img::ImgFlt, kernel::Array{Float32, 3})
	c, w, h = size(img)
	kh, kw, d = size(kernel)
	x_off = fld(kw - 1, 2)
	y_off = fld(kh - 1, 2)
	cpyImg::ImgFlt = ImgFlt(undef, (3, w, h))
	acc::Array{Float32, 1} = Float32[0, 0, 0]
	kerInd = (- x_off : x_off, - y_off : y_off)
	for x in 1 : w
		for y in 1 : h
			acc .= 0
			for kx in kerInd[1]
				for ky in kerInd[2]
					acc[1] += img[1, clamp(kx + x, 1, w), clamp(ky + y, 1, h)] * kernel[ky + y_off + 1, kx + x_off + 1, 1]
					acc[2] += img[2, clamp(kx + x, 1, w), clamp(ky + y, 1, h)] * kernel[ky + y_off + 1, kx + x_off + 1, 2]
					acc[3] += img[3, clamp(kx + x, 1, w), clamp(ky + y, 1, h)] * kernel[ky + y_off + 1, kx + x_off + 1, 3]
				end
			end
			cpyImg[:, x, y] .= acc
		end
	end
	img .= cpyImg
end


text:
"""
Abstract type for forecasts that are stored in a system.
Users never create them or get access to them.
Stores references to time series data, so a disk read may be required for access.
"""
abstract type ForecastInternal <: InfrastructureSystemsType end

"""
Abstract type for forecasts supplied to users. They are not stored in a system. Instead,
they are generated on demand for the user.
Users can create them. The system will convert them to a subtype of ForecastInternal for
storage.
Time series data is stored as a field, so reads will always be from memory.
"""
abstract type Forecast <: Any end

get_horizon(forecast::Forecast) = length(get_data(forecast))
get_initial_time(forecast::Forecast) = TimeSeries.timestamp(get_data(forecast))[1]
get_scaling_factors(forecast::Forecast) = get_data(forecast)

function get_resolution(forecast::Forecast)
    data = get_data(forecast)
    return TimeSeries.timestamp(data)[2] - TimeSeries.timestamp(data)[1]
end

struct ForecastKey
    forecast_type::Type{<:ForecastInternal}
    initial_time::Dates.DateTime
    label::String
end

const ForecastsByType = Dict{ForecastKey, ForecastInternal}

"""
Forecast container for a component.
"""
mutable struct Forecasts
    data::ForecastsByType
    time_series_storage::Union{Nothing, TimeSeriesStorage}
end

function Forecasts()
    return Forecasts(ForecastsByType(), nothing)
end

Base.length(forecasts::Forecasts) = length(forecasts.data)
Base.isempty(forecasts::Forecasts) = isempty(forecasts.data)

function set_time_series_storage!(
    forecasts::Forecasts,
    storage::Union{Nothing, TimeSeriesStorage},
)
    if !isnothing(forecasts.time_series_storage) && !isnothing(storage)
        @show forecasts.time_series_storage
        throw(ArgumentError(
            "The time_series_storage reference is already set. Is this component being " *
            "added to multiple systems?",
        ))
    end

    forecasts.time_series_storage = storage
end

function add_forecast!(
    forecasts::Forecasts,
    forecast::T;
    skip_if_present = false,
) where {T <: ForecastInternal}
    key = ForecastKey(T, get_initial_time(forecast), get_label(forecast))
    if haskey(forecasts.data, key)
        if skip_if_present
            @warn "forecast $key is already present, skipping overwrite"
        else
            throw(ArgumentError("forecast $key is already stored"))
        end
    else
        forecasts.data[key] = forecast
    end
end

function remove_forecast!(
    ::Type{T},
    forecasts::Forecasts,
    initial_time::Dates.DateTime,
    label::AbstractString,
) where {T <: ForecastInternal}
    key = ForecastKey(T, initial_time, label)
    if !haskey(forecasts.data, key)
        throw(ArgumentError("forecast $key is not stored"))
    end

    pop!(forecasts.data, key)
end

function clear_forecasts!(forecasts::Forecasts)
    empty!(forecasts.data)
end

function get_forecast(
    ::Type{T},
    forecasts::Forecasts,
    initial_time::Dates.DateTime,
    label::AbstractString,
) where {T <: ForecastInternal}
    key = ForecastKey(T, initial_time, label)
    if !haskey(forecasts.data, key)
        throw(ArgumentError("forecast $key is not stored"))
    end

    return forecasts.data[key]
end

function get_forecast_initial_times(forecasts::Forecasts)::Vector{Dates.DateTime}
    initial_times = Set{Dates.DateTime}()
    for key in keys(forecasts.data)
        push!(initial_times, key.initial_time)
    end

    return sort!(Vector{Dates.DateTime}(collect(initial_times)))
end

function get_forecast_initial_times(
    ::Type{T},
    forecasts::Forecasts,
) where {T <: ForecastInternal}
    initial_times = Set{Dates.DateTime}()
    for key in keys(forecasts.data)
        if key.forecast_type <: T
            push!(initial_times, key.initial_time)
        end
    end

    return sort!(Vector{Dates.DateTime}(collect(initial_times)))
end

function get_forecast_initial_times(
    ::Type{T},
    forecasts::Forecasts,
    label::AbstractString,
) where {T <: ForecastInternal}
    initial_times = Set{Dates.DateTime}()
    for key in keys(forecasts.data)
        if key.forecast_type <: T && key.label == label
            push!(initial_times, key.initial_time)
        end
    end

    return sort!(Vector{Dates.DateTime}(collect(initial_times)))
end

function get_forecast_initial_times!(
    initial_times::Set{Dates.DateTime},
    forecasts::Forecasts,
)
    for key in keys(forecasts.data)
        push!(initial_times, key.initial_time)
    end
end

function get_forecast_labels(
    ::Type{T},
    forecasts::Forecasts,
    initial_time::Dates.DateTime,
) where {T <: ForecastInternal}
    labels = Set{String}()
    for key in keys(forecasts.data)
        if key.forecast_type <: T && key.initial_time == initial_time
            push!(labels, key.label)
        end
    end

    return Vector{String}(collect(labels))
end

struct ForecastSerializationWrapper
    forecast::ForecastInternal
    type::DataType
end

function JSON2.write(io::IO, forecasts::Forecasts)
    return JSON2.write(io, encode_for_json(forecasts))
end

function JSON2.write(forecasts::Forecasts)
    return JSON2.write(encode_for_json(forecasts))
end

function encode_for_json(forecasts::Forecasts)
    # Store a flat array of forecasts. Deserialization can unwind it.
    data = Vector{ForecastSerializationWrapper}()
    for (key, forecast) in forecasts.data
        push!(data, ForecastSerializationWrapper(forecast, key.forecast_type))
    end

    return data
end

function JSON2.read(io::IO, ::Type{Forecasts})
    forecasts = Forecasts()
    for raw_forecast in JSON2.read(io)
        forecast_type = getfield(
            InfrastructureSystems,
            Symbol(strip_module_name(string(raw_forecast.type))),
        )
        forecast = JSON2.read(JSON2.write(raw_forecast.forecast), forecast_type)
        add_forecast!(forecasts, forecast)
    end

    return forecasts
end

function Base.length(forecast::Forecast)
    return get_horizon(forecast)
end

function Base.getindex(forecast::Forecast, args...)
    return _split_forecast(forecast, getindex(get_data(forecast), args...))
end

Base.first(forecast::Forecast) = head(forecast, 1)

Base.last(forecast::Forecast) = tail(forecast, 1)

Base.firstindex(forecast::Forecast) = firstindex(get_data(forecast))

Base.lastindex(forecast::Forecast) = lastindex(get_data(forecast))

Base.lastindex(forecast::Forecast, d) = lastindex(get_data(forecast), d)

Base.eachindex(forecast::Forecast) = eachindex(get_data(forecast))

Base.iterate(forecast::Forecast, n = 1) = iterate(get_data(forecast), n)

"""
Refer to TimeSeries.when(). Underlying data is copied.
"""
function when(forecast::Forecast, period::Function, t::Integer)
    new = _split_forecast(forecast, TimeSeries.when(get_data(forecast), period, t))

end

"""
Return a forecast truncated starting with timestamp.
"""
function from(forecast::T, timestamp) where {T <: Forecast}
    return T(get_label(forecast), TimeSeries.from(get_data(forecast), timestamp))
end

"""
Return a forecast truncated after timestamp.
"""
function to(forecast::T, timestamp) where {T <: Forecast}
    return T(get_label(forecast), TimeSeries.to(get_data(forecast), timestamp))
end

"""
Return a forecast with only the first num values.
"""
function head(forecast::Forecast)
    return _split_forecast(forecast, TimeSeries.head(get_data(forecast)))
end

function head(forecast::Forecast, num)
    return _split_forecast(forecast, TimeSeries.head(get_data(forecast), num))
end

"""
Return a forecast with only the ending num values.
"""
function tail(forecast::Forecast)
    return _split_forecast(forecast, TimeSeries.tail(get_data(forecast)))
end

function tail(forecast::Forecast, num)
    return _split_forecast(forecast, TimeSeries.tail(get_data(forecast), num))
end

"""
Creates a new forecast from an existing forecast with a split TimeArray.
"""
function

_split_forecast(forecast::T, data::TimeSeries.TimeArray;) where {T <: Forecast}
    vals = []
    for (fname, ftype) in zip(fieldnames(T), fieldtypes(T))
        if ftype <: TimeSeries.TimeArray
            val = data
        elseif ftype <: InfrastructureSystemsInternal
            # Need to create a new UUID.
            continue
        else
            val = getfield(forecast, fname)
        end

        push!(vals, val)
    end

    return T(vals...)
end

function get_resolution(ts::TimeSeries.TimeArray)
    tstamps = TimeSeries.timestamp(ts)
    timediffs = unique([tstamps[ix] - tstamps[ix - 1] for ix in 2:length(tstamps)])

    res = []

    for timediff in timediffs
        if mod(timediff, Dates.Millisecond(Dates.Day(1))) == Dates.Millisecond(0)
            push!(res, Dates.Day(timediff / Dates.Millisecond(Dates.Day(1))))
        elseif mod(timediff, Dates.Millisecond(Dates.Hour(1))) == Dates.Millisecond(0)
            push!(res, Dates.Hour(timediff / Dates.Millisecond(Dates.Hour(1))))
        elseif mod(timediff, Dates.Millisecond(Dates.Minute(1))) == Dates.Millisecond(0)
            push!(res, Dates.Minute(timediff / Dates.Millisecond(Dates.Minute(1))))
        elseif mod(timediff, Dates.Millisecond(Dates.Second(1))) == Dates.Millisecond(0)
            push!(res, Dates.Second(timediff / Dates.Millisecond(Dates.Second(1))))
        else
            throw(DataFormatError("cannot understand the resolution of the timeseries"))
        end
    end

    if length(res) > 1
        throw(DataFormatError("timeseries has non-uniform resolution: this is currently not supported"))
    end

    return res[1]
end

get_columns(::Type{<:ForecastInternal}, ta::TimeSeries.TimeArray) = nothing


text:
using FIRLSFilterDesign
using Test

function generate_bands_vec(N; n_duplicates = 0, fs = 1.)
    bands_vec = rand(N)
    bands_vec = vcat(bands_vec, bands_vec[1:n_duplicates])
    sort!(bands_vec)
    bands_vec .-= bands_vec[1]
    bands_vec ./= bands_vec[end]
    bands_vec .*= fs/2
    return bands_vec
end

vec2mat(vec) = @views hcat(vec[1:end-1], vec[2:end])

function all_perms_q(f,W)
    M = 5
    q_odd = FIRLSFilterDesign.get_q(M, f, W, FIRLSFilterDesign.FIR_I())
    q_even = FIRLSFilterDesign.get_q(M, f, W, FIRLSFilterDesign.FIR_II())
    return q_odd, q_even
end

function all_perms_b(f,W,D)
    M = 5
    b_FIR_I     = FIRLSFilterDesign.get_b(M, f, D, W, FIRLSFilterDesign.FIR_I())
    b_FIR_II    = FIRLSFilterDesign.get_b(M, f, D, W, FIRLSFilterDesign.FIR_II())
    b_FIR_III   = FIRLSFilterDesign.get_b(M, f, D, W, FIRLSFilterDesign.FIR_III())
    b_FIR_IV    = FIRLSFilterDesign.get_b(M, f, D, W, FIRLSFilterDesign.FIR_IV())
    return b_FIR_I, b_FIR_II, b_FIR_III, b_FIR_IV
end

function all_perms_firlsdesign(args...; kwargs...)
    order_odd = 10
    order_even = order_odd + 1
    h_FIR_I     = firls_design(order_odd,  args..., false; kwargs...)
    h_FIR_II    = firls_design(order_even, args..., false; kwargs...)
    h_FIR_III   = firls_design(order_odd,  args..., true;  kwargs...)
    h_FIR_IV    = firls_design(order_even, args..., true;  kwargs...)
    return h_FIR_I, h_FIR_II, h_FIR_III, h_FIR_IV
end

function all_perms_toimpulseresponse(a)
    h_FIR_I     = FIRLSFilterDesign._to_impulse_response(a, FIRLSFilterDesign.FIR_I())
    h_FIR_II    = FIRLSFilterDesign._to_impulse_response(a, FIRLSFilterDesign.FIR_II())
    h_FIR_III   = FIRLSFilterDesign._to_impulse_response(a, FIRLSFilterDesign.FIR_III())
    h_FIR_IV    = FIRLSFilterDesign._to_impulse_response(a, FIRLSFilterDesign.FIR_IV())
    return h_FIR_I, h_FIR_II, h_FIR_III, h_FIR_IV
end

# Randomized inputs
fs_rand = 1 + rand() * 10
bands_D_mat_rand = vec2mat(generate_bands_vec(100, n_duplicates = 10, fs = fs_rand))
D_mat_rand = rand(size(bands_D_mat_rand)...)
bands_W_mat_rand = vec2mat(generate_bands_vec(100, n_duplicates = 10, fs = fs_rand))
W_mat_rand = rand(size(bands_W_mat_rand)...)
# Fixed inputs
fs = 1.
bands_D_mat = [0. .25; 0.25 0.25; .25 .5]
D_mat = [2. 1.; 10. 10.; 0. 0.]
bands_W_mat = [0. .25; 0.25 0.25; .25 .5]
W_mat = [1. 2.; 11. 11.; 2. 1.]

hs_rand1 = all_perms_firlsdesign(bands_D_mat_rand, D_mat_rand; fs = fs_rand)
hs_rand2 = all_perms_firlsdesign(bands_D_mat_rand, D_mat_rand, fill(rand(), size(D_mat_rand)); fs = fs_rand)
# hs_rand3 = all_perms_firlsdesign(bands_D_mat_rand, D_mat_rand, bands_D_mat_rand, ones(size(D_mat_rand)); fs = fs_rand)

qs_fixed = all_perms_q(bands_D_mat, W_mat)
qs_fixed_data = [
    [
        1.5
        0.0
        -0.2026423672846756
        0.0
        0.0
        0.0
        -0.02251581858718621
        0.0
        0.0
        0.0
        -0.008105694691387062
   ],
   [
        1.5
        0.0
        -0.2026423672846756
        0.0
        0.0
        0.0
        -0.02251581858718621
        0.0
        0.0
        0.0
        -0.008105694691387062
        0.0
   ]
]

bs_fixed = all_perms_b(bands_D_mat, D_mat, W_mat)
bs_fixed_data = [
    [
        1.0833333333333333
        0.6919896805485019
        0.0
        -0.24427841957880944
        -0.025330295910584485
        0.12128235798585402
    ],
    [
        0.9761255156564024
        0.32778347331995833
       -0.2012359779620491
       -0.16111118889304135
        0.08495418615790015
        0.080981441609415
    ],
    [
        NaN
        0.6919896805485019
        0.7011228412339804
        0.2442784195788095
        0.0
        0.12128235798585399
    ],
    [
        0.40432442716331324
        0.7913393068108028
        0.48582662723339226
        0.06673443948955127
        0.03518917608697089
        0.19550649463397424
    ]
]

hs_fixed = all_perms_firlsdesign(bands_D_mat, D_mat, W_mat)
hs_fixed_data = [
    [
        0.07788241325853798
        -0.0017268410299129148
        -0.08309128924173821
        0.1009958144655533
        0.5217570192015153
        0.7495102634610803
        0.5217570192015153
        0.1009958144655533
        -0.08309128924173821
        -0.0017268410299129148
        0.07788241325853798
    ],
    [
        0.056631088208928274
        0.05686190577259725
        -0.05915886754220385
        -0.034409300873168014
        0.30433805233189765
        0.6880664081058022
        0.6880664081058022
        0.30433805233189765
        -0.034409300873168014
        -0.05915886754220385
        0.05686190577259725
        0.056631088208928274
    ],
    [
        0.10513740773549088
        0.056950501912317146
        0.2320211851531876
        0.4742540923933572
        0.432645024065513
        0.0
       -0.432645024065513
       -0.4742540923933572
       -0.2320211851531876
       -0.056950501912317146
       -0.10513740773549088
    ],
    [
        0.14347815111515985
        0.06421192516629874
        0.12735250873084986
        0.3641022891441058
        0.5103455604394119
        0.24763921096456235
        -0.24763921096456235
        -0.5103455604394119
        -0.3641022891441058
        -0.12735250873084986
        -0.06421192516629874
        -0.14347815111515985
    ]
]

function get_half(h)
    N = length(h)
    M = N÷2
    if iseven(N)
        @views h[1:M], h[end:-1:M+1]
    else
        @views h[1:M+1], h[end:-1:M+1]
    end
end

results = all_perms_toimpulseresponse(rand(5))
halfs = get_half.(results)

all_approx(x,y) = all(@. abs(x - y) < 1e16)

@testset "FIRLSFilterDesign.jl" begin
    # q vector test
    @test all_approx(qs_fixed[1], qs_fixed[2][1:end-1])
    @test all_approx(qs_fixed[1], qs_fixed_data[1]) # odd filter length
    @test all_approx(qs_fixed[2], qs_fixed_data[2]) # even filter length
    # b vector test
    @test all_approx(bs_fixed[1],        bs_fixed_data[1])        # odd  filter length, symmetric
    @test all_approx(bs_fixed[2],        bs_fixed_data[2])        # even filter length, symmetric
    @test all_approx(bs_fixed[3][2:end], bs_fixed_data[3][2:end]) # odd  filter length, antisymmetric
    @test all_approx(bs_fixed[4],        bs_fixed_data[4])        # even filter length, antisymmetric
    # Symmetry of filter coefficients
    @test all(halfs[1][1] .==  halfs[1][2]) # odd  filter length, symmetric
    @test all(halfs[2][1] .==  halfs[2][2]) # even filter length, symmetric
    @test all(halfs[3][1] .== -halfs[3][2]) # odd  filter length, antisymmetric
    @test all(halfs[4][1] .== -halfs[4][2]) # even filter length, antisymmetric
    # Tests that show equivalence between passing no weights and passing all equal weights
    @test all_approx(hs_rand1[1], hs_rand2[1]) # odd  filter length, symmetric
    @test all_approx(hs_rand1[2], hs_rand2[2]) # even filter length, symmetric
    @test all_approx(hs_rand1[3], hs_rand2[3]) # odd  filter length, antisymmetric
    @test all_approx(hs_rand1[4], hs_rand2[4]) # even filter length, antisymmetric
    # Fixed inputs
    @test all_approx(hs_fixed[1], hs_fixed_data[1]) # odd  filter length, symmetric
    @test all_approx(hs_fixed[2], hs_fixed_data[2]) # even filter length, symmetric
    @test all_approx(hs_fixed[3], hs_fixed_data[3]) # odd  filter length, antisymmetric
    @test all_approx(hs_fixed[4], hs_fixed_data[4]) # even filter length, antisymmetric
end


text:
export LoggerPipelineStage

"""
    LoggerPipelineStage <: AbstractPipelineStage

Logs all updates from `stream` into `output`

# Arguments
- `output`: (optional), an output stream used to print log statements

See also: [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref), [`EmptyPipelineStage`](@ref), [`CompositePipelineStage`](@ref)
"""
struct LoggerPipelineStage{T} <: AbstractPipelineStage
    output :: T
    prefix :: String
end

LoggerPipelineStage()               = LoggerPipelineStage(Core.stdout, "Log")
LoggerPipelineStage(output::IO)     = LoggerPipelineStage(output, "Log")
LoggerPipelineStage(prefix::String) = LoggerPipelineStage(Core.stdout, prefix)

apply_pipeline_stage(stage::LoggerPipelineStage, factornode, tag::Type{ <: Val{ T } },    stream) where { T } = stream |> tap((v) -> Core.println(stage.output, "[$(stage.prefix)][$(functionalform(factornode))][$(T)]: $v"))
apply_pipeline_stage(stage::LoggerPipelineStage, factornode, tag::Tuple{ Val{ T }, Int }, stream) where { T } = stream |> tap((v) -> Core.println(stage.output, "[$(stage.prefix)][$(functionalform(factornode))][$(T):$(tag[2])]: $v"))


text:
recurs_sym_type(ex::Any) =
  (ex==None || typeof(ex)==Symbol || length(ex.args)==1) ? eval(ex) : Expr(ex.head, ex.args[1], recurs_sym_type(ex.args[2]))
macro c(ret_type, func, arg_types, lib)
  local _arg_types = Expr(:tuple, [recurs_sym_type(a) for a in arg_types.args]...)
  local _ret_type = recurs_sym_type(ret_type)
  local _args_in = Any[ symbol(string('a',x)) for x in 1:length(_arg_types.args) ]
  local _lib = eval(lib)
  quote
    $(esc(func))($(_args_in...)) = ccall( ($(string(func)), $(Expr(:quote, _lib)) ), $_ret_type, $_arg_types, $(_args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  real_t = recurs_sym_type(real_t)
  quote
    typealias $fake_t $real_t
  end
end

# header: /usr/local/include/sundials/sundials_band.h
@ctypedef DlsMat Ptr{:Void}
@c DlsMat NewDenseMat (:Int64,:Int64) shlib
@c DlsMat NewBandMat (:Int64,:Int64,:Int64,:Int64) shlib
@c None DestroyMat (:DlsMat,) shlib
@c Ptr{:Int32} NewIntArray (:Int32,) shlib
@c Ptr{:Int64} NewLintArray (:Int64,) shlib
@c Ptr{:realtype} NewRealArray (:Int64,) shlib
@c None DestroyArray (Ptr{:None},) shlib
@c None AddIdentity (:DlsMat,) shlib
@c None SetToZero (:DlsMat,) shlib
@c None PrintMat (:DlsMat,) shlib
@c Ptr{Ptr{:realtype}} newDenseMat (:Int64,:Int64) shlib
@c Ptr{Ptr{:realtype}} newBandMat (:Int64,:Int64,:Int64) shlib
@c None destroyMat (Ptr{Ptr{:realtype}},) shlib
@c Ptr{:Int32} newIntArray (:Int32,) shlib
@c Ptr{:Int64} newLintArray (:Int64,) shlib
@c Ptr{:realtype} newRealArray (:Int64,) shlib
@c None destroyArray (Ptr{:None},) shlib
@c Int64 BandGBTRF (:DlsMat,Ptr{:Int64}) shlib
@c Int64 bandGBTRF (Ptr{Ptr{:realtype}},:Int64,:Int64,:Int64,:Int64,Ptr{:Int64}) shlib
@c None BandGBTRS (:DlsMat,Ptr{:Int64},Ptr{:realtype}) shlib
@c None bandGBTRS (Ptr{Ptr{:realtype}},:Int64,:Int64,:Int64,Ptr{:Int64},Ptr{:realtype}) shlib
@c None BandCopy (:DlsMat,:DlsMat,:Int64,:Int64) shlib
@c None bandCopy (Ptr{Ptr{:realtype}},Ptr{Ptr{:realtype}},:Int64,:Int64,:Int64,:Int64,:Int64) shlib
@c None BandScale (:realtype,:DlsMat) shlib
@c None bandScale (:realtype,Ptr{Ptr{:realtype}},:Int64,:Int64,:Int64,:Int64) shlib
@c None bandAddIdentity (Ptr{Ptr{:realtype}},:Int64,:Int64) shlib

# header: /usr/local/include/sundials/sundials_config.h

# header: /usr/local/include/sundials/sundials_dense.h
@c Int64 DenseGETRF (:DlsMat,Ptr{:Int64}) shlib
@c None DenseGETRS (:DlsMat,Ptr{:Int64},Ptr{:realtype}) shlib
@c Int64 denseGETRF (Ptr{Ptr{:realtype}},:Int64,:Int64,Ptr{:Int64}) shlib
@c None denseGETRS (Ptr{Ptr{:realtype}},:Int64,Ptr{:Int64},Ptr{:realtype}) shlib
@c Int64 DensePOTRF (:DlsMat,) shlib
@c None DensePOTRS (:DlsMat,Ptr{:realtype}) shlib
@c Int64 densePOTRF (Ptr{Ptr{:realtype}},:Int64) shlib
@c None densePOTRS (Ptr{Ptr{:realtype}},:Int64,Ptr{:realtype}) shlib
@c Int32 DenseGEQRF (:DlsMat,Ptr{:realtype},Ptr{:realtype}) shlib
@c Int32 DenseORMQR (:DlsMat,Ptr{:realtype},Ptr{:realtype},Ptr{:realtype},Ptr{:realtype}) shlib
@c Int32 denseGEQRF (Ptr{Ptr{:realtype}},:Int64,:Int64,Ptr{:realtype},Ptr{:realtype}) shlib
@c Int32 denseORMQR (Ptr{Ptr{:realtype}},:Int64,:Int64,Ptr{:realtype},Ptr{:realtype},Ptr{:realtype},Ptr{:realtype}) shlib
@c None DenseCopy (:DlsMat,:DlsMat) shlib
@c None denseCopy (Ptr{Ptr{:realtype}},Ptr{Ptr{:realtype}},:Int64,:Int64) shlib
@c None DenseScale (:realtype,:DlsMat) shlib
@c None denseScale (:realtype,Ptr{Ptr{:realtype}},:Int64,:Int64) shlib
@c None denseAddIdentity (Ptr{Ptr{:realtype}},:Int64) shlib

# header: /usr/local/include/sundials/sundials_direct.h

# header: /usr/local/include/sundials/sundials_fnvector.h

# header: /usr/local/include/sundials/sundials_iterative.h
# enum ANONYMOUS
@ctypedef ANONYMOUS Uint32
const PREC_NONE = 0
const PREC_LEFT = 1
const PREC_RIGHT = 2
const PREC_BOTH = 3
# end
# enum ATimesFn
@ctypedef ATimesFn Uint32
const MODIFIED_GS = 1
const CLASSICAL_GS = 2
# end
@ctypedef PSolveFn Ptr{:Void}
@c Int32 ModifiedGS (Ptr{:N_Vector},Ptr{Ptr{:realtype}},:Int32,:Int32,Ptr{:realtype}) shlib
@c Int32 ClassicalGS (Ptr{:N_Vector},Ptr{Ptr{:realtype}},:Int32,:Int32,Ptr{:realtype},:N_Vector,Ptr{:realtype}) shlib
@c Int32 QRfact (:Int32,Ptr{Ptr{:realtype}},Ptr{:realtype},:Int32) shlib
@c Int32 QRsol (:Int32,Ptr{Ptr{:realtype}},Ptr{:realtype},Ptr{:realtype}) shlib

# header: /usr/local/include/sundials/sundials_math.h
@c realtype RPowerI (:realtype,:Int32) shlib
@c realtype RPowerR (:realtype,:realtype) shlib
@c realtype RSqrt (:realtype,) shlib
@c realtype RAbs (:realtype,) shlib
@c realtype RExp (:realtype,) shlib

# header: /usr/local/include/sundials/sundials_nvector.h

# header: /usr/local/include/sundials/sundials_spbcgs.h
@ctypedef SpbcgMemRec Void
@ctypedef SpbcgMem Ptr{:Void}
@c SpbcgMem SpbcgMalloc (:Int32,:N_Vector) shlib
@c Int32 SpbcgSolve (:SpbcgMem,Ptr{:None},:N_Vector,:N_Vector,:Int32,:realtype,Ptr{:None},:N_Vector,:N_Vector,:ATimesFn,:PSolveFn,Ptr{:realtype},Ptr{:Int32},Ptr{:Int32}) shlib
@c None SpbcgFree (:SpbcgMem,) shlib

# header: /usr/local/include/sundials/sundials_spgmr.h
@ctypedef SpgmrMemRec Void
@ctypedef SpgmrMem Ptr{:Void}
@c SpgmrMem SpgmrMalloc (:Int32,:N_Vector) shlib
@c Int32 SpgmrSolve (:SpgmrMem,Ptr{:None},:N_Vector,:N_Vector,:Int32,:Int32,:realtype,:Int32,Ptr{:None},:N_Vector,:N_Vector,:ATimesFn,:PSolveFn,Ptr{:realtype},Ptr{:Int32},Ptr{:Int32}) shlib
@c None SpgmrFree (:SpgmrMem,) shlib

# header: /usr/local/include/sundials/sundials_sptfqmr.h
@ctypedef SptfqmrMemRec Void
@ctypedef SptfqmrMem Ptr{:Void}
@c SptfqmrMem SptfqmrMalloc (:Int32,:N_Vector) shlib
@c Int32 SptfqmrSolve (:SptfqmrMem,Ptr{:None},:N_Vector,:N_Vector,:Int32,:realtype,Ptr{:None},:N_Vector,:N_Vector,:ATimesFn,:PSolveFn,Ptr{:realtype},Ptr{:Int32},Ptr{:Int32}) shlib
@c None SptfqmrFree (:SptfqmrMem,) shlib

# header: /usr/local/include/sundials/sundials_types.h


text:
# This file was generated, do not modify it. # hide
cmach = machine(creg, X, y)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials