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 = "", bdd_and_unbdd = false,
    bdd_and_bdd::Bool = false,
    plotroot = figurespath(m, \"forecast\"), titles = [],
    plot_handles = fill(plot(), length(vars)), verbose = :low,
    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`
- `bdd_and_unbdd::Bool`: if true, then unbounded means and bounded bands are plotted
- `bdd_and_bdd::Bool`: if true, then bounded means and bounded bands are plotted
- `untrans::Bool`: whether to plot untransformed (model units) history and forecast
- `fourquarter::Bool`: whether to plot four-quarter history and forecast
- `plotroot::String`: 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`

See `?histforecast` 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 = "",
                                   bdd_and_unbdd::Bool = false, bdd_and_bdd::Bool = false,
                                   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,
                                   kwargs...)
    # Determine output_vars
    if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    elseif bdd_and_bdd
        if bdd_and_unbdd
            error("Only one of bdd_and_unbdd and bdd_and_bdd can be true")
        end
        hist_prod  = :hist
        fcast_prod = :bddforecast
    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,
                    bdd_and_unbdd = bdd_and_unbdd, 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, 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,
                                                fileformat = plot_extension())
            @show output_file
            DSGE.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

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

@recipe function f(hf::HistForecast;
                   start_date = hf.args[2].means[1, :date],
                   end_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,
                   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]

    # 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"

        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

    # 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
end


text:
using StringViews, Test

b = Vector{UInt8}("foobar")
s = StringView(b)
ss = SubString(s, 2, 5) # "ooba"
abc = StringView(0x61:0x63)
invalid = StringView([0x8b, 0x52, 0x9b, 0x8d])
su = StringView("föôẞαr")

@testset "construction/conversion" begin
    @test StringView(s) === s
    @test Vector{UInt8}(s) === Array{UInt8}(s) === codeunits(s) === b
    @test Vector{UInt8}(StringView(@view b[1:3])) == b[1:3]
    @test codeunits(String(s)) == s.data
    @test Vector{UInt8}(abc) == collect(0x61:0x63)
    @test Symbol(s) == :foobar
    @test Symbol(abc) == :abc

    @test StringView("foo") isa StringView{Base.CodeUnits{UInt8,String}}

    @test s isa StringViews.DenseStringView
    @test StringView(@view b[1:3]) isa StringViews.DenseStringView
    @test StringView("foo") isa StringViews.DenseStringView
    @test StringView(@view codeunits("foobar")[1:3]) isa StringViews.DenseStringView

    @test pointer(s) == pointer(b) == Base.unsafe_convert(Ptr{UInt8}, s)
    @test Base.unsafe_convert(Ptr{Int8}, s) == Ptr{Int8}(pointer(s))
    @test pointer(s, 3) == pointer(b, 3)
    @test_throws MethodError pointer(abc)

    @test ncodeunits(s) == sizeof(s) == length(b)
    @test codeunit(s) == UInt8
    @test codeunit(s,3) == b[3]

    @test Base.print_to_string(s) == "foobar"
    @test Base.print_to_string(abc) == "abc"
end

@testset "substrings" begin
    @test Vector{UInt8}(ss) == Array{UInt8}(ss) == codeunits(ss) == b[2:5]
    @test codeunits(ss) isa Base.FastContiguousSubArray
    @test Symbol(ss) == :ooba

    @test pointer(ss) == pointer(b) + 1 == Base.unsafe_convert(Ptr{UInt8}, ss)
    @test ncodeunits(ss) == sizeof(ss) == length(b)-2
    @test codeunit(ss) == UInt8
    @test codeunit(ss,3) == b[4]

    @test Base.print_to_string(ss) == "ooba"

    @test cmp("foobar","bar") == cmp(ss,"bar") == -cmp("bar",ss) == cmp(ss,StringView("bar"))
    @test ss == StringView("ooba") == "ooba" == ss == "ooba"
    @test isvalid(ss)
end

@testset "iteration" begin
    for str in (s, ss, abc, invalid, su)
        sS = String(str)
        @test sS == str
        @test length(sS) == length(str)
        @test collect(sS) == collect(str) ==
              getindex.(sS, eachindex(sS)) == getindex.(str, eachindex(sS))
        @test collect(eachindex(sS)) == collect(eachindex(str))
        @test sS[1:end] == str[1:end]
        @test sS[nextind(sS,1):prevind(sS,end)] == str[nextind(str,1):prevind(str,end)]
    end
end

@testset "regular expressions" begin
    for str in (s,ss)
        @test [m.match for m in collect(eachmatch(r"[aeiou]+", str))] == ["oo", "a"]
        @test occursin(r"o+", str) && !occursin(r"z+", str)
        @test startswith(str, r"o+") == (str[1:2] == "oo")
        @test startswith(str, r"f+") == (str[1:2] == "fo")
        @test endswith(str, r"[aeiou]") == (str[end] == 'a')
        @test endswith(str, r"[q-z]") == (str[end] == 'r')
        @test findnext(r"o+", str, 4) === nothing
    end
    @test findnext(r"[aeiou]+", s, 1) == 2:3
    @test findnext(r"[aeiou]+", ss, 1) == 1:2
end

@testset "parsing" begin
    for val in (true, 1234, 1234.5, 1234.5f0, 4.5+3.25im)
        sval = string(val)
        for str in (StringView(sval), SubString("foo"*sval*"bar", 4, 3+length(sval)))
            @test parse(typeof(val), str) === val
        end
    end
end

@testset "searching" begin
    for str in (s, ss, abc)
        sS, n = String(str), lastindex(str)
        @test startswith(str, "foo") == startswith(sS, "foo")
        @test endswith(str, "bar") == endswith(sS, "bar")
        @test replace(str, 'o'=>"xy") == replace(sS, 'o'=>"xy")
        @test replace(str, ('o','a')=>'x') == replace(sS, ('o','a')=>'x')
        @test findnext(==('b'), str, 1) === findnext(==('b'), sS, 1)
        @test findprev(==('b'), str, n) === findprev(==('b'), sS, n)
        @test findnext(==("ba"), str, 1) === findnext(==("ba"), sS, 1)
        @test findprev(==("ba"), str, n) === findprev(==("ba"), sS, n)
    end
    @test chomp(StringView("foo\n")) == "foo"

    # issue #5
    let v = [0x32, 0x30, 0x32, 0x31, 0x2d, 0x31, 0x31, 0x2d, 0x31, 0x30, 0x20, 0x32, 0x31, 0x3a, 0x34, 0x32, 0x3a, 0x30, 0x35, 0x2e, 0x31, 0x31, 0x35, 0x38, 0x30, 0x37],
        pat = r"(\.[\d]{3})\d+" => s"\g<1>"
        @test replace(String(copy(v)), pat) == replace(StringView(v), pat)
    end
end

@testset "miscellaneous" begin
    @test cmp("foobar","bar") == cmp(s,"bar") == -cmp("bar",s) == cmp(s,StringView("bar"))
    @test s == StringView("foobar") == "foobar" == s == "foobar" != StringView("bar")
    @test cmp(abc, "bar") == cmp("abc","bar")

    @test Base.typemin(s) isa StringView{Vector{UInt8}}
    @test Base.typemin(s) == ""

    @test isascii(s)
    @test !isascii(StringView("fööbār"))

    @test isvalid(s)
    @test isvalid(abc)
    @test !isvalid(invalid)
    @test !invoke(isvalid, Tuple{StringView}, invalid)

    for str in (s, abc, invalid, ss, su), f in (hash, reverse)
        @test f(str) == f(String(str))
    end
end

ext: jl


text:
#=
Julia script to get keywords in a target corpus
Earl K. Brown, ekbrown@byu.edu
based on some logic in Python code written by Adam Davies
=#

using Formatting, DataFrames

function get_num_wds_freqs(pathway)
    # Helper function to retrieve the number of words and the frequencies of those words in a directory
    cd(pathway)
    files = [file for file in readdir() if occursin(r"\.txt"i, file)]
    freqs = Dict{String, Int64}()
    num_wds = 0
    for file in files
        open(file) do infile
            whole_file = read(infile, String)
            wds = split(whole_file, r"[^-'a-záéíóúüñ]+"i, keepempty = false)
            wds = [uppercase(wd) for wd in wds]
            num_wds += length(wds)
            for wd in wds
                freqs[wd] = get(freqs, wd, 0) + 1
            end
        end
    end
    println("There are " * format(num_wds, commas = true) * " words in $pathway")
    return (num_wds, freqs)
end

function get_keywords(target_dir, ref_dir, num_keywords, min_freq)
    #=
    Get keywords in .txt files within a target directory, comparing them with words in .txt files within a reference directory.
    param: target_dir - the directory with the target corpus in .txt files
    param: ref_dir - the directory with the reference corpus in .txt files
    param: num_keywords - number of keywords desired
    param: min_freq - the minimum frequency of keywords in the target corpus
    return value: a DataFrame with three columns: (1) keyword, (2) frequency in target corpus, (3) keyness score
    =#

    # get number of words and freqs in target and reference directories
    target_num_wds, target_freqs = get_num_wds_freqs(target_dir)
    ref_num_wds, ref_freqs = get_num_wds_freqs(ref_dir)

    # calculate frequency ratio between target corpus and reference corpus
    rel_freq = target_num_wds / ref_num_wds

    # calculate keyness
    keywords = Dict{String, Float64}()
    for wd in [i[1] for i in sort(collect(target_freqs), by=x->x[2], rev = true)]
        if target_freqs[wd] >= min_freq
            if haskey(ref_freqs, wd)
                keywords[wd] = log2((target_freqs[wd] * rel_freq) / ref_freqs[wd])
            end
        end
    end

    # sort keywords and limit to number desired by user
    top_keywords = sort(collect(keywords), by = x -> x[2], rev = true)[1:num_keywords]

    # push keywords to DataFrame
    df = DataFrame(keyword = String[], freq = Int64[], keyness = Float64[])
    # df = DataFrame(A = Int64[], B = Int64[])
    for kw in top_keywords
        push!(df, [kw[1], target_freqs[kw[1]], round(kw[2], digits = 2)])
    end

    return df
end

### test the function
target_dir = "/Users/ekb5/Documents/LING_580R/gen_conf/"
ref_dir = "/Users/ekb5/Corpora/Brown/"
num_keywords = 10
min_freq = 3

@time println(get_keywords(target_dir, ref_dir, num_keywords, min_freq))


text:
TecatorStatement = "Statement of permission from Tecator (the original data source).These data are recorded on a Tecator" *
"\n Infratec Food and Feed Analyzer working in the wavelength range 850 - 1050 nm by the Near Infrared" *
"\n Transmission (NIT) principle. Each sample contains finely chopped pure meat with different moisture, fat" *
"\n and protein contents.If results from these data are used in a publication we want you to mention the" *
"\n instrument and company name (Tecator) in the publication. In addition, please send a preprint of your " *
"\n article to Karin Thente, Tecator AB, Box 70, S-263 21 Hoganas, Sweden. The data are available in the " *
"\n public domain with no responsability from the original data source. The data can be redistributed as long " *
"\n as this permission note is attached. For more information about the instrument - call Perstorp Analytical's" *
"\n representative in your area."

datapath = joinpath(@__DIR__, "..", "data")

"""
    ChemometricsToolsDatasets()

Displays a list of the available datasets in ChemometricsTools.jl .

"""
function ChemometricsToolsDatasets()
    dircontents = readdir(datapath)
    dircontents = [ f for f in dircontents if f != "Readme.md" ]
    return Dict( (1:length(dircontents)) .=> dircontents )
end

"""
    ChemometricsToolsDataset(file::Int)

Loads in a dataset included in the ChemometricsTools.jl package.

"""
function ChemometricsToolsDataset(filename::String)
    if filename == "tecator.csv"
        println(TecatorStatement)
    end
    if filename != "Readme.md"
        read( Base.joinpath( datapath, filename ) )
    else
        println("Don't load the markdown Readme as a csv... You're better than this.")
    end
end

"""
    ChemometricsToolsDataset(file::Int)

Loads in a dataset included in the ChemometricsTools.jl package.

"""
function ChemometricsToolsDataset(file::Int)
    if readdir(datapath)[file] == "tecator.csv"
        println(TecatorStatement)
    end
    read( Base.joinpath( datapath, readdir(datapath)[file] ) )
