text:
@testset "FBM" begin
    rng = MersenneTwister(42)
    h = 0.3
    k = FBMKernel(h = h)
    v1 = rand(rng, 3); v2 = rand(rng, 3)
    @test k(v1,v2) ≈ (sqeuclidean(v1, zero(v1))^h + sqeuclidean(v2, zero(v2))^h - sqeuclidean(v1-v2, zero(v1-v2))^h)/2 atol=1e-5

    # kernelmatrix tests
    m1 = rand(rng, 3, 3)
    m2 = rand(rng, 3, 3)
    Kref = kernelmatrix(k, m1, m1)
    @test kernelmatrix(k, m1) ≈ Kref atol=1e-5
    K = zeros(3, 3)
    kernelmatrix!(K, k, m1, m1)
    @test K ≈ Kref atol=1e-5
    fill!(K, 0)
    kernelmatrix!(K, k, m1)
    @test K ≈ Kref atol=1e-5

    x1 = rand(rng)
    x2 = rand(rng)
    @test kernelmatrix(k, x1*ones(1,1), x2*ones(1,1))[1] ≈ k(x1, x2) atol=1e-5

    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"
    test_ADs(FBMKernel, ADs = [:ReverseDiff])
    @test_broken "Tests failing for kernelmatrix(k, x) for ForwardDiff and Zygote"
end


text:
abstract type AbstractStandard{T} end

function Simplexe(A::Array{T, 2}, b::Array{S, 1}, c::Array{W, 1}; method::Type{METHODE} = StandardSimplexe,
        verbose::Bool = false) where {METHODE <: AbstractStandard, T, S, W}
    nT = promote_type(T, S, W)
    !(supertype(nT) <: Real) && (nT = Float64)
    return method(Array{nT, 2}(A), Array{nT, 1}(b), Array{nT, 1}(c), verbose = verbose)
end

include("BaseEnumeration.jl")
include("Simplexe.jl")

ext: jl