end

"""
    FindCommonVariables( DFArray::Array{DataFrame,1} )

Finds exact matches across variable names over an array of DataFrames.
Variables that are not mutual to each dataframe are discarded.

For intelligent/soft matching please see ExploreCommonVariables().

"""
function FindCommonVariables( DFArray::Array{DataFrame,1} )
    common = []
    for df in DFArray
        curvars = Set( names( df ) )
        common = ( length( common ) > 0) ? curvars : intersect(common, curvars)
    end
    return common
end

# function ExploreCommonVariables( DFArray::Array{DataFrame,1}, allowed_distance = 1 )
#     common_names = []
#     common_types = []
#     for df in DFArray
#         cur_names = names( df )
#         cur_types = types( df )
#         if length( common ) > 0
#             common_names = cur_names
#         else
#             common_names = intersect(common, cur_names)
#         end
#     end
#     return common_names, common_types
# end


text:
# to do: alt, ...

wm(w::WordType) = Maton(length(w) + 1,
                        Set([1]), Set([length(w) + 1]),
                        [(k,k + 1) for k in 1:length(w)],
                        [c for c in w])

botm() = Maton(0, Set{Int}(), Set{Int}(),
               Pair{Int,Int}[],
               eltype(ABET)[])

topm() = Maton(1, Set([1]), Set([1]),
               [(1,1) for _ in ABET],
               [c for c in ABET])

# alt requires a ren (for rename)


text:
using LinearAlgebra, GenericLinearAlgebra, Arblib

if !isdefined(@__MODULE__, :ColumnNorm)
    # Make ColumnNorm work on old julia versions, it is introduced in Julia 1.7
    ColumnNorm() = Val(true)
end

function approximate_fekete(initial_points, basis; show_det = false, s = 3,alg=[:BigFloat, :Arb, :high_prec][2])
    # P contains the basis transformation from the monomial basis to the new basis.
    V = [pol(point...) for point in initial_points, pol in basis]
    if alg == :BigFloat
        V, P, point_indices = approximate_fekete_bigfloat(V,s=s,show_det=show_det)
    elseif alg == :Arb
        V, P, point_indices = approximate_fekete_arb(V,s=s,show_det=show_det)
    elseif alg == :high_prec
        V, P, point_indices = approximate_fekete_highprec(V,s=s,show_det=show_det)
    end
    return V, P, initial_points[point_indices]
end