text:
function add_scal(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  A = randn(n,n)
  B = randn(n,n)

  Benchmarker.add_data(csv, "add", Benchmarker.measure(reps, add, A, B))
  Benchmarker.add_data(csv, "scal", Benchmarker.measure(reps, scal, A))

end

function add(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  K = copy(A)
  X = copy(B)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  K = X + K

  finish = time_ns()
  return ((finish-start)*1e-9, K)
end

function scal(A::DenseArray{Float64, 2})
  K = copy(A)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  K = 3.0 * K

  finish = time_ns()
  return ((finish-start)*1e-9, K)
end


text:
using FANN
using Base.Test

# desired input
X = [ 1.0 3.0 5.0 7.0;
	  2.0 4.0 6.0 8.0]

# desired ouputs
Y = [0.0 1.0 1.0 0.0]

# definitions
M = size(X, 2)
nin = size(X, 1)
nout = size(Y, 1)

# create dataset
dset = DataSet(X, Y)

# 1) try to load input data from the struct
inputs = pointer_to_array(unsafe_load(dset.data).input, M)
outputs = pointer_to_array(unsafe_load(dset.data).output, M)
for i = 1:M
	@test vec(X[:, i]) == pointer_to_array(inputs[i], nin)
	@test vec(Y[:, i]) == pointer_to_array(outputs[i], nout)
end

# 2)  try to save it and the read file
savedset(dset, "here")
data = readdlm("here")
rm("here")

# first row is header data
header = data[1, :]
@test header[1] == M
@test header[2] == nin
@test header[3] == nout

# then read rest
for i = 1:M
	actual = Float64[data[2*i, 1:nin]...] # very ugly
	expected = X[:, i]
	@test norm(actual - expected) < 1e-5

	actual = Float64[data[2*i + 1, 1:nout]...] # very ugly
	expected = Y[:, i]
	@test norm(actual - expected) < 1e-5
end


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

using Markdown
using InteractiveUtils

# ╔═╡ cd0b0c44-4457-4ab6-ac00-3519ccdf1388
import Pkg; Pkg.add("Plots")

# ╔═╡ b7a5404c-9570-11eb-0562-47e5f0701cd0
using MolecularIntegrals, Profile, ProfileSVG

# ╔═╡ 39a6c055-e9e6-4736-afe5-e46d51c97139
using MolecularIntegrals:Fgamma

# ╔═╡ 1eb1804e-2894-40b7-86dc-98ec4db2b3d7
using Plots

# ╔═╡ f49f47dd-cad0-49e8-9dd4-04ae57a8fb28
md"""
# MolecularIntegrals.jl

## TODO/Tasks:
- Profile code
- Benchmark timings
- Replace as many hand-written routines with library functions
  - dist2 with hypot^2?
  - [X] gammainc, gser, gcf in nuclear attraction ints with Special Functions gamma_inc
  - Similar replacement of beta functions?
- Treat constants more intelligently, and make them `const`
- Move ProfileSVG, Test, Plots dependencies to tools and tests subdirectories
"""

# ╔═╡ b07aa4a1-4372-4bc7-9432-2c8aa29b6485
md"## Fgamma investigations"

# ╔═╡ 20dede5e-5b86-4176-8f4a-52609f73971b
begin
	xs = 0:0.01:2
	plot(xs,Fgamma.(1,xs),label="m=0")
	plot!(xs,Fgamma.(2,xs),label="m=1")
	plot!(xs,Fgamma.(3,xs),label="m=2")
end

# ╔═╡ 6f3be0b3-856d-45d8-a870-187f09587195
md"## Profiling and Timing"

# ╔═╡ 3b737810-5ac7-4c9e-9307-6e04763ca3d1
bfs = build_basis(ethane,"6-31G")

# ╔═╡ 11bb6095-c88e-41b1-a569-bc705005db32
length(bfs)

# ╔═╡ 5d3b0917-97d6-4588-8d0c-921288051428
MolecularIntegrals.all_twoe_ints(bfs)

# ╔═╡ aa89c6af-49d9-41db-94a6-3c87a92a4298
md"The Huzinaga and HGP methods showed 7.1 and 44.1 seconds, respectively. The HGP time is very surprising. This was once relatively fast. I guess it shouldn't be too hard to find efficiencies!"

# ╔═╡ 135dae13-7bba-4936-9006-942d2798e158
@profview MolecularIntegrals.all_twoe_ints(bfs)

# ╔═╡ Cell order:
# ╠═b7a5404c-9570-11eb-0562-47e5f0701cd0
# ╟─f49f47dd-cad0-49e8-9dd4-04ae57a8fb28
# ╠═b07aa4a1-4372-4bc7-9432-2c8aa29b6485
# ╠═39a6c055-e9e6-4736-afe5-e46d51c97139
# ╠═cd0b0c44-4457-4ab6-ac00-3519ccdf1388
# ╠═1eb1804e-2894-40b7-86dc-98ec4db2b3d7
# ╠═20dede5e-5b86-4176-8f4a-52609f73971b
# ╟─6f3be0b3-856d-45d8-a870-187f09587195
# ╠═3b737810-5ac7-4c9e-9307-6e04763ca3d1
# ╠═11bb6095-c88e-41b1-a569-bc705005db32
# ╠═5d3b0917-97d6-4588-8d0c-921288051428
# ╠═aa89c6af-49d9-41db-94a6-3c87a92a4298
# ╠═135dae13-7bba-4936-9006-942d2798e158


text:
type Account
    id::Int64
    acctype::Symbol
    accname::Compat.UTF8String
end

id(acc::Account) = acc.id
id!(acc::Account, o::Int64) = (acc.id = o)

name(acc::Account) = acc.accname
accounttype(acc::Account) = acc.acctype


text:
function planeFit(X)
    n = size(X)[1]
    X′ = zeros(n, 3)
    Xmean = mean(X, 1)
    for i in 1:n
        X′[i, :] = X[i, :] - Xmean
    end
    A = zeros(3, 3)
    B = zeros(3)
    for i in 1:n
        A[1, 1] += X′[i, 1] ^ 2
        A[1, 2] += X′[i, 1] * X′[i, 2]
        A[1, 3] += X′[i, 1]
        A[2, 2] += X′[i, 2] ^ 2
        A[2, 3] += X′[i, 2]
        A[3, 3] += 1.0
        B[1] += X′[i, 1] * X′[i, 3]
        B[2] += X′[i, 2] * X′[i, 3]
        B[3] += X′[i, 3]
    end
    A[2, 1] = A[1, 2]
    A[3, 1] = A[1, 3]
    A[3, 2] = A[2, 3]

    Y = \(A, B)
    Y[3] = Y[3] - Y[1] * Xmean[1] - Y[2] * Xmean[2] + Xmean[3]
    return Y, Xmean
end

function sphereFit(X, iternum)
    n = size(X)[1]
    Xmean = mean(X, 1)

    function update(a, b, c, Xmean)
        L̄ = 0.0
        L̄a = 0.0
        L̄b = 0.0
        L̄c = 0.0
        function L(x, y, z, a, b, c)
            return sqrt((x - a)^2 + (y - b)^2 + (z - c)^2)
        end
        for i in 1:n
            L̄ += L(X[i, 1], X[i, 2], X[i, 3], a, b, c)
            L̄a += (a - X[i, 1]) / L(X[i, 1], X[i, 2], X[i, 3], a, b, c)
            L̄b += (b - X[i, 2]) / L(X[i, 1], X[i, 2], X[i, 3], a, b, c)
            L̄c += (c - X[i, 3]) / L(X[i, 1], X[i, 2], X[i, 3], a, b, c)
        end
        L̄ /= n
        L̄a /= n
        L̄b /= n
        L̄c /= n

        return Xmean[1] + L̄ * L̄a, Xmean[2] + L̄ * L̄b, Xmean[3] + L̄ * L̄c, L̄
    end

    a = Xmean[1]
    b = Xmean[2]
    c = Xmean[3]
    r = 0.0
    for i in 1:iternum
        a, b, c, r = update(a, b, c, Xmean)
    end
    return a, b, c, r
end

type param
    θ₀::Float64
    θ₁::Float64
    C::Vector{Float64}
    sampleSize::Int64
end

function cylinderFit(X)

    function G(X, W, rSqr, params)
        S = zeros(3, 3)
        P = eye(3, 3)
        sinθ₀ = sin(params.θ₀)
        cosθ₀ = cos(params.θ₀)
        sinθ₁ = sin(params.θ₁)
        cosθ₁ = cos(params.θ₁)
        W[1] = cosθ₀ * sinθ₁
        W[2] = sinθ₀ * sinθ₁
        W[3] = cosθ₁
        S[1, 1] = 0.0
        S[1, 2] = -W[3]
        S[1, 3] = W[2]
        S[2, 1] = W[3]
        S[2, 2] = 0.0
        S[2, 3] = W[1]
        S[3, 1] = -W[2]
        S[3, 2] = W[1]
        S[3, 3] = 0.0
        P -= W * W'
        A = zeros(3, 3)
        B = zeros(3)
        Y = zeros(params.sampleSize, 3)
        averageSqrLength = 0.0
        sqrLength = zeros(params.sampleSize)
        for i in 1:params.sampleSize
            Y[i, :] = P * vec(X[i, :])
            sqrLength[i] = dot(vec(Y[i, :]), vec(Y[i, :]))
            A += vec(Y[i, :]) * vec(Y[i, :])'
            B += sqrLength[i] * vec(Y[i, :])
            averageSqrLength += sqrLength[i]
        end
        A /= params.sampleSize
        B /= params.sampleSize
        averageSqrLength /= params.sampleSize
        Ahat = - S * A * S
        params.C = (Ahat * B) / trace(Ahat * A)
        error = 0.0
        rSqr = 0.0
        for i in 1:params.sampleSize
            term = sqrLength[i] - averageSqrLength - 2 * dot(vec(Y[i, :]), params.C)
            error += term * term
            diff = params.C - vec(Y[i, :])
            rSqr += dot(diff, diff)
        end
        error /= params.sampleSize
        rSqr /= params.sampleSize
        return  error, W, rSqr, params
    end

    params = param(0, 0, zeros(3), size(X)[1])
    average = mean(X, 1)
    for i in 1:params.sampleSize
        X[i, :] -= average
    end
    jMax = 30
    iMax = 120
    W = zeros(3)
    C = zeros(3)
    θ₀ = 0.0
    θ₁ = 0.0
    rSqr = 0.0
    minError = Inf
    for i in 1:iMax
        params.θ₀ = 2 * pi * (i / iMax)
        for j in 1:jMax
            params.θ₁ = pi / 2.0 * (j / jMax)
            currentW = zeros(3)
            currentRSqr = 0.0
            error, currentW, currentRSqr, params = G(X, currentW, currentRSqr, params)
            if error < minError
                minError = error
                W =  currentW
                C = params.C
                θ₀ = params.θ₀
                θ₁ = params.θ₁
                rSqr = currentRSqr
            end
        end
    end
    return  W, C + vec(average), sqrt(rSqr), θ₀, θ₁
end


text:
using FileIO
using JLD2
using Distances
using PyPlot
using StatsBase
using Statistics

function read_evecs(dataset::String, k::Int64)
    data1 = load("results/$dataset-$k.jld2")
    data2 = load("results/$dataset-$k-random.jld2")

    # Get mode Z-eigenvector
    z_evecs = data2["z_evecs"]
    dists = pairwise(Cityblock(), z_evecs, dims=2)
    cutoff = 1e-5
    n = size(z_evecs, 2)
    # group vectors by distance
    components = zeros(Int64, n)
    for i = 1:n
        if components[i] == 0
            components[i] = maximum(components) + 1
        end
        for j = 1:n
            if dists[i, j] < cutoff && components[j] == 0
                components[j] = components[i]
            end
        end
    end
    # find mode vector index
    biggest = argmax(counts(components))
    z_mode_ind = 0
    for i = 1:n
        if components[i] == biggest
            z_mode_ind = i
            break
        end
    end

    return data1["cec_c"], data1["hec_c"], z_evecs[:, z_mode_ind]
end

function summary_statistics(dataset::String, k::Int64)
    data = load("results/$dataset-$k.jld2")
    nnodes, nnz = data["nnodes"], data["nnz"]
    println("# nodes: $nnodes")
    println("nnz (ignoring symmetries): $nnz")
end

# Top ranked items
function top_ranked(dataset::String, numtop::Int64, singlerow::Bool=false)
    label_mat = Array{String}(undef, numtop, 9)
    col_ind = 1
    for k in [3, 4, 5]
        data = load("results/$dataset-$(k).jld2")
        lcc_labels = data["lcc_labels"]
        cec_c, hec_c, zec_c = read_evecs(dataset, k)

        get_top_labels(c::Vector{Float64}) =
            [string(s) for s in lcc_labels[sortperm(c, rev=true)[1:numtop]]]

        label_mat[:, col_ind] = get_top_labels(cec_c); col_ind += 1
        label_mat[:, col_ind] = get_top_labels(zec_c); col_ind += 1
        label_mat[:, col_ind] = get_top_labels(hec_c); col_ind += 1

        local_label_mat = Array{String}(undef, numtop, 3)
        local_label_mat[:, 1] = get_top_labels(cec_c)
        local_label_mat[:, 2] = get_top_labels(zec_c)
        local_label_mat[:, 3] = get_top_labels(hec_c)

        if !singlerow
            for i in 1:numtop
                #str = [join(split(s)[2:end], " ") for s in collect(local_label_mat[i, :])]
                #str = join([lowercase(s) for s in str], " & ")
                str = join([lowercase(s) for s in collect(local_label_mat[i, :])], " & ")
                str = "& $i & $(str) \\\\"
                println(str)
            end
            println("--------------------")
        end
    end

    if singlerow
        for i in 1:numtop
            str = join([lowercase(s) for s in collect(label_mat[i, :])], " & ")
            str = "$i & $(str) \\\\"
            println(str)
        end
    end
end

# Rank correlation plot
function rank_corr(dataset::String, k::Int64, index_start::Int64)
    close()
    cec_c, hec_c, zec_c = read_evecs(dataset, k)
    @show length(cec_c), length(hec_c), length(zec_c)

    start = log10(index_start - 1)
    finish = log10(length(cec_c) - 1)
    ran = [convert(Int64, round(v)) for v in 10 .^ range(start, stop=finish, length=500)]

    sp = sortperm(cec_c)
    cec_c, hec_c, zec_c = cec_c[sp], hec_c[sp], zec_c[sp]
    cec_hec = [corspearman(cec_c[(end-s):end], hec_c[(end-s):end]) for s in ran]
    cec_zec = [corspearman(cec_c[(end-s):end], zec_c[(end-s):end]) for s in ran]
    sp = sortperm(zec_c)
    cec_c, hec_c, zec_c = cec_c[sp], hec_c[sp], zec_c[sp]
    zec_cec = [corspearman(zec_c[(end-s):end], cec_c[(end-s):end]) for s in ran]
    zec_hec = [corspearman(zec_c[(end-s):end], hec_c[(end-s):end]) for s in ran]
    sp = sortperm(hec_c)
    cec_c, hec_c, zec_c = cec_c[sp], hec_c[sp], zec_c[sp]
    hec_cec = [corspearman(hec_c[(end-s):end], cec_c[(end-s):end]) for s in ran]
    hec_zec = [corspearman(hec_c[(end-s):end], zec_c[(end-s):end]) for s in ran]

    xs = collect(ran) .+ 1
    semilogx(xs, cec_hec, linestyle="-",  lw=1,    label="CEC-HEC")
    semilogx(xs, cec_zec, linestyle="-",  lw=3.5,  label="CEC-ZEC")
    semilogx(xs, zec_cec, linestyle="--", lw=2.25, label="ZEC-CEC")
    semilogx(xs, zec_hec, linestyle="--", lw=1,    label="ZEC-HEC")
    semilogx(xs, hec_cec, linestyle=":",  lw=2.25, label="HEC-CEC")
    semilogx(xs, hec_zec, linestyle=":",  lw=3.5,  label="HEC-ZEC")

    fsz = 20
    ax = gca()
    ax[:tick_params](axis="both", labelsize=fsz-2, length=7, width=1.5)
    ax[:tick_params](axis="x",    which="minor", length=4, width=1)
    if k == 3; legend(fontsize=fsz-4, frameon=false); end
    title("$(dataset) ($(k)-uniform)", fontsize=fsz)
    xlabel("Number of top-ranked elements", fontsize=fsz)
    ylabel("Spearman's rank corr. coeff.", fontsize=fsz)
    savefig("$(dataset)-$(k)-rankcorr.eps", bbox_inches="tight")
end

function get_rank(dataset::String, k::Int, target_label::String)
    cec_c, hec_c, zec_c = read_evecs(dataset, k)
    n = length(cec_c)
    sp_cec = zeros(Int64, n); sp_cec[sortperm(cec_c, rev=true)] = collect(1:n)
    sp_zec = zeros(Int64, n); sp_zec[sortperm(zec_c, rev=true)] = collect(1:n)
    sp_hec = zeros(Int64, n); sp_hec[sortperm(hec_c, rev=true)] = collect(1:n)
    data = load("results/$dataset-$k.jld2")
    labels = data["lcc_labels"]

    for (label, cec_rk, zec_rk, hec_rk) in zip(labels, sp_cec, sp_zec, sp_hec)
        if occursin(target_label, string(label))
            println("$label $(cec_rk) $(zec_rk) $(hec_rk)")
        end
    end
end
;

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
#==============================================================================#
# CostExplorer.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/ce-2017-10-25.normal.json
#==============================================================================#

__precompile__()

module CostExplorer

using AWSCore

"""
    using AWSSDK.CostExplorer.get_cost_and_usage
    get_cost_and_usage([::AWSConfig], arguments::Dict)
    get_cost_and_usage([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.ce
    ce([::AWSConfig], "GetCostAndUsage", arguments::Dict)
    ce([::AWSConfig], "GetCostAndUsage", <keyword arguments>)

# GetCostAndUsage Operation

Retrieves cost and usage metrics for your account. You can specify which cost and usage-related metric, such as `BlendedCosts` or `UsageQuantity`, that you want the request to return. You can also filter and group your data by various dimensions, such as `SERVICE` or `AZ`, in a specific time range. For a complete list of valid dimensions, see the `[GetDimensionValues](http://docs.aws.amazon.com/aws-cost-management/latest/APIReference/API_GetDimensionValues.html)` operation. Master accounts in an organization in AWS Organizations have access to all member accounts.

# Arguments

## `TimePeriod = [ ... ]`
Sets the start and end dates for retrieving AWS costs. The start date is inclusive, but the end date is exclusive. For example, if `start` is `2017-01-01` and `end` is `2017-05-01`, then the cost and usage data is retrieved from `2017-01-01` up to and including `2017-04-30` but not including `2017-05-01`.
```
 TimePeriod = [
        "Start" => <required> ::String,
        "End" => <required> ::String
    ]
```

## `Granularity = "DAILY" or "MONTHLY"`
Sets the AWS cost granularity to `MONTHLY` or `DAILY`. If `Granularity` isn't set, the response object doesn't include the `Granularity`, either `MONTHLY` or `DAILY`.

## `Filter = [ ... ]`
Filters AWS costs by different dimensions. For example, you can specify `SERVICE` and `LINKED_ACCOUNT` and get the costs that are associated with that account's usage of that service. You can nest `Expression` objects to define any combination of dimension filters. For more information, see [Expression](http://docs.aws.amazon.com/aws-cost-management/latest/APIReference/API_Expression.html).
```
 Filter = [
        "Or" =>  [structure, ...],
        "And" =>  [structure, ...],
        "Not" =>  structure,
        "Dimensions" =>  [
            "Key" =>  "AZ", "INSTANCE_TYPE", "LINKED_ACCOUNT", "OPERATION", "PURCHASE_TYPE", "REGION", "SERVICE", "USAGE_TYPE", "USAGE_TYPE_GROUP", "RECORD_TYPE", "OPERATING_SYSTEM", "TENANCY", "SCOPE", "PLATFORM", "SUBSCRIPTION_ID", "LEGAL_ENTITY_NAME", "DEPLOYMENT_OPTION", "DATABASE_ENGINE", "CACHE_ENGINE" or "INSTANCE_TYPE_FAMILY",
            "Values" =>  [::String, ...]
        ],
        "Tags" =>  [
            "Key" =>  ::String,
            "Values" =>  [::String, ...]
        ]
    ]
```

## `Metrics = [::String, ...]`
Which metrics are returned in the query. For more information about blended and unblended rates, see [Why does the "blended" annotation appear on some line items in my bill?](https://aws.amazon.com/premiumsupport/knowledge-center/blended-rates-intro/).

Valid values are `AmortizedCost`, `BlendedCost`, `UnblendedCost`, and `UsageQuantity`.

**Note**
> If you return the `UsageQuantity` metric, the service aggregates all usage numbers without taking into account the units. For example, if you aggregate `usageQuantity` across all of EC2, the results aren't meaningful because EC2 compute hours and data transfer are measured in different units (for example, hours vs. GB). To get more meaningful `UsageQuantity` metrics, filter by `UsageType` or `UsageTypeGroups`.

`Metrics` is required for `GetCostAndUsage` requests.

## `GroupBy = [[ ... ], ...]`
You can group AWS costs using up to two different groups, either dimensions, tag keys, or both.

When you group by tag key, you get all tag values, including empty strings.

Valid values are `AZ`, `INSTANCE_TYPE`, `LEGAL_ENTITY_NAME`, `LINKED_ACCOUNT`, `OPERATION`, `PLATFORM`, `PURCHASE_TYPE`, `SERVICE`, `TAGS`, `TENANCY`, and `USAGE_TYPE`.
```
 GroupBy = [[
        "Type" =>  "DIMENSION" or "TAG",
        "Key" =>  ::String
    ], ...]
```

## `NextPageToken = ::String`
The token to retrieve the next set of results. AWS provides the token when the response from a previous call has more results than the maximum page size.

# Returns

`GetCostAndUsageResponse`

# Exceptions

`LimitExceededException`, `BillExpirationException`, `DataUnavailableException`, `InvalidNextTokenException` or `RequestChangedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/ce-2017-10-25/GetCostAndUsage)
"""
@inline get_cost_and_usage(aws::AWSConfig=default_aws_config(); args...) = get_cost_and_usage(aws, args)

@inline get_cost_and_usage(aws::AWSConfig, args) = AWSCore.Services.ce(aws, "GetCostAndUsage", args)

@inline get_cost_and_usage(args) = get_cost_and_usage(default_aws_config(), args)

"""
    using AWSSDK.CostExplorer.get_dimension_values
    get_dimension_values([::AWSConfig], arguments::Dict)
    get_dimension_values([::AWSConfig]; TimePeriod=, Dimension=, <keyword arguments>)

    using AWSCore.Services.ce
    ce([::AWSConfig], "GetDimensionValues", arguments::Dict)
    ce([::AWSConfig], "GetDimensionValues", TimePeriod=, Dimension=, <keyword arguments>)

# GetDimensionValues Operation

Retrieves all available filter values for a specified filter over a period of time. You can search the dimension values for an arbitrary string.

# Arguments

## `SearchString = ::String`
The value that you want to search the filter values for.

## `TimePeriod = [ ... ]` -- *Required*
The start and end dates for retrieving the dimension values. The start date is inclusive, but the end date is exclusive. For example, if `start` is `2017-01-01` and `end` is `2017-05-01`, then the cost and usage data is retrieved from `2017-01-01` up to and including `2017-04-30` but not including `2017-05-01`.
```
 TimePeriod = [
        "Start" => <required> ::String,
        "End" => <required> ::String
    ]
```

## `Dimension = "AZ", "INSTANCE_TYPE", "LINKED_ACCOUNT", "OPERATION", "PURCHASE_TYPE", "REGION", "SERVICE", "USAGE_TYPE", "USAGE_TYPE_GROUP", "RECORD_TYPE", "OPERATING_SYSTEM", "TENANCY", "SCOPE", "PLATFORM", "SUBSCRIPTION_ID", "LEGAL_ENTITY_NAME", "DEPLOYMENT_OPTION", "DATABASE_ENGINE", "CACHE_ENGINE" or "INSTANCE_TYPE_FAMILY"` -- *Required*
The name of the dimension. Each `Dimension` is available for different a `Context`. For more information, see `Context`.

## `Context = "COST_AND_USAGE" or "RESERVATIONS"`
The context for the call to `GetDimensionValues`. This can be `RESERVATIONS` or `COST_AND_USAGE`. The default value is `COST_AND_USAGE`. If the context is set to `RESERVATIONS`, the resulting dimension values can be used in the `GetReservationUtilization` operation. If the context is set to `COST_AND_USAGE` the resulting dimension values can be used in the `GetCostAndUsage` operation.

If you set the context to `COST_AND_USAGE`, you can use the following dimensions for searching:

*   AZ - The Availability Zone. An example is `us-east-1a`.

*   DATABASE_ENGINE - The Amazon Relational Database Service database. Examples are Aurora or MySQL.

*   INSTANCE_TYPE - The type of EC2 instance. An example is `m4.xlarge`.

*   LEGAL_ENTITY_NAME - The name of the organization that sells you AWS services, such as Amazon Web Services.

*   LINKED_ACCOUNT - The description in the attribute map that includes the full name of the member account. The value field contains the AWS ID of the member account.

*   OPERATING_SYSTEM - The operating system. Examples are Windows or Linux.

*   OPERATION - The action performed. Examples include `RunInstance` and `CreateBucket`.

*   PLATFORM - The EC2 operating system. Examples are Windows or Linux.

*   PURCHASE_TYPE - The reservation type of the purchase to which this usage is related. Examples include On-Demand Instances and Standard Reserved Instances.

*   SERVICE - The AWS service such as Amazon DynamoDB.

*   USAGE_TYPE - The type of usage. An example is DataTransfer-In-Bytes. The response for the `GetDimensionValues` operation includes a unit attribute. Examples include GB and Hrs.

*   USAGE_TYPE_GROUP - The grouping of common usage types. An example is EC2: CloudWatch – Alarms. The response for this operation includes a unit attribute.

*   RECORD_TYPE - The different types of charges such as RI fees, usage costs, tax refunds, and credits.

If you set the context to `RESERVATIONS`, you can use the following dimensions for searching:

*   AZ - The Availability Zone. An example is `us-east-1a`.

*   CACHE_ENGINE - The Amazon ElastiCache operating system. Examples are Windows or Linux.

*   DEPLOYMENT_OPTION - The scope of Amazon Relational Database Service deployments. Valid values are `SingleAZ` and `MultiAZ`.

*   INSTANCE_TYPE - The type of EC2 instance. An example is `m4.xlarge`.

*   LINKED_ACCOUNT - The description in the attribute map that includes the full name of the member account. The value field contains the AWS ID of the member account.

*   PLATFORM - The EC2 operating system. Examples are Windows or Linux.

*   REGION - The AWS Region.

*   SCOPE (Utilization only) - The scope of a Reserved Instance (RI). Values are regional or a single Availability Zone.

*   TAG (Coverage only) - The tags that are associated with a Reserved Instance (RI).

*   TENANCY - The tenancy of a resource. Examples are shared or dedicated.

## `NextPageToken = ::String`
The token to retrieve the next set of results. AWS provides the token when the response from a previous call has more results than the maximum page size.

# Returns

`GetDimensionValuesResponse`

# Exceptions

`LimitExceededException`, `BillExpirationException`, `DataUnavailableException`, `InvalidNextTokenException` or `RequestChangedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/ce-2017-10-25/GetDimensionValues)
"""
@inline get_dimension_values(aws::AWSConfig=default_aws_config(); args...) = get_dimension_values(aws, args)

@inline get_dimension_values(aws::AWSConfig, args) = AWSCore.Services.ce(aws, "GetDimensionValues", args)

@inline get_dimension_values(args) = get_dimension_values(default_aws_config(), args)

"""
    using AWSSDK.CostExplorer.get_reservation_coverage
    get_reservation_coverage([::AWSConfig], arguments::Dict)
    get_reservation_coverage([::AWSConfig]; TimePeriod=, <keyword arguments>)

    using AWSCore.Services.ce
    ce([::AWSConfig], "GetReservationCoverage", arguments::Dict)
    ce([::AWSConfig], "GetReservationCoverage", TimePeriod=, <keyword arguments>)

# GetReservationCoverage Operation

Retrieves the reservation coverage for your account. This allows you to see how much of your Amazon Elastic Compute Cloud, Amazon ElastiCache, Amazon Relational Database Service, or Amazon Redshift usage is covered by a reservation. An organization's master account can see the coverage of the associated member accounts. For any time period, you can filter data about reservation usage by the following dimensions:

*   AZ

*   CACHE_ENGINE

*   DATABASE_ENGINE

*   DEPLOYMENT_OPTION

*   INSTANCE_TYPE

*   LINKED_ACCOUNT

*   OPERATING_SYSTEM

*   PLATFORM

*   REGION

*   SERVICE

*   TAG

*   TENANCY

To determine valid values for a dimension, use the `GetDimensionValues` operation.

# Arguments

## `TimePeriod = [ ... ]` -- *Required*
The start and end dates of the period for which you want to retrieve data about reservation coverage. You can retrieve data for a maximum of 13 months: the last 12 months and the current month. The start date is inclusive, but the end date is exclusive. For example, if `start` is `2017-01-01` and `end` is `2017-05-01`, then the cost and usage data is retrieved from `2017-01-01` up to and including `2017-04-30` but not including `2017-05-01`.
```
 TimePeriod = [
        "Start" => <required> ::String,
        "End" => <required> ::String
    ]
```

## `GroupBy = [[ ... ], ...]`
You can group the data by the following attributes:

*   AZ

*   CACHE_ENGINE

*   DATABASE_ENGINE

*   DEPLOYMENT_OPTION

*   INSTANCE_TYPE

*   LINKED_ACCOUNT

*   OPERATING_SYSTEM

*   PLATFORM

*   REGION

*   TAG

*   TENANCY
```
 GroupBy = [[
        "Type" =>  "DIMENSION" or "TAG",
        "Key" =>  ::String
    ], ...]
```

## `Granularity = "DAILY" or "MONTHLY"`
The granularity of the AWS cost data for the reservation. Valid values are `MONTHLY` and `DAILY`.

If `GroupBy` is set, `Granularity` can't be set. If `Granularity` isn't set, the response object doesn't include `Granularity`, either `MONTHLY` or `DAILY`.

## `Filter = [ ... ]`
Filters utilization data by dimensions. You can filter by the following dimensions:

*   AZ

*   CACHE_ENGINE

*   DATABASE_ENGINE

*   DEPLOYMENT_OPTION

*   INSTANCE_TYPE

*   LINKED_ACCOUNT

*   OPERATING_SYSTEM

*   PLATFORM

*   REGION

*   SERVICE

*   TAG

*   TENANCY

`GetReservationCoverage` uses the same `[Expression](http://docs.aws.amazon.com/aws-cost-management/latest/APIReference/API_Expression.html)` object as the other operations, but only `AND` is supported among each dimension. You can nest only one level deep. If there are multiple values for a dimension, they are OR'd together.

If you don't provide a `SERVICE` filter, Cost Explorer defaults to EC2.
```
 Filter = [
        "Or" =>  [structure, ...],
        "And" =>  [structure, ...],
        "Not" =>  structure,
        "Dimensions" =>  [
            "Key" =>  "AZ", "INSTANCE_TYPE", "LINKED_ACCOUNT", "OPERATION", "PURCHASE_TYPE", "REGION", "SERVICE", "USAGE_TYPE", "USAGE_TYPE_GROUP", "RECORD_TYPE", "OPERATING_SYSTEM", "TENANCY", "SCOPE", "PLATFORM", "SUBSCRIPTION_ID", "LEGAL_ENTITY_NAME", "DEPLOYMENT_OPTION", "DATABASE_ENGINE", "CACHE_ENGINE" or "INSTANCE_TYPE_FAMILY",
            "Values" =>  [::String, ...]
        ],
        "Tags" =>  [
            "Key" =>  ::String,
            "Values" =>  [::String, ...]
        ]
    ]
```

## `NextPageToken = ::String`
The token to retrieve the next set of results. AWS provides the token when the response from a previous call has more results than the maximum page size.

# Returns

`GetReservationCoverageResponse`

# Exceptions

`LimitExceededException`, `DataUnavailableException` or `InvalidNextTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/ce-2017-10-25/GetReservationCoverage)
"""
@inline get_reservation_coverage(aws::AWSConfig=default_aws_config(); args...) = get_reservation_coverage(aws, args)

@inline get_reservation_coverage(aws::AWSConfig, args) = AWSCore.Services.ce(aws, "GetReservationCoverage", args)

@inline get_reservation_coverage(args) = get_reservation_coverage(default_aws_config(), args)

"""
    using AWSSDK.CostExplorer.get_reservation_purchase_recommendation
    get_reservation_purchase_recommendation([::AWSConfig], arguments::Dict)
    get_reservation_purchase_recommendation([::AWSConfig]; Service=, <keyword arguments>)

    using AWSCore.Services.ce
    ce([::AWSConfig], "GetReservationPurchaseRecommendation", arguments::Dict)
    ce([::AWSConfig], "GetReservationPurchaseRecommendation", Service=, <keyword arguments>)

# GetReservationPurchaseRecommendation Operation

Gets recommendations for which reservations to purchase. These recommendations could help you reduce your costs. Reservations provide a discounted hourly rate (up to 75%) compared to On-Demand pricing.

AWS generates your recommendations by identifying your On-Demand usage during a specific time period and collecting your usage into categories that are eligible for a reservation. After AWS has these categories, it simulates every combination of reservations in each category of usage to identify the best number of each type of RI to purchase to maximize your estimated savings.

For example, AWS automatically aggregates your EC2 Linux, shared tenancy, and c4 family usage in the US West (Oregon) Region and recommends that you buy size-flexible regional reservations to apply to the c4 family usage. AWS recommends the smallest size instance in an instance family. This makes it easier to purchase a size-flexible RI. AWS also shows the equal number of normalized units so that you can purchase any instance size that you want. For this example, your RI recommendation would be for `c4.large`, because that is the smallest size instance in the c4 instance family.

# Arguments

## `AccountId = ::String`
The account ID that is associated with the recommendation.

## `Service = ::String` -- *Required*
The specific service that you want recommendations for.

## `AccountScope = "PAYER" or "LINKED"`
The account scope that you want recommendations for. `PAYER` means that AWS includes the master account and any member accounts when it calculates its recommendations. `LINKED` means that AWS includes only member accounts when it calculates its recommendations.

Valid values are `PAYER` and `LINKED`.

## `LookbackPeriodInDays = "SEVEN_DAYS", "THIRTY_DAYS" or "SIXTY_DAYS"`
The number of previous days that you want AWS to consider when it calculates your recommendations.

## `TermInYears = "ONE_YEAR" or "THREE_YEARS"`
The reservation term that you want recommendations for.

## `PaymentOption = "NO_UPFRONT", "PARTIAL_UPFRONT", "ALL_UPFRONT", "LIGHT_UTILIZATION", "MEDIUM_UTILIZATION" or "HEAVY_UTILIZATION"`
The reservation purchase option that you want recommendations for.

## `ServiceSpecification = ["EC2Specification" =>  ["OfferingClass" =>  "STANDARD" or "CONVERTIBLE"]]`
The hardware specifications for the service instances that you want recommendations for, such as standard or convertible EC2 instances.

## `PageSize = ::Int`
The number of recommendations that you want returned in a single response object.

## `NextPageToken = ::String`
The pagination token that indicates the next set of results that you want to retrieve.

# Returns

`GetReservationPurchaseRecommendationResponse`

# Exceptions

`LimitExceededException`, `DataUnavailableException` or `InvalidNextTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/ce-2017-10-25/GetReservationPurchaseRecommendation)
"""
@inline get_reservation_purchase_recommendation(aws::AWSConfig=default_aws_config(); args...) = get_reservation_purchase_recommendation(aws, args)

@inline get_reservation_purchase_recommendation(aws::AWSConfig, args) = AWSCore.Services.ce(aws, "GetReservationPurchaseRecommendation", args)

@inline get_reservation_purchase_recommendation(args) = get_reservation_purchase_recommendation(default_aws_config(), args)

"""
    using AWSSDK.CostExplorer.get_reservation_utilization
    get_reservation_utilization([::AWSConfig], arguments::Dict)
    get_reservation_utilization([::AWSConfig]; TimePeriod=, <keyword arguments>)

    using AWSCore.Services.ce
    ce([::AWSConfig], "GetReservationUtilization", arguments::Dict)
    ce([::AWSConfig], "GetReservationUtilization", TimePeriod=, <keyword arguments>)

# GetReservationUtilization Operation

Retrieves the reservation utilization for your account. Master accounts in an organization have access to member accounts. You can filter data by dimensions in a time period. You can use `GetDimensionValues` to determine the possible dimension values. Currently, you can group only by `SUBSCRIPTION_ID`.

# Arguments

## `TimePeriod = [ ... ]` -- *Required*
Sets the start and end dates for retrieving Reserved Instance (RI) utilization. The start date is inclusive, but the end date is exclusive. For example, if `start` is `2017-01-01` and `end` is `2017-05-01`, then the cost and usage data is retrieved from `2017-01-01` up to and including `2017-04-30` but not including `2017-05-01`.
```
 TimePeriod = [
        "Start" => <required> ::String,
        "End" => <required> ::String
    ]
```

## `GroupBy = [[ ... ], ...]`
Groups only by `SUBSCRIPTION_ID`. Metadata is included.
```
 GroupBy = [[
        "Type" =>  "DIMENSION" or "TAG",
        "Key" =>  ::String
    ], ...]
```

## `Granularity = "DAILY" or "MONTHLY"`
If `GroupBy` is set, `Granularity` can't be set. If `Granularity` isn't set, the response object doesn't include `Granularity`, either `MONTHLY` or `DAILY`. If both `GroupBy` and `Granularity` aren't set, `GetReservationUtilization` defaults to `DAILY`.

## `Filter = [ ... ]`
Filters utilization data by dimensions. You can filter by the following dimensions:

*   AZ

*   CACHE_ENGINE

*   DATABASE_ENGINE

*   DEPLOYMENT_OPTION

*   INSTANCE_TYPE

*   LINKED_ACCOUNT

*   OPERATING_SYSTEM

*   PLATFORM

*   REGION

*   SERVICE

*   SCOPE

*   TENANCY

`GetReservationUtilization` uses the same `[Expression](http://docs.aws.amazon.com/aws-cost-management/latest/APIReference/API_Expression.html)` object as the other operations, but only `AND` is supported among each dimension, and nesting is supported up to only one level deep. If there are multiple values for a dimension, they are OR'd together.
```
 Filter = [
        "Or" =>  [structure, ...],
        "And" =>  [structure, ...],
        "Not" =>  structure,
        "Dimensions" =>  [
            "Key" =>  "AZ", "INSTANCE_TYPE", "LINKED_ACCOUNT", "OPERATION", "PURCHASE_TYPE", "REGION", "SERVICE", "USAGE_TYPE", "USAGE_TYPE_GROUP", "RECORD_TYPE", "OPERATING_SYSTEM", "TENANCY", "SCOPE", "PLATFORM", "SUBSCRIPTION_ID", "LEGAL_ENTITY_NAME", "DEPLOYMENT_OPTION", "DATABASE_ENGINE", "CACHE_ENGINE" or "INSTANCE_TYPE_FAMILY",
            "Values" =>  [::String, ...]
        ],
        "Tags" =>  [
            "Key" =>  ::String,
            "Values" =>  [::String, ...]
        ]
    ]
```

## `NextPageToken = ::String`
The token to retrieve the next set of results. AWS provides the token when the response from a previous call has more results than the maximum page size.

# Returns

`GetReservationUtilizationResponse`

# Exceptions

`LimitExceededException`, `DataUnavailableException` or `InvalidNextTokenException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/ce-2017-10-25/GetReservationUtilization)
"""
@inline get_reservation_utilization(aws::AWSConfig=default_aws_config(); args...) = get_reservation_utilization(aws, args)

@inline get_reservation_utilization(aws::AWSConfig, args) = AWSCore.Services.ce(aws, "GetReservationUtilization", args)

@inline get_reservation_utilization(args) = get_reservation_utilization(default_aws_config(), args)

"""
    using AWSSDK.CostExplorer.get_tags
    get_tags([::AWSConfig], arguments::Dict)
    get_tags([::AWSConfig]; TimePeriod=, <keyword arguments>)

    using AWSCore.Services.ce
    ce([::AWSConfig], "GetTags", arguments::Dict)
    ce([::AWSConfig], "GetTags", TimePeriod=, <keyword arguments>)

# GetTags Operation

Queries for available tag keys and tag values for a specified period. You can search the tag values for an arbitrary string.

# Arguments

## `SearchString = ::String`
The value that you want to search for.

## `TimePeriod = [ ... ]` -- *Required*
The start and end dates for retrieving the dimension values. The start date is inclusive, but the end date is exclusive. For example, if `start` is `2017-01-01` and `end` is `2017-05-01`, then the cost and usage data is retrieved from `2017-01-01` up to and including `2017-04-30` but not including `2017-05-01`.
```
 TimePeriod = [
        "Start" => <required> ::String,
        "End" => <required> ::String
    ]
```

## `TagKey = ::String`
The key of the tag that you want to return values for.

## `NextPageToken = ::String`
The token to retrieve the next set of results. AWS provides the token when the response from a previous call has more results than the maximum page size.

# Returns

`GetTagsResponse`

# Exceptions

`LimitExceededException`, `BillExpirationException`, `DataUnavailableException`, `InvalidNextTokenException` or `RequestChangedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/ce-2017-10-25/GetTags)
"""
@inline get_tags(aws::AWSConfig=default_aws_config(); args...) = get_tags(aws, args)

@inline get_tags(aws::AWSConfig, args) = AWSCore.Services.ce(aws, "GetTags", args)

@inline get_tags(args) = get_tags(default_aws_config(), args)

end # module CostExplorer

#==============================================================================#
# End of file
#==============================================================================#


text:
function normalize_theta!(scales::AbstractArray, θ::AbstractArray)
    @assert length(scales) == size(θ, 1)
    @inbounds for (i, ti) in enumerate(eachrow(θ))
        scales[i] = norm(ti, 2)
        normalize!(ti, 2)
    end
    return
end

function rescale_xi!(Ξ::AbstractArray, scales::AbstractArray)
    @assert length(scales) == size(Ξ, 1)
    @inbounds for (si, ti) in zip(scales, eachrow(Ξ))
        ti .= ti / si
    end
    return
end

function rescale_theta!(θ::AbstractArray, scales::AbstractArray)
    @assert length(scales) == size(θ, 1)
    @inbounds for (i, ti) in enumerate(eachrow(θ))
        ti .= ti * scales[i]
    end
    return
end

"""
    sparse_regression(X, Y, basis, p, t, maxiter, opt, denoise, normalize, convergence_error)
    sparse_regression!(Xi, X, Y, basis, p, t, maxiter, opt, denoise, normalize, convergence_error)
    spares_regression!(Xi, Theta, Y, maxiter, opt, denoise, normalize, convergence_error)

Performs a sparse regression via the algorithm `opt <: AbstractOptimizer`. `maxiter` specifies the upper bound of the iterations
of the optimizer, `convergence_error` the breaking condition due to convergence.
`denoise` defines if the matrix holding candidate trajectories should be thresholded via the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
`normalize` normalizes the matrix holding candidate trajectories via the L2-Norm over each function.

If the data matrices `X`, `Y` are given with a `Basis` `basis` and the additional information for parameters `p` and time points of the
measurements `t`, it returns the coefficient matrix `Xi` and the iterations taken.
This function is also available in place, which returns just the iterations.

If `Xi`, `Theta` and `Y` are given, the sparse regression will find the coefficients `Xi`, which minimize the objective and return the iterations needed.

# Example

```julia
opt = STRRidge()
maxiter = 10
c_error = 1e-3

Xi, iters = sparse_regression(X, Y, basis, [], [], maxiter, opt, false, false, c_error)

iters = sparse_regression!(Xi,X, Y, basis, [], [], maxiter, opt, false, false, c_error)

Xi2 = zeros(size(Y, 1), size(X, 1))
iters = sparse_regression!(Xi2, X, Y, maxiter, opt, false, false, c_error)
```
"""
function sparse_regression(X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis, p::AbstractArray, t::AbstractVector , maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool, convergence_error) where T <: Optimize.AbstractOptimizer
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)

    Ξ = zeros(eltype(X), length(Ψ), ny)
    scales = ones(eltype(X), length(Ψ))
    θ = Ψ(X, p, t)

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimize.init!(Ξ, opt, θ', Ẋ')
    iters = Optimize.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter, convergence_error = convergence_error)

    normalize ? rescale_xi!(Ξ, scales) : nothing

    return Ξ, iters