function approximate_fekete_bigfloat(V; s = 3, show_det = false)
    # we do the QR factorizations in floating point precision, and the basis changes in high precision (BigFloats)
    P = Matrix{BigFloat}(I, size(V, 2), size(V, 2)) # to keep track of the basis change
    for k = 1:s
        F = qr(Float64.(V))
        U = BigFloat.(I / F.R)
        V = V * U # should be approximately equal to Q
        P = P * U #keep track of the basis change
    end
    # Find the good points:
    F2 = qr(Float64.(V'), ColumnNorm()) # this is deprecated in 1.7, but the new thing (ColumnNorm()) doesn't work in 1.6
    point_indices = F2.p[1:size(V,2)]
    # Do a last basis change to get a good basis for these points
    F = qr(Float64.(V[point_indices,:]))
    U = BigFloat.(I / F.R)
    V = V[point_indices,:] * U # the new V is square
    P = P * U
    if show_det
        println(det(V))
    end
    return V, P, point_indices
end

function approximate_fekete_arb(V::Matrix{T}; s = 3,show_det = false) where T
    # We do the QR factorizations in floating point precision, and the basis changes in high precision (Arb)
    prec = precision(eltype(V))
    V = ArbMatrix(V,prec=prec)
    P = ArbMatrix(size(V, 2), size(V, 2),prec=prec) # to keep track of the basis change
    Arblib.one!(P)
    for k = 1:s
        F = qr(Float64.(V))
        U = ArbMatrix(I / F.R,prec=prec)
        Arblib.approx_mul!(V,V,U)
        Arblib.approx_mul!(P,P,U)
    end
    # Find the good points:
    F2 = qr(Float64.(V'), ColumnNorm())
    point_indices = F2.p[1:size(V,2)]
    # Do a last basis change to get a good basis for these points
    V = V[point_indices,:]
    F = qr(Float64.(V))
    U = ArbMatrix(I / F.R,prec=prec)
    Arblib.approx_mul!(V,V,U)
    Arblib.approx_mul!(P,P,U)
    if show_det
        det_V = Arb(prec=prec)
        Arblib.det!(det_V,V)
        println(det_V)
    end
    return T.(V), T.(P), point_indices
end

function approximate_fekete_highprec(V::Matrix{T}; s = 3,show_det = false) where T
    # We do one basis change, and for that we compute the qr factorization in high precision with BigFloats
    # Then we do the basis change with Arb for speed. Note that especially the QR factorization can take long for large matrices (say 200+ rows/cols)
    # Find the basis change matrix (from the QR factorization)
    prec=precision(V[1,1])
    F = qr(BigFloat.(V))
    U = ArbMatrix(F.R,prec=prec)
    Arblib.approx_inv!(U,U)

    #Do the basis change for V
    V = ArbMatrix(V,prec=prec)
    Arblib.approx_mul!(V,V,U)

    # Do a QR factorization of V^T to get a subset of good sample points
    F2 = qr(BigFloat.(V'), ColumnNorm())
    point_indices = F2.p[1:size(V,2)]
    # Do a last basis change to get a good basis for these points
    V = V[point_indices,:]

    F = qr(BigFloat.(V))
    U2 = ArbMatrix(F.R,prec=prec)
    Arblib.approx_inv!(U2,U2)
    Arblib.approx_mul!(V,V,U2)
    Arblib.approx_mul!(U,U,U2)

    if show_det
        det_V = Arb(prec=prec)
        Arblib.det!(det_V,V)
        println(det_V)
    end
    return T.(V), T.(U), point_indices
end

# end #of module


text:
#using FixedSizeArrays
using CompScienceMeshes
using Test

for T in [Float32, Float64]
    local p = point(T, 0.0, 1.0)
    q = point(T, 2.0, 1.0)
    r = point(T, 2.0, 2.0)
    local a = point(T, 1.0, 0.0)
    b = point(T, 2.0, 0.0)
    c = point(T, 1.0, 2.0)

    @test CompScienceMeshes.leftof(p, a,b) == true
    @test CompScienceMeshes.intersectlines(p,r, a,c) == T.([1.0, 1.5])

    α = CompScienceMeshes.sutherlandhodgman2d([p,q,r], [a,b,c])
    @test length(α) == 4
    @test α[1] == T.([1.2, 1.6])
    @test α[2] == T.([1.0, 1.5])
    @test α[3] == T.([1.0, 1.0])
    @test α[4] == T.([1.5, 1.0])

    β = CompScienceMeshes.sutherlandhodgman([p,q,r], [a,b,c])
    @test length(β) == length(α)
    for i in eachindex(β) @test β[i] == α[i] end

    local P = simplex(p, q, r)
    Q = simplex(a, b, c)
    R = intersection(P,Q)

    @test length(R) == 2
    @test R[1].vertices == [α[1], α[2], α[3]]
    @test R[2].vertices == [α[1], α[3], α[4]]

    t = point(T, 1.0, 2.0, 3.0)
    p = point(T, 0.0, 1.0, 0.0) + t
    q = point(T, 2.0, 1.0, 0.0) + t
    r = point(T, 2.0, 2.0, 0.0) + t
    a = point(T, 1.0, 0.0, 0.0) + t
    b = point(T, 2.0, 0.0, 0.0) + t
    c = point(T, 1.0, 2.0, 0.0) + t

    P = simplex(p,q,r)
    Q = simplex(a,b,c)
    S = intersection(P,Q)

    # Test whether translation and intersection commute
    @test length(S) == 2
    for i in 1:2
        for j in 1:3
            for k in 1:2
                @test S[i].vertices[j][k] ≈ R[i].vertices[j][k]+t[k]
            end
            @test S[i].vertices[j][3] ≈ t[3]
        end
    end
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
struct Wildcard
    type ::Symbol
end

JSON2.write(io::IO, wild::Wildcard; kwargs...) = Base.write(io, "[\"", wild.type, "\"]")

const MATCH = Wildcard(:MATCH)
const ALL = Wildcard(:ALL)
const ALLSMALLER = Wildcard(:ALLSMALLER)

is_wild(a::T) where {T<:Wildcard} = true
is_wild(a) = false

struct TraitInput end
struct TraitOutput end
struct TraitState end

struct Dependency{Trait, IdT <: Union{String, NamedTuple}}
    id ::IdT
    property ::String
    Dependency{Trait}(id::T, property::String) where {Trait, T} = new{Trait, T}(id, property)
end

dep_id_string(dep::Dependency{Trait, String}) where {Trait} = dep.id
dep_id_string(dep::Dependency{Trait, <:NamedTuple}) where {Trait} = sorted_json(dep.id)

dependency_tuple(dep::Dependency) = (id = dep_id_string(dep), property = dep.property)

const Input = Dependency{TraitInput}
const State = Dependency{TraitState}
const Output = Dependency{TraitOutput}

"""
    Base.==(a::Dependency, b::Dependency)

We use "==" to denote two deps that refer to the same prop on
the same component. In the case of wildcard deps, this means
the same prop on *at least one* of the same components.
"""
function Base.:(==)(a::Dependency, b::Dependency)
    (a.property == b.property) && is_id_matches(a, b)
end

function Base.isequal(a::Dependency, b::Dependency)
    return a == b
end

#The regular unique works using Set and will not work correctly, because in our case "equal" dependencies have different hashes
#This implementation is not algorithmically efficient (O(n^2)), but it is quite suitable for checking the uniqueness of outputs
function check_unique(deps::Vector{<:Dependency})
    tmp = Dependency[]
    for dep in deps
        dep in tmp && return false
        push!(tmp, dep)
    end
    return true
end

is_id_matches(a::Dependency{Trait1, String}, b::Dependency{Trait2, String}) where {Trait1, Trait2} = a.id == b.id
is_id_matches(a::Dependency{Trait1, String}, b::Dependency{Trait2, <:NamedTuple}) where {Trait1, Trait2} = false
is_id_matches(a::Dependency{Trait1, <:NamedTuple}, b::Dependency{Trait2, String}) where {Trait1, Trait2} = false
function is_id_matches(a::Dependency{Trait1, <:NamedTuple}, b::Dependency{Trait2, <:NamedTuple}) where {Trait1, Trait2}
    (Set(keys(a.id)) != Set(keys(b.id))) && return false

    for key in keys(a.id)
        a_value = a.id[key]
        b_value = b.id[key]

        (a_value == b_value) && continue

        a_wild = is_wild(a_value)
        b_wild = is_wild(b_value)

        (!a_wild && !b_wild) && return false #Both not wild
        !(a_wild && b_wild) && continue #One wild, one not
        ((a_value == ALL) || (b_value == ALL)) && continue #at least one is ALL
        ((a_value == MATCH) || (b_value == MATCH)) && return false #one is MATCH and one is ALLSMALLER

    end
    return true
end

struct CallbackDeps
    output ::Vector{<:Output}
    input ::Vector{<:Input}
    state ::Vector{<:State}
    multi_out::Bool
    CallbackDeps(output, input, state, multi_out) = new(output, input, state, multi_out)
    CallbackDeps(output::Output, input, state = State[]) = new(output, input, state, false)
    CallbackDeps(output::Vector{<:Output}, input, state = State[]) = new(output, input, state, true)
end

Base.convert(::Type{Vector{<:T}}, v::T) where {T<:Dependency} = [v]

struct ClientsideFunction
    namespace ::String
    function_name ::String
end
struct Callback
    func ::Union{Function, ClientsideFunction}
    dependencies ::CallbackDeps
end

is_multi_out(cb::Callback) = cb.dependencies.multi_out == true
get_output(cb::Callback) = cb.dependencies.output
get_output(cb::Callback, i) = cb.dependencies.output[i]
first_output(cb::Callback) = first(cb.dependencies.output)

struct PreventUpdate <: Exception

end

struct NoUpdate
end

no_update() = NoUpdate()

const ExternalSrcType = Union{String, Dict{String, String}}


text:
import GeometryBasics: AbstractPoint, Point, Point2, AbstractPolygon, area, Point3, Polygon

import GeneralPolygonClipper as gpc

export ConvexPolygon, area, centroid, normal, coordinates, simple2poly

struct ConvexPolygon{Dim,T,P<:AbstractPoint{Dim,T},L<:AbstractVector{P}} <: AbstractPolygon{Dim,T}
    contour::L
end

ConvexPolygon(x::AbstractVector, y::AbstractVector) =
    ConvexPolygon([Point2(xx, yy) for (xx,yy) in zip(x,y)])

ConvexPolygon(x::AbstractVector, y::AbstractVector, z::AbstractVector) =
    ConvexPolygon([Point3(xx, yy,zz) for (xx,yy,zz) in zip(x,y,z)])

#ConvexPolygon(pts::E) where {E<:AbstractVector{P}} where {P<:AbstractPoint{Dim,T}} where {Dim,T} = ConvexPolygon{Dim,T,P,E}(pts)

import GeometryBasics.coordinates
coordinates(p::ConvexPolygon) = p.contour

area(p::ConvexPolygon{2,T}) where {T} = area(coordinates(p))

normal(p::ConvexPolygon{2,T}) where {T} = area(coordinages(p))

crossprod(u::Point3, v::Point3) = (u[2]*v[3] - u[3]*v[2],
                                   u[3]*v[1] - u[1]*v[3],
                                   u[1]*v[2] - u[2]*v[1])

function normal(p::ConvexPolygon{3,T}) where {T}
    v = coordinates(p)
    nv = length(v)
    x,y,z = crossprod(v[end], v[begin])

    for i in 2:nv
        prd = crossprod(v[i-1], v[i])
        x += prd[1]
        y += prd[2]
        z += prd[3]
    end

    return Point3(x/2, y/2, z/2)

end

area(p::ConvexPolygon{3}) = hypot(normal(p)...)

"""
`centroid(p)`

Computes the centroid of a [`ConvexPolygon`](@ref).
"""
function centroid(p::ConvexPolygon{2,T}) where {T}

    A = area(p)

    v = coordinates(p)
    nv = length(v)

    # xᵢyᵢ₊₁ - xᵢ₊₁yᵢ  First index = last index
    tmp = v[end][1]*v[begin][2] - v[begin][1]*v[end][2]

    Cx = (v[end][1] + v[begin][1]) * tmp
    Cy = (v[end][2] + v[begin][2]) * tmp

    for i in firstindex(v):lastindex(v)-1
        tmp = v[i][1]*v[i+1][2] - v[i+1][1]*v[i][2]
        Cx += (v[i][1] + v[i+1][1]) * tmp
        Cy += (v[i][2] + v[i+1][2]) * tmp
    end

    return Point2{T}(Cx/6A, Cy/6A)
end

function centroid(p::ConvexPolygon{3,T}) where {T}
    v = coordinates(p)
    nv = length(v)
    nrm = normal(p)
    A = hypot(nrm...)
    n⃗ = nrm ./ A

    C = Point3{T}(0,0,0)
    v₀ = v[begin]
    for i = (firstindex(v)+1):(lastindex(v)-1)
        centr = (v₀ + v[i] + v[i+1]) ./ 3
        Aᵢ⃗ = crossprod(v[i]-v₀, v[i+1]-v₀) ./ 2
        s = sign(sum(Aᵢ⃗ .* n⃗))
        C += s * centr .* hypot(Aᵢ⃗...)
    end
    return C

end

import GeneralPolygonClipper: Vertex
import Base.union

"""
`gpc_operations(op::GPCOperation, p1, p2)`

Performs boolean operations on 2D polygons. First the polygons
are converted to [`GPCPolygon`](@ref). Then, performs the boolean
operation specified by parameter `op` (see enum [`GPCOperation`](@ref)).

Finally, convert the results to [`ConvexPolygon`](@ref). There can be
more than one resulting polygon depending on the exact nature of `p1`, `p2`
and `op`.

"""
function gpc_operations(op::gpc.GPCOperation, p1::ConvexPolygon{2,Float64},
                        p2::ConvexPolygon{2,Float64})

    gpc1 = gpc.GPCPolygon([false], [Vertex(p[1], p[2]) for p in coordinates(p1)])
    gpc2 = gpc.GPCPolygon([false], [Vertex(p[1], p[2]) for p in coordinates(p2)])

    gpcout = gpc.gpc_polygon_clip(op, gpc1, gpc2)

    # If both polygons were convex, the output would be a single polygon
    # In other cases there can be more than one output polygon.

    return ConvexPolygon.([[Point2{Float64}(v.x, v.y) for v in g] for g in gpcout.contours])

end

import Base.union
union(p1::ConvexPolygon{2,Float64}, p2::ConvexPolygon{2,Float64}) =
    gpc_operations(gpc.GPC_UNION, p1, p2)
import Base.intersect
intersect(p1::ConvexPolygon{2,Float64}, p2::ConvexPolygon{2,Float64}) =
    gpc_operations(gpc.GPC_INT, p1, p2)
import Base.(-)
(-)(p1::ConvexPolygon{2,Float64}, p2::ConvexPolygon{2,Float64}) =
    gpc_operations(gpc.GPC_DIFF, p1, p2)

import Base.xor
xor(p1::ConvexPolygon{2,Float64}, p2::ConvexPolygon{2,Float64}) =
    gpc_operations(gpc.GPC_XOR, p1, p2)

simple2poly(p::SimplePolygon{2,Float64}) = Polygon(coordinates(p))


text:
using Graphs, SimpleWeightedGraphs
include("debruijnnets.jl")

"""
Sample from a directed GNP ensemble.
"""
function gnp_graph(N, p)
    edgedict = Dict{Tuple, Integer}()
    for u in range(1, N)
        for v in range(1, N)
            if u == v
                continue
            end
            if rand() < p
                edgedict[(u, v)] = 1.
            end
        end
    end
    src, dest, weights = vectors_from_edgedict(edgedict)
    G = SimpleWeightedDiGraph(src, dest, weights)
    return G
end

"""
Sample from a directed Preferential Attachment ensemble
where a back link is sampled with uniform probability.
"""
function pa_graph(N, m)
    edgedict = Dict{Tuple, Integer}()
    inneighbors = Dict(node=> 0 for node in range(1, N))
    # Start with a clique on m nodes
    for i in range(1, m)
        for j in range(1, m)
            if i !=j
                edgedict[(i,j)] = 1
                inneighbors[j] += 1
            end
        end
    end

    # Grow the network up to N nodes
    for n in range(m+1, N)
        # Compute degree weights of nodes so far
        nodes = [u for u in range(1, n-1)]
        node_weights = [Float64(inneighbors[ne]) for ne in nodes]
        node_weights ./= sum(node_weights)
        node_weights = ProbabilityWeights(node_weights)
        # Sample neighbors for node n
        out_neighbors = sample(nodes, node_weights, m)
        for ne in out_neighbors
            edgedict[(n, ne)] = 1
            inneighbors[ne] += 1
            # With uniform probability, add a back link
            if rand() < (1. / length(out_neighbors))
                edgedict[(ne,n)] = 1
                inneighbors[n] += 1
            end
        end
    end
    src, dest, weights = vectors_from_edgedict(edgedict)
    G = SimpleWeightedDiGraph(src, dest, weights)
    return G
end


text:
using PhantomSegmentation
using Test
using TestSetExtensions

using DICOM
using DICOMUtils


text:
# functions that are used to compute financial frictions
# steady-state values from parameter values
@inline function ζ_spb_fn(z, σ, spr)
    zetaratio = ζ_bω_fn(z, σ, spr)/ζ_zω_fn(z, σ, spr)
    nk = nk_fn(z, σ, spr)
    return -zetaratio/(1-zetaratio)*nk/(1-nk)
end

@inline function ζ_bω_fn(z, σ, spr)
    nk          = nk_fn(z, σ, spr)
    μstar       = μ_fn(z, σ, spr)
    ω_star       = ω_fn(z, σ)
    Γstar       = Γ_fn(z, σ)
    Gstar       = G_fn(z, σ)
    dΓ_dω_star   = dΓ_dω_fn(z)
    dG_dω_star   = dG_dω_fn(z, σ)
    d2Γ_dω2star = d2Γ_dω2_fn(z, σ)
    d2G_dω2star = d2G_dω2_fn(z, σ)
    return ω_star*μstar*nk*(d2Γ_dω2star*dG_dω_star - d2G_dω2star*dΓ_dω_star)/
        (dΓ_dω_star - μstar*dG_dω_star)^2/spr/(1 - Γstar + dΓ_dω_star*(Γstar - μstar*Gstar)/
            (dΓ_dω_star - μstar*dG_dω_star))
end

@inline function ζ_zω_fn(z, σ, spr)
    μstar = μ_fn(z, σ, spr)
    return ω_fn(z, σ)*(dΓ_dω_fn(z) - μstar*dG_dω_fn(z, σ))/
        (Γ_fn(z, σ) - μstar*G_fn(z, σ))
end

nk_fn(z, σ, spr) = 1 - (Γ_fn(z, σ) - μ_fn(z, σ, spr)*G_fn(z, σ))*spr
μ_fn(z, σ, spr)  =
    (1 - 1/spr)/(dG_dω_fn(z, σ)/dΓ_dω_fn(z)*(1 - Γ_fn(z, σ)) + G_fn(z, σ))
ω_fn(z, σ)        = exp(σ*z - σ^2/2)
G_fn(z, σ)        = cdf(Normal(), z-σ)
Γ_fn(z, σ)        = ω_fn(z, σ)*(1 - cdf(Normal(), z)) + cdf(Normal(), z-σ)
dG_dω_fn(z, σ)    = pdf(Normal(), z)/σ
d2G_dω2_fn(z, σ)  = -z*pdf(Normal(), z)/ω_fn(z, σ)/σ^2
dΓ_dω_fn(z)       = 1 - cdf(Normal(), z)
d2Γ_dω2_fn(z, σ)  = -pdf(Normal(), z)/ω_fn(z, σ)/σ
dG_dσ_fn(z, σ)    = -z*pdf(Normal(), z-σ)/σ
d2G_dωdσ_fn(z, σ) = -pdf(Normal(), z)*(1 - z*(z-σ))/σ^2
dΓ_dσ_fn(z, σ)    = -pdf(Normal(), z-σ)
d2Γ_dωdσ_fn(z, σ) = (z/σ-1)*pdf(Normal(), z)


text:
using DDF

using StaticArrays
using Test

@testset "Fun D=$D P=$P R=$R" for D in 0:Dmax, P in (Pr, Dl), R in 0:D
    S = Float64
    T = Rational{Int64}
    mfd = hypercube_manifold(Val(D), S)

    z = zero(Fun{D,P,R,D,S,T}, mfd)
    e = id(Fun{D,P,0,D,S,SVector{D,S}}, mfd)
    f = rand(Fun{D,P,R,D,S,T}, mfd)
    g = rand(Fun{D,P,R,D,S,T}, mfd)
    h = rand(Fun{D,P,R,D,S,T}, mfd)
    a = rand(T)
    b = rand(T)

    # Vector space

    @test +f == f
    @test (f + g) + h == f + (g + h)

    @test z + f == f
    @test f + z == f

    @test f + (-f) == z
    @test (-f) + f == z
    @test f - g == f + (-g)

    @test (a * b) * f == a * (b * f)
    @test f * (a * b) == (f * a) * b

    @test one(T) * a == a
    @test one(T) * f == f
    @test f * one(T) == f

    @test zero(T) * a == zero(T)
    @test zero(T) * f == z
    @test f * zero(T) == z

    if a != zero(T)
        @test a * inv(a) == one(T)

        @test inv(a) * (a * f) == f
        @test (f * a) * inv(a) == f

        @test inv(a) * f == a \ f
        @test f * inv(a) == f / a
    end

    @test (a + b) * f == a * f + b * f
    @test f * (a + b) == f * a + f * b

    @test a * (f + g) == a * f + a * g
    @test (f + g) * a == f * a + g * a

    # Nonlinear transformations

    e2 = map(x -> 2x, e)
    @test [2x for x in e] == [x for x in e2]
    @test map(+, f, g) == f + g
end


text:
using Distributions

export value, delta, gammamatrix, Δrtk, rΔrtk, rtk, srtk,
       aggrstress, UP, DOWN

## Constructors -------------------------------------------------
function RiskFactor(σ::Vector{Float64},
                    corr::Array{Float64, 2},
                    x0::Vector{Float64},
                    h::Vector{Float64},
                    add::Vector{Bool})
  return RiskFactor((σ * σ') .* corr, x0, h, add)
end

## Valuation ----------------------------------------------------
"""
`value(t::Int, zb::ZeroBond, x::Vector{Float64}, rf::RiskFactor,
  cap_mkt::SSTCapMkt)`

Calculates the value of a zero bond (after sensitivity `x`)
"""
function value(t::Int,
               zb::ZeroBond,
               x::Vector{Float64},
               rf::RiskFactor,
               cap_mkt::SSTCapMkt)
  x_spot = cap_mkt.spot[zb.index] + x[zb.index]
  val = zb.nom / (1 + x_spot)^zb.τ
  if t != 0
    val *= (1 + cap_mkt.spot[t] + x[t])^t
  end
  return val
end

"""
`value(t::Int, si::StockIndex,x::Vector{Float64}, rf::RiskFactor,
  cap_mkt::SSTCapMkt)`

Calculates the value of a stock (after sensitivity `x`)
"""
function value(t::Int,
               si::StockIndex,
               x::Vector{Float64},
               rf::RiskFactor,
               cap_mkt::SSTCapMkt)
  return x[si.index] * si.nom * (1 + cap_mkt.stock_increase)^t
end

"""
`value(t::Int, assts::Vector{Asset}, x::Vector{Float64},
  rf::RiskFactor, cap_mkt::SSTCapMkt)`

Calculates the value of the asset portfolio
(after sensitivity `x`)
"""
function value(t::Int,
               assts::Vector{Asset},
               x::Vector{Float64},
               rf::RiskFactor,
               cap_mkt::SSTCapMkt)
  val = 0.0
  for asset in assts
    val += value(t, asset, x, rf, cap_mkt)
  end
  return val
end

"""
`value(t::Int, liabs::Liabilities, x::Vector{Float64},
  rf::RiskFactor, cap_mkt::SSTCapMkt)`

Calculates the value of the liability portfolio
(after sensitivity `x`)
"""
function value(t::Int,
               liabs::Liabilities,
               x::Vector{Float64},
               rf::RiskFactor,
               cap_mkt::SSTCapMkt)
  T = length(cap_mkt.spot)
  x_spot = cap_mkt.spot[1:T] + x[1:T]
  x_mort = x[liabs.index_mort]
  val = 0.0
  for τ in 1:length(liabs.B_PX)
    val +=
      prod(1 .- liabs.qx[1:τ] .* x_mort) *
      liabs.B_PX[τ] / (1 + x_spot[τ])^τ
  end
  if t != 0
    val *= (1 + cap_mkt.spot[t] + x[t])^t
  end
  return val
end

"""
`rtk(t::Int, assets::Vector{Asset}, liabs::Liabilities,
  x::Vector{Float64}, rf::RiskFactor, cap_mkt::SSTCapMkt)`

Calculates the risk bearing capital (after sensitivity `x`)
"""
rtk(t::Int,
    assets::Vector{Asset},
    liabs::Liabilities,
    x::Vector{Float64},
    rf::RiskFactor,
    cap_mkt::SSTCapMkt) =
  value(t, assets, x, rf, cap_mkt) -
  value(t,liabs, x, rf, cap_mkt)

## capital calculation ------------------------------------------
const UP, DOWN = 1, -1

"""
`srtk(shift::Int, assets::Vector{Asset}, liabs::Liabilities,
  rf::RiskFactor, cap_mkt::SSTCapMkt)`

Calculates (linear) sensitivities for rtk
"""
function srtk(shift::Int,
             assets::Vector{Asset},
             liabs::Liabilities,
             rf::RiskFactor,
             cap_mkt::SSTCapMkt)
  x = deepcopy(rf.x0)
  n = length(x)
  rtk_shift = Array(Float64, n)
  for i = 1:n
    x[i] += shift * rf.h[i]
    rtk_shift[i] = rtk(1, assets, liabs, x, rf, cap_mkt)
    x[i] -= shift * rf.h[i]  ## restore old value for y[i]
  end
  return rtk_shift
end

"""
`srtk(shift_1::Int, shift_2::Int, assets::Vector{Asset},
  liabs::Liabilities, rf::RiskFactor, cap_mkt::SSTCapMkt)`

Calculates quadratic sensitivities for rtk
"""
function srtk(shift_1::Int,
             shift_2::Int,
             assets::Vector{Asset},
             liabs::Liabilities,
             rf::RiskFactor,
             cap_mkt::SSTCapMkt)
  x = deepcopy(rf.x0)
  n = length(x)
  rtk_shift_shift = Array(Float64, n, n)
  for i = 1:n
    for k = 1:n
      x[i] += shift_1 * rf.h[i]
      x[k] += shift_2 * rf.h[k]
      rtk_shift_shift[i, k] =
        rtk(1, assets, liabs, x, rf, cap_mkt)
      x[i] -= shift_1 * rf.h[i]  ## restore old value for x[i]
      x[k] -= shift_2 * rf.h[k]  ## restore old value for x[i]
    end
  end
  return rtk_shift_shift
end

Δ(rf::RiskFactor) =
  Float64[rf.add[i] ?  rf.h[i] : rf.x0[i] * rf.h[i]
          for i = 1:length(rf.x0)]

"""
`delta(assets::Vector{Asset}, liabs::Liabilities, rf::RiskFactor,
  cap_mkt::SSTCapMkt)`

Calculates the δ-vector
"""
delta(assets::Vector{Asset},
      liabs::Liabilities,
      rf::RiskFactor,
      cap_mkt::SSTCapMkt) =
  (srtk(UP, assets, liabs, rf, cap_mkt) -
     srtk(DOWN, assets, liabs, rf, cap_mkt)) ./ (2Δ(rf))

"""
`gammamatrix(assets::Vector{Asset}, liabs::Liabilities,
  rf::RiskFactor, cap_mkt::SSTCapMkt)`

Calculates the Γ-matrix
"""
function gammamatrix(assets::Vector{Asset},
               liabs::Liabilities,
               rf::RiskFactor,
               cap_mkt::SSTCapMkt)
  Δx = Δ(rf)
  rtk_uu = srtk(UP, UP, assets, liabs, rf, cap_mkt)
  rtk_ud = srtk(UP, DOWN, assets, liabs, rf, cap_mkt)
  rtk_du = srtk(DOWN, UP, assets, liabs, rf, cap_mkt)
  rtk_dd = srtk(DOWN, DOWN, assets, liabs, rf, cap_mkt)
  Γ_diag =
    (srtk(UP, assets, liabs, rf, cap_mkt) +
       srtk(DOWN, assets, liabs, rf, cap_mkt) -
       2rtk(1,assets, liabs, rf.x0, rf, cap_mkt)) ./ (Δx .* Δx)
  Γ = Array(Float64, length(rf.x0), length(rf.x0))
  for i = 1:length(rf.x0)
    for k = 1:(i-1)
      Γ[i,k] = (rtk_uu[i,k] -
                  rtk_ud[i,k] -
                  rtk_du[i,k] +
                  rtk_dd[i,k]) / (4 * Δx[i] * Δx[k])
      Γ[k,i] = Γ[i,k]
    end
    Γ[i,i] = Γ_diag[i]
  end
  return Γ
end

"""
`Δrtk(Δx::Vector{Float64}, δ::Vector{Float64},
  Γ::Matrix{Float64})`

Calculates the shocked rtk based on shocks Δx
"""
Δrtk(Δx::Vector{Float64},
     δ::Vector{Float64},
     Γ::Matrix{Float64}) = (Δx ⋅ δ + 0.5 * Δx' * Γ * Δx)[1]

## random values for Δrtk
function rΔrtk(n_scen::Int,
               assets::Vector{Asset},
               liabs::Liabilities,
               rf::RiskFactor,
               cap_mkt::SSTCapMkt,
               x_index::Vector{Int},
               )
  r_Δx = rand(MvNormal(zeros(Float64, length(x_index)),
                       rf.Σ[x_index, x_index]),
              n_scen)
  r_Δrtk = Array(Float64, n_scen)
  δ = delta(assets, liabs, rf, cap_mkt)[x_index]
  Γ = gammamatrix(assets, liabs, rf, cap_mkt)[x_index, x_index]
  for mc in 1:n_scen
    r_Δrtk[mc] = Δrtk(r_Δx[:, mc], δ, Γ)
  end
  return r_Δrtk
end

##  Aggregate stress-scenarios to randomly generated values
##    r_Δrtk_no_stress.
##  We use the same approximation as in the shift method,
##    namely, two different stress scenarios cannot happen
##    in the same year.
function aggrstress(stress::Stress, r_Δrtk_no_stress)
  r_Δrtk = deepcopy(r_Δrtk_no_stress)
  n_scen = length(r_Δrtk)
  i = 0
  for scen = 1:stress.n
    if stress.target[scen]
      n_adj = floor(Integer, n_scen * stress.prob[scen])
      for j = 1:n_adj
        r_Δrtk[i + j] += min(0, stress.Δrtk[scen])
      end
      i += n_adj
    end
  end
  return r_Δrtk
end


text:
using .Networks

using Dates
using JuMP
using Printf

@with_kw mutable struct BalanceMarket <: AbstractMarket
    configs = EnergyMarketConfigs(problem_name = "BalanceMarket",
                                CONSIDER_TSOACTIONS_LIMITATIONS=true,
                                CONSIDER_TSOACTIONS_IMPOSITIONS=true,
                                REF_SCHEDULE_TYPE=PSCOPF.TSO()  #This can be the
                                                                # TSO if TSOBilevel updates the tso_schedule according to the preceding market decided values
                                                                # or the market too cause DECIDED values refer to the preceding ECH and it is the market who launched last at the preceding ech
                                )
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using {{PKG}}
using Documenter

DocMeta.setdocmeta!({{PKG}}, :DocTestSetup, :(using {{PKG}}); recursive=true)

makedocs(;
    modules=[{{PKG}}],
    authors="xKDR Forum",
    repo="https://github.com/xKDR/{{PKG}}.jl/blob/{commit}{path}#{line}",
    sitename="${{PKG}}.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xKDR.github.io/{{PKG}}.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/xKDR/{{PKG}}.jl",
    target = "build",
    devbranch="main"
)


text:
clenshawcurtis_plan(μ) = length(μ) > 1 ? FFTW.plan_r2r!(μ, FFTW.REDFT00) : ones(μ)'

"""
Compute nodes and weights of the Clenshaw—Curtis quadrature rule with a Jacobi weight.
"""
clenshawcurtis{T<:AbstractFloat}(N::Int,α::T,β::T) = clenshawcurtis(N,α,β,clenshawcurtis_plan(zeros(T,N)))
clenshawcurtis{T<:AbstractFloat}(N::Int,α::T,β::T,plan) = T[cospi(k/(N-one(T))) for k=0:N-1],clenshawcurtisweights(N,α,β,plan)

"""
Compute weights of the Clenshaw—Curtis quadrature rule with a Jacobi weight.
"""
clenshawcurtisweights{T<:AbstractFloat}(N::Int,α::T,β::T) = clenshawcurtisweights(N,α,β,clenshawcurtis_plan(zeros(T,N)))
function clenshawcurtisweights{T<:AbstractFloat}(N::Int,α::T,β::T,plan)
    μ = chebyshevjacobimoments1(N,α,β)
    scale!(μ,inv(N-one(T)))
    plan*μ
    μ[1]/=2;μ[N]/=2
    return μ
end

# Chebyshev-T coefficients to values at Clenshaw-Curtis nodes

applyTN_plan(x) = length(x) > 1 ? FFTW.plan_r2r!(x, FFTW.REDFT00) : ones(x)'

applyTN!{T<:AbstractFloat}(x::Vector{T}) = applyTN!(x,applyTN_plan(x))
function applyTN!{T<:AbstractFloat}(x::Vector{T},plan)
    x[1] *= 2; x[end] *=2
    plan*x
    scale!(x,half(T))
end
applyTN{T<:AbstractFloat}(x::Vector{T},plan) = applyTN!(copy(x),plan)
applyTN{T<:AbstractFloat}(x::Vector{T}) = applyTN!(copy(x))

# Values at Clenshaw-Curtis nodes to Chebyshev-T coefficients

applyTNinv_plan(x) = length(x) > 1 ? FFTW.plan_r2r!(x, FFTW.REDFT00) : ones(x)'

applyTNinv!{T<:AbstractFloat}(x::Vector{T}) = applyTNinv!(x,applyTNinv_plan(x))
function applyTNinv!{T<:AbstractFloat}(x::Vector{T},plan)
    plan*x
    x[1] /= 2;x[end] /= 2
    scale!(x,inv(length(x)-one(T)))
end
applyTNinv{T<:AbstractFloat}(x::Vector{T},plan) = applyTNinv!(copy(x),plan)
applyTNinv{T<:AbstractFloat}(x::Vector{T}) = applyTNinv!(copy(x))

# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1

applyUN_plan(x) = length(x) > 0 ? FFTW.plan_r2r!(x, FFTW.RODFT00) : ones(x)'

applyUN!{T<:AbstractFloat}(x::AbstractVector{T}) = applyUN!(x,applyUN_plan(x))
function applyUN!{T<:AbstractFloat}(x::AbstractVector{T},plan)
    plan*x
    scale!(x,half(T))
end
applyUN{T<:AbstractFloat}(x::AbstractVector{T},plan) = applyUN!(copy(x),plan)
applyUN{T<:AbstractFloat}(x::AbstractVector{T}) = applyUN!(copy(x))


text:
"""
    RuleConfig{T}

The configuration for what rules to use.
`T`: **traits**. This should be a `Union` of all special traits needed for rules to be
allowed to be defined for your AD. If nothing special this should be set to `Union{}`.

**AD authors** should define a subtype of `RuleConfig` to use when calling `frule`/`rrule`.

**Rule authors** can dispatch on this config when defining rules.
For example:
```julia
# only define rrule for `pop!` on AD systems where mutation is supported.
rrule(::RuleConfig{>:SupportsMutation}, typeof(pop!), ::Vector) = ...

# this definition of map is for any AD that defines a forwards mode
rrule(conf::RuleConfig{>:HasForwardsMode}, typeof(map), ::Vector) = ...

# this definition of map is for any AD that only defines a reverse mode.
# It is not as good as the rrule that can be used if the AD defines a forward-mode as well.
rrule(conf::RuleConfig{>:Union{NoForwardsMode, HasReverseMode}}, typeof(map), ::Vector) = ...
```

For more details see [rule configurations and calling back into AD](@ref config).
"""
abstract type RuleConfig{T} end

# Broadcast like a scalar
Base.Broadcast.broadcastable(config::RuleConfig) = Ref(config)

abstract type ReverseModeCapability end

"""
HasReverseMode

This trait indicates that a `RuleConfig{>:HasReverseMode}` can perform reverse mode AD.
If it is set then [`rrule_via_ad`](@ref) must be implemented.
"""
struct HasReverseMode <: ReverseModeCapability end

"""
NoReverseMode

This is the complement to [`HasReverseMode`](@ref). To avoid ambiguities [`RuleConfig`]s
that do not support performing reverse mode AD should be `RuleConfig{>:NoReverseMode}`.
"""
struct NoReverseMode <: ReverseModeCapability end

abstract type ForwardsModeCapability end

"""
HasForwardsMode

This trait indicates that a `RuleConfig{>:HasForwardsMode}` can perform forward mode AD.
If it is set then [`frule_via_ad`](@ref) must be implemented.
"""
struct HasForwardsMode <: ForwardsModeCapability end

"""
NoForwardsMode

This is the complement to [`HasForwardsMode`](@ref). To avoid ambiguities [`RuleConfig`]s
that do not support performing forwards mode AD should be `RuleConfig{>:NoForwardsMode}`.
"""
struct NoForwardsMode <: ForwardsModeCapability end

"""
    frule_via_ad(::RuleConfig{>:HasForwardsMode}, ȧrgs, f, args...; kwargs...)

This function has the same API as [`frule`](@ref), but operates via performing forwards mode
automatic differentiation.
Any `RuleConfig` subtype that supports the [`HasForwardsMode`](@ref) special feature must
provide an implementation of it.

See also: [`rrule_via_ad`](@ref), [`RuleConfig`](@ref) and the documentation on
[rule configurations and calling back into AD](@ref config)
"""
function frule_via_ad end

"""
    rrule_via_ad(::RuleConfig{>:HasReverseMode}, f, args...; kwargs...)

This function has the same API as [`rrule`](@ref), but operates via performing reverse mode
automatic differentiation.
Any `RuleConfig` subtype that supports the [`HasReverseMode`](@ref) special feature must
provide an implementation of it.

See also: [`frule_via_ad`](@ref), [`RuleConfig`](@ref) and the documentation on
[rule configurations and calling back into AD](@ref config)
"""
function rrule_via_ad end


text:
module TestTSML
using Test
using TSML.System

include("test_valdate.jl")
include("test_decisiontree.jl")
include("test_readerwriter.jl")
include("test_statifier.jl")
include("test_monotonicer.jl")
include("test_cliwrapper.jl")
include("test_tsclassifier.jl")
include("test_outliernicer.jl")

if LIB_SKL_AVAILABLE
    include("test_scikitlearn.jl")
else
    @info "Skipping scikit-learn tests."
end
if LIB_CRT_AVAILABLE
    include("test_caret.jl")
else
    @info "Skipping CARET tests."
end

end


text:
using HMMs
using Random
using Brobdingnag

import HMMs: GenTDist

Random.seed!(1)
(mu1,sigma1,nu1) = HMMs.randomparameters( GenTDist )
t1 = GenTDist( mu1, sigma1, nu1 )
y1 = rand( t1, 1_000_000 );
w1 = ones(length(y1));

parameters = HMMs.randomparameters( GenTDist )
HMMs.fit_mle!( GenTDist, parameters, y1, w1, Dict{Symbol,Any}() )
@assert( maximum(abs.(parameters - [mu1,sigma1,nu1])) < 1e-2 )

(mu2,sigma2,nu2) = HMMs.randomparameters( GenTDist )
t2 = GenTDist( mu2, sigma2, nu2 )
y2 = fill(NaN,1_000_000);
w2 = fill(NaN,length(y2));
for i = 1:length(y2)
    w2[i] = rand()
    y2[i] = rand() < w2[i] ? rand( t1 ) : rand( t2 )
end

parameters = HMMs.randomparameters( GenTDist )
HMMs.fit_mle!( GenTDist, parameters, y2, w2, Dict{Symbol,Any}() )
# this doesn't seem to work but I'm not sure it's well-founded
#@assert( maximum(abs.(parameters - [mu1,sigma1,nu1])) < 1e-2 )

graph = HMMs.fullyconnected(1)
hmm1 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=1 )
y1 = rand( hmm1, 10_000 );

hmm2 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y1 );
HMMs.em( hmm2, debug=2 )

graph = HMMs.fullyconnected(2)
hmm1 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=1 )
y1 = rand( hmm1, 100_000 );

hmm2 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y1 );
HMMs.em( hmm2, debug=2 )

hmm1 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=1 )
y1 = rand( hmm1, 10_000 );

hmm2 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y1 );
HMMs.em( hmm2, debug=2 )

hmm2 = HMMs.randomhmm( graph, dist=GenTDist, calc=Brob, seed=3 )
HMMs.setobservations( hmm2, y1 );
HMMs.em( hmm2, debug=2 )


text:
module DEAC
include("./DEACfitness.jl")
include("./DEACisf_fitness.jl")
include("./DEACpopulation.jl")
include("./DEACreject.jl")
include("./DEACfrequency.jl")

using DelimitedFiles
using UUIDs
using Random
using Statistics
using LinearAlgebra

export deac

function mutant_indexes(rng::MersenneTwister,populationSize::Int64,idx0::Int64)
    idx1 = idx0;
    idx2 = idx0;
    idx3 = idx0;
    while idx1 == idx0
        idx1 = rand(rng,1:populationSize);
    end

    while idx2 == idx0 || idx2 == idx1
        idx2 = rand(rng,1:populationSize);
    end
    while idx3 == idx0 || idx3 == idx1 || idx3 == idx2
        idx3 = rand(rng,1:populationSize);
    end
    idx1,idx2,idx3
end

function set_differential_weights(rng::MersenneTwister,new_differential_weights::Array{Float64,1},
                              differential_weights::Array{Float64,1},SAdifferentialWeightProb::Float64,
                              SAdifferentialWeightShift::Float64,SAdifferentialWeight::Float64)
    @inbounds for i = eachindex(new_differential_weights)
        if rand(rng) < SAdifferentialWeightProb
            #new_differential_weights[i] = SAdifferentialWeightShift + SAdifferentialWeight*rand(rng);
            new_differential_weights[i] = 2.0*rand(rng);
        else
            new_differential_weights[i] = differential_weights[i];
        end
    end
    nothing
end

function set_crossover_probs(rng::MersenneTwister,new_crossover_probs::Array{Float64,1},
                                   crossover_probs::Array{Float64,1},SAcrossoverProb::Float64)
    @inbounds for i = eachindex(new_crossover_probs)
        if rand(rng) < SAcrossoverProb
            new_crossover_probs[i] = rand(rng);
        else
            new_crossover_probs[i] = crossover_probs[i];
        end
    end
    nothing
end

function set_mutateInd(rng::MersenneTwister,mutateInd::BitArray{2},new_crossover_probs::Array{Float64,1})
    @inbounds for j = 1:size(mutateInd,2)
        cR = new_crossover_probs[j];
        for i = 1:size(mutateInd,1)
            mutateInd[i,j] = (rand(rng) < cR);
        end
    end
    nothing
end

function mutate(rng::MersenneTwister,mP::Array{Float64,2},P::Array{Float64,2},
                 mIdx::BitArray{2},new_differential_weights::Array{Float64,1})
    @inbounds for j = 1:size(mP,2)
        F = new_differential_weights[j];
        mIdx1,mIdx2,mIdx3 = mutant_indexes(rng,size(mP,2),j);
        for i = 1:size(mP,1)
            if mIdx[i,j]
                mP[i,j] = abs(P[i,mIdx1] + (F * (P[i,mIdx2]-P[i,mIdx3])));
            else
                mP[i,j] = P[i,j];
            end
        end
    end
    nothing
end

function generate_pop(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1},gpfn::String)
    gpfunc = getfield(DEACpopulation, Symbol(gpfn));
    gpfunc(rng,genomeSize,populationSize,omega)
end

function replace_pop(p::Array{Float64,2},c::Array{Float64,2},rInd::BitArray{1})
    @inbounds for i = 1:size(rInd,1)
        idx = rInd[i];
        if idx
            for j = 1:size(p,1)
                p[j,i] = c[j,i];
            end
        end
    end
    nothing
end

function replace_fitness(pFitness::Array{Float64,1},cFitness::Array{Float64,1},rInd::BitArray{1})
    @inbounds for i = 1:size(rInd,1)
        idx = rInd[i];
        if idx
            pFitness[i] = cFitness[i];
        end
    end
    nothing
end

function replace_crossover_probs(crossover_probs::Array{Float64,1},new_crossover_probs::Array{Float64,1},rInd::BitArray{1})
    replace_fitness(crossover_probs,new_crossover_probs,rInd);
end

function replace_differential_weights(differential_weights::Array{Float64,1},new_differential_weights::Array{Float64,1},rInd::BitArray{1})
    replace_fitness(differential_weights,new_differential_weights,rInd);
end

function set_minP(minP::Array{Float64,1},P::Array{Float64,2},minidx::Int64)
    @inbounds for i = 1:size(P,1)
        minP[i] = P[i,minidx];
    end
    nothing
end

function set_isf_term(imaginary_time::Array{Float64,1},frequency_bins::Array{Float64,1},beta::Float64)
    b = beta;
    isf_term = Array{Float64,2}(undef,size(imaginary_time,1),size(frequency_bins,1));
    for i in 1:size(imaginary_time,1)
        t = imaginary_time[i];
        bo2mt = b/2 - t;
        for j in 1:size(frequency_bins,1)
            f = frequency_bins[j];
            isf_term[i,j] = cosh(bo2mt*f);
        end
    end
    isf_term
end

function beam_distance(d::Array{Float64,1},frequency::Array{Float64,1},P::Array{Float64,2})
    #FIXME this method is slow
    for j in 1:(size(P,2))
        d[j] = 1.0;
        for i in 1:(size(frequency,1) - 2)
            d[j] += abs((P[i+2,j] - P[i,j])*frequency[i+1] - (frequency[i+2] - frequency[i])*P[i+1,j] + frequency[i+2]*P[i,j] - frequency[i]*P[i+2,j])/sqrt((frequency[i+2] - frequency[i])^2 + (P[i+2,j] - P[i,j])^2);
        end
        nothing
    end
end

function deac(  imaginary_time::Array{Float64,1},
                isf::Array{Float64,1},
                isf_error::Array{Float64,1},
                frequency::Array{Float64,1};
                smooth::Bool = false,
                normalize::Bool = false,
                use_inverse_first_moment::Bool = false,
                first_moment::Float64 = -1.0,
                third_moment::Float64 = -1.0,
                third_moment_error::Float64 = 0.01,
                temperature::Float64 = 1.2,
                number_of_generations::Int64 = 1000,
                population_size::Int64 = 512,
                crossoverProb::Float64=0.9,
                SAcrossoverProb::Float64=0.1,
                differentialWeight::Float64=0.9,
                SAdifferentialWeightProb::Float64=0.1,
                SAdifferentialWeightShift::Float64=0.1,
                SAdifferentialWeight::Float64=0.9,
                rejectFunc::String="smallerFit",
                stop_minimum_fitness::Float64 = 1.0e-8,
                seed::Int64 = 1,
                track_stats::Bool = false,
                number_of_blas_threads::Int64 = 0)

    #_bts = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
    #println("Number of BLAS threads: $(_bts)");
    if number_of_blas_threads > 0
        LinearAlgebra.BLAS.set_num_threads(number_of_blas_threads)
    end

    reject = getfield(DEACreject, Symbol(rejectFunc));

    #Set uuid
    u4 = uuid4();

    #Set RNG
    rng = MersenneTwister(seed);

    #Set constants
    ntau = size(imaginary_time,1);
    ndsf = size(frequency,1);
    npop = population_size;
    isf_r = 1.0 ./ isf;
    beta = 1/temperature;

    moment0=isf[1];
    isf_term = set_isf_term(imaginary_time,frequency,beta);
    isf_term2 = copy(isf_term);

    df = frequency[2:end] .- frequency[1:size(frequency,1) - 1];
    dfrequency = zeros(size(frequency,1));
    dfrequency2 = zeros(size(frequency,1));
    for i in 1:(size(frequency,1) - 1)
        dfrequency[i] = df[i]/2
        dfrequency2[i+1] = df[i]/2
    end
    dfrequency3 = dfrequency .+ dfrequency2;

    isf_term .*= dfrequency';
    isf_term2 .*= dfrequency2';
    isf_term .+= isf_term2;

    isf_m = Array{Float64}(undef,ntau,npop);

    #Generate population and set initial fitness

    P = Array{Float64,2}(undef,ndsf,npop);
    P = rand!(rng,P);

    # Normalize population
    if normalize
        normalization_factor = cosh.((beta/2) .* frequency);
        normalization_term = normalization_factor .* dfrequency3;
        normalization = Array{Float64,1}(undef,npop);
        mul!(normalization,P',normalization_term);
        normalization ./= moment0;
        P ./= normalization';
    end
    P_new = Array{Float64,2}(undef,ndsf,npop);
    fitness = zeros(npop);
    fitness_new = Array{Float64,1}(undef,npop);

    if first_moment > 0.0
        first_moments_factor = frequency .* sinh.((beta/2) .* frequency);
        first_moments_term = first_moments_factor .* dfrequency;
        first_moments_term2 = first_moments_factor .* dfrequency2;
        first_moments_term .+= first_moments_term2;

        first_moments = Array{Float64,1}(undef,npop);

        mul!(first_moments', first_moments_term', P);
        #FIXME see if faster
        #mul!(first_moments, P', first_moments_term);
    end

    if third_moment > 0.0
        third_moments_factor = (frequency .^ 3) .* sinh.((beta/2) .* frequency);
        third_moments_term = third_moments_factor .* dfrequency;
        third_moments_term2 = third_moments_factor .* dfrequency2;
        third_moments_term .+= third_moments_term2;

        third_moments = Array{Float64,1}(undef,npop);

        mul!(third_moments', third_moments_term', P);
        #FIXME see if faster
        #mul!(third_moments, P', third_moments_term);
    end

    #set isf_m and calculate fitness
    mul!(isf_m,isf_term,P);

    if use_inverse_first_moment
        dt = imaginary_time[2:end] .- imaginary_time[1:size(imaginary_time,1) - 1];
        dimaginary_time = zeros(size(imaginary_time,1));
        dimaginary_time2 = zeros(size(imaginary_time,1));
        for i in 1:(size(imaginary_time,1) - 1)
            dimaginary_time[i] = dt[i]/2
            dimaginary_time2[i+1] = dt[i]/2
        end
        dimaginary_time .+= dimaginary_time2;
        inverse_first_moment = dot(isf,dimaginary_time);
        inverse_first_moment_error = sqrt(dot(isf_error.^2,dimaginary_time.^2));
        inverse_first_moments = Array{Float64,1}(undef,npop);
        mul!(inverse_first_moments,isf_m',dimaginary_time);
    end

    broadcast!((x,y,z)->(((x-y)/z)^2),isf_m,isf,isf_m,isf_error);
    mean!(fitness',isf_m);

    if use_inverse_first_moment
        broadcast!((x,y,z)->(((x-y)/z)^2),inverse_first_moments,inverse_first_moment,inverse_first_moments,inverse_first_moment_error);
        fitness .+= inverse_first_moments;
    end

    if first_moment > 0.0
        fitness .+= (first_moments .- first_moment).^2;
    end
    if third_moment > 0.0
        broadcast!((x,y,z)->(((x-y)/z)^2),third_moments,third_moment,third_moments,third_moment_error);
        fitness .+= third_moments;
    end

    if smooth
        total_beam_distances = Array{Float64,1}(undef,npop);
        beam_distance(total_beam_distances,frequency,P);
        fitness .*= total_beam_distances.^2;
    end

    crossover_probs = ones(npop) .* crossoverProb;
    differential_weights = ones(npop) .* differentialWeight;
    new_crossover_probs = Array{Float64,1}(undef,npop);
    new_differential_weights = Array{Float64,1}(undef,npop);

    #Initialize statistics arrays

    if track_stats
        avgFitness = zeros(number_of_generations);
        minFitness = zeros(number_of_generations);
        stdFitness = zeros(number_of_generations);
    end

    #Crude Preallocation
    mIdx = trues(ndsf,npop);
    generation = 1;

    rInd = falses(npop); #rejection indexes

    minP = P[:,argmin(fitness)];
    minFit = minimum(fitness);
    avgFit = mean(fitness);
    stdFit = std(fitness,mean=avgFit);

    #Begin loop
    @inbounds for outer generation = 1:number_of_generations
        avgFit = mean(fitness);
        minFit = minimum(fitness);
        stdFit = std(fitness,mean=avgFit);

        #Get Statistics
        if track_stats
            avgFitness[generation] = avgFit;
            minFitness[generation] = minFit;
            stdFitness[generation] = stdFit;
        end

        #Stopping criteria
        if minFit <= stop_minimum_fitness
            break
        end

        set_crossover_probs(rng,new_crossover_probs,crossover_probs,SAcrossoverProb);
        set_differential_weights(rng,new_differential_weights,differential_weights,SAdifferentialWeightProb,
                                 SAdifferentialWeightShift,SAdifferentialWeight);

        #Set mutant population FIXME crossover built-in here do I want this?
        set_mutateInd(rng,mIdx,new_crossover_probs);
        mutate(rng,P_new,P,mIdx,new_differential_weights);

        # Normalization
        if normalize
            mul!(normalization,P_new',normalization_term);
            normalization ./= moment0;
            P_new ./= normalization';
        end

        #Rejection
        if first_moment > 0.0
            mul!(first_moments', first_moments_term', P_new);
        end

        if third_moment > 0.0
            mul!(third_moments', third_moments_term', P_new);
        end

        mul!(isf_m,isf_term,P_new);

        if use_inverse_first_moment
            mul!(inverse_first_moments,isf_m',dimaginary_time);
        end

        broadcast!((x,y,z)->(((x-y)/z)^2),isf_m,isf,isf_m,isf_error);
        mean!(fitness_new',isf_m);

        if use_inverse_first_moment
            broadcast!((x,y,z)->(((x-y)/z)^2),inverse_first_moments,inverse_first_moment,inverse_first_moments,inverse_first_moment_error);
            fitness_new .+= inverse_first_moments;
        end

        if first_moment > 0.0
            fitness_new .+= (first_moments .- first_moment).^2;
        end

        if third_moment > 0.0
            broadcast!((x,y,z)->(((x-y)/z)^2),third_moments,third_moment,third_moments,third_moment_error);
            fitness_new .+= third_moments;
        end

        if smooth
            beam_distance(total_beam_distances,frequency,P_new);
            fitness_new .*= total_beam_distances.^2;
        end

        reject(rng,rInd,fitness_new,fitness);
        replace_pop(P,P_new,rInd);
        replace_fitness(fitness,fitness_new,rInd);
        replace_crossover_probs(crossover_probs,new_crossover_probs,rInd);
        replace_differential_weights(differential_weights,new_differential_weights,rInd);

    end
    #Set best candidate solution from final generation
    minidx = argmin(fitness);
    set_minP(minP,P,minidx);

    minS = minP .* exp.((beta/2) .* frequency) ./ 2;
    println(minimum(fitness));
    if track_stats
        return u4,frequency,minS,minP,minimum(fitness),generation,avgFitness[1:generation],minFitness[1:generation],stdFitness[1:generation],P,fitness,rng,crossover_probs,differential_weights
    end
    return u4,frequency,minS,minP,minimum(fitness),generation,zeros(2),zeros(2),zeros(2),P,fitness,rng,crossover_probs,differential_weights
end
end


text:
# Needleman-Wunsch Algorithm
# ==========================
#
# The Needleman-Wunsch algorithm.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioAlignments.jl/blob/master/LICENSE.md

mutable struct NeedlemanWunsch{T<:Union{Signed,AbstractFloat}}
    trace::Matrix{Trace}
    H::Vector{T}
    E::Vector{T}

    function NeedlemanWunsch{T}(m::Integer, n::Integer) where T
        trace = Matrix{Trace}(undef, (m + 1, n + 1))
        fill!(trace, 0xff)
        H = Vector{T}(undef, m + 1)
        E = Vector{T}(undef, m)
        return new{T}(trace, H, E)
    end

    function NeedlemanWunsch{T}() where T
        return NeedlemanWunsch{T}(0, 0)
    end
end

function ensureroom!(nw::NeedlemanWunsch, m, n)
    if size(nw.trace, 1) < m + 1 || size(nw.trace, 2) < n + 1
        # TODO: resize the trace matrix if possible
        nw.trace = Matrix{Trace}(undef, (m + 1, n + 1))
        nw.H = resize!(nw.H, m + 1)
        nw.E = resize!(nw.E, m)
    end
    return nw
end

# Generate a code block that updates the current score and trace.
macro update()
    esc(quote
        e = E[i]
        g = h_diag + submat[a[i],b_j]
        h = max(e, f, g)
        h_diag = H[i+1]
        H[i+1] = h
        t = trace[i+1,j+1] | ft
        e == h && (t |= TRACE_DELETE)
        f == h && (t |= TRACE_INSERT)
        g == h && (t |= TRACE_MATCH)
        trace[i+1,j+1] = t
    end)
end

macro update_nextcol(gap_init, gap_extend)
    esc(quote
        e′ = e + $gap_extend
        e = max(e′, h + $gap_init)
        E[i] = e
        et = TRACE_NONE
        e == e′ && (et |= TRACE_EXTDEL)
        trace[i+1,j+2] = et
    end)
end

macro update_nextrow(gap_init, gap_extend)
    esc(quote
        f′ = f + $gap_extend
        f = max(f′, h + $gap_init)
        ft = TRACE_NONE
        f == f′ && (ft |= TRACE_EXTINS)
    end)
end

#    start   middle  end gaps of a
#    |       |       |
# a: ----ACCA---ATGTG---
# b: AAATACGATTGATGAGGGT

function run!(
        nw::NeedlemanWunsch{T},
        a, b,
        submat::AbstractSubstitutionMatrix{T},
        # a
        start_gap_open_a::T,
        start_gap_extend_a::T,
        middle_gap_open_a::T,
        middle_gap_extend_a::T,
        end_gap_open_a::T,
        end_gap_extend_a::T,
        # b
        start_gap_open_b::T,
        start_gap_extend_b::T,
        middle_gap_open_b::T,
        middle_gap_extend_b::T,
        end_gap_open_b::T,
        end_gap_extend_b::T) where T

    m = length(a)
    n = length(b)
    ensureroom!(nw, m, n)
    H = nw.H
    E = nw.E
    trace = nw.trace

    start_gap_init_a = start_gap_open_a + start_gap_extend_a
    middle_gap_init_a = middle_gap_open_a + middle_gap_extend_a
    end_gap_init_a = end_gap_open_a + end_gap_extend_a

    start_gap_init_b = start_gap_open_b + start_gap_extend_b
    middle_gap_init_b = middle_gap_open_b + middle_gap_extend_b
    end_gap_init_b = end_gap_open_b + end_gap_extend_b

    # ↓: insertion
    # →: deletion
    # ↘: match/mismatch
    H[1] = T(0)
    trace[1,1] = TRACE_NONE
    if m == 0 && n == 0
        return H[1]
    elseif m == 0
        for j in 1:n
            trace[1,j+1] = TRACE_DELETE
        end
        return start_gap_open_a + start_gap_extend_a * n
    elseif n == 0
        for i in 1:m
            trace[i+1,1] = TRACE_INSERT
        end
        return start_gap_open_b + start_gap_extend_b * m
    else
        # initialize the first column of H, E, and trace
        @inbounds for i in 1:m
            H[i+1] = start_gap_open_b + start_gap_extend_b * i
            E[i] = H[i+1] + (i == m ? end_gap_init_a : middle_gap_init_a)
            trace[i+1,1] = TRACE_INSERT
            trace[i+1,2] = TRACE_NONE
        end

        # run dynamic programming column by column (except the last column)
        @inbounds for j in 1:n-1
            b_j = b[j]
            h_diag = H[1]
            H[1] = start_gap_open_a + start_gap_extend_a * j
            f = H[1] + middle_gap_init_b
            ft = TRACE_NONE
            trace[1,j+1] = TRACE_DELETE
            # inner loop along the sequence a
            for i in 1:m-1
                @update
                @update_nextcol middle_gap_init_a middle_gap_extend_a
                @update_nextrow middle_gap_init_b middle_gap_extend_b
            end
            # fill the final row
            let i = m
                @update
                @update_nextcol end_gap_init_a end_gap_extend_a
            end
        end

        # fill the last column
        let j = n
            b_j = b[j]
            h_diag = H[1]
            H[1] = start_gap_open_a + start_gap_extend_a * j
            f = H[1] + end_gap_init_b
            ft = TRACE_NONE
            trace[1,j+1] = TRACE_DELETE
            # inner loop along the sequence a
            @inbounds for i in 1:m-1
                @update
                @update_nextrow end_gap_init_b end_gap_extend_b
            end
            # finally, fill the last element of the DP matrices
            let i = m
                @update
            end
        end

        return H[m+1]
    end
end

function run!(nw::NeedlemanWunsch, a, b, submat, gap_open, gap_extend)
    return run!(nw, a, b, submat,
        # a
        gap_open,
        gap_extend,
        gap_open,
        gap_extend,
        gap_open,
        gap_extend,
        # b
        gap_open,
        gap_extend,
        gap_open,
        gap_extend,
        gap_open,
        gap_extend
    )
end

function traceback(nw::NeedlemanWunsch, a, b, endpos)
    anchors = Vector{AlignmentAnchor}()
    i, j = endpos
    was_extins = false
    was_extdel = false
    @start_traceback
    while i ≥ 1 || j ≥ 1
        t = nw.trace[i+1,j+1]
        @assert !(was_extins && was_extdel)
        if was_extins
            was_extins = t & TRACE_EXTINS > 0
            @anchor OP_INSERT
        elseif was_extdel
            was_extdel = t & TRACE_EXTDEL > 0
            @anchor OP_DELETE
        elseif t & TRACE_INSERT > 0
            was_extins = t & TRACE_EXTINS > 0
            @anchor OP_INSERT
        elseif t & TRACE_DELETE > 0
            was_extdel = t & TRACE_EXTDEL > 0
            @anchor OP_DELETE
        elseif t & TRACE_MATCH > 0
            if a[i] == b[j]
                @anchor OP_SEQ_MATCH
            else
                @anchor OP_SEQ_MISMATCH
            end
        else
            error("failed to trace back")
        end
    end
    @finish_traceback
    return AlignedSequence(a, anchors)
end


text:
function maph_moments_and_rand_test(;N=1000)
    #QQQQ create an MAPH Distributions
    Λ₄, λ45, λ54, Λ₅ = 7.0, 2.0, 7.0, 14.0
    μ41, μ42, μ43, μ51, μ52, μ53 = 2.0, 1.0, 2.0, 1.0, 4.0, 2.0
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]

    maph = MAPHDist([0.3,0.7]', T_example, T0_example)

    data = [rand(maph) for _ in 1:N]
    π_est = [count((x)->x.a == i, data)/length(data) for i in 1:3]

    μ_est = [mean(first.(filter((x)->x.a == i, data))) for i in 1:3]

    σ2_est = [var(first.(filter((x)->x.a == i, data))) for i in 1:3]

    @show mean(maph), μ_est
    @show var(maph),σ2_est

    scv_est = [scv(first.(filter((x)->x.a == i, data))) for i in 1:3]

    @show scv(maph), scv_est

    # first_moment = mean(maph)
    # variance = var(maph)

    # (norm(first_moment .- μ_est) < 10/sqrt(N)) || (return false)
    # (norm(variance .- σ2_est) < 10/sqrt(N)) || (return false)
    # (norm(scv(maph) .- scv_est) < 20/sqrt(N)) || (return false)
    # (norm(absorption_probs(maph) .- π_est) < 5/sqrt(N)) || (return false)

    return true
end


text:
using MomentClosure, ModelingToolkit
using Documenter

makedocs(;
    modules=[MomentClosure, ModelingToolkit, Symbolics],
    authors="Augustinas Sukys",
    repo="https://github.com/augustinas1/MomentClosure.jl/blob/{commit}{path}#L{line}",
    sitename="MomentClosure.jl",
    doctest=false,
    format=Documenter.HTML(;
        mathengine=Documenter.Writers.HTMLWriter.MathJax2(),
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://augustinas1.github.io/MomentClosure.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Theory" => Any[
            "theory/moment_expansion.md",
            "theory/moment_closure_approximations.md"
        ],
        "Tutorials" => Any[
            "tutorials/using_momentclosure.md",
            "tutorials/common_issues.md",
            "tutorials/geometric_reactions+conditional_closures.md",
            "tutorials/P53_system_example.md",
            "tutorials/derivative_matching_example.md"
        ],
        "API" => Any[
            "api/momentclosure_api.md"
        ]
    ],

)

deploydocs(;
    repo="github.com/augustinas1/MomentClosure.jl",
    devbranch="main"
)


text:
function test_iterator_takewhile()
    @test collect(takewhile(x -> x^2 < 1, 1:10)) == []
    @test collect(takewhile(x -> x^2 < 25, 1:10)) == [1, 2, 3, 4]
    @test collect(takewhile(x -> x^2 < 1000, 1:10)) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    @test eltype(takewhile(_ -> true, 1:10)) == typeof(1)

    @test collect(takewhile(1:10) do x; x^2 < 25; end) == [1, 2, 3, 4]
end

function test_iterator_dropwhile()
    @test collect(dropwhile(x -> x^2 < 1, 1:10)) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    @test collect(dropwhile(x -> x^2 < 25, 1:10)) == [5, 6, 7, 8, 9, 10]
    @test collect(dropwhile(x -> x^2 < 1000, 1:10)) == []
    @test eltype(dropwhile(_ -> true, 1:10)) == typeof(1)

    @test collect(dropwhile(1:10) do x; x^2 < 25; end) == [5, 6, 7, 8, 9, 10]
end

function test_iterator_tmap()
    @test collect(tmap(x -> x + 1, Int, 1:10)) == [2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
    @test collect(tmap(x -> x + 1, BigInt, 1:10)) == [2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
    @test eltype(tmap(x -> x + 1, Int, 1:10)) == Int
    @test eltype(tmap(x -> x + 1, BigInt, 1:10)) == BigInt
    @test length(tmap(x -> x + 1, Int, 1:10)) == 10

    @test collect(tmap(+, Int, 1:5, 11:20)) == [12, 14, 16, 18, 20]
    @test length(tmap(+, Int, 1:5, 11:20)) == 5
end

create(it::DefaultNestedIterator{Array{Int,1}}, value) = begin
    f = filter(x -> x != value, it.source)
    length(f) > 1 ? nested(f, it.level + 1) : f
end
combine(::DefaultNestedIterator{Array{Int,1}}, outer, inner) = "$outer$inner"

Base.eltype(it::DefaultNestedIterator{Array{Int,1}}) = eltype(typeof(it))
Base.eltype(::Type{DefaultNestedIterator{Array{Int,1}}}) = String

function test_iterator_nested_iterator()
    @test collect(nested([big(1), big(2), big(3)])) == []

    @test collect(nested(collect(1:3))) == ["123","132","213","231","312","321"]
    @test eltype(nested(collect(1:3))) == String
end

struct DigitPermutationIterator <: NestedIterator
    source::Array{Int,1}
    maxlen::Int
end

digitperm(s, l) = l == 1 ? s : DigitPermutationIterator(s, l)

create(it::DigitPermutationIterator, value) = begin
    f = filter(x -> x != value, it.source)
    (it.maxlen > 2 && length(f) > 1) ? digitperm(f, it.maxlen - 1) : f
end
combine(::DigitPermutationIterator, outer, inner) = "$outer$inner"

Base.eltype(it::DigitPermutationIterator) = eltype(typeof(it))
Base.eltype(::Type{DigitPermutationIterator}) = String

function test_iterator_nested_iterator2()
    @test collect(digitperm(collect(1:4), 2)) == ["12","13","14","21","23","24","31","32","34","41","42","43"]
    @test eltype(digitperm(collect(1:4), 2)) == String
end

function test_iterator_all()
    print(rpad("DataStructure.Iterator...", 50, ' '))

    test_iterator_takewhile()
    test_iterator_dropwhile()
    test_iterator_tmap()
    test_iterator_nested_iterator()
    test_iterator_nested_iterator2()

    println("PASS")
end


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

### Gauss quadrature rules for triangular elements

""" Gauss-Legendre quadrature, 1 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI1}})
    weights = (0.5, )
    points = ((1.0/3.0, 1.0/3.0), )
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI1}})
    return 1
end

""" Gauss-Legendre quadrature, 3 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI3}})
    weights = (1.0/6.0, 1.0/6.0, 1.0/6.0)
    points = (
              (2.0/3.0, 1.0/6.0),
              (1.0/6.0, 2.0/3.0),
              (1.0/6.0, 1.0/6.0)
             )
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI3}})
    return 2
end

""" Gauss-Legendre quadrature, 3 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI3B}})
    weights = (1.0/6.0, 1.0/6.0, 1.0/6.0)
    points = (
              (0.0, 1.0/2.0),
              (1.0/2.0, 0.0),
              (1.0/2.0, 1.0/2.0)
             )
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI3B}})
    return 2
end

""" Gauss-Legendre quadrature, 4 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI4}})
    weights = (
        1.5902069087198858469718450103758e-01,
        9.0979309128011415302815498962418e-02,
        1.5902069087198858469718450103758e-01,
        9.0979309128011415302815498962418e-02)
    points = (
        (1.5505102572168219018027159252941e-01, 1.7855872826361642311703513337422e-01),
        (6.4494897427831780981972840747059e-01, 7.5031110222608118177475598324603e-02),
        (1.5505102572168219018027159252941e-01, 6.6639024601470138670269327409637e-01),
        (6.4494897427831780981972840747059e-01, 2.8001991549907407200279599420481e-01))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI4}})
    return 3