end

function sparse_regression!(Ξ::AbstractArray, X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis, p::AbstractArray , t::AbstractVector, maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool, convergence_error) where T <: Optimize.AbstractOptimizer
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)
    @assert size(Ξ) == (length(Ψ), ny)

    scales = ones(eltype(X), length(Ψ))
    θ = Ψ(X, p, t)

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimize.init!(Ξ, opt, θ', Ẋ')
    iters = Optimize.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter, convergence_error = convergence_error)

    normalize ? rescale_xi!(Ξ, scales) : nothing

    return iters
end

# For pareto
function sparse_regression!(Ξ::AbstractArray, θ::AbstractArray, Ẋ::AbstractArray, maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool, convergence_error) where T <: Optimize.AbstractOptimizer

    scales = ones(eltype(Ξ), size(θ, 1))

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimize.init!(Ξ, opt, θ', Ẋ')'
    iters = Optimize.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter, convergence_error = convergence_error)

    normalize ? rescale_xi!(Ξ, scales) : nothing
    normalize ? rescale_theta!(θ, scales) : nothing

    return iters
end

# One Variable on multiple derivatives
function SInDy(X::AbstractArray{S, 1}, Ẋ::AbstractArray, Ψ::Basis; kwargs...) where S <: Number
    return SInDy(X', Ẋ, Ψ; kwargs...)
end

# Multiple on one
function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 1}, Ψ::Basis; kwargs...) where S <: Number
    return SInDy(X, Ẋ', Ψ; kwargs...)
end

# General
"""
    SInDy(X, Y, basis; p, t, opt, maxiter, convergence_error, denoise, normalize)
    SInDy(X, Y, basis, lambdas; weights, f_target, p, t, opt, maxiter, convergence_error, denoise, normalize)

Performs Sparse Identification of Nonlinear Dynamics given the data matrices `X` and `Y` via the `AbstractBasis` `basis.`
Keyworded arguments include the parameter (values) of the basis `p` and the timepoints `t` which are passed in optionally.
`opt` is an `AbstractOptimizer` useable for sparse regression, `maxiter` the maximum iterations to perform and `convergence_error` the
bound which causes the optimizer to stop.
`denoise` defines if the matrix holding candidate trajectories should be thresholded via the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
`normalize` normalizes the matrix holding candidate trajectories via the L2-Norm over each function.

If `SInDy` is called with an additional array of thresholds contained in `lambdas`, it performs a multi objective optimization over all thresholds.
The best candidate is determined via the `AbstractScalarizationMethod` given in `alg`. The evaluation has two fields, the sparsity of the coefficients and the L2-Norm error at index 1 and 2 respectively.

Returns a `SInDyResult`. If the pareto optimization is used, the result combines the best candidate for each row of `Y`.
"""
function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 2}, Ψ::Basis; p::AbstractArray = [], t::AbstractVector = [], maxiter::Int64 = 10, opt::T = Optimize.STRRidge(), denoise::Bool = false, normalize::Bool = true, convergence_error = eps()) where {T <: Optimize.AbstractOptimizer, S <: Number}
    Ξ, iters = sparse_regression(X, Ẋ, Ψ, p, t, maxiter, opt, denoise, normalize, convergence_error)
    convergence = iters < maxiter
    SparseIdentificationResult(Ξ, Ψ, iters, opt, convergence, Ẋ, X, p = p)