end

""" Gauss-Legendre quadrature, 4 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI4B}})
    weights = (-27.0/96.0, 25.0/96.0, 25.0/96.0, 25.0/96.0)
    points = (
              (1.0/3.0, 1.0/3.0),
              (1.0/5.0, 1.0/5.0),
              (1.0/5.0, 3.0/5.0),
              (3.0/5.0, 1.0/5.0))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI4B}})
    return 3
end

""" Gauss-Legendre quadrature, 6 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI6}})
    P1 = 0.11169079483905
    P2 = 0.0549758718227661
    A = 0.445948490915965
    B = 0.091576213509771
    weights = (P2, P2, P2, P1, P1, P1)
    points = ((B, B), (1.0-2.0*B, B), (B, 1.0-2.0*B),
              (A, 1.0-2*A), (A, A), (1.0 - 2.0*A, A))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI6}})
    return 4
end

""" Gauss-Legendre quadrature, 7 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI7}})
    A = 0.470142064105115
    B = 0.101286507323456
    P1 = 0.066197076394253
    P2 = 0.062969590272413
    weights = (9/80, P1, P1, P1, P2, P2, P2)
    points = ((1/3, 1/3), (A, A), (1.0-2.0*A, A), (A, 1.0-2.0*A),
              (B, B), (1.0-2.0*B, B), (B, 1.0-2.0*B))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI7}})
    return 5
end

""" Gauss-Legendre quadrature, 12 point rule on triangle. """
function get_quadrature_points(::Type{Val{:GLTRI12}})
    A = 0.063089014491502
    B = 0.249286745170910
    C = 0.310352451033785
    D = 0.053145049844816
    P1 = 0.025422453185103
    P2 = 0.058393137863189
    P3 = 0.041425537809187
    weights = (P1, P1, P1, P2, P2, P2, P3, P3, P3, P3, P3, P3)
    points = ((A, A), (1.0-2.0*A, A), (A, 1.0-2.0*A), (B, B), (1.0-2.0*B, B),
              (B, 1.0-2.0*B), (C, D), (D, C), (1.0-C-D, C), (1.0-C-D, D),
              (C, 1.0-C-D), (D, 1.0-C-D))
    return zip(weights, points)
end

function get_order(::Type{Val{:GLTRI12}})
    return 6
end


text:
# Test multivariate optimization
let
    rosenbrock = Optim.UnconstrainedProblems.examples["Rosenbrock"]
    f = rosenbrock.f
    g! = rosenbrock.g!
    h! = rosenbrock.h!
    initial_x = rosenbrock.initial_x

    d1 = DifferentiableFunction(f)
    d2 = DifferentiableFunction(f, g!)
    d3 = TwiceDifferentiableFunction(f, g!, h!)

    Optim.optimize(f, initial_x, BFGS())
    Optim.optimize(f, g!, initial_x, BFGS())
    Optim.optimize(f, g!, h!, initial_x, BFGS())
    Optim.optimize(d2, initial_x, BFGS())
    Optim.optimize(d3, initial_x, BFGS())

    Optim.optimize(f, initial_x, BFGS(), Optim.Options())
    Optim.optimize(f, g!, initial_x, BFGS(), Optim.Options())
    Optim.optimize(f, g!, h!, initial_x, BFGS(), Optim.Options())
    Optim.optimize(d2, initial_x, BFGS(), Optim.Options())
    Optim.optimize(d3, initial_x, BFGS(), Optim.Options())

    Optim.optimize(d1, initial_x, BFGS())
    Optim.optimize(d2, initial_x, BFGS())

    Optim.optimize(d1, initial_x, GradientDescent())
    Optim.optimize(d2, initial_x, GradientDescent())

    Optim.optimize(d1, initial_x, LBFGS())
    Optim.optimize(d2, initial_x, LBFGS())

    Optim.optimize(f, initial_x, NelderMead())

    Optim.optimize(d3, initial_x, Newton())

    Optim.optimize(f, initial_x, SimulatedAnnealing())

    optimize(f, initial_x, BFGS())
    optimize(f, g!, initial_x, BFGS())
    optimize(f, g!, h!, initial_x, BFGS())

    optimize(f, initial_x, GradientDescent())
    optimize(f, g!, initial_x, GradientDescent())
    optimize(f, g!, h!, initial_x, GradientDescent())

    optimize(f, initial_x, LBFGS())
    optimize(f, g!, initial_x, LBFGS())
    optimize(f, g!, h!, initial_x, LBFGS())

    optimize(f, initial_x, NelderMead())
    optimize(f, g!, initial_x, NelderMead())
    optimize(f, g!, h!, initial_x, NelderMead())

    optimize(f, g!, h!, initial_x, Newton())

    optimize(f, initial_x, SimulatedAnnealing())
    optimize(f, g!, initial_x, SimulatedAnnealing())
    optimize(f, g!, h!, initial_x, SimulatedAnnealing())

    options = Optim.Options(g_tol = 1e-12, iterations = 10,
                                  store_trace = true, show_trace = false)
    res = optimize(f, g!, h!,
    	           initial_x,
    	           BFGS(),
    	           options)

    options_g = Optim.Options(g_tol = 1e-12, iterations = 10,
                                  store_trace = true, show_trace = false)
    options_f = Optim.Options(g_tol = 1e-12, iterations = 10,
                                  store_trace = true, show_trace = false)
    res = optimize(f, g!, h!,
    	           initial_x,
    	           GradientDescent(),
    	           options_g)

    res = optimize(f, g!, h!,
    	           initial_x,
    	           LBFGS(),
    	           options_g)

    res = optimize(f, g!, h!,
    	           initial_x,
    	           NelderMead(),
                   options_f)

    res = optimize(f, g!, h!,
    	           initial_x,
    	           Newton(),
                   options_g)
    options_sa = Optim.Options(iterations = 10, store_trace = true,
                                     show_trace = false)
    res = optimize(f, g!, h!,
    	           initial_x,
    	           SimulatedAnnealing(),
    	           options_sa)

    res = optimize(f, g!, h!,
    	           initial_x,
    	           BFGS(),
    	           options_g)
    options_ext = Optim.Options(g_tol = 1e-12, iterations = 10,
                                      store_trace = true, show_trace = false,
                                      extended_trace = true)
    res_ext = optimize(f, g!, h!,
                       initial_x,
                       BFGS(),
                       options_ext)

   @test Optim.method(res) == "BFGS"
   @test Optim.minimum(res) ≈ 0.055119582904897345
   @test Optim.minimizer(res) ≈ [0.7731690866149542; 0.5917345966396391]
   @test Optim.iterations(res) == 10
   @test Optim.f_calls(res) == 48
   @test Optim.g_calls(res) == 48
   @test Optim.converged(res) == false
   @test Optim.x_converged(res) == false
   @test Optim.f_converged(res) == false
   @test Optim.g_converged(res) == false
   @test Optim.x_tol(res) == 1e-32
   @test Optim.f_tol(res) == 1e-32
   @test Optim.g_tol(res) == 1e-12
   @test Optim.iteration_limit_reached(res) == true
   @test Optim.initial_state(res) == [0.0; 0.0]
   @test haskey(Optim.trace(res_ext)[1].metadata,"x")

   # just testing if it runs
   Optim.trace(res)
   Optim.f_trace(res)
   Optim.g_norm_trace(res)
   @test_throws ErrorException Optim.x_trace(res)
   @test_throws ErrorException Optim.x_lower_trace(res)
   @test_throws ErrorException Optim.x_upper_trace(res)
   @test_throws ErrorException Optim.lower_bound(res)
   @test_throws ErrorException Optim.upper_bound(res)
   @test_throws ErrorException Optim.rel_tol(res)
   @test_throws ErrorException Optim.abs_tol(res)
   options_extended = Optim.Options(store_trace = true, extended_trace = true)
   res_extended = Optim.optimize(f, g!, initial_x, BFGS(), options_extended)
   @test haskey(Optim.trace(res_extended)[1].metadata,"~inv(H)")
   @test haskey(Optim.trace(res_extended)[1].metadata,"g(x)")
   @test haskey(Optim.trace(res_extended)[1].metadata,"x")
   options_extended_nm = Optim.Options(store_trace = true, extended_trace = true)
   res_extended_nm = Optim.optimize(f, g!, initial_x, NelderMead(), options_extended_nm)
   @test haskey(Optim.trace(res_extended_nm)[1].metadata,"centroid")
   @test haskey(Optim.trace(res_extended_nm)[1].metadata,"step_type")
end

# Test univariate API
let
    f(x) = 2x^2+3x+1
    res = optimize(f, -2.0, 1.0, GoldenSection())
    @test Optim.method(res) == "Golden Section Search"
    @test Optim.minimum(res) ≈ -0.125
    @test Optim.minimizer(res) ≈ -0.749999994377939
    @test Optim.iterations(res) == 38
    @test Optim.iteration_limit_reached(res) == false
    @test_throws ErrorException Optim.trace(res)
    @test_throws ErrorException Optim.x_trace(res)
    @test_throws ErrorException Optim.x_lower_trace(res)
    @test_throws ErrorException Optim.x_upper_trace(res)
    @test_throws ErrorException Optim.f_trace(res)
    @test Optim.lower_bound(res) == -2.0
    @test Optim.upper_bound(res) == 1.0
    @test Optim.rel_tol(res) ≈ 1.4901161193847656e-8
    @test Optim.abs_tol(res) ≈ 2.220446049250313e-16
    @test_throws ErrorException Optim.initial_state(res)
    @test_throws ErrorException Optim.g_norm_trace(res)
    @test_throws ErrorException Optim.g_calls(res)
    @test_throws ErrorException Optim.x_converged(res)
    @test_throws ErrorException Optim.f_converged(res)
    @test_throws ErrorException Optim.g_converged(res)
    @test_throws ErrorException Optim.x_tol(res)
    @test_throws ErrorException Optim.f_tol(res)
    @test_throws ErrorException Optim.g_tol(res)
    options =
    res = optimize(f, -2.0, 1.0, GoldenSection(), store_trace = true, extended_trace = true)

    # Right now, these just "test" if they run
    Optim.x_trace(res)
    Optim.x_lower_trace(res)
    Optim.x_upper_trace(res)
end


text:
# LIMIT clause.

mutable struct LimitClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}
    offset::Union{Int, Nothing}
    limit::Union{Int, Nothing}
    with_ties::Bool

    LimitClause(;
                over = nothing,
                offset = nothing,
                limit = nothing,
                with_ties = false) =
        new(over, offset, limit, with_ties)
end

LimitClause(limit; over = nothing, offset = nothing, with_ties = false) =
    LimitClause(over = over, offset = offset, limit = limit, with_ties = with_ties)

LimitClause(offset, limit; over = nothing, with_ties = false) =
    LimitClause(over = over, offset = offset, limit = limit, with_ties = with_ties)

LimitClause(range::UnitRange, over = nothing, with_ties = false) =
    LimitClause(over = over, offset = first(range) - 1, limit = length(range), with_ties = with_ties)

"""
    LIMIT(; over = nothing, offset = nothing, limit = nothing, with_ties = false)
    LIMIT(limit; over = nothing, offset = nothing, with_ties = false)
    LIMIT(offset, limit; over = nothing, with_ties = false)
    LIMIT(start:stop; over = nothing, with_ties = false)

A `LIMIT` clause.

# Examples

```jldoctest
julia> c = FROM(:person) |>
           LIMIT(1) |>
           SELECT(:person_id);

julia> print(render(c))
SELECT "person_id"
FROM "person"
FETCH FIRST 1 ROW ONLY
```
"""
LIMIT(args...; kws...) =
    LimitClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(LIMIT), pats::Vector{Any}) =
    dissect(scr, LimitClause, pats)

function PrettyPrinting.quoteof(c::LimitClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(LIMIT))
    if c.offset !== nothing
        push!(ex.args, c.offset)
    end
    push!(ex.args, c.limit)
    if c.with_ties
        push!(ex.args, Expr(:kw, :with_ties, c.with_ties))
    end
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::LimitClause, c′) =
    LimitClause(over = rebase(c.over, c′), offset = c.offset, limit = c.limit, with_ties = c.with_ties)