end

function SInDy(X::AbstractArray{S, 1}, Ẋ::AbstractArray, Ψ::Basis, thresholds::AbstractArray; kwargs...) where S <: Number
    return SInDy(X', Ẋ, Ψ, thresholds; kwargs...)
end

function SInDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 1}, Ψ::Basis, thresholds::AbstractArray; kwargs...) where S <: Number
    return SInDy(X, Ẋ', Ψ, thresholds; kwargs...)
end

function SInDy(X::AbstractArray{S, 2}, DX::AbstractArray{S, 2}, Ψ::Basis, thresholds::AbstractArray ; alg::Optimize.AbstractScalarizationMethod = WeightedSum(), p::AbstractArray = [], t::AbstractVector = [], maxiter::Int64 = 10, opt::T = Optimize.STRRidge(),denoise::Bool = false, normalize::Bool = true, convergence_error = eps()) where {T <: Optimize.AbstractOptimizer, S <: Number}
    @assert size(X)[end] == size(DX)[end]
    nx, nm = size(X)
    ny, nm = size(DX)

    θ = Ψ(X, p, t)

    scales = ones(eltype(X), length(Ψ))

    ξ = zeros(eltype(X), length(Ψ), ny)
    Ξ = zeros(eltype(X), length(Ψ), ny)

    iters = 0

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? DataDrivenDiffEq.normalize_theta!(scales, θ) : nothing

    # Set two paretofronts
    opt_front = ParetoFront(ny, scalarization = alg)
    tmp_front = ParetoFront(ny, scalarization = alg)

    for (j, threshold) in enumerate(thresholds)
        set_threshold!(opt, threshold)
        iters = sparse_regression!(ξ, θ, DX, maxiter, opt, false, false, convergence_error)
        normalize ? DataDrivenDiffEq.rescale_xi!(ξ, scales) : nothing

        if j < 2
            for (i, ξi) in enumerate(eachcol(ξ))
                set_candidate!(tmp_front, i, [norm(ξi, 0); norm(DX[i, :] .- θ'*ξi)], ξi, iters, threshold)
                set_candidate!(opt_front, i, [norm(ξi, 0); norm(DX[i, :] .- θ'*ξi)], ξi, iters, threshold)
            end
        else
            @inbounds for (i, ξi) in enumerate(eachcol(ξ))
                set_candidate!(tmp_front, i, [norm(ξi, 0); norm(DX[i, :] .- θ'*ξi)], ξi, iters, threshold)
            end
            conditional_add!(opt_front, tmp_front)
        end
    end

    for i in 1:ny
        Ξ[:, i] .= parameter(opt_front[i])
    end

    set_threshold!(opt, threshold(opt_front))

    return SparseIdentificationResult(Ξ, Ψ, iter(opt_front), opt, iter(opt_front) < maxiter, DX, X, p = p)
end


text:
 using Plots
 pdf_fixed(v) = pdf(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
 cdf_fixed(v) = cdf(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
 cdf_gradient_fixed(v) = cdf_gradient(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
 cdf_hessian_fixed(v) = cdf_hessian(v[2:end], linear_polynomial_basis(v[2:end]), v[1])

pdf_vstar = y -> pdf_fixed(vcat(y, vstar[2:end]))
cdf_vstar = y -> cdf_fixed(vcat(y, vstar[2:end]))
#plt(cdf_vstar, .0001, 1000.0)
#plt(pdf_vstar, .0001, 50)

function vcdfplot(vstar; lower = 0.0, upper=1000)
    cdf_vstar = y -> cdf_fixed(vcat(y, vstar[2:end]))
    plt(cdf_vstar, lower, upper)
end

"""
Multiple graphs per plot
"""
function vcdfplot_sequence(vstars; lower = 0.0, upper = 1000, res = nothing)
    cdf_vstar = y -> cdf_fixed(vcat(y, vstar[2:end]))
    Plots.plot()
    for i = 1:length(vstars)
        vstar = vstars[i]
        cdf_vstar = y -> cdf_fixed(vcat(y, vstar[2:end]))
        if res!=nothing
            #label = abs(res[i]-0.25) < .01 ? "converged" : "did not converge"
            color = abs(res[i]-0.25) < .01 ? "blue" : "red"
            plt!(cdf_vstar, lower, upper; color = color)
        end
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
struct Normalize{S <: PreMetric} <: PreMetric
    dist::S
end
"""
   normalize(dist::PreMetric)

   Normalize a metric, so that `evaluate` always return a Float64 between 0 and 1 (or a `missing` if one element is missing)
"""
function normalize(dist::PreMetric)
    isnormalized(dist) ? dist : Normalize{typeof(dist)}(dist)
end
isnormalized(dist::Normalize) = true

function evaluate(dist::Normalize{<: Union{Levenshtein, DamerauLevenshtein}}, s1, s2, max_dist = 1.0)
    (ismissing(s1) | ismissing(s2)) && return missing
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 1.0
    d = evaluate(dist.dist, s1, s2, ceil(Int, len2 * max_dist))
    out = d / len2
    out > max_dist ? 1.0 : out
end

function evaluate(dist::Normalize{<: QGramDistance}, s1, s2, max_dist = 1.0)
    (ismissing(s1) | ismissing(s2)) && return missing
    # When string length < q for qgram distance, returns s1 == s2
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 <= dist.dist.q - 1 && return convert(Float64, s1 != s2)
    if typeof(dist.dist) <: QGram
        evaluate(dist.dist, s1, s2) / (len1 + len2 - 2 * dist.dist.q + 2)
    else
        evaluate(dist.dist, s1, s2)
    end
end

"""
   Winkler(dist; p::Real = 0.1, threshold::Real = 0.7, maxlength::Integer = 4)

Creates the `Winkler{dist, p, threshold, maxlength}` distance

`Winkler{dist, p, threshold, length)` modifies the string distance `dist` to decrease the
distance between  two strings, when their original distance is below some `threshold`.
The boost is equal to `min(l,  maxlength) * p * dist` where `l` denotes the
length of their common prefix and `dist` denotes the original distance
"""
struct Winkler{S <: PreMetric} <: PreMetric
    dist::S
    p::Float64          # scaling factor. Default to 0.1
    threshold::Float64  # boost threshold. Default to 0.7
    maxlength::Integer      # max length of common prefix. Default to 4
    Winkler{S}(dist::S, p, threshold, maxlength) where {S <: PreMetric} = new(dist, p, threshold, maxlength)
end

function Winkler(dist::PreMetric; p = 0.1, threshold = 0.7, maxlength = 4)
    p * maxlength <= 1 || throw("scaling factor times maxlength of common prefix must be lower than one")
    Winkler{typeof(normalize(dist))}(normalize(dist), 0.1, 0.7, 4)
end
isnormalized(dist::Winkler) = true

function evaluate(dist::Winkler, s1, s2, max_dist = 1.0)
    # cannot do min_score because of boosting threshold
    score = evaluate(dist.dist, s1, s2)
    if score <= 1 - dist.threshold
        l = common_prefix(s1, s2)[1]
        score -= min(l, dist.maxlength) * dist.p * score
    end
    return score
end

"""
   Partial(dist)

Creates the `Partial{dist}` distance

`Partial{dist}` modifies the string distance `dist` to return the
minimum distance  between the shorter string and substrings of the longer string

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta Braves"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> evaluate(Partial(RatcliffObershelp()), s1, s2)
0.5483870967741935
```
"""
struct Partial{S <: PreMetric} <: PreMetric
    dist::S
    Partial{S}(dist::S) where {S <: PreMetric} = new(dist)
end
Partial(dist::PreMetric) = Partial{typeof(normalize(dist))}(normalize(dist))
isnormalized(dist::Partial) = true

function evaluate(dist::Partial, s1, s2, max_dist = 1.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 == len2 && return evaluate(dist.dist, s1, s2, max_dist)
    len1 == 0 && return 0.0
    out = 1.0
    for x in qgrams(s2, len1)
        curr = evaluate(dist.dist, s1, x, max_dist)
        out = min(out, curr)
        max_dist = min(out, max_dist)
    end
    return out
end

function evaluate(dist::Partial{RatcliffObershelp}, s1, s2, max_dist = 1.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 == len2 && return evaluate(dist.dist, s1, s2)
    out = 1.0
    for r in matching_blocks(s1, s2)
        # Make sure the substring of s2 has length len1
        s2_start = r[2] - r[1] + 1
        s2_end = s2_start + len1 - 1
        if s2_start <= 0
            s2_end += 1 - s2_start
            s2_start += 1 - s2_start
        elseif s2_end > len2
            s2_start += len2 - s2_end
            s2_end += len2 - s2_end
        end
        curr = evaluate(dist.dist, s1, _slice(s2, s2_start - 1, s2_end))

        out = min(out, curr)
    end
    return out
end

"""
   TokenSort(dist)

Creates the `TokenSort{dist}` distance

`TokenSort{dist}` modifies the string distance `dist` to adjust for differences
in word orders by reording words alphabetically.

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta Braves"
julia> s1 = "New York Mets vs Atlanta Braves"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> evaluate(TokenSort(RatcliffObershelp()), s1, s2)
0.0
```
"""
struct TokenSort{S <: PreMetric} <: PreMetric
    dist::S
    TokenSort{S}(dist::S) where {S <: PreMetric} = new(dist)
end
TokenSort(dist::PreMetric) = TokenSort{typeof(normalize(dist))}(normalize(dist))
isnormalized(dist::TokenSort) = true

# http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
function evaluate(dist::TokenSort, s1::AbstractString, s2::AbstractString, max_dist = 1.0)
    s1 = join(sort!(split(s1)), " ")
    s2 = join(sort!(split(s2)), " ")
    evaluate(dist.dist, s1, s2, max_dist)
end

"""
   TokenSet(dist)

Creates the `TokenSet{dist}` distance

`TokenSet{dist}` modifies the string distance `dist` to adjust for differences
in word orders and word numbers by comparing the intersection of two strings with each string.

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> evaluate(TokenSet(RatcliffObershelp()), s1, s2)
0.0
```
"""
struct TokenSet{S <: PreMetric} <: PreMetric
    dist::S
    TokenSet{S}(dist::S) where {S <: PreMetric} = new(dist)
end
TokenSet(dist::PreMetric) = TokenSet{typeof(normalize(dist))}(normalize(dist))
isnormalized(dist::TokenSet) = true

# http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
function evaluate(dist::TokenSet, s1::AbstractString, s2::AbstractString, max_dist = 1.0)
    v1 = unique!(sort!(split(s1)))
    v2 = unique!(sort!(split(s2)))
    v0 = intersect(v1, v2)
    s0 = join(v0, " ")
    s1 = join(v1, " ")
    s2 = join(v2, " ")
    isempty(s0) && return evaluate(dist.dist, s1, s2, max_dist)
    score_01 = evaluate(dist.dist, s0, s1, max_dist)
    max_dist = min(max_dist, score_01)
    score_02 = evaluate(dist.dist, s0, s2, max_dist)
    max_dist = min(max_dist, score_02)
    score_12 = evaluate(dist.dist, s1, s2, max_dist)
    min(score_01, score_02, score_12)
end

"""
   TokenMax(dist)

Creates the `TokenMax{dist}` distance

`TokenMax{dist}` is the minimum of the base distance `dist`,
its [`Partial`](@ref) modifier, its [`TokenSort`](@ref) modifier, and its
[`TokenSet`](@ref) modifier, with penalty terms depending on string lengths.

### Examples
```julia-repl
julia> s1 = "New York Mets vs Atlanta"
julia> s2 = "Atlanta Braves vs New York Mets"
julia> evaluate(TokenMax(RatcliffObershelp()), s1, s2)
0.05
```
"""
struct TokenMax{S <: PreMetric} <: PreMetric
    dist::S
    TokenMax{S}(dist::S) where {S <: PreMetric} = new(dist)
end

TokenMax(dist::PreMetric) = TokenMax{typeof(normalize(dist))}(normalize(dist))
isnormalized(dist::TokenMax) = true

function evaluate(dist::TokenMax, s1::AbstractString, s2::AbstractString, max_dist = 1.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    score = evaluate(dist.dist, s1, s2, max_dist)
    min_score = min(max_dist, score)
    unbase_scale = 0.95
    # if one string is much shorter than the other, use partial
    if length(s2) >= 1.5 * length(s1)
        partial_scale = length(s2) > (8 * length(s1)) ? 0.6 : 0.9
        score_partial = 1 - partial_scale * (1 - evaluate(Partial(dist.dist), s1, s2, 1 - (1 - max_dist) / partial_scale))
        min_score = min(max_dist, score_partial)
        score_sort = 1 - unbase_scale * partial_scale *
                (1 - evaluate(TokenSort(Partial(dist.dist)), s1, s2, 1 - (1 - max_dist) / (unbase_scale * partial_scale)))
        max_dist = min(max_dist, score_sort)
        score_set = 1 - unbase_scale * partial_scale *
                (1 - evaluate(TokenSet(Partial(dist.dist)), s1, s2, 1 - (1 - max_dist) / (unbase_scale * partial_scale)))
        return min(score, score_partial, score_sort, score_set)
    else
        score_sort = 1 - unbase_scale *
                (1 - evaluate(TokenSort(dist.dist), s1, s2, 1 - (1 - max_dist) / unbase_scale))
        max_dist = min(max_dist, score_sort)
        score_set = 1 - unbase_scale *
                (1 - evaluate(TokenSet(dist.dist), s1, s2, 1 - (1 - max_dist) / unbase_scale))
        return min(score, score_sort, score_set)
    end
end


text:
"""
	gmtsplit(cmd0::String="", arg1=nothing; kwargs...)

Reads a series of (x,y[,z]) records [or optionally (x,y,z,d,h)] and splits this into separate lists
of (x,y[,z]) series, such that each series has a nearly constant azimuth through the x,y plane.

Full option list at [`gmtsplit`]($(GMTdoc)gmtsplit.html)

Parameters
----------

- **A** | **azim_tol** :: [Type => Str | Array]

    Write out only those segments which are within +/- tolerance degrees of azimuth in heading,
    measured clockwise from North, [0 - 360].
    ($(GMTdoc)gmtsplit.html#a)
- **C** | **course_change** :: [Type => Number]

    Terminate a segment when a course change exceeding course_change degrees of heading is detected.
    ($(GMTdoc)gmtsplit.html#c)
- **D** | **min_dist** | **min_distance** :: [Type => Number]

    Do not write a segment out unless it is at least minimum_distance units long.
    ($(GMTdoc)gmtsplit.html#d)
- **F** | **filter** :: [Type => Str | Array]

    Filter the z values and/or the x,y values, assuming these are functions of d coordinate.
    xy_filter and z_filter are filter widths in distance units.
    ($(GMTdoc)gmtsplit.html#f)
- **Q** | **xyzdh** :: [Type => Str]

    Specify your desired output using any combination of xyzdh, in any order.
    ($(GMTdoc)gmtsplit.html#q)
- **S** | **dh** | **dist_head** :: [Type => Bool]

    Both d and h are supplied. In this case, input contains x,y,z,d,h. [Default expects (x,y,z) input,
    and d,h are computed from delta x, delta y.
    ($(GMTdoc)gmtsplit.html#s)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_bi)
- $(GMT.opt_bo)
- $(GMT.opt_di)
- $(GMT.opt_do)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_swap_xy)
"""
function gmtsplit(cmd0::String="", arg1=nothing; kwargs...)

    proggy = (GMTver > v"6.1.1") ? "gmtsplit" : "splitxyz"
	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic(proggy, cmd0, arg1)

	d = init_module(false, kwargs...)[1]		    # Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:V_params :bi :bo :di :do :e :f :g :h :i :yx])
	cmd  = parse_these_opts(cmd, d, [[:A :azim_tol], [:C :course_change], [:D :min_dist :min_distance], [:F :filter],
	                                 [:Q :xyzdh], [:S :dh :dist_head]])
	common_grd(d, cmd0, cmd, proggy * " ", arg1)	# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtsplit(arg1; kw...) = gmtsplit("", arg1; kw...)
const splitxyz = gmtsplit			# Alias


text:
using RecyclePermanent
using JLD2

include("../models/gkModel.jl")

@load "data/gkData.jld" ys θ
ys = ys[1:100]
model = makeGKModel(θ, ys)

N = 80 * model.n

println("Single-threaded")
rpio = RPIO{model.Tx}(N, model.n, 1)
@time simpleEstimate(model, rpio)
@time simpleEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimate(model, rpio)

println("Multi-threaded")
rpio = RPIO{model.Tx}(N, model.n, Threads.nthreads())
@time simpleEstimate(model, rpio)
@time simpleEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time biasedEstimate(model, rpio)
@time rpEstimate(model, rpio)
@time rpEstimate(model, rpio)


text:
# This file is a part of GERDADeepLearning.jl, licensed under the MIT License (MIT).

using StatsBase, LsqFit

export EfficiencyCurve
type EfficiencyCurve
  name :: Symbol # Tl_DEP, Bi_FEP, Tl_SEP, Tl_FEP, :Phy2v
  cut_values :: Vector{Float64} # in ascending order
  efficiencies :: Vector{Float64}
  std :: Vector{Float64} # zeros if not used
  # count_before :: Float64
  # counts_after :: Vector{Float64}

  EfficiencyCurve(name::Symbol, cut_values::Vector{Float64}) = new(name, cut_values, zeros(Float64, length(cut_values)), zeros(Float64, length(cut_values)))
  EfficiencyCurve(name :: Symbol, cut_values :: Vector{Float64}, efficiencies :: Vector{Float64}, std :: Vector{Float64}) = new(name, cut_values, efficiencies, std)
end

export EfficiencyCollection
type EfficiencyCollection
  curves :: Array{EfficiencyCurve}
  name :: AbstractString
  title :: AbstractString
end

Base.getindex(curve::EfficiencyCurve, range) = EfficiencyCurve(curve.name, curve.cut_values[range], curve.efficiencies[range], curve.std[range])

function Base.getindex(curves::EfficiencyCollection, name::Symbol)
  indices = find(curve -> curve.name == name, curves.curves)
  if length(indices) != 1
    throw(ArgumentError("Efficiency curve not contained: $name"))
  end
  return curves.curves[indices[1]]
end

Base.getindex(curves::EfficiencyCollection, range) = EfficiencyCollection([curve[range] for curve in curves.curves], curves.name, curves.title)

export save_effs
function save_effs(env::DLEnv, effs::EfficiencyCollection)
  dir = resolvepath(env, "models", "efficiencies-$(effs.name)")
  isdir(dir) || mkdir(dir)
  for curve in effs.curves
    save_eff_txt(curve, joinpath(dir, "$(curve.name).txt"))
  end
end

export load_effs
function load_effs(env::DLEnv, name::AbstractString)
  dir = resolvepath(env, "models", "efficiencies-$(name)")
  curves = EfficiencyCurve[]
  for filename in readdir(dir)
    cname = filename
    if endswith(cname, ".txt")
      cname = cname[1:end-4]
    end
    curve = load_eff_txt(joinpath(dir, filename), Symbol(cname))
    push!(curves, curve)
  end
  return sort_peaks!(EfficiencyCollection(curves, name, name))
end

function sort_peaks!(curves::EfficiencyCollection)
  order = [:Tl_DEP, :Bi_FEP, :Tl_SEP, :Tl_FEP]
  list = EfficiencyCurve[]
  for key in order
    push!(list, curves[key])
  end
  curves.curves = list
  return curves
end

export save_eff_txt
function save_eff_txt(eff::EfficiencyCurve, filename::AbstractString)
  if eff.std != nothing
    writedlm(filename, hcat(eff.cut_values, eff.efficiencies, eff.std))
  else
    writedlm(filename, hcat(eff.cut_values, eff.efficiencies))
  end
end

export load_effs_txt
function load_effs_txt(filenames::Vector{String}, names::Vector{Symbol}, name::AbstractString)
  curves = [load_eff_txt(filenames[i], names[i]) for i in 1:length(filenames)]
  return EfficiencyCollection(curves, name, name)
end

export load_eff_txt
function load_eff_txt(filename::AbstractString, name::Symbol)
  data = readdlm(filename)
  cut_values = data[:,1]

  efficiencies = data[:,2]
  # Efficiencies may be in percent
  if maximum(efficiencies) > 10
    efficiencies /= 100
  end

  # Std if available
  if size(data, 2) > 2
    std = data[:,3]
  else
    std = zeros(Float64, length(cut_values))
  end
  # Order cut values so they are ascending
  if cut_values[1] > cut_values[end]
    cut_values = cut_values[end:-1:1]
    efficiencies = efficiencies[end:-1:1]
    std = std[end:-1:1]
  end
  return EfficiencyCurve(name, cut_values, efficiencies, std)
end

export random_guessing
function random_guessing()
  sig = EfficiencyCurve(:Tl_DEP, [0.0, 1], [0.0, 1], [0.0, 0])
  bkg = EfficiencyCurve(:Bi_FEP, [0.0, 1], [0.0, 1], [0.0, 0])
  return EfficiencyCollection([sig, bkg], "rand", "Random guessing")
end

export get_peak_efficiencies
function get_peak_efficiencies(env::DLEnv, lib::EventLibrary, name::AbstractString)
  return get_peak_efficiencies(env, lib, name, name)
end
function get_peak_efficiencies(env::DLEnv, lib::EventLibrary, name::AbstractString, title::AbstractString)
  if isdir(resolvepath(env, "models", "efficiencies-$(name)"))
    effs = load_effs(env, name)
  else
    effs = peak_efficiency_curves(lib, name)
    effs = cut_efficiency_fluctuations(effs)
    save_effs(env, effs)
  end
  effs.title = title
  return effs
end

export peak_efficiency_curves
function peak_efficiency_curves(lib::EventLibrary, name::AbstractString; cut_values=100, peak_names::Vector{Symbol}=[:Tl_DEP, :Bi_FEP, :Tl_SEP, :Tl_FEP], peak_centers::Vector{Float64}=[1592.5, 1620.7, 2103.5, 2614.5], half_window::Float64=10.0, bin_width::Float64=0.5, psd_key=:psd)
  if isa(cut_values, Integer)
    cut_values = equal_counts_cut_values(lib, cut_values; psd_key=psd_key)
  else
        cut_values = collect(cut_values)
  end
  curves = [peak_efficiency_curve(peak_names[i], lib, peak_centers[i], cut_values; half_window=half_window, bin_width=bin_width, psd_key=psd_key) for i in 1:length(peak_centers)]
  return EfficiencyCollection(curves, name, name)
end

export peak_efficiency_curve
"""
Cut values must be in ascending order.
"""
function peak_efficiency_curve(curve_name::Symbol, lib::EventLibrary, peak_center::Float64, cut_values::Array{Float64}; half_window::Float64=10.0, bin_width::Float64=0.5, psd_key=:psd)
  curve = EfficiencyCurve(curve_name, cut_values)

  start_energy = peak_center-half_window-bin_width
  end_energy = peak_center+half_window+bin_width
  nbins = Int64(round((end_energy-start_energy)/bin_width))

  bin_edges = linspace(start_energy, end_energy, nbins)
  bin_width = bin_edges[2] - bin_edges[1]

  lib = filter(lib, :E, E -> (E>start_energy) && (E<end_energy))

  energy_hist = fit(Histogram, lib[:E], bin_edges)
  amplitude, error = fit_peak(bin_edges, energy_hist.weights, peak_center, half_window)
  amplitude_before = max(0.1, amplitude)

  for (i,cut_value) in enumerate(cut_values)
    filter!(lib, psd_key, psd -> psd > cut_value)
    energy_hist = fit(Histogram, lib[:E], bin_edges)
    amplitude, error = fit_peak(bin_edges, energy_hist.weights, peak_center, half_window)
    amplitude = max(0, amplitude)
    curve.efficiencies[i] = amplitude / amplitude_before
    curve.std[i] = sqrt(amplitude * (1- min(1,amplitude/amplitude_before))) / amplitude_before
  end

  return curve
end

"""
Hitogram fit: Fit a gaussian with variable width + simple background model to every peak in peak_centers. The background is taken from the half_window.
Points outside half_window+margin are thrown away.

Returns the best fit amplitude with respective error of the peak height in absolut event counts.
"""
function fit_peak(total_x_axis, total_energy_hist::Array, peak_center, half_window)
  # Define peak function (Gaussian)
  # Parameters: Constant, Linear, Gauss amplitude, Gauss sigma
  peak(e, p) = p[1] + p[2]*(e-peak_center) + p[3] * exp(- (e-peak_center).^2 ./ (2*p[4]^2))

  indices = find(x -> (x >= peak_center-half_window) && (x <= peak_center+half_window), total_x_axis)
  x_axis = total_x_axis[indices]
  energy_hist = total_energy_hist[indices]
  energy_hist = convert(Array{Float64}, energy_hist)

  mean_count = mean(energy_hist)

  fit_result = curve_fit(peak, x_axis, energy_hist, [mean_count*0.8, 0.0, mean_count*0.2, 3.0])
  fit_amplitude = fit_result.param[3]
  try
    fit_error = estimate_errors(fit_result)
    fit_amplitude_error = fit_error[3]
    return fit_amplitude, fit_amplitude_error
  catch exc
    info("Fit Errors could not be calculated $exc")
    return fit_amplitude, Inf
  end
end

export equal_counts_cut_values
"""
Returns an array of count cut values in descending order the first of which is max_cut_value (all events excluded). The cut values are derived from the classifier distribution so that the total event count decreases with constant speed.
"""
function equal_counts_cut_values(events::EventLibrary, count; psd_key=:psd)
  psd = events[psd_key]
  bin_edges = linspace(Float64(minimum(psd)), maximum(psd), count*10)
  psd_hist = fit(Histogram, psd, bin_edges)
  psd_counts = psd_hist.weights

  target_sizes = linspace(0, eventcount(events), count)
  cut_values = [bin_edges[_integration_index(psd_counts, target_sizes[i])] for i in 1:length(target_sizes)]
  return cut_values
end
function _integration_index(hist, target)
  sum = 0
  for (i,val) in enumerate(hist)
    if sum >= target return i end
    sum += val
  end
  return length(hist)
end

export cut_efficiency_fluctuations
function cut_efficiency_fluctuations(effs::EfficiencyCollection)
  for i in 1:length(effs.curves[1].cut_values)
    for curve in effs.curves
      if curve.std[i] > curve.efficiencies[i]
        info("Cut efficiencies at index $i due to high statistical error.")
        return effs[1:i]
      end
    end
  end
  return effs
end

export background_rejection_at
function background_rejection_at(signal_eff::AbstractFloat, effs::EfficiencyCollection; signal_peak=:Tl_DEP, bkg_peak=:Bi_FEP)
  sig = effs[signal_peak]
  bkg = effs[bkg_peak]

  for i in length(sig.cut_values):-1:1
    if sig.efficiencies[i] > signal_eff
      return i, sig.cut_values[i], sig.efficiencies[i], 1 - bkg.efficiencies[i]
    end
  end
  info("Signal efficiency $signal_eff is never reached!")
end

export background_rejection_std_at
function background_rejection_std_at(signal_eff::AbstractFloat, effs::EfficiencyCollection; signal_peak=:Tl_DEP, bkg_peak=:Bi_FEP)
    sig = effs[signal_peak]
    bkg = effs[bkg_peak]
    i, cut, eff, rej = background_rejection_at(signal_eff, effs; signal_peak=signal_peak, bkg_peak=bkg_peak)
    if (sig.std[i] == 0) || (bkg.std[i] == 0)
        return -1
    else
        return sqrt(sig.std[i]^2 + bkg.std[i]^2)
    end
end

export equal_rejection_and_efficiency
function equal_rejection_and_efficiency(effs::EfficiencyCollection; signal_peak=:Tl_DEP, bkg_peak=:Bi_FEP)
  sig = effs[signal_peak]
  bkg = effs[bkg_peak]

  for i in length(sig.cut_values):-1:1
    if sig.efficiencies[i] > 1-bkg.efficiencies[i]
      return i, sig.cut_values[i], sig.efficiencies[i], 1 - bkg.efficiencies[i]
    end
  end
  info("Signal efficiency is never better than background rejection!")
end

export find_cut_value
function find_cut_value(effs::EfficiencyCollection, cut_value)
  if isa(cut_value, Real)
    return cut_value
  elseif cut_value == "90% DEP"
    cut_value_result = background_rejection_at(0.9, effs)[2]
  elseif cut_value == "Equal"
    cut_value_result = equal_rejection_and_efficiency(effs)[2]
  else
    info("Illegal cut_value argument: $cut_value")
    return cut_value
  end
  info("Determined cut value $cut_value = $cut_value_result")
  cut_value_result
end

export curve
curve(eff::EfficiencyCurve) = (eff.cut_values, eff.efficiencies)

export roc_curve
""" Complete: :left, :right, true, false
"""
function roc_curve(effs::EfficiencyCollection; signal_effs=linspace(0.01, 0.99, 400).^0.1, signal_peak=:Tl_DEP, bkg_peak=:Bi_FEP, complete=true)

  real_signal_effs = Float64[]
  real_bkg_rej = Float64[]

  if complete == true || complete == :left
    push!(real_signal_effs, 0)
    push!(real_bkg_rej, 1)
  end

  for signal_eff in signal_effs
    values = background_rejection_at(signal_eff, effs; signal_peak=signal_peak, bkg_peak=bkg_peak)
    if values != nothing
      push!(real_signal_effs, values[3])
      push!(real_bkg_rej, values[4])
    else break
    end
  end

  if complete == true || complete == :right
    push!(real_signal_effs, 1)
    push!(real_bkg_rej, 0)
  end

  return real_signal_effs, real_bkg_rej
end

export load_current_effs
function load_current_effs(detector_name::AbstractString)
  target_peaks = [:Tl_DEP, :Bi_FEP, :Tl_SEP, :Tl_FEP]
  current_eff = nothing
  try
    if !contains(detector_name, "GD") # coax
      detector_id = get_detector_index(detector_name)
      files = ["/home/iwsatlas1/pholl/workspace/Andreas Efficiencies/$det-$(lpad(detector_id,2,0))-$pname.txt" for pname in ["208Tl DEP", "212Bi FEP", "208Tl SEP", "208Tl FEP"]]
      return load_effs_txt(files, target_peaks, "Current implementation")
    else # BEGe
      files = ["/home/iwsatlas1/pholl/workspace/AoE Efficiencies/efficiency_$(pname)_$(detector_name).txt" for pname in ["DEP", "FEP", "SEP", "FEP2"]]
      return load_effs_txt(files, target_peaks, "A/E")
    end
  catch err
    info("Could not load current efficiencies for detector $detector_name.")
  end
end


text:
function is_known_state(memory, history, max_index)
  for i = 1:max_index
    if history[i, :] == memory
      return true
    end
  end

  return false
end

function reallocate!(memory)
  # Find where the most blocks are hiding and remove them
  max_blocks = maximum(memory)
  index = minimum(findin(memory, max_blocks));
  memory[index] = 0

  # Redistribute the blocks on the rest
  for i = 1:max_blocks
    len = length(memory)
    memory[mod(index + i - 1, len) + 1] += 1
  end
end

function solve_part_1!(memory)
  len = length(memory)
  history = zeros(10000, len);

  step_count = 0
  while !is_known_state(memory, history, step_count)
    # Store current state in memory
    history[step_count + 1, :] = memory
    step_count += 1

    # Reallocate memory
    reallocate!(memory)
  end

  return step_count
end

function solve_part_2(memory)
  # Copy the original state of the memory
  original = copy(memory)

  # Reallocate the memory
  step_count = 0
  while true
    reallocate!(memory)
    step_count += 1

    # Until it becomes the same as the original
    memory == original && break
  end

  return step_count
end

input = [4, 1, 15, 12, 0, 9, 9, 5, 5, 8, 7, 3, 14, 5, 12, 3];
input_test_1 = [0, 2, 7, 0]
input_test_2 = [2, 4, 1, 2]

@printf("Steps for part 1 test: %d\n", solve_part_1!(input_test_1))
@printf("Steps for part 1: %d\n", solve_part_1!(input))

@printf("Steps for part 2 test: %d\n", solve_part_2(input_test_2))
@printf("Steps for part 2: %d\n", solve_part_2(input))


text:
using Zygote
using Random
using Flux
using UnPack
using Functors

include("operators.jl")
include("diffrules.jl")
Random.seed!(1)

mutable struct EquivConv
    in::Props
    out::Props
    paths::AbstractVector{NamedTuple}
    pathsmap::AbstractVector
    dx::Real
    rmax::Real
end
@functor EquivConv

function Flux.trainable(m::EquivConv)
    [x.op.radfunc for x in m.paths]
end
function sep(x, nranks)
    [body]
end
"""
    function EquivConv(
        name,
        in.ranks,
        out.ranks;
        dims = 3,
        dx = 1.0,
        rmax = 1.0,
        rank_max = max(1, in.ranks..., out.ranks...),
        σ = identity
    )
"""
function EquivConv(
    in_::Props,
    out::Props,
    dx,
    rmax;
    dims = 3,
    rank_max = max(1, length(in_.nranks), length(out.nranks)) - 1,
)
    Zygote.ignore() do
        grid = Grid(dx, rmax; rank_max)
        nin = length(in_.ranks)
        nout = length(out.ranks)

        paths = []
        pathsmap = [[] for i = 1:nout]
        # iterate all tensor product paths
        # if isempty(paths)

        for lf in (rmax == 0 ? [0] : 0:rank_max)
            for (i, li, si) in zip(1:nin, in_.ranks, in_.parities)
                loutrange = abs(li - lf):min(li + lf, rank_max)
                parity = si * (-1)^lf
                for (o, lo) in enumerate(out.ranks)
                    if (lo in loutrange) && out.parities[o] == parity
                        ranks = (li, lf, lo)
                        op = LinearOperator(:neural; grid, ranks, rmax)
                        path = (; i, o, ranks, op)
                        push!(paths, path)
                    end
                    # end
                end
            end
        end

        for (i, path) in enumerate(paths)
            push!(pathsmap[path.o], i)
        end

        return EquivConv(in_, out, paths, pathsmap, dx, rmax)
    end
end

"""
"""
function (f::EquivConv)(X::AbstractArray; remake = true)
    in = out = paths = pathsmap = dx = rmax = 0
    Zygote.ignore() do
        @unpack in, out, paths, pathsmap, dx, rmax = f
    end
    pathfunc = i -> paths[i].op(get(X, in, paths[i].i))
    if remake
        for x in paths
            remake!(x.op)
        end
    end

    res = cat(
        [
            sum([pathfunc(i) for i in pathsmap[o]]) for o in eachindex(pathsmap)
        ]...,
        dims = 4,
    )

end

struct EquivProd
    in::Props
    out::Props
    paths::AbstractVector{NamedTuple}
    spectra::Bool
end
function EquivProd(
    in::Props;
    rank_max = max(1, length(in.ranks)),
    spectra = false,
)
    Zygote.ignore() do
        paths = []
        nranks_out = zeros(Int, rank_max + 1)
        for (i1, li1, si1) in zip(1:length(in.ranks), in.ranks, in.parities)
            for (i2, li2, si2) in
                zip(0:i1, [0, in.ranks[1:i1]...], [1, in.parities[1:i1]...])
                loutrange = abs(li1 - li2):min(li1 + li2, rank_max)
                so = si1 * si2
                for lo in loutrange
                    prod = FieldProd(li1, li2, lo)
                    if prod !== nothing && (i1 != i2 || prod != prod111)
                        nranks_out[lo+1] += 1
                        path = (; i1, i2, lo, so, prod)
                        push!(paths, path)
                    end
                end
            end
        end
        if spectra
            nranks_out = (sum(nranks_out),)
        end
        out = Props(nranks_out)
        sort!(paths, by = x -> x.lo)
        EquivProd(in, out, paths, spectra)
    end
end
function (f::EquivProd)(X::AbstractArray)
    in, out, paths, spectra = f.in, f.out, f.paths, f.spectra
    res = [
        path.i2 === 0 ? get(X, in, path.i1) :
            path.prod(get(X, in, path.i1), get(X, in, path.i2))
        for path in paths
    ]
    if spectra
        res = [fieldnorm(x) for x in res]
        # res = fieldnorm.(res)
    end
    cat(res..., dims = 4)
end
struct LocalDense
    m::Dense
end
@functor LocalDense
function LocalDense(args...)
    Random.seed!(1)
    LocalDense(Dense(args...))
end

function Flux.trainable(m::LocalDense)
    return [m.m]
    # return vcat(m.W, [m.b])
end

function (m::LocalDense)(x::AbstractArray)
    @unpack W, b, σ = m.m
    out, in = size(W)
    res =
        σ.(cat(
            [
                b[i] .+ sum([W[i, j] * x[:, :, :, j:j] for j = 1:in]) for i = 1:out
            ]...,
            dims = 4,
        ))
end
struct EquivDense
    in::Props
    out::Props
    W::AbstractVector
    b::AbstractVector
    σ::Any
end
function Flux.trainable(m::EquivDense)
    return [m.W, m.b]
    # return vcat(m.W, [m.b])
end

function EquivDense(in, out; σ = identity)
    Zygote.ignore() do
        W = [ones(b, a) / a for (a, b) in zip(in.nranks, out.nranks)]
        b = zeros(out.nranks[1])
        EquivDense(in, out, W, b, σ)
    end
end
function (f::EquivDense)(X::AbstractArray)
    @unpack in, out, W, b, σ = f
    res = [
        [
            sum([
                wi[j] * X[:, :, :, in.grouped_slices[a][j]]
                for j in eachindex(in.grouped_slices[a])
            ]) for wi in eachrow(W[a])
        ] for a in eachindex(in.grouped_slices)
    ]
    # res = [w * [X[:,:,:,s] for s in S] for (w, S) in zip(W, in.grouped_slices)]
    res1 = [σ.(res[1][i] .+ b[i]) for i in eachindex(b)]
    # res1 = [σ.(x.+bi) for (x,bi) in zip(res[1], b)]
    # res[1] = res1
    res = [res1]
    cat(vcat(res...)..., dims = 4)
end

function nae(yhat, y; sumy = sum(abs.(y)))
    if sumy == 0
        error()
    end
    sum(abs.(yhat .- y)) / sumy
end

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
module TrustRegionMethods

export trust_region_solver, TrustRegionResult, ForwardDiff_wrapper, SolverStoppingCriterion,
    Dogleg, GeneralizedEigenSolver

using ArgCheck: @argcheck
import DiffResults
using DocStringExtensions: FIELDS, FUNCTIONNAME, SIGNATURES, TYPEDEF
import ForwardDiff
using KrylovKit: eigsolve
using LinearAlgebra: diag, Diagonal, dot, I, issuccess, lu, norm, Symmetric, UniformScaling
using UnPack: @unpack

####
#### building blocks
####

"""
$(TYPEDEF)

A *linear model* for system of nonlinear equations, relative to the origin.

The L2 norm induces the objective function

```math
m(p) = p' J' r + 1/2 p' J' J p
```
approximating some `f(x) = \\| r(x) \\|^2_2` as

```math
1/2 \\| f(x + p) \\|_2^2 \approx 1/2 \\| r + J p \\|_2^2 = \\| r \\|_2^2 + p'J'r + 1/2 p' J' J p
```
"""
struct ResidualModel{TR,TJ}
    r::TR
    J::TJ
    function ResidualModel(r::TR, J::TJ) where {TR,TJ}
        @argcheck all(isfinite, r) "Non-finite residuals $(r)."
        @argcheck all(isfinite, J) "Non-finite Jacobian $(J)."
        n = length(r)
        @argcheck size(J) ≡ (n, n) "Non-conformable residual and Jacobian."
        new{TR,TJ}(r, J)
    end
end

"""
$(TYPEDEF)

A *quadratic model* for a minimization problem, relative to the origin, with the objective
function

```math
m(p) = p' g + 1/2 g' A g
```
"""
struct MinimizationModel{TG,TA}
    g::TG
    A::TA
    # FIXME constructor enforces symmetry of `A`, document, but do we actually rely on it?
end

"""
$(SIGNATURES)

Convert a model for a residual to a minimization model using the L2 norm.
"""
function MinimizationModel(model::ResidualModel)
    @unpack r, J = model
    MinimizationModel(J' * r, Symmetric(J' * J))
end

"""
$(SIGNATURES)

Find the Cauchy point (the minimum of the linear part, subject to ``\\| p \\| ≤ Δ``) of the problem.

Return three values:

1. the Cauchy point vector `pC`,

2. the (Euclidean) *norm* of `pC`

3. a boolean indicating whether the constraint was binding.

Caller guarantees non-zero gradient.
"""
function cauchy_point(Δ::Real, model::ResidualModel)
    @unpack r, J = model
    g = J' * r
    q = g' * (J' * J) * g
    g_norm = norm(g, 2)
    @argcheck g_norm > 0
    τ = if q ≤ 0              # practically 0 (semi-definite form) but allow for float error
        one(q)
    else
        min(one(q), g_norm^3 / (Δ * q))
    end
    (-τ * Δ / g_norm) .* g, τ * Δ, τ ≥ 1
end

"""
$(SIGNATURES)

Finds `\tau` such that

```math
\\| p_C + τ D \\|_2 = Δ
```

`pC_norm` is `\\| p_C \\|_2^2`, and can be provided when available.

Caller guarantees that

1. `pC_norm ≤ Δ`,

2. `norm(pC + D, 2) ≥ Δ`.

3. `dot(pC, D) ≥ 0`.

These ensure a meaningful solution `0 ≤ τ ≤ 1`. None of them are checked explicitly.
"""
function dogleg_boundary(Δ, D, pC, pC_norm = norm(pC, 2))
    a = sum(abs2, D)
    b = 2 * dot(D, pC)
    c = abs2(pC_norm) - abs2(Δ)
    (-b + √(abs2(b) - 4 * a * c)) / (2 * a)
end

"""
$(SIGNATURES)

Return `(F, is_valid_F)` where

1. `F` is a form (usually a factorization) of the argument that supports `F \\ r` for
   vectors `r`,

2. `is_valid_F` is a boolean indicating whether `F` can be used in this form (eg `false` for
   a singular factorization).
"""
function _factorize(J::AbstractMatrix)
    LU = lu(J; check = false)
    LU, issuccess(LU)
end

function _factorize(J::Diagonal)
    J, all(x -> x > 0, diag(J))
end

"""
$(SIGNATURES)

Calculate the unconstrained optimum of the model, return its norm as the second value.

When the second value is *infinite*, the unconstrained optimum should not be used as this
indicates a singular problem.
"""
function unconstrained_optimum(model::ResidualModel)
    @unpack r, J = model
    F, is_valid_F = _factorize(J)
    if is_valid_F
        pU = -(F \ r)
        pU, norm(pU, 2)
    else
        ∞ = oftype(one(eltype(r)) / one(eltype(F)), Inf)
        fill(∞, length(r)), ∞
    end
end

####
#### solvers
####

"""
`$(FUNCTIONNAME)(method, Δ, model)`

Optimize the `model` with the given constraint `Δ` using `method`. Returns the following
values:

- the optimum,
- the norm of the optimum,
- a boolean indicating whether the constraint binds.
"""
function solve_model end

"""
The dogleg method.

A simple and efficient heuristic for solving local trust region problems. Does not
necessarily find the optimum, but improves on the Cauchy point.
"""
struct Dogleg end

"""
$(SIGNATURES)

Implementats the branch of the dogleg method where it is already determined that the
unconstrained optimum is lies outside the `Δ` ball. Returns the same kind of results as
[`solve_model`](@ref).
"""
function dogleg_implementation(Δ, model, pU, pU_norm)
    pC, pC_norm, on_boundary = cauchy_point(Δ, model)
    if on_boundary || isinf(pU_norm)
        pC, pC_norm, on_boundary
    else
        D = pU .- pC
        τ = dogleg_boundary(Δ, D, pC, pC_norm)
        pC .+ D .* τ, Δ, true
    end
end

function solve_model(::Dogleg, Δ, model::ResidualModel)
    pU, pU_norm = unconstrained_optimum(model)
    if pU_norm ≤ Δ
        pU, pU_norm, false
    else
        dogleg_implementation(Δ, model, pU, pU_norm)
    end
end

####
#### Generalized eigenvalue solver (Adachi et al 2017).
####
#### NOTE:
####
#### - notation mostly follows the paper
####
#### - implementation has ellipsoidal norm in some parts, for future generalizations,
####   currently unused
####
#### - “hard case” is not implemented yet, falls back to dogleg

"""
Generalized eigenvalue solver (Adachi et al 2017).
"""
struct GeneralizedEigenSolver end

"""
$(SIGNATURES)

Kernel for the generalized eigenvalue solver (methods specialize on the ellipsoidal norm
matrix `B`). Solves the `M̃(λ)` pencil in Adachi et al (2017), eg Algorithm 5.1.

Returns

- `λ`, the largest eigenvalue,

- the gap to the next eigenvalue (see note below)

- `y1` and `y2`, the two parts of the corresponding generalized eigenvalue,

All values are real, methods are type stable.

When theorerical assumptions are violated, `gap` will be non-finite and a debug statement is
emitted. No other values should be used in this case.
"""
function ges_kernel(Δ, model::MinimizationModel, B::UniformScaling)
    @unpack g, A = model
    n = length(g)
    G = ((g * g') ./ abs2(Δ))
    M = [-A G; B -A]
    λs, vs, info = eigsolve(M, 2, :LR)
    @argcheck info.converged ≥ 2 "Eigensolver did not converge."
    T = promote_type(Float64, eltype(M))
    ϵ = √eps(T)
    is_practically_real(z) = imag(z) ≤ n*(ϵ*abs(z) + ϵ)
    λ1, λ2 = λs
    v1 = first(vs)
    if is_practically_real(λ1) && all(is_practically_real, v1)
        λ = real(λ1)::T
        gap = abs(λ - λ2)::T
        y = real.(v1)::Vector{T}
        y1 = y[1:n]
        y2 = y[(n + 1):end]
        λ, gap, y1, y2
    else
        @debug "rightmost eigenvalue not real" M λs vs
        ∞ = T(Inf)
        y∞ = fill(∞, n)
        ∞, ∞, y∞, y∞
    end
end

"""
$(SIGNATURES)

Ellipsoidal norm ``\\| x \\|_B = x'Bx``.
"""
ellipsoidal_norm(x, ::UniformScaling) = norm(x, 2)

function solve_model(::GeneralizedEigenSolver, Δ, model::ResidualModel)
    B = I                       # FIXME: we hardcode Euclidean norm for now
    pU, pU_norm = unconstrained_optimum(model)
    if pU_norm < Δ
        pU, pU_norm, false
    else
        model′ = MinimizationModel(model)
        λ, gap, y1, y2 = ges_kernel(Δ, model′, B)
        τ = √(eps(typeof(λ)) / gap)
        if isfinite(gap) && norm(y1, 2) > τ
            # “easy” case, we can generate a candidate
            @unpack g = model′
            p = (-sign(dot(g, y2)) * Δ / ellipsoidal_norm(y1, B)) .* y1
            p, Δ, true
        else
            # FIXME hard case, this needs to be implemented
            # here we bail and fall back to dogleg
            # NOTE this is also the fallback for eigensolver gone wrong
            @debug "hard case not implemented, falling back to dogleg"
            dogleg_implementation(Δ, model, pU, pU_norm)
        end
    end
end

####
#### trust region steps
####

struct TrustRegionParameters{T}
    η::T
    Δ̄::T
    function TrustRegionParameters(η::T, Δ̄::T) where {T <: Real}
        @argcheck 0 < η < 0.25
        @argcheck Δ̄ > 0
        new{T}(η, Δ̄)
    end
end

TrustRegionParameters(η, Δ̄) = TrustRegionParameterS(promote(η, Δ̄)...)

TrustRegionParameters(; η = 0.125, Δ̄ = Inf) = TrustRegionParameters(η, Δ̄)

"""
$(SIGNATURES)

Ratio between predicted (using `model`, at `p`) and actual reduction (taken from the
residual value `r′`). Will return an arbitrary negative number for infeasible coordinates.
"""
function reduction_ratio(model::ResidualModel, p, r′)
    @unpack r, J = model
    r2 = sum(abs2, r)
    r′2 = sum(abs2, r′)
    ρ = (r2 - r′2) / (r2 - sum(abs2, r .+ J * p))
    isfinite(ρ) ? ρ : -one(ρ)
end

"""
$(SIGNATURES)

Checks residual and Jacobian and throws an appropriate error if they are not both finite.

Internal, each evaluation of `f` should call this.
"""
function _check_residual_Jacobian(x, fx)
    finite_value = all(isfinite, fx.residual)
    finite_Jacobian = all(isfinite, fx.Jacobian)
    if !finite_value || !finite_Jacobian
        if finite_value == finite_Jacobian == false
            msg = "residual and Jacobian"
        else
            msg = !finite_value ? "residual" : "Jacobian"
        end
        DomainError(x, "Non-finite values in $(msg)")
    end
end

"""
$(SIGNATURES)

Take a trust region step using `local_method`.

`f` is the function that returns the residual and the Jacobian (see
[`trust_region_solver`](@ref)).

`Δ` is the trust region radius, `x` is the position, `fx = f(x)`. Caller ensures that the
latter is feasible.
"""
function trust_region_step(parameters::TrustRegionParameters, local_method, f, Δ, x, fx)
    @unpack η, Δ̄ = parameters
    model = ResidualModel(fx.residual, fx.Jacobian)
    p, p_norm, on_boundary = solve_model(local_method, Δ, model)
    x′ = x .+ p
    fx′ = f(x′)
    _check_residual_Jacobian(x′, fx′)
    ρ = reduction_ratio(model, p, fx′.residual)
    Δ′ =
        if ρ < 0.25
            p_norm / 4
        elseif ρ > 0.75 && on_boundary
            min(2 * Δ, Δ̄)
        else
            Δ
        end
    if ρ ≥ η
        Δ′, x′, fx′            # reasonable reduction, use new position
    else
        Δ′, x, fx              # very small reduction, try again from original
    end
end

struct SolverStoppingCriterion{T <: Real}
    residual_norm_tolerance::T
    function SolverStoppingCriterion(residual_norm_tolerance::T) where {T}
        @argcheck residual_norm_tolerance > 0
        new{T}(residual_norm_tolerance)
    end
end

function SolverStoppingCriterion(; residual_norm_tolerance = √eps())
    SolverStoppingCriterion(residual_norm_tolerance)
end

function check_stopping_criterion(nsc::SolverStoppingCriterion, fx)
    r_norm = norm(fx.residual, 2)
    converged = r_norm ≤ nsc.residual_norm_tolerance
    converged, (converged = converged , residual_norm = r_norm)
end

"""
$(TYPEDEF)

# Fields

$(FIELDS)
"""
struct TrustRegionResult{T<:Real,TX<:AbstractVector{T},TFX}
    "The final trust region radius."
    Δ::T
    "The last value (the root only when converged)."
    x::TX
    "`f(x)` at the last `x`."
    fx::TFX
    "The Euclidean norm of the residual at `x`."
    residual_norm::T
    "A boolean indicating convergence."
    converged::Bool
    "Number of iterations (≈ number of function evaluations)."
    iterations::Int
end

function TrustRegionResult(Δ::T1, x::AbstractVector{T2}, fx, residual_norm::T3, converged,
                           iterations) where {T1 <: Real, T2 <: Real, T3 <: Real}
    T = promote_type(T1, T2, T3)
    TrustRegionResult(T(Δ), T.(x), fx, T(residual_norm), converged, iterations)
end

function Base.show(io::IO, ::MIME"text/plain", trr::TrustRegionResult)
    @unpack Δ, x, fx, residual_norm, converged, iterations = trr
    _sig(x) = round(x; sigdigits = 3)
    print(io, "Nonlinear solver using trust region method ")
    if converged
        printstyled(io, "converged"; color = :green)
    else
        printstyled(io, "didn't converge"; color = :red)
    end
    println(io, " after $(iterations) steps")
    print(io, "  with ")
    printstyled(io, "‖x‖₂ = $(_sig(residual_norm)), Δ = $(_sig(Δ))\n"; color = :blue)
    println(io, "  x = ", _sig.(x))
    println(io, "  r = ", _sig.(fx.residual))
end

"""
$(SIGNATURES)

Solve `f ≈ 0` using trust region methods, starting from `x`.

`f` is a callable (function) that

1. takes vectors real numbers of the same length as `x`,

2. returns a value with properties `residual` and `Jacobian`, containing a vector and a
   conformable matrix, with finite values for both or a non-finite residual (for infeasible
   points; Jacobian is ignored). A `NamedTuple` can be used for this, but any structure with
   these properties will be accepted. Importantly, this is treated as a single object and
   can be used to return extra information (a “payload”), which will be ignored by this
   function.

Returns a [`TrustRegionResult`](@ref) object.
"""
function trust_region_solver(f, x;
                             parameters = TrustRegionParameters(),
                             local_method = Dogleg(),
                             stopping_criterion = SolverStoppingCriterion(),
                             maximum_iterations = 500,
                             Δ = 1.0)
    @argcheck Δ > 0
    fx = f(x)
    _check_residual_Jacobian(x, fx)
    iterations = 1
    while true
        Δ, x, fx = trust_region_step(parameters, local_method, f, Δ, x, fx)
        iterations += 1
        do_stop, convergence_statistics = check_stopping_criterion(stopping_criterion, fx)
        reached_max_iter = iterations ≥ maximum_iterations
        if do_stop || reached_max_iter
            @unpack converged, residual_norm = convergence_statistics
            return TrustRegionResult(Δ, x, fx, residual_norm, converged, iterations)
        end
    end
end

####
#### AD shims
####

###
### ForwardDiff
###

"""
$(TYPEDEF)

A buffer for wrapping a (nonlinear) mapping `f` to calculate the Jacobian with
`ForwardDiff.jacobian`.
"""
struct ForwardDiffBuffer{TX,TF,TR,TC}
    "Buffer for inputs."
    x::TX
    "A function that maps vectors to residual vectors."
    f::TF
    "`DiffResults` buffer for Jacobians."
    result::TR
    "Gradient configuration."
    cfg::TC
end

"""
$(SIGNATURES)

Wrap an ``ℝⁿ`` function `f` in a callable that can be used in [`trust_region_solver`](@ref)
directly. Remaining parameters are passed on to `ForwardDiff.JacobianConfig`, and can be
used to set eg chunk size.

Non-finite residuals will be treated as infeasible.

```jldoctest
julia> f(x) = [1.0 2; 3 4] * x - ones(2)
f (generic function with 1 method)

julia> ff = ForwardDiff_wrapper(f, 2)
AD wrapper via ForwardDiff for ℝⁿ→ℝⁿ function, n = 2

julia> ff(ones(2))
(r = [2.0, 6.0], J = [1.0 2.0; 3.0 4.0])

julia> trust_region_solver(ff, [100.0, 100.0])
Nonlinear solver using trust region method converged after 9 steps
  with ‖x‖₂ = 3.97e-15, Δ = 128.0
  x = [-1.0, 1.0]
  r = [-1.78e-15, 3.55e-15]
```
"""
function ForwardDiff_wrapper(f, n, jacobian_config...)
    x = zeros(n)
    result = DiffResults.JacobianResult(x)
    cfg = ForwardDiff.JacobianConfig(f, x, jacobian_config...)
    ForwardDiffBuffer(x, f, result, cfg)
end

function Base.show(io::IO, ::MIME"text/plain", fdb::ForwardDiffBuffer)
    n = length(DiffResults.value(fdb.result))
    print(io, "AD wrapper via ForwardDiff for ℝⁿ→ℝⁿ function, n = $(n)")
end

function (fdb::ForwardDiffBuffer)(y)
    @unpack x, f, result, cfg = fdb
    # copy to our own buffer to work with types other than Float64
    ForwardDiff.jacobian!(result, f, copy!(x, y), cfg)
    residual = copy(DiffResults.value(result))
    (residual = residual, Jacobian = copy(DiffResults.jacobian(result)))
end

end # module