text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""NetworkPolicySpec provides the specification of a NetworkPolicy

    IoK8sApiNetworkingV1NetworkPolicySpec(;
        egress=nothing,
        ingress=nothing,
        podSelector=nothing,
        policyTypes=nothing,
    )

    - egress::Vector{IoK8sApiNetworkingV1NetworkPolicyEgressRule} : List of egress rules to be applied to the selected pods. Outgoing traffic is allowed if there are no NetworkPolicies selecting the pod (and cluster policy otherwise allows the traffic), OR if the traffic matches at least one egress rule across all of the NetworkPolicy objects whose podSelector matches the pod. If this field is empty then this NetworkPolicy limits all outgoing traffic (and serves solely to ensure that the pods it selects are isolated by default). This field is beta-level in 1.8
    - ingress::Vector{IoK8sApiNetworkingV1NetworkPolicyIngressRule} : List of ingress rules to be applied to the selected pods. Traffic is allowed to a pod if there are no NetworkPolicies selecting the pod (and cluster policy otherwise allows the traffic), OR if the traffic source is the pod's local node, OR if the traffic matches at least one ingress rule across all of the NetworkPolicy objects whose podSelector matches the pod. If this field is empty then this NetworkPolicy does not allow any traffic (and serves solely to ensure that the pods it selects are isolated by default)
    - podSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : Selects the pods to which this NetworkPolicy object applies. The array of ingress rules is applied to any pods selected by this field. Multiple network policies can select the same set of pods. In this case, the ingress rules for each are combined additively. This field is NOT optional and follows standard label selector semantics. An empty podSelector matches all pods in this namespace.
    - policyTypes::Vector{String} : List of rule types that the NetworkPolicy relates to. Valid options are \"Ingress\", \"Egress\", or \"Ingress,Egress\". If this field is not specified, it will default based on the existence of Ingress or Egress rules; policies that contain an Egress section are assumed to affect Egress, and all policies (whether or not they contain an Ingress section) are assumed to affect Ingress. If you want to write an egress-only policy, you must explicitly specify policyTypes [ \"Egress\" ]. Likewise, if you want to write a policy that specifies that no egress is allowed, you must specify a policyTypes value that include \"Egress\" (since such a policy would not include an Egress section and would otherwise default to just [ \"Ingress\" ]). This field is beta-level in 1.8
"""
mutable struct IoK8sApiNetworkingV1NetworkPolicySpec <: SwaggerModel
    egress::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1NetworkPolicyEgressRule} } # spec name: egress
    ingress::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1NetworkPolicyIngressRule} } # spec name: ingress
    podSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: podSelector
    policyTypes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: policyTypes

    function IoK8sApiNetworkingV1NetworkPolicySpec(;egress=nothing, ingress=nothing, podSelector=nothing, policyTypes=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("egress"), egress)
        setfield!(o, Symbol("egress"), egress)
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("ingress"), ingress)
        setfield!(o, Symbol("ingress"), ingress)
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("podSelector"), podSelector)
        setfield!(o, Symbol("podSelector"), podSelector)
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("policyTypes"), policyTypes)
        setfield!(o, Symbol("policyTypes"), policyTypes)
        o
    end
end # type IoK8sApiNetworkingV1NetworkPolicySpec

const _property_map_IoK8sApiNetworkingV1NetworkPolicySpec = Dict{Symbol,Symbol}(Symbol("egress")=>Symbol("egress"), Symbol("ingress")=>Symbol("ingress"), Symbol("podSelector")=>Symbol("podSelector"), Symbol("policyTypes")=>Symbol("policyTypes"))
const _property_types_IoK8sApiNetworkingV1NetworkPolicySpec = Dict{Symbol,String}(Symbol("egress")=>"Vector{IoK8sApiNetworkingV1NetworkPolicyEgressRule}", Symbol("ingress")=>"Vector{IoK8sApiNetworkingV1NetworkPolicyIngressRule}", Symbol("podSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("policyTypes")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }) = collect(keys(_property_map_IoK8sApiNetworkingV1NetworkPolicySpec))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1NetworkPolicySpec[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1NetworkPolicySpec[property_name]

function check_required(o::IoK8sApiNetworkingV1NetworkPolicySpec)
    (getproperty(o, Symbol("podSelector")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }, name::Symbol, val)
end


text:
#!/usr/bin/env julia
# File              : History.jl
# Author            : Denis Dumoulin <denis.dumoulin@uclouvain.be>
# Date              : 30.08.2021
# Last Modified Date: 01.09.2021
"""
    History

A structure allowing to record useful data.

# Fields
 - `t`: timelapse.
 - `X`: position `[X, Y, α]`.
 - `Ẋ`: velocity `[U, V, ω]`.
 - `acp`: aerodynamic coefficients obtained with pressure integration.
 - `acn`: aerodynamic coefficients obtained with control volume approach.
 - `aci`: aerodynamic coefficients obtained with momentum conservation.
 - `ϕs`: velocity potential on body surface.
 - `cps`: pressure coefficients on body surface.
 - `P`: impulse `[px, py, pm]`.
 - `Γ`: total body circulation.
 - `θ`: shedding angle (0 is the bissector).
"""
mutable struct History
    t::Array{Float64, 1}
    X::Array{Float64, 2}
    Ẋ::Array{Float64, 2}
    acp::Array{Float64, 2}
    acn::Array{Float64, 2}
    aci::Array{Float64, 2}
    ϕs::Array{Float64, 2}
    cps::Array{Float64, 2}
    P::Array{Float64, 2}
    NP::Array{Float64, 2}
    Γ::Array{Float64, 1}
    θ::Array{Float64, 1}
end

function History(N, dt, T)
    n = 2 + Int(T/dt)
    return History(zeros(n), zeros(n, 3), zeros(n, 3), zeros(n, 3), zeros(n, 3), zeros(n, 3),
                   zeros(n, N+1), zeros(n, N+1), zeros(n, 3), zeros(n, 3), zeros(n), zeros(n))
end

"""
    updatehistory!(h::History, i, t, X, Ẋ, acp, acn, aci, ϕs, cps, Γ, θ)

Update history of the airfoil.
"""
function updatehistory!(h::History, i, t, X, Ẋ, Γ, θ)
    h.t[i] = t
    h.X[i, :] = X
    h.Ẋ[i, :] = Ẋ
    h.Γ[i] = Γ
    h.θ[i] = θ
end

"""
    updatehistory!(h::History, sym::Symbol, i, val)

Update field `sym` of the history.
"""
function updatehistory!(h::History, sym::Symbol, i, val)
    v = getfield(h, sym)
    v[idx, :] = val
    setfield!(h, sym, v)
end

"""
    getlastvalues(h::History, sym::Symbol, i, n)

Return the `n` last values of `h.sym` up to index `i`.
Values with negative indices are replaced with zeros.
"""
function getlastvalues(h::History, sym::Symbol, i, n)
    v = getfield(h, sym)
    i < n &&  return [circshift(v, (-i, 0))[end-n+1:end, :][j, :] for j in n:-1:1]
    return [v[i-n+1:i, :][j, :] for j in n:-1:1]
end

"""
    plotaero(h::history, save=false)

Plot the aerodynamic coefficients obtained with all methods.

# Arguments
 - `h`: the `History` to read from.
 - `save=false`: whether to save the plots as a .png image.
"""
function plotaero(h::History, save=false)
    p1 = plot(h.t[6:end], h.acn[6:end, 1])
    plot!(p1, h.t[6:end], h.aci[6:end, 1])
    plot!(p1, h.t[6:end], h.acp[6:end, 1],
          title = "Cd", legend = false)

    p2 = plot(h.t[6:end], h.acn[6:end, 2])
    plot!(p2, h.t[6:end], h.aci[6:end, 2])
    plot!(p2, h.t[6:end], h.acp[6:end, 2],
          title = "Cl", legend = false)

    p3 = plot(h.t[6:end], h.acn[6:end, 3], label="Control volume")
    plot!(p3, h.t[6:end], h.aci[6:end, 3], label="Impulse conservation")
    plot!(p3, h.t[6:end], h.acp[6:end, 3],
          title = "Cm", label="Cp integration")

    p = plot(p1, p2, p3, layout = (1, 3), size=(900, 300))
    save && savefig("aero.png")
    return p
end

"""
    plotcps(h::History, t, save=false)

Plot the pressure coefficients at time `t`.

# Arguments
 - `h`: the `History` to read from.
 - `t`: time at which to evaluate cps.
 - `save=false`: whether to save the plots as a .png image.
"""
function plotcps(h::History, t, save=false)
    idx = findfirst(x -> x≥t, h.t)
    N = 1 + div(size(h.cps, 2) - 1, 2)
    p1 = plot(h.cps[idx, N:-1:1], label="Suction Side", size=(600,600))
    plot!(p1, h.cps[idx, N:end], label="Pressure Side", legend=:topleft)

    save && savefig("cps.png")
    return p1
end

"""
    plotϕs(h::History, t, save=false)

Plot the velocity potential distribution at time `t`.

# Arguments
 - `h`: the `History` to read from.
 - `t`: time at which to evaluate ϕs.
 - `save=false`: whether to save the plots as a .png image.
"""
function plotϕs(h::History, t, save=false)
    idx = findfirst(x -> x≥t, h.t)
    N = 1 + div(size(h.ϕs, 2) - 1, 2)
    p1 = plot(h.ϕs[idx, N:-1:1], label="Suction Side", size=(600,600))
    plot!(p1, h.ϕs[idx, N:end], label="Pressure Side", legend=:topleft)

    save && savefig("phis.png")
    return p1
end

ext: jl


text:
module SingularityAppBuilder

using Pkg

export build_recipe, build_app

function header(bootstrap="docker", source="julia")
    return "Bootstrap: $bootstrap\nFrom: $source"
end

newline!(recipe) = recipe*"\n\n"

function environment(readonly, depot_path="/opt/julia")
    if readonly
        return "%environment\n\texport JULIA_DEPOT_PATH=:$depot_path"
    else
        return "%environment\n\texport JULIA_DEPOT_PATH=$depot_path"
    end
end

function post(project, depot_path, readonly)
    section = "%post\n"
    section * "\t# Julia packages\n\texport JULIA_DEPOT_PATH=/opt/julia\n\t"*
    "export PATH=/usr/local/julia/bin:\$PATH\n"
    section * post(project)
    if readonly
        # fix permissions
        section * "\t# Permissions\n\tchmod -R 645 $depot_path"
    end
    return section
end

function post(path::AbstractString)
    # TODO
    # return "julia -e 'using Pkg; pkg"activate $path; instantiate"'"
end

function post(packages::Pkg.Types.PackageSpec)
    # TODO
end

function build_recipe(project, path::AbstractString; bootstrap="docker", source="julia",
        readonly=true, depot_path="/opt/julia")
    recipe = header(bootstrap, source)
    newline!(recipe)
    recipe*environment(readonly, depot_path)
    newline!(recipe)

    return path
end

end # module


text:
const AnteriorCingulum = BilateralStructure("anterior cingulum")

const AnteriorCingulateCortex = BilateralStructure("anterior cingulate cortex")

const PosteriorCingulateCortex = BilateralStructure("posterior cingulate cortex")

const PosteriorCingulum = BilateralStructure("posterior cingulum")

const CingulateCingulum = BilateralStructure("cingulate cingulum")
const ParahippocampalCingulum = BilateralStructure("parahippocampal cingulum")

#=
    children = (
        AnteriorCingulateCortex,
        AnteriorCingulum,
        PosteriorCingulateCortex,
        PosteriorCingulum
   )
=#
const CingulateGyrus = BilateralStructure("cingulate gyrus")


text:
abstract type WeightMatrix end

mutable struct BlockBootstrapWeightMatrix <: WeightMatrix
    seed::Int64
    obs::Array{Float64, 1}
    block_size::Int64
    n_bootstrap::Int64
    weights::Array{Float64, 2}
end
function BlockBootstrapWeightMatrix(seed::Int64, log_returns::Array{Float64, 1},
    block_size::Int64, n_bootstrap::Int64)
    # Step 1: Apply a Moving Block Bootstrap to the Log Return Path
    Random.seed!(seed)
    n = size(log_returns, 1)
    block_ind = 1:n-block_size+1
    b_samples = Array{Float64,2}(undef, n, n_bootstrap)
    for i in 1:n_bootstrap
        b_samples[:,i] = log_returns[block_bootstrap_index(block_ind, n, block_size)]
    end

    # Step 2: Calculate Summary Statistics for the Bootstrapped Samples
    dist = get_summary_stats(b_samples, log_returns)

    # Step 3: Calculate Inverse Weight Matrix
    W = inv(cov(dist))
    return BlockBootstrapWeightMatrix(seed, log_returns, block_size, n_bootstrap, W)
end

function block_bootstrap_index(block_ind, n, b)
    rand_blocks = sample(block_ind, floor(Int,n/b))
    sample_ind = transpose(repeat(rand_blocks, 1, b))
    sample_ind = sample_ind[:]
    addition_vec = repeat(0:b-1,floor(Int,n/b))
    sample_ind = sample_ind + addition_vec
    return sample_ind
end

function select_moments(obs)
    sample_mean = mean(obs)
    sample_std = std(obs)
    sample_kurt = normal_kurtosis(obs)
    sample_ks_stat = 0.0
    sample_hurst = generalized_hurts_exp(obs)
    sample_acf_lag1 = autocor(obs, [1])[1]
    return [sample_mean, sample_std, sample_kurt, sample_ks_stat, sample_hurst]
end


text:
using Octrees
using Test
using LinearAlgebra
using Random
using Plots
include("../src/misc/plotstuff.jl")
include("../src/misc/helper.jl")    # import nearest neighbors

@testset "Octrees.jl" begin
    # getindex
    xloc = randstring("01", 10)
    yloc = randstring("01", 10)
    zloc = randstring("01", 10)

    # constructor
    code = LocationalCode(xloc, yloc, zloc)
    @test code[1:5] == LocationalCode(xloc[1:5], yloc[1:5], zloc[1:5])
    @test code[1:end] == LocationalCode(xloc, yloc, zloc)
    @test code[3:end] == LocationalCode(xloc[3:end], yloc[3:end], zloc[3:end])

    # ops
    @test code * [1,0,0] == LocationalCode(xloc * "1", yloc * "0", zloc * "0")
    @test code * "100" == LocationalCode(xloc * "1", yloc * "0", zloc * "0")

    # Test locational code functionality
    for N in [100, 150, 200]
        X = randn(3, N)
        octree = Octree(X)
        for i in rand(1:N, 10)
            @test i in code_to_cell(octree, LocationalCode(octree, X[:, i])).data.ids
        end
    end

    # Check if brute force checking for nearest neighbors is equivalent to
    # finding nearest neighbors via octree.
    for iteration in 1:10
        n = 100
        X = randn(3, n)
        p = X[:, 1]
        r = 0.5
        octree = Octree(X)
        octree_neighbors = knn(octree, p, r)
        actual_neighbors = nearest_neighbors(X, p, r) |> Set
        @test octree_neighbors == actual_neighbors
    end

    # Sample from sphere
    for iteration in 1:10
        n = 300
        X = randn(3, n)
        X ./= mapslices(norm, X, dims=1)
        p = X[:, 1]
        r = 0.2
        octree = Octree(X)
        octree_neighbors = knn(octree, p, r)
        actual_neighbors = nearest_neighbors(X, p, r) |> Set
        @test octree_neighbors == actual_neighbors
    end

    # Plot octree of data projected onto xy-plane.
    X = randn(3,100)
    X[3,:] .= 0
    plt = plot!(Octree(X))
    savefig(plt, "testoutput.png")
end


text:
using PkgTemplates

lib_template = Template(
    dir="..",
    plugins=[
        PackageCompilerLib(lib_name="cg")
    ]
)
lib_template("CG")

max_stars_repo_name: simonbyrne/libcg

max_issues_repo_name: simonbyrne/libcg

max_forks_repo_name: simonbyrne/libcg


text:
################################################################################
# MISCELLANOUS
################################################################################
"Given the output of `airfoil_props()` it generates a file in aerodyn format"
function generate_aerodyn_file(file_name, re, AoA_d, clmin, clmax,
                                AoA_zl_d, lift_slope, cldata, cdfdata,
                                cdpdata, cmdata; numPoints=100)

    ncrit = "-"

    # Generates splines
    Top_Xtr = zeros(AoA_d) # Dummy arrays
    Bot_Xtr = zeros(AoA_d)
    (Cl_spl, Cd_spl, Cdp_spl, Cm_spl, Top_Xtr_spl, Bot_Xtr_spl, aoaZeroCl,
        aoaMinDrag, cdmin, aoamin, aoamax) = _makeSpline(AoA_d, cldata, cdfdata, cdpdata,
                                                            cmdata, Top_Xtr, Bot_Xtr)

    # Creates file
    _file_name = file_name[end-3:end]==".dat" ? file_name : file_name*".dat"

    write(_file_name,"    See filename for airfoil ID
    Cl and Cd values uncorrected
    line
    1        Number of airfoil tables in this file
    $(re/1e6)     Reynolds numbers in millions
    0.0      Control setting
    $(ncrit)     Stall angle (deg)
    $aoaZeroCl   Zero Cn angle of attack (deg)
    0.0   Cn slope for zero lift (dimensionless)
    0.0   Cn extrapolated to value at positive stall angle of attack
    0.0   Cn at stall value for negative angle of attack
    $aoaMinDrag     Angle of attack for minimum CD (deg)
    $cdmin   Minimum CD value\n")

    f = open(_file_name, "a")
    open(_file_name, "a") do x
        for aoa in range(aoamin, aoamax, length=numPoints)
            write(x, "$(aoa)\t $(Cl_spl(aoa))\t $(Cd_spl(aoa))\t $(Cm_spl(aoa)) \n")
        end
        write(x, "EOT")
    end

end

"Plots the contour of an airfoil given in x,y"
function plot_airfoil(x::Array{Float64,1}, y::Array{Float64,1};
                      label="", style="-k", figfactor=1.0,
                      title_str="Airfoil geometry", fig_id="airfoil_geometry",
                      alpha=1.0)
  # Sizes the figure
  figsize = [7*1.5,5*0.5]*figfactor
  xmin, xmax = -0.05, 1.05
  yrange = (xmax-xmin)/figsize[1] * figsize[2]
  ymin, ymax = -yrange/2, yrange/2
  fig = figure(fig_id, figsize=(figsize[1], figsize[2]))
  ax = gca()
  ax.set_xlim([xmin, xmax])
  ax.set_ylim([ymin, ymax])

  ax.plot(x,y, style, label=label, alpha=alpha)
  ax.set_xlabel("x")
  ax.set_ylabel("y")
  ax.grid(true, color="0.8", linestyle="--")
  ax.title.set_text(title_str)

  # if label!=""; legend(loc="best"); end;
end

"Receives a .dat file as pulled from airfoiltools.com containing the x and y
contour coordinates of an airfoil, and returns arrays x and y."
function readcontour(file_name; header_len=1, delim=" ", path="", output="arrays")

    x, y = Float64[], Float64[]

    open(joinpath(path, file_name)) do f
        for (i, line) in enumerate(eachline(f))

            # Ignores header
            if i <= header_len
                nothing
                # Parses each line
            else
                this_x, this_y = split(line, delim; keepempty=false)
                push!(x, Base.parse(Float64, this_x))
                push!(y, Base.parse(Float64, this_y))
            end

        end
    end

    if output=="arrays"
        return x, y
    elseif output=="matrix"
        return hcat(x, y)
    else
        error("Invalid `output` argument $(output).")
    end
end

"""
  Receives an airfoil contour and splits it up in upper and lower surfaces as
  divided by the chord line. It returns `(upper, lower)` with `upper=(x,y)` the
  points of the upper surface, ditto for `lower`. Both `upper` and `lower` are
  given in increasing order in x (i.e., from leading to trailing edge).
"""
function splitcontour(x, y)
    # ERROR CASES
    if !(x[1] in [0.0, 1.0])
        error("Invalid contour. x[1] must be either 0.0 or 1.0, got $(x[1]).")
    end

    # Flag indicating whether the contour start at the trailing or leading edge
    start_TE = x[1]==1.0

    # Find the opposite end of the contour
    end_i = -1
    for (i, xi) in enumerate(x)
        if i == 1
            nothing
          # Case of starting from the trailing edge
        elseif start_TE && xi > x[i-1]
            end_i = i - 1
            break
        # Case of leading edge
        elseif !start_TE  && xi < x[i-1]
            end_i = i - 1
            break
        end
    end

    # ERROR CASE
    if end_i == -1
        error("Logic error! End of airfoil not found!")
    end

    # Splits them up
    x_sec1, y_sec1 = x[1:end_i], y[1:end_i]
    x_sec2, y_sec2 = x[end_i:end], y[end_i:end]

    # Sorts them from LE to TE
    if x_sec1[1] > 0.5; reverse!(x_sec1); reverse!(y_sec1); end;
    if x_sec2[1] > 0.5; reverse!(x_sec2); reverse!(y_sec2); end;

    # Determines upper and lower surfaces
    if mean(y_sec1) > mean(y_sec2)
        upper = [x_sec1, y_sec1]
        lower = [x_sec2, y_sec2]
    else
        upper = [x_sec2, y_sec2]
        lower = [x_sec1, y_sec1]
    end

    return upper, lower
end

### INTERNAL FUNCTIONS #########################################################
"Receives the output of `airfoil_props()` and returns splines of them"
function _makeSpline(Alpha, Cl, Cd, Cdp, Cm, Top_Xtr, Bot_Xtr)

  # Splines
  Cl_spl = Dierckx.Spline1D(Alpha, Cl)
  Cd_spl = Dierckx.Spline1D(Alpha, Cd+Cdp)
  Cdp_spl = Dierckx.Spline1D(Alpha, Cdp)
  Cm_spl = Dierckx.Spline1D(Alpha, Cm)
  Top_Xtr_spl = Dierckx.Spline1D(Alpha, Top_Xtr)
  Bot_Xtr_spl = Dierckx.Spline1D(Alpha, Bot_Xtr)

  # Finds the zero-lift AOA
  zerolift(aoa) = Cl_spl(aoa)
  aoaZeroCl = Roots.fzero(zerolift, -8,10)

  # Finds AOA of minimum drag
  step = 1e-6
  zeroMinDrag(aoa) = (Cd_spl(aoa) - Cd_spl(aoa+step)) / step
  aoaMinDrag = Roots.fzero(zeroMinDrag, -8.0, 10.0)
  cdmin = Cd_spl(aoaMinDrag)

  aoamin = minimum(Alpha)
  aoamax = maximum(Alpha)

  return (Cl_spl, Cd_spl, Cdp_spl, Cm_spl, Top_Xtr_spl, Bot_Xtr_spl,
          aoaZeroCl, aoaMinDrag, cdmin, aoamin, aoamax)
end
### END OF TOOLS ###############################################################


text:
#==============================================================================#
# EKS.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/eks-2017-11-01.normal.json
#==============================================================================#

__precompile__()

module EKS

using AWSCore

"""
    using AWSSDK.EKS.create_cluster
    create_cluster([::AWSConfig], arguments::Dict)
    create_cluster([::AWSConfig]; name=, roleArn=, resourcesVpcConfig=, <keyword arguments>)

    using AWSCore.Services.eks
    eks([::AWSConfig], "POST", "/clusters", arguments::Dict)
    eks([::AWSConfig], "POST", "/clusters", name=, roleArn=, resourcesVpcConfig=, <keyword arguments>)

# CreateCluster Operation

Creates an Amazon EKS control plane.

The Amazon EKS control plane consists of control plane instances that run the Kubernetes software, like `etcd` and the API server. The control plane runs in an account managed by AWS, and the Kubernetes API is exposed via the Amazon EKS API server endpoint.

Amazon EKS worker nodes run in your AWS account and connect to your cluster's control plane via the Kubernetes API server endpoint and a certificate file that is created for your cluster.

The cluster control plane is provisioned across multiple Availability Zones and fronted by an Elastic Load Balancing Network Load Balancer. Amazon EKS also provisions elastic network interfaces in your VPC subnets to provide connectivity from the control plane instances to the worker nodes (for example, to support `kubectl exec`, `logs`, and `proxy` data flows).

After you create an Amazon EKS cluster, you must configure your Kubernetes tooling to communicate with the API server and launch worker nodes into your cluster. For more information, see [Managing Cluster Authentication](http://docs.aws.amazon.com/eks/latest/userguide/managing-auth.html) and [Launching Amazon EKS Worker Nodes](http://docs.aws.amazon.com/eks/latest/userguide/launch-workers.html)in the *Amazon EKS User Guide*.

# Arguments

## `name = ::String` -- *Required*
The unique name to give to your cluster.

## `version = ::String`
The desired Kubernetes version for your cluster. If you do not specify a value here, the latest version available in Amazon EKS is used.

## `roleArn = ::String` -- *Required*
The Amazon Resource Name (ARN) of the IAM role that provides permissions for Amazon EKS to make calls to other AWS API operations on your behalf. For more information, see [Amazon EKS Service IAM Role](http://docs.aws.amazon.com/eks/latest/userguide/service_IAM_role.html) in the **Amazon EKS User Guide**

## `resourcesVpcConfig = [ ... ]` -- *Required*
The VPC subnets and security groups used by the cluster control plane. Amazon EKS VPC resources have specific requirements to work properly with Kubernetes. For more information, see [Cluster VPC Considerations](http://docs.aws.amazon.com/eks/latest/userguide/network_reqs.html) and [Cluster Security Group Considerations](http://docs.aws.amazon.com/eks/latest/userguide/sec-group-reqs.html) in the *Amazon EKS User Guide*.
```
 resourcesVpcConfig = [
        "subnetIds" => <required> [::String, ...],
        "securityGroupIds" =>  [::String, ...]
    ]
```

## `clientRequestToken = ::String`
Unique, case-sensitive identifier you provide to ensure the idempotency of the request.

# Returns

`CreateClusterResponse`

# Exceptions

`ResourceInUseException`, `ResourceLimitExceededException`, `InvalidParameterException`, `ClientException`, `ServerException`, `ServiceUnavailableException` or `UnsupportedAvailabilityZoneException`.

# Example: To create a new cluster

The following example creates an Amazon EKS cluster called prod.

Input:
```
[
    "version" => "1.10",
    "name" => "prod",
    "clientRequestToken" => "1d2129a1-3d38-460a-9756-e5b91fddb951",
    "resourcesVpcConfig" => [
        "securityGroupIds" => [
            "sg-6979fe18"
        ],
        "subnetIds" => [
            "subnet-6782e71e",
            "subnet-e7e761ac"
        ]
    ],
    "roleArn" => "arn:aws:iam::012345678910:role/eks-service-role-AWSServiceRoleForAmazonEKS-J7ONKE3BQ4PI"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/eks-2017-11-01/CreateCluster)
"""
@inline create_cluster(aws::AWSConfig=default_aws_config(); args...) = create_cluster(aws, args)

@inline create_cluster(aws::AWSConfig, args) = AWSCore.Services.eks(aws, "POST", "/clusters", args)

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

"""
    using AWSSDK.EKS.delete_cluster
    delete_cluster([::AWSConfig], arguments::Dict)
    delete_cluster([::AWSConfig]; name=)

    using AWSCore.Services.eks
    eks([::AWSConfig], "DELETE", "/clusters/{name}", arguments::Dict)
    eks([::AWSConfig], "DELETE", "/clusters/{name}", name=)

# DeleteCluster Operation

Deletes the Amazon EKS cluster control plane.

**Note**
> If you have active services in your cluster that are associated with a load balancer, you must delete those services before deleting the cluster so that the load balancers are deleted properly. Otherwise, you can have orphaned resources in your VPC that prevent you from being able to delete the VPC. For more information, see [Deleting a Cluster](http://docs.aws.amazon.com/eks/latest/userguide/delete-cluster.html) in the *Amazon EKS User Guide*.

# Arguments

## `name = ::String` -- *Required*
The name of the cluster to delete.

# Returns

`DeleteClusterResponse`

# Exceptions

`ResourceInUseException`, `ResourceNotFoundException`, `ClientException`, `ServerException` or `ServiceUnavailableException`.

# Example: To delete a cluster

This example command deletes a cluster named `devel` in your default region.

Input:
```
[
    "name" => "devel"
]
```

Output:
```
Dict(

)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/eks-2017-11-01/DeleteCluster)
"""
@inline delete_cluster(aws::AWSConfig=default_aws_config(); args...) = delete_cluster(aws, args)

@inline delete_cluster(aws::AWSConfig, args) = AWSCore.Services.eks(aws, "DELETE", "/clusters/{name}", args)

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

"""
    using AWSSDK.EKS.describe_cluster
    describe_cluster([::AWSConfig], arguments::Dict)
    describe_cluster([::AWSConfig]; name=)

    using AWSCore.Services.eks
    eks([::AWSConfig], "GET", "/clusters/{name}", arguments::Dict)
    eks([::AWSConfig], "GET", "/clusters/{name}", name=)

# DescribeCluster Operation

Returns descriptive information about an Amazon EKS cluster.

The API server endpoint and certificate authority data returned by this operation are required for `kubelet` and `kubectl` to communicate with your Kubernetes API server. For more information, see [Create a kubeconfig for Amazon EKS](http://docs.aws.amazon.com/eks/latest/userguide/create-kubeconfig.html).

**Note**
> The API server endpoint and certificate authority data are not available until the cluster reaches the `ACTIVE` state.

# Arguments

## `name = ::String` -- *Required*
The name of the cluster to describe.

# Returns

`DescribeClusterResponse`

# Exceptions

`ResourceNotFoundException`, `ClientException`, `ServerException` or `ServiceUnavailableException`.

# Example: To describe a cluster

This example command provides a description of the specified cluster in your default region.

Input:
```
[
    "name" => "devel"
]
```

Output:
```
Dict(
    "cluster" => Dict(
        "version" => "1.10",
        "name" => "devel",
        "arn" => "arn:aws:eks:us-west-2:012345678910:cluster/devel",
        "certificateAuthority" => Dict(
            "data" => "LS0tLS1CRUdJTiBDRVJUSUZJQ0FURS0tLS0tCk1JSUN5RENDQWJDZ0F3SUJBZ0lCQURBTkJna3Foa2lHOXcwQkFRc0ZBREFWTVJNd0VRWURWUVFERXdwcmRXSmwKY201bGRHVnpNQjRYRFRFNE1EVXpNVEl6TVRFek1Wb1hEVEk0TURVeU9ESXpNVEV6TVZvd0ZURVRNQkVHQTFVRQpBeE1LYTNWaVpYSnVaWFJsY3pDQ0FTSXdEUVlKS29aSWh2Y05BUUVCQlFBRGdnRVBBRENDQVFvQ2dnRUJBTTZWCjVUaG4rdFcySm9Xa2hQMzRlVUZMNitaRXJOZGIvWVdrTmtDdWNGS2RaaXl2TjlMVmdvUmV2MjlFVFZlN1ZGbSsKUTJ3ZURyRXJiQyt0dVlibkFuN1ZLYmE3ay9hb1BHekZMdmVnb0t6b0M1N2NUdGVwZzRIazRlK2tIWHNaME10MApyb3NzcjhFM1ROeExETnNJTThGL1cwdjhsTGNCbWRPcjQyV2VuTjFHZXJnaDNSZ2wzR3JIazBnNTU0SjFWenJZCm9hTi8zODFUczlOTFF2QTBXb0xIcjBFRlZpTFdSZEoyZ3lXaC9ybDVyOFNDOHZaQXg1YW1BU0hVd01aTFpWRC8KTDBpOW4wRVM0MkpVdzQyQmxHOEdpd3NhTkJWV3lUTHZKclNhRXlDSHFtVVZaUTFDZkFXUjl0L3JleVVOVXM3TApWV1FqM3BFbk9RMitMSWJrc0RzQ0F3RUFBYU1qTUNFd0RnWURWUjBQQVFIL0JBUURBZ0trTUE4R0ExVWRFd0VCCi93UUZNQU1CQWY4d0RRWUpLb1pJaHZjTkFRRUxCUUFEZ2dFQkFNZ3RsQ1dIQ2U2YzVHMXl2YlFTS0Q4K2hUalkKSm1NSG56L2EvRGt0WG9YUjFVQzIrZUgzT1BZWmVjRVZZZHVaSlZCckNNQ2VWR0ZkeWdBYlNLc1FxWDg0S2RXbAp1MU5QaERDSmEyRHliN2pVMUV6VThTQjFGZUZ5ZFE3a0hNS1E1blpBRVFQOTY4S01hSGUrSm0yQ2x1UFJWbEJVCjF4WlhTS1gzTVZ0K1Q0SU1EV2d6c3JRSjVuQkRjdEtLcUZtM3pKdVVubHo5ZEpVckdscEltMjVJWXJDckxYUFgKWkUwRUtRNWEzMHhkVWNrTHRGQkQrOEtBdFdqSS9yZUZPNzM1YnBMdVoyOTBaNm42QlF3elRrS0p4cnhVc3QvOAppNGsxcnlsaUdWMm5SSjBUYjNORkczNHgrYWdzYTRoSTFPbU90TFM0TmgvRXJxT3lIUXNDc2hEQUtKUT0KLS0tLS1FTkQgQ0VSVElGSUNBVEUtLS0tLQo="
        ),
        "createdAt" => 1.527807879988e9,
        "endpoint" => "https://A0DCCD80A04F01705DD065655C30CC3D.yl4.us-west-2.eks.amazonaws.com",
        "resourcesVpcConfig" => Dict(
            "securityGroupIds" => [
                "sg-6979fe18"
            ],
            "subnetIds" => [
                "subnet-6782e71e",
                "subnet-e7e761ac"
            ],
            "vpcId" => "vpc-950809ec"
        ),
        "roleArn" => "arn:aws:iam::012345678910:role/eks-service-role-AWSServiceRoleForAmazonEKS-J7ONKE3BQ4PI",
        "status" => "ACTIVE"
    )
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/eks-2017-11-01/DescribeCluster)
"""
@inline describe_cluster(aws::AWSConfig=default_aws_config(); args...) = describe_cluster(aws, args)

@inline describe_cluster(aws::AWSConfig, args) = AWSCore.Services.eks(aws, "GET", "/clusters/{name}", args)

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

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

    using AWSCore.Services.eks
    eks([::AWSConfig], "GET", "/clusters", arguments::Dict)
    eks([::AWSConfig], "GET", "/clusters", <keyword arguments>)

# ListClusters Operation

Lists the Amazon EKS clusters in your AWS account in the specified region.

# Arguments

## `maxResults = ::Int`
The maximum number of cluster results returned by `ListClusters` in paginated output. When this parameter is used, `ListClusters` only returns `maxResults` results in a single page along with a `nextToken` response element. The remaining results of the initial request can be seen by sending another `ListClusters` request with the returned `nextToken` value. This value can be between 1 and 100. If this parameter is not used, then `ListClusters` returns up to 100 results and a `nextToken` value if applicable.

## `nextToken = ::String`
The `nextToken` value returned from a previous paginated `ListClusters` request where `maxResults` was used and the results exceeded the value of that parameter. Pagination continues from the end of the previous results that returned the `nextToken` value.

**Note**
> This token should be treated as an opaque identifier that is only used to retrieve the next items in a list and not for other programmatic purposes.

# Returns

`ListClustersResponse`

# Exceptions

`InvalidParameterException`, `ClientException`, `ServerException` or `ServiceUnavailableException`.

# Example: To list your available clusters

This example command lists all of your available clusters in your default region.

Input:
```
[

]
```

Output:
```
Dict(
    "clusters" => [
        "devel",
        "prod"
    ]
)
```

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/eks-2017-11-01/ListClusters)
"""
@inline list_clusters(aws::AWSConfig=default_aws_config(); args...) = list_clusters(aws, args)

@inline list_clusters(aws::AWSConfig, args) = AWSCore.Services.eks(aws, "GET", "/clusters", args)

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

end # module EKS

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


text:
#!/usr/bin/env julia --project

include("../../../experiments/OceanBoxGCM/simple_box.jl")
ClimateMachine.init()

const FT = Float64

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    include("../refvals/test_windstress_refvals.jl")

    # simulation time
    timestart = FT(0)    # s
    timeend = FT(3600)  # s
    timespan = (timestart, timeend)

    # DG polynomial order
    N = Int(4)

    # Domain resolution
    Nˣ = Int(5)
    Nʸ = Int(5)
    Nᶻ = Int(5)
    resolution = (N, Nˣ, Nʸ, Nᶻ)

    # Domain size
    Lˣ = 1e6    # m
    Lʸ = 1e6    # m
    H = 400   # m
    dimensions = (Lˣ, Lʸ, H)

    BC = (
        OceanBC(Impenetrable(NoSlip()), Insulating()),
        OceanBC(Impenetrable(FreeSlip()), Insulating()),
        OceanBC(Penetrable(KinematicStress()), Insulating()),
    )

    run_simple_box(
        "test_windstress_short_imex",
        resolution,
        dimensions,
        timespan,
        HomogeneousBox,
        imex = true,
        BC = BC,
        Δt = 60,
        refDat = refVals.imex,
    )

    run_simple_box(
        "test_windstress_short_explicit",
        resolution,
        dimensions,
        timespan,
        HomogeneousBox,
        imex = false,
        BC = BC,
        Δt = 180,
        refDat = refVals.explicit,
    )
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)

fibTailRecursive(a, b=0, c=1) = a === 0 ? b : fibTailRecursive(a - 1, c, b + c)

fibFormula(n) = round(((5 ^ .5 + 1) / 2) ^ n / (5 ^ .5))

function fibLinear(n)
  x,y = (0,1)
  for i = 1:n
    x,y = (y, x+y)
  end
  return x
end

function fibMatrix(n)
  prevMatrix = BigInt[1 1; 1 0]
  matrix = prevMatrix ^ n
  matrix[2, 1]
end

input = length(ARGS) === 0 ? 29 : parse(Int, ARGS[1])
println(fibMatrix(input))


text:
using MOEA
using Test
import Random: seed!
seed!(42)

for tests in [
    "DE.jl"
]
    include(tests)
end


text:
using Flux.Tracker: gradient

struct ODEAdjointSensitivityFunction{dgType,rateType,uType,F,J,PJ,UF,PF,G,JC,GC,A,SType,DG,MM,TJ,PJT,PJC,CP,INT} <: SensitivityFunction
  f::F
  jac::J
  paramjac::PJ
  uf::UF
  pf::PF
  g::G
  J::TJ
  pJ::PJT
  dg_val::dgType
  jac_config::JC
  g_grad_config::GC
  paramjac_config::PJC
  alg::A
  numparams::Int
  numindvar::Int
  f_cache::rateType
  discrete::Bool
  y::uType
  sol::SType
  dg::DG
  mass_matrix::MM
  checkpoints::CP
  integrator::INT
end

function ODEAdjointSensitivityFunction(f,jac,paramjac,uf,pf,g,u0,
                                      jac_config,g_grad_config,paramjac_config,
                                      p,f_cache,alg,discrete,y,sol,dg,mm,checkpoints)
  numparams = length(p)
  numindvar = length(u0)
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(alg)
  J = isautojacvec ? nothing : similar(sol.prob.u0, numindvar, numindvar)
  pJ = if !isquad(alg)
    isautojacvec ? nothing : similar(sol.prob.u0, numindvar, numparams)
  else
    nothing
  end
  integrator = if ischeckpointing(alg)
    integ = init(sol.prob, sol.alg, save_on=false)
    integ.u = y
    integ
  else
    nothing
  end
  dg_val = similar(u0, numindvar) # number of funcs size
  ODEAdjointSensitivityFunction(f,jac,paramjac,uf,pf,g,J,pJ,dg_val,
                               jac_config,g_grad_config,paramjac_config,
                               alg,numparams,numindvar,f_cache,
                               discrete,y,sol,dg,mm,checkpoints,integrator)
end

# u = λ'
function (S::ODEAdjointSensitivityFunction)(du,u,p,t)
  idx = length(S.y)
  y = S.y
  isautojacvec = DiffEqBase.has_jac(S.f) ? false : get_jacvec(S.alg)

  if isbcksol(S.alg)
    λ     = @view u[1:idx]
    dλ    = @view du[1:idx]
    grad  = @view u[idx+1:end-idx]
    dgrad = @view du[idx+1:end-idx]
    _y    = @view u[end-idx+1:end]
    dy    = @view du[end-idx+1:end]
    copyto!(y, _y)
    isautojacvec || S.sol.prob.f(dy, _y, p, t)
  else
    if ischeckpointing(S.alg)
      # assuming that in the forward direction `t0` < `t1`, and the
      # `checkpoints` vector is sorted with respect to the forward direction
      tidx = findlast(x->x <= t, S.checkpoints)
      t0 = S.checkpoints[tidx]
      dt = t-t0
      integrator = S.integrator
      if abs(dt) > integrator.opts.dtmin
        S.sol(integrator.u, t0)
        copyto!(integrator.uprev, integrator.u)
        integrator.t = t0
        # set `iter` to some arbitrary integer so that there won't be max maxiters error
        integrator.iter=100
        u_modified!(integrator, true)
        step!(integrator, dt, true)
        # `integrator.u` is aliased to `y`
      else
        S.sol(y,t)
      end
    else
      S.sol(y,t)
    end
    if isquad(S.alg)
      λ     = u
      dλ    = du
    else
      λ     = @view u[1:idx]
      dλ    = @view du[1:idx]
      grad  = @view u[idx+1:end]
      dgrad = @view du[idx+1:end]
    end
  end

  if !isautojacvec
    if DiffEqBase.has_jac(S.f)
      S.f.jac(S.J,y,p,t) # Calculate the Jacobian into J
    else
      S.uf.t = t
      jacobian!(S.J, S.uf, y, S.f_cache, S.alg, S.jac_config)
    end
    mul!(dλ',λ',S.J)
  elseif isquad(S.alg)
    _dy, back = Tracker.forward(y) do u
      out_ = map(zero, u)
      S.f(out_, u, p, t)
      Tracker.collect(out_)
    end
    dλ[:] = Tracker.data(back(λ)[1])
    isbcksol(S.alg) && (dy[:] = Tracker.data(_dy))
  else
    _dy, back = Tracker.forward(y, S.sol.prob.p) do u, p
      out_ = map(zero, u)
      S.f(out_, u, p, t)
      Tracker.collect(out_)
    end
    dλ[:], dgrad[:] = map(Tracker.data, back(λ))
    isbcksol(S.alg) && (dy[:] = Tracker.data(_dy))
  end

  dλ .*= -one(eltype(λ))

  if !S.discrete
    if S.dg != nothing
      S.dg(S.dg_val,y,p,t)
    else
      S.g.t = t
      gradient!(S.dg_val, S.g, y, S.alg, S.g_grad_config)
    end
    dλ .+= S.dg_val
  end

  if !isquad(S.alg) && !isautojacvec
    if DiffEqBase.has_paramjac(S.f)
      S.f.paramjac(S.pJ,y,S.sol.prob.p,t) # Calculate the parameter Jacobian into pJ
    else
      jacobian!(S.pJ, S.pf, S.sol.prob.p, S.f_cache, S.alg, S.paramjac_config)
    end
    mul!(dgrad',λ',S.pJ)
  end
  nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
function ODEAdjointProblem(sol,g,t=nothing,dg=nothing,
                           alg=SensitivityAlg();
                           checkpoints=sol.t,
                           callback=CallbackSet(),mass_matrix=I)

  f = sol.prob.f
  tspan = (sol.prob.tspan[2],sol.prob.tspan[1])
  t != nothing && (tspan = (t[end],t[1]))
  discrete = t != nothing

  isinplace = DiffEqBase.isinplace(sol.prob)
  p = sol.prob.p
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  isautojacvec = DiffEqBase.has_jac(f) ? false : get_jacvec(alg)
  p == nothing && error("You must have parameters to use parameter sensitivity calculations!")
  uf = DiffEqDiffTools.UJacobianWrapper(f,tspan[2],p)
  pg = UGradientWrapper(g,tspan[2],p)

  u0 = zero(sol.prob.u0)

  if DiffEqBase.has_jac(f) || isautojacvec
    jac_config = nothing
  else
    jac_config = build_jac_config(alg,uf,u0)
  end

  if !discrete
    if dg != nothing || isautojacvec
      pg_config = nothing
    else
      pg_config = build_grad_config(alg,pg,u0,p)
    end
  else
    pg_config = nothing
  end

  y = copy(sol(tspan[1])) # TODO: Has to start at interpolation value!
  paramjac_config = nothing
  pf = nothing
  if !isquad(alg)
    if DiffEqBase.has_paramjac(f) || isautojacvec
      paramjac_config = nothing
    else
      pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],y)
      paramjac_config = build_param_jac_config(alg,pf,y,p)
    end
  end

  len = isquad(alg) ? length(u0) : length(u0)+length(p)
  λ = similar(u0, len)
  sense = ODEAdjointSensitivityFunction(f,f.jac,f.paramjac,
                                       uf,pf,pg,u0,jac_config,pg_config,paramjac_config,
                                       p,deepcopy(u0),alg,discrete,
                                       y,sol,dg,mass_matrix,checkpoints)

  if discrete
    cur_time = Ref(length(t))
    function time_choice(integrator)
      cur_time[] > 0 ? t[cur_time[]] : nothing
    end
    affect! = let isq = isquad(alg), λ=λ, t=t, y=y, cur_time=cur_time, idx=length(u0)
      function (integrator)
        p, u = integrator.p, integrator.u
        λ  = isq ? λ : @view(λ[1:idx])
        g(λ,y,p,t[cur_time[]],cur_time[])
        if isq
          u .+= λ
        else
          u = @view u[1:idx]
          u .= λ .+ @view integrator.u[1:idx]
        end
        u_modified!(integrator,true)
        cur_time[] -= 1
      end
    end
    cb = IterativeCallback(time_choice,affect!,eltype(tspan);initial_affect=true)

    _cb = CallbackSet(cb,callback)
  else
    _cb = callback
  end

  z0 = isbcksol(alg) ? [vec(zero(λ)); vec(y)] : vec(zero(λ))
  ODEProblem(sense,z0,tspan,p,callback=_cb)
end

struct AdjointSensitivityIntegrand{pType,uType,rateType,S,AS,F,PF,PJC,A,PJT}
  sol::S
  adj_sol::AS
  f::F
  p::pType
  y::uType
  λ::uType
  pf::PF
  f_cache::rateType
  pJ::PJT
  paramjac_config::PJC
  alg::A
end

function AdjointSensitivityIntegrand(sol,adj_sol,alg=SensitivityAlg())
  f = sol.prob.f
  tspan = sol.prob.tspan
  p = sol.prob.p
  # we need to copy here, because later, we will call
  # `ReverseDiff.compile(ReverseDiff.GradientTape(pf′, (y, p)))`
  y = similar(sol.prob.u0)
  λ = similar(adj_sol.prob.u0)
  # we need to alias `y`
  pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],y)
  f_cache = similar(y)
  isautojacvec = DiffEqBase.has_paramjac(f) ? false : get_jacvec(alg)
  pJ = isautojacvec ? nothing : Matrix{eltype(sol.prob.u0)}(undef,length(sol.prob.u0),length(p))

  if DiffEqBase.has_paramjac(f) || isautojacvec
    paramjac_config = nothing
  else
    paramjac_config = build_param_jac_config(alg,pf,y,p)
  end
  AdjointSensitivityIntegrand(sol,adj_sol,f,p,y,λ,pf,f_cache,pJ,paramjac_config,alg)
end

function (S::AdjointSensitivityIntegrand)(out,t)
  y = S.y
  S.sol(y,t)
  λ = S.λ
  S.adj_sol(λ,t)
  λ .*= -one(eltype(λ))
  isautojacvec = DiffEqBase.has_paramjac(S.f) ? false : get_jacvec(S.alg)
  # y is aliased
  S.pf.t = t

  if !isautojacvec
    if DiffEqBase.has_paramjac(S.f)
      S.f.paramjac(S.pJ,y,S.p,t) # Calculate the parameter Jacobian into pJ
    else
      jacobian!(S.pJ, S.pf, S.p, S.f_cache, S.alg, S.paramjac_config)
    end
    mul!(out',λ',S.pJ)
  else
    _, back = Tracker.forward(y, S.p) do u, p
      out_ = map(zero, u)
      S.f(out_, u, p, t)
      Tracker.collect(out_)
    end
    out[:] = Tracker.data(back(λ)[2])
  end
  out'
end

function (S::AdjointSensitivityIntegrand)(t)
  out = similar(S.p)
  S(out,t)
end

function adjoint_sensitivities_u0(sol,alg,g,t=nothing,dg=nothing;
                                  abstol=1e-6,reltol=1e-3,
                                  iabstol=abstol, ireltol=reltol,sensealg=SensitivityAlg(checkpointing=!sol.dense,quad=false),
                                  checkpoints=sol.t,
                                  kwargs...)
  isquad(sensealg) && error("Can't get sensitivities of u0 with quadrature.")
  adj_prob = ODEAdjointProblem(sol,g,t,dg,sensealg,checkpoints=checkpoints)
  adj_sol = solve(adj_prob,alg;abstol=abstol,reltol=reltol,save_everystep=false,kwargs...)
  -adj_sol[end][1:length(sol.prob.u0)],
    adj_sol[end][(1:length(sol.prob.p)) .+ length(sol.prob.u0)]'
end

function adjoint_sensitivities(sol,alg,g,t=nothing,dg=nothing;
                               abstol=1e-6,reltol=1e-3,
                               iabstol=abstol, ireltol=reltol,sensealg=SensitivityAlg(checkpointing=!sol.dense),
                               checkpoints=sol.t,
                               kwargs...)
  adj_prob = ODEAdjointProblem(sol,g,t,dg,sensealg,checkpoints=checkpoints)
  isq = isquad(sensealg)
  adj_sol = solve(adj_prob,alg;abstol=abstol,reltol=reltol,save_everystep=isq,kwargs...)
  !isq && return adj_sol[end][(1:length(sol.prob.p)) .+ length(sol.prob.u0)]'
  integrand = AdjointSensitivityIntegrand(sol,adj_sol)

  if t == nothing
    res,err = quadgk(integrand,sol.prob.tspan[1],sol.prob.tspan[2],
                   atol=iabstol,rtol=ireltol)
  else
    res = zero(integrand.p)'
    for i in 1:length(t)-1
      res .+= quadgk(integrand,t[i],t[i+1],
                     atol=iabstol,rtol=ireltol)[1]
    end
  end
  res
end


text:
# Overloads
export acim, linearresponse, correlationsum, birkhoffcov, birkhoffvar
export lyapunov, covariancefunction, MA_process, AR_process

"""
    acim(L)

Output a Fun object giving the acim of the associated map"""
acim(K::SolutionInvWrapper) = K.op\K.u
acim(L::Operator,u::Fun=uniform(domainspace(L))) = SolutionInv(L,u)\u
acim(M::AbstractMarkovMap) = acim(Transfer(M))

"""
    linearresponse(L, X::Fun)

Output a Fun object giving the first-order change in the acim of the map under perturbation X"""
function linearresponse(K::SolutionInvWrapper,X::Fun)
  @assert all(isapprox.(X.(∂(domain(X))),0;atol=sqrt(eps(cfstype(X))*arclength(domain(X)))))
  K\(-acim(K)*X)'
end

"""
    zero_to(A::Fun, ρ::Fun=uniform(space(A)))

Output ρ*(A - ∫A dρ)"""
function zero_to(A::Fun,ρ::Fun=uniform(space(A)))
  ρA = ρ*A
  ρA -= ρ*sum(ρA)/sum(ρ)
end

"""
    correlationsum(L, A::Fun, ρ=acim(L))

Output resolvent of transfer operator applied to ρ*(A - ∫A dρ)"""
correlationsum(K::SolutionInvWrapper,A::Fun,ρ=acim(K)) = K \ zero_to(A,ρ)

"""
    birkhoffcov(L, A::Fun, B::Fun)

Output covariance of CLT-normalised Birkhoff sums of A and B under the map"""
function birkhoffcov(K::SolutionInvWrapper,A::Fun,B::Fun)
  ρ = acim(K); Az = zero_to(A,ρ); Bz = zero_to(B,ρ)
  sum(B*(K\Az)) + sum(A*(K\Bz)) - sum(B*Az)
end

"""
    birkhoffvar(L, A::Fun)

Output diffusion coefficient of A under the map using Green-Kubo formula"""
function birkhoffvar(K::SolutionInvWrapper,A::Fun)
  Az = zero_to(A,acim(K))
  2sum(A*(K\Az)) - sum(A*Az)
end

for OP in (:linearresponse,:correlationsum,:birkhoffvar)
  @eval $OP(L::Operator,X::Fun) = $OP(SolutionInv(L),X)
  @eval $OP(M::AbstractMarkovMap,X::Fun) = $OP(Transfer(M),X)
  @eval $OP(obj::Operator,X) = $OP(obj,Fun(X,domainspace(obj)))
  @eval $OP(obj::AbstractMarkovMap,X) = $OP(obj,Fun(X,domain(obj)))
end
birkhoffcov(L::Operator,X::Fun,Y::Fun) = birkhoffcov(SolutionInv(L),X,Y)
birkhoffcov(M::AbstractMarkovMap,X::Fun,Y::Fun) = birkhoffcov(Transfer(M),X,Y)

# lag-correlation function
for TYP=(:AbstractMarkovMap,:SolutionInvWrapper)
  @eval covariancefunction(K::$TYP,A::Fun,n::Int) = covariancefunction(Transfer(K),A,n)
  @eval covariancefunction(K::$TYP,A::Fun,B::Fun,n::Int) = covariancefunction(Transfer(K),A,B,n)
  @eval covariancefunction(K::$TYP,A::Fun;tol=eps(norm(A)^2)) = covariancefunction(Transfer(K),A;tol=tol)
  @eval covariancefunction(K::$TYP,A::Fun,B::Fun;tol=eps(norm(A)*norm(B))) = covariancefunction(Transfer(K),A,B;tol=tol)
end
  """
      cfA = covariancefunction(L, A, n)

  Compute the lag-covariance function against transfer operator L for observable A for n time steps
  in each direction. Specifically, the (k+1)th entry is the expectation of A A∘f^(k).
  """
function covariancefunction(L::Operator,A::Fun,n::Int;r=acim(L))
  @compat cf = Array{cfstype(A)}(undef,n+1) # TODO: update to coefficienttype
  Az = zero_to(A,r)
  cf[1] = sum(A*Az)
  for i = 1:n
    Az = L*Az
    Az -= r*sum(Az)
    chop!(Az)
    cf[i+1] = sum(A*Az)
  end
  cf
end

"""
    cf = covariancefunction(L, A; tol = eps(A))

Choose n so that the covariance declines to a given tolerance
"""
function covariancefunction(L::Operator,A::Fun;r=acim(L),tol=eps(norm(A.coefficients)^2))
  n = 16
  @compat cf = Array{cfstype(A)}(undef,n+1)
  Az = zero_to(A,r)
  cf[1] = sum(A*Az)
  for i = 1:n
    Az = L*Az
    Az -= r*sum(Az)
    chop!(Az)
    cf[i+1] = sum(A*Az)
  end
  while maximum(abs.(cf[div(2n,3):end]))>tol && n <= 2^20
    n_old = n
    n = 2n
    pad!(cf,n+1)
    for i = n_old+1:n
      Az = L*Az
      Az -= r*sum(Az)
      chop!(Az)
      cf[i+1] = sum(A*Az)
    end
    # chop!(Az)
  end

  if n==2^20
    warn("Maximum number of coefficients reached (n=$n)")
    cf
  else
    chop!(cf)
  end
end

"""
    cfA, cfB = covariancefunction(L, A, B, n)

Compute the lag-covariance function against transfer operator L between observables A and B for n time steps
in each direction. Specifically, the (k+1)th entry of cfA is the expectation of B A∘f^(k) and vice versa for cfB.
"""
function covariancefunction(L::Operator,A::Fun,B::Fun,n::Int;r=acim(L))
  @compat cfA = Array{cfstype(A)}(undef,n+1)  # TODO: update to coefficienttype
  @compat cfB = Array{cfstype(A)}(undef,n+1)
  Az = zero_to(A,r); Bz = zero_to(B,r)
  cfA[1] = sum(A*Bz)
  cfB[1] = cfA[1]
  for i = 1:n
    Az = L*Az;Bz = L*Bz
    Az -= r*sum(Az); Bz -= r*sum(Bz)
    chop!(Az); chop!(Bz)
    cfA[i+1] = sum(A*Bz); cfB[i+1] = sum(B*Az)
  end
  cfA,cfB
end

"""
    cfA, cfB = covariancefunction(L, A, B; tol=eps(A*B))

Chooses n so that the covariance declines to a given tolerance.
"""
function covariancefunction(L::Operator,A::Fun,B::Fun;r=acim(L),tol=eps(norm(coefficients(A))*norm(coefficients(B))))
  n = 16
  cfA,cfB = covariancefunction(L,A,B,n;r=r)
  while max(maximum(abs.(cfA[div(2n,3):end])),maximum(abs.(cfB[div(2n,3):end])))>tol && n <= 2^20
    n = 2n
    cfA,cfB = covariancefunction(L,A,B,n;r=r)   # chop!(Az)
  end

  if n==2^20
    warn("Maximum number of coefficients reached (n=$n)")
    cfA,cfB
  else
    chop!(cfA),chop!(cfB)
  end
end

# Lyapunov exponent
"""
    lyapunov(f, r=acim(f), sp=Space(rangedomain(f)))

Calculate Lyapunov exponent associated with f"""
function lyapunov(f::AbstractMarkovMap,r=acim(f),sp=Space(rangedomain(f)))
  sum(Fun(x->transferfunction(x,f,x->log(abs(f'(x)))*r(x)),sp))
end
function lyapunov(K::SolutionInvWrapper)
  L = Transfer(K)
  lyapunov(markovmap(L),acim(K),rangespace(L))
end
lyapunov(m::Operator) = lyapunov(markovmap(m),acim(m),rangespace(m))

struct LyapContainer{M<:AbstractMarkovMap,rr}
  m::M
  tr::rr
end
(lc::LyapContainer)(x) = log(abs(lc.m'(x)))*lc.tr(x)

function lyapunov(f::ComposedMarkovMap,r=acim(f))#,sp=Space(rangedomain(f)))
  tr = copy(r)
  lyap = sum(Fun(x->transferfunction(x,f.maps[end],LyapContainer(f.maps[end],tr)),rangedomain(f.maps[end])))
  for k = complength(f)-1:-1:1
    tr = Fun(x->transferfunction(x,f.maps[k+1],tr),rangedomain(f.maps[k+1]))
    lyap += sum(Fun(x->transferfunction(x,f.maps[k],LyapContainer(f.maps[k],tr)),rangedomain(f.maps[k])))
  end
  lyap
end

# Gaussian process parametrisation
function logMA_process(L,A::Fun)
  @assert isreal(A)
  corf = covariancefunction(L,A)
  # d = Circle(); sp = Taylor(d)
  Corfn = Fun(CosSpace(),[corf[1];2corf[2:end]])
  # lCf = log(Corfn)/2
  # SlogMA = Fun(SinSpace(),coefficients(lCf)[2:end])
  # MA = cos(SlogMA)*exp(lCf)
  #  MA = 1+expm1(Fun(sp,chop!(coefficients(log(Corfn)/2))))
  #  # TODO: because the special functions in ApproxFun are buggy
  # chop!(real.(coefficients(MA)[1:2:end]))

  ## OR
  # Corfn = Fun(Laurent(d),[corf[1];repeat(corf[2:end],inner=2)])
  # logMA = Fun(sp,
  # coefficients(Fun(x->log(abs(Corfn(e^(im*x)))),CosSpace())))/2

  logMA = Fun(Taylor(),coefficients(log(Corfn))/2)
end

"""
    MA_process(f,A::Fun)

Calculate coefficients of Gaussian MA process with the same autocorrelation function as A under the map f.

See Appendix A.2 of Wormell, C.L. & Gottwald, G.A. 'On the Validity of Linear Response Theory in High-Dimensional Deterministic Dynamical Systems' (2018).
"""
function MA_process(L,A::Fun)
  logMA = logMA_process(L,A)
  chop!(real.(coefficients(Fun(x->exp(logMA(x)),space(logMA)))))
end

"""
   AR_process(f,A::Fun)

Calculate coefficients of Gaussian AR process with the same autocorrelation function as A under the map f.

See Appendix A.2 of Wormell, C.L. & Gottwald, G.A. 'On the Validity of Linear Response Theory in High-Dimensional Deterministic Dynamical Systems' (2018).
"""
function AR_process(L,A::Fun)
  logMA = logMA_process(L,A)
  chop!(real.(coefficients(Fun(x->exp(-logMA(x)),space(logMA)))))
end


text:
runDenoising(dataArray,dataMask,dataVar,alpha=0.05,maskSize=nothing;connectedMask=false)=runDenoising(CPU1(),dataArray,dataMask,dataVar,alpha,maskSize,connectedMask)

"""

    runDenoising(dataArray,dataMask,dataVar,alpha=0.05,maskSize=nothing)

Runs the denoising step.
"""
function runDenoising(resource::CPU1,dataArray,dataMask,dataVar,alpha=0.05,maskSize=nothing,connectedMask=false)
    ballSize = cumsum([1;4;8;16;36;92;212;477])

    dim = size(dataArray)
    ndim = ndims(dataArray)
    ntime = dim[ndim]
    coord = dim[1:(ndim-1)]
    #coord = size(dataArray)[1:ndim-1]
    fullLength = prod(coord)
    dataMaskInd = (LinearIndices(dataMask))[findall(dataMask)]
    nvox = length(dataMaskInd)
    dataArray = copy(transpose(reshape(dataArray,(fullLength,ntime))[dataMaskInd,:]))
    iter = floor(Integer,log2(ntime))-1
    Dmax = exp2(iter)
    from = round(Int,exp2(iter+1))
    quant = alpha/(iter+1)
    realIter = floor(Integer,log2(from-1))-1
    thrs = [quantile.(Chisq(i),1-quant) for i=exp2.(0:(realIter-1))]

    if maskSize==nothing
        p = (1000/fullLength)^(1/length(coord))
        maskSize =ceil.(Integer,p.*[coord...])
        println("Mask size : $maskSize")
    elseif maskSize=="full"
        maskSize = [coord...]
    end

    ### Define the projection matrix
    dataProj = Array{Float64}(undef,(from-1,nvox))
    ## Inialize with the finest partition
    ## locations of time indexes in the finest partition
    loc = ceil.(Integer,collect(1:ntime)*Dmax/ntime)
    ## lengths of the finest partition intervals // returned as a vector
    num = counts(loc,maximum(loc))
    ## storage locations
    to = from-1
    from = to-length(num)+1
    ## projections
    for i=1:nvox
        dataProj[from:to,i]=counts(loc,maximum(loc),Weights(dataArray[:,i]))
    end

    ## loop from finer to thicker partitions
    for idx=(iter-1):-1:0
        ## new lengths of the partition intervals
        numNew = num[1:2:length(num)]+num[2:2:length(num)]
        ## new storage locations
        toNew = from-1
        fromNew = toNew-length(numNew)+1
        ## new projections
        dataProj[fromNew:toNew,:] = dataProj[from:2:to,:]+dataProj[from+1:2:to,:]
        ## normalize old projections
        dataProj[from:to,:] = dataProj[from:to,:]./sqrt.(num)
        ## update
        num = numNew
        to = toNew
        from = fromNew
    end
    ## last normalization is not done in the loop
    dataProj[1,:] = dataProj[1,:]./sqrt.(num)

    ## ########################### Build geometry #############################
    ## matrix of all the 3D coordinates at column indexes of data.matrix corresponds to row indexes in dataCoord

    dataCoord = CartesianIndices(coord)
    dataMask = dataMask[dataCoord]
    idOne = one(dataCoord[1])
    idBox = CartesianIndex(maskSize...)
    idMax = dataCoord[end]
    #dataCoord = reshape(dataCoord,fullLength,ndim-1)
    #arrCoord = reshape(1:fullLength,coord)
    arrCoord = LinearIndices(coord)
    ## transforms the data into a matrix of projections
    nproj = size(dataProj,1)

    ### Analysis
    println(nvox," voxels")
    #prog = Progress(nvox,1)
    resVisited = pmap(1:nvox) do pixIdx
        ## Find the neighbors of the current voxel pixIdx
        @info "$(pixIdx) voxels out of $(nvox)"
        pixInd = dataMaskInd[pixIdx]

        dCo = dataCoord[pixInd]

        inf = max(dCo-idBox,idOne)
        sup = min(idMax,dCo+idBox)
        #mask = BitSet(arrCoord[UnitRange.(inf.I,sup.I)...])

        mask = dataMask[UnitRange.(inf.I,sup.I)...]
        if connectedMask
            mask = label_components(mask)
            mask = OffsetArray(mask,UnitRange.((inf-dCo).I,(sup-dCo).I))
            mask = (mask .== mask[zero(dCo)])
        else
            mask = OffsetArray(mask,UnitRange.((inf-dCo).I,(sup-dCo).I))
        end

        mask = arrCoord[findall(mask) .+ dCo]

        #intersect!(mask,BitSet(dataMaskInd))
        maskIdx = [findfirst(isequal(x),dataMaskInd) for x=mask]

        goodPix = multitestH0(dataProj[:,maskIdx].-dataProj[:,pixIdx],dataVar[maskIdx].+dataVar[pixIdx],thrs)[:]
        neighborsInd = collect(mask)[goodPix]
        dist = vec([sum(abs2.(c.I)) for c in dataCoord[neighborsInd] .- dCo])
        neighborsInd = neighborsInd[sortperm(dist)]

        neighborsIdx = map(neighborsInd) do x
            findfirst(isequal(x),dataMaskInd)
        end
        ## projection of the dynamics of the neighbors
        neighborsProj = dataProj[:,neighborsIdx]

        ## Denoising procedure
        ## number of possible balls
        nV = sum(ballSize.<=length(neighborsIdx))

        ## Iv is the matrix of the projection estimates build over the successive balls
        iv = Array{Float64}(undef,(nproj,nV))
        ## datavarIv vector of the associated variances
        datavarIv = Array{Float64}(undef,nV)

        ## Initialize
        limits = kV = 1
        iv[:,kV] = neighborsProj[:,1]
        datavarIv[kV] = dataVar[pixIdx]

        while kV<nV
            kV+=1
            limitsNew=ballSize[kV]
            ringIdx=(limits+1):limitsNew
            jvTemp=mean(neighborsProj[:,ringIdx],dims=2)
            dataVarJvKv = mean(dataVar[ringIdx])/length(ringIdx)
            ## test thresholds with Bonferroni correction adapted to both partition number and interior balls

            thrs = [quantile(Chisq(i),1-quant/(kV-1)) for i=exp2.(0:(realIter-1))]
            ## test time coherence
            testcoh = multitestH0(iv[:,1:(kV-1)] .- jvTemp,datavarIv[1:(kV-1)] .+ dataVarJvKv,thrs)

            ## if no time coherence with previous estimates
            if !all(testcoh)
                kV = kV-1
                break
            end
            ## otherwise update projection estimates
            limits=limitsNew
            iv[:,kV] = mean(neighborsProj[:,1:limits],dims=2)
            datavarIv[kV] = dataVar[pixIdx]/limits
        end
        ## the denoised dynamics rescaled

        ix = mean(dataArray[:,neighborsIdx[1:limits]],dims=2)

        #### returns a Dict containing:
        #### 'Vx' a vector containing all the neighbors indexes used to build the denoised dynamic
        #### 'Ix' a vector containing the denoised dynamic
        #### 'Px' a vector containing the denoised projection
        #### 'Lx' a matrix containing the original coordinates of the neighbors in data.array
        #### 'Cx' a vector containing the original coordinates of the center
        Dict("Lx" => dataCoord[neighborsInd[1:limits]],"Cx"=>dataCoord[pixInd],"Px" => vec(iv[:,kV]),"Ix" =>vec(ix),"Vx"=>vec(neighborsIdx[1:limits]))
    end
    Dict("infoDen" => resVisited,"dataProj" => dataProj, "var"=> dataVar)
end

function getDenoisingResults(dataArray,denoisDicArr)
    for i in 1:length(denoisDicArr["infoDen"])
        tmp = denoisDicArr["infoDen"][i]
        ## get the spatial coordinates of the center
        mid = tmp["Cx"]
        ## replace the dynamics in dataArray by its denoised version
        dataArray[mid...,:] = tmp["Ix"]
    end
    dataArray
end

function makeDataProj(dataArray,dataMask)
    dim = size(dataArray)
    ndim = ndims(dataArray)
    ntime = dim[ndim]
    iter = floor(Integer,log2(ntime))-1
    Dmax = exp2(iter)
    from = round(Int,exp2(iter+1))
    coord = dim[1:(ndim-1)]
    fullLength = prod(coord)
    dataMaskInd = (LinearIndices(dataMask))[findall(dataMask)]
    nvox = length(dataMaskInd)
    dataArray = copy(transpose(reshape(dataArray,(fullLength,ntime))[dataMaskInd,:]))

     ### Define the projection matrix
    dataProj = Array{Float64}(undef,(from-1,nvox))
    ## Inialize with the finest partition
    ## locations of time indexes in the finest partition
    loc = ceil.(Integer,collect(1:ntime)*Dmax/ntime)
    ## lengths of the finest partition intervals // returned as a vector
    num = counts(loc,maximum(loc))
    ## storage locations
    to = from-1
    from = to-length(num)+1
    ## projections
    for i=1:nvox
        dataProj[from:to,i]=counts(loc,maximum(loc),Weights(dataArray[:,i]))
    end

    ## loop from finer to thicker partitions
    for idx=(iter-1):-1:0
        ## new lengths of the partition intervals
        numNew = num[1:2:length(num)]+num[2:2:length(num)]
        ## new storage locations
        toNew = from-1
        fromNew = toNew-length(numNew)+1
        ## new projections
        dataProj[fromNew:toNew,:] = dataProj[from:2:to,:]+dataProj[from+1:2:to,:]
        ## normalize old projections
        dataProj[from:to,:] = dataProj[from:to,:]./sqrt.(num)
        ## update
        num = numNew
        to = toNew
        from = fromNew
    end
    ## last normalization is not done in the loop
    dataProj[1,:] = dataProj[1,:]./sqrt.(num)
    dataProj
end

function getDataArraySmall(dataArray,dataMask)
    dim = size(dataArray)
    ndim = ndims(dataArray)
    ntime = dim[ndim]
    iter = floor(Integer,log2(ntime))-1
    Dmax = exp2(iter)
    from = round(Int,exp2(iter+1))
    coord = dim[1:(ndim-1)]
    fullLength = prod(coord)
    dataMaskInd = (LinearIndices(dataMask))[findall(dataMask)]
    nvox = length(dataMaskInd)
    dataArray = copy(transpose(reshape(dataArray,(fullLength,ntime))[dataMaskInd,:]))
end


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

using BAT
using Test

@testset "test_whiten" begin
end


text:
struct RawCor{T <: AbstractFloat}
    name::String
    n::Int
    m₁::Int
    m₂::Int
    vars::OrderedDict{Col, Int}
    rows::OrderedDict{Row, Tuple{Int,Int,Symbol}}
    cols::OrderedDict{Col, Vector{Pair{Symbol,T}}}
    rhs::Dict{Row, T}
    ranges::Dict{Row, T}
    bounds::OrderedDict{Col, Vector{Pair{Symbol,T}}}
    objgiven::Bool
    objname::Symbol
    rhsname::Symbol

    function RawCor(n::Integer, m₁::Integer, m₂::Integer,
                    vars::OrderedDict{Col, Int}, rows::OrderedDict{Row, Tuple{Int,Int,Symbol}},
                    cols::OrderedDict{Col, Vector{Pair{Symbol,T}}}, rhs::Dict{Row, T},
                    ranges::Dict{Row, T}, bounds::OrderedDict{Col, Vector{Pair{Symbol,T}}};
                    objgiven::Bool = false,
                    objname::Symbol = OBJ,
                    rhsname::Symbol = RHS,
                    name::String = "LP") where T <: AbstractFloat
        return new{T}(name,
                      n, m₁, m₂,
                      vars, rows, cols,
                      rhs, ranges, bounds,
                      objgiven, objname, rhsname)
    end
end

function parse_cor(::Type{T}, filename::AbstractString) where T <: AbstractFloat
    # Initialize auxiliary variables
    name       = "SLP"
    mode       = NAME
    rowidx     = 1
    eqrowidx   = 0
    ineqrowidx = 0
    varidx     = 1
    objgiven   = false
    objname    = OBJ
    rhsname    = RHS
    # Define sections
    vars    = OrderedDict{Col, Int}()
    rows    = OrderedDict{Row, Tuple{Int,Int,Symbol}}()
    cols    = OrderedDict{Col, Vector{Pair{Symbol,T}}}()
    rhs     = Dict{Row, T}()
    ranges  = Dict{Row, T}()
    bounds  = OrderedDict{Col, Vector{Pair{Symbol,T}}}()
    # Parse the file
    open(filename) do io
        firstline = split(readline(io))
        if Symbol(firstline[1]) == NAME
            name = join(firstline[2:end], " ")
        else
            throw(ArgumentError("`NAME` field is expected on the first line."))
        end
        for line in eachline(io)
            if mode == END
                # Parse finished
                break
            end
            words = split(line)
            first_word = Symbol(words[1])
            if first_word in COR_MODES && first_word != mode
                mode = first_word
                continue
            end
            if mode == ROWS
                rowsym = Symbol(words[2])
                rows[rowsym] =
                    words[1] == "N" ? (objgiven = true; objname = rowsym; (rowidx, 0, OBJ)) :
                    words[1] == "L" ? (ineqrowidx += 1; (rowidx, ineqrowidx, LEQ))            :
                    words[1] == "G" ? (ineqrowidx += 1; (rowidx, ineqrowidx, GEQ))            :
                    (eqrowidx   += 1; (rowidx, eqrowidx, EQ))
                rowidx += 1
            elseif mode == COLUMNS
                var = Symbol(words[1])
                if get!(vars, var, 0) == 0
                    vars[var] = varidx
                    varidx    += 1
                end
                for idx = 2:2:length(words)
                    push!(get!(cols, var, Pair{Symbol,T}[]),
                          Pair(Symbol(words[idx]), convert(T, parse(Float64, words[idx+1]))))
                end
            elseif mode == RHS
                _rhsname = Symbol(words[1])
                if _rhsname  != RHS && rhsname == RHS
                    rhsname = _rhsname
                end
                rhsname == _rhsname || error("Multiple RHS names given.")
                for idx = 2:2:length(words)
                    rhs[Symbol(words[idx])] = convert(T, parse(Float64, words[idx+1]))
                end
            elseif mode == RANGES
                for idx = 2:2:length(words)
                    ranges[Symbol(words[idx])] = convert(T, parse(Float64, words[idx+1]))
                end
            elseif mode == BOUNDS
                var = Symbol(words[3])
                bnd = words[1] == "LO" ? LOWER :
                      words[1] == "UP" ? UPPER :
                      words[1] == "BV" ? BINARY :
                      words[1] == "I"  ? INTEGER :
                      words[1] == "LI" ? INTEGER_LOWER :
                      words[1] == "UI" ? INTEGER_UPPER :
                      words[1] == "FR" ? FREE  : FIXED
                push!(get!(bounds, var, Pair{Symbol,T}[]),
                      Pair(bnd, convert(T, parse(Float64, bnd in [FREE, BINARY, INTEGER] ? "0" : words[4]))))
            else
                throw(ArgumentError("$(mode) is not a valid cor file mode."))
            end
        end
    end
    # Return raw data
    return RawCor(varidx - 1,
                  eqrowidx,
                  ineqrowidx,
                  vars,
                  rows,
                  cols,
                  rhs,
                  ranges,
                  bounds;
                  objgiven = objgiven,
                  objname = objname,
                  rhsname = rhsname,
                  name = name)
end
parse_cor(filename::AbstractString) = parse_cor(Float64, filename)

function sparsity(cor::RawCor)
    return 1 - length(cor.cols) / (cor.n * (cor.m₁ + cor.m₂))
end


text:
# Seperate file to make it easier to include separately from REPL for running pieces
using VectorizationBase, OffsetArrays
using VectorizationBase: data
using Test

const W64S = VectorizationBase.pick_vector_width(Float64)
const W64 = Int(VectorizationBase.register_size() ÷ sizeof(Float64))
const W32 = Int(VectorizationBase.register_size() ÷ sizeof(Float32))
const VE = Core.VecElement
randnvec(N = Val{W64}()) = Vec(ntuple(_ -> Core.VecElement(randn()), N))
function tovector(u::VectorizationBase.VecUnroll{_N,W,_T}) where {_N,W,_T}
    T = _T === VectorizationBase.Bit ? Bool : _T
    N = _N + 1; i = 0
    x = Vector{T}(undef, N * W)
    for n ∈ 1:N
        v = VectorizationBase.data(u)[n]
        for w ∈ 0:W-1
            x[(i += 1)] = VectorizationBase.extractelement(v, w)
        end
    end
    x
end
tovector(v::VectorizationBase.AbstractSIMDVector{W}) where {W} = [VectorizationBase.extractelement(v,w) for w ∈ 0:W-1]
tovector(v::VectorizationBase.LazyMulAdd) = tovector(VectorizationBase._materialize(v))
tovector(x) = x
tovector(i::MM{W,X}) where {W,X} = collect(range(data(i), step = X, length = W))
tovector(i::MM{W,X,I}) where {W,X,I<:Union{Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64}} = collect(range(data(i), step = I(X), length = I(W)))
A = randn(13, 17); L = length(A); M, N = size(A);

trunc_int(x::Integer, ::Type{T}) where {T} = x % T
trunc_int(x, ::Type{T}) where {T} = x
size_trunc_int(x::Signed, ::Type{T}) where {T} = signed(x % T)
size_trunc_int(x::Unsigned, ::Type{T}) where {T} = unsigned(x % T)
size_trunc_int(x, ::Type{T}) where {T} = x

check_within_limits(x, y) = @test x ≈ y
function check_within_limits(x::Vector{T}, y) where {T <: Integer}
    if Bool(VectorizationBase.has_feature(Val(:x86_64_avx512dq)))
        return @test x ≈ y
    end
    r = typemin(Int32) .≤ y .≤ typemax(Int32)
    xs = x[r]; ys = y[r]
    @test xs ≈ ys
end

maxi(a,b) = max(a,b)
mini(a,b) = min(a,b)
function maxi(a::T1,b::T2) where {T1<:Base.BitInteger,T2<:Base.BitInteger}
    T = promote_type(T1,T2)
    T(a > b ? a : b)
end
function mini(a::T1,b::T2) where {T1<:Base.BitInteger,T2<:Base.BitInteger}
    _T = promote_type(T1,T2)
    T = if T1 <: Signed || T2 <: Signed
        signed(_T)
    else
        _T
    end
    T(a < b ? a : b)
end
maxi_fast(a,b) = Base.FastMath.max_fast(a,b)
mini_fast(a,b) = Base.FastMath.min_fast(a,b)
maxi_fast(a::Base.BitInteger, b::Base.BitInteger) = maxi(a, b)
mini_fast(a::Base.BitInteger, b::Base.BitInteger) = mini(a, b)

include("accuracy.jl")

nothing


text:
@testset "read" begin
    testsized1 = "'Twas brillig, and the slithy toves\n" *
        "      Did gyre and gimble in the wabe:\n" *
        "All mimsy were the borogoves,\n" *
        "      And the mome raths outgrabe." *
        "An some extra nonsense at the end, just for good measure"

    io = IOBuffer(testsized1)
    s = SizedStream(io, 139)

    data = read(s)

    expected = "'Twas brillig, and the slithy toves\n" *
        "      Did gyre and gimble in the wabe:\n" *
        "All mimsy were the borogoves,\n" *
        "      And the mome raths outgrabe."
    @test String(data) == expected

    close(s)
end


text:
using TrajectoryOptimization
using LinearAlgebra
using SparseArrays

# Simple keplarian dynamics with thrust acceleration ODE function
function keplar!(out, x, u, p, t)
    # Earth gravitational parameter (likely would never be a static parameter,
    # but using here for testing purposes...)
    μ = p[1]

    # Requirements
    @views r    = norm(x[1:3])

    # Derivatives
    @views out[1:3] .= x[4:6]
    @views out[4:6] .= -(μ / r^3).*x[1:3] .+ u
    return nothing
end

# Analytical Jacobians
function stateJac!(out, x, u, p, t)
    μ = p[1]
    out .= 0.0

    # ddr/dv
    for i in 1:3
        out[i, 3 + i] = 1.0
    end

    # ddv/dr
    @views r = norm(x[1:3])
    @views mul!(out[4:6, 1:3], x[1:3], transpose(x[1:3]))
    @views out[4:6, 1:3] .*= (3.0 / r^2)
    for i in 1:3; out[3 + i, i] -= 1.0; end
    @views out[4:6, 1:3] .*= (μ/r^3)
    return nothing
end
function controlJac!(out, x, u, p, t)
    out .= 0.0
    for i in 1:3
        out[3 + i, i] = 1.0
    end
    return nothing
end
function staticJac!(out, x, u, p, t)
    out .= 0.0
    @views r = norm(x[1:3])
    @views out[4:6] .= -(1/r^3).*x[1:3]
end

# ===== Construct Path Functions
nFuncs      = 6
nStates     = 6
nControls   = 3
nStatic     = 1
stateSP     = sparse([1, 2, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6],
                     [4, 5, 6, 1, 2, 3, 1, 2, 3, 1, 2, 3],
                     [true for _ in 1:12])
controlSP   = sparse([4,5,6],[1,2,3],[true for _ in 1:3])
staticSP    = sparse([4,5,6],[1,1,1],[true for _ in 1:3])
adpf        = PathFunction(Dynamics(), keplar!, nFuncs, nStates, nControls, nStatic)
anpf        = PathFunction(Dynamics(), keplar!, stateJac!, controlJac!, staticJac!, nothing,
                nFuncs, nStates, nControls, nStatic, stateSP, controlSP, staticSP, [])

# ===== Test function evaluation
# Create inputs
adout       = Vector{Float64}(undef,
                TrajectoryOptimization.GetNumberOfFunctions(adpf))
anout       = Vector{Float64}(undef,
                TrajectoryOptimization.GetNumberOfFunctions(anpf))
states      = rand(TrajectoryOptimization.GetNumberOfStates(adpf))
controls    = rand(TrajectoryOptimization.GetNumberOfControls(adpf))
static      = rand(TrajectoryOptimization.GetNumberOfStatics(adpf))
time        = 0.0

# Test that nuber of parameters are correct
@test length(adout)     == nFuncs
@test length(states)    == nStates
@test length(controls)  == nControls
@test length(static)    == nStatic

# Evaluate functions with method
TrajectoryOptimization.EvaluateFunction(adpf, adout, states, controls, static, time)
TrajectoryOptimization.EvaluateFunction(anpf, anout, states, controls, static, time)

# Evaluate wrapped function
outWrapped  = Vector{Float64}(undef, nFuncs)
keplar!(outWrapped, states, controls, static, time)

# Test for equality
for i in 1:nFuncs
    @test adout[i] == outWrapped[i];
    @test anout[i] == outWrapped[i];
end

# ===== Test Sparsity Detection and Jacobian evaluation
# State Jacobian
(rs,cs,_)       = findnz(stateSP)
adstateJac      = zeros(6,6)
anstateJac      = zeros(6,6)
adstateJacSP    = sparse(rs,cs,ones(length(rs)))
anstateJacSP    = sparse(rs,cs,ones(length(rs)))
stateJacWrapped = zeros(6,6)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.State(),
    adpf, adstateJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.State(),
    anpf, anstateJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.State(),
    adpf, adstateJacSP, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.State(),
    anpf, anstateJacSP, states, controls, static, time)
stateJac!(stateJacWrapped, states, controls, static, time)
for col in 1:6
    for row in 1:6
        @test adpf.stateSP[row,col] == stateSP[row,col]
        @test adstateJac[row, col] ≈ stateJacWrapped[row, col]
        @test anstateJac[row, col] ≈ stateJacWrapped[row, col]
        @test adstateJacSP[row, col] ≈ stateJacWrapped[row, col]
        @test anstateJacSP[row, col] ≈ stateJacWrapped[row, col]
    end
end

# Control Jacobian
(rs,cs,_)       = findnz(controlSP)
adcontrolJac    = zeros(6,3)
ancontrolJac    = zeros(6,3)
adcontrolJacSP  = sparse(rs,cs,ones(length(rs)))
ancontrolJacSP  = sparse(rs,cs,ones(length(rs)))
controlJacWrap  = zeros(6,3)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Control(),
    adpf, adcontrolJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Control(),
    anpf, ancontrolJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Control(),
    adpf, adcontrolJacSP, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Control(),
    anpf, ancontrolJacSP, states, controls, static, time)
controlJac!(controlJacWrap, states,  controls, static, time)
for col in 1:3
    for row in 1:6
        @test adpf.controlSP[row,col] == controlSP[row,col]
        @test adcontrolJac[row, col] ≈ controlJacWrap[row, col]
        @test ancontrolJac[row, col] ≈ controlJacWrap[row, col]
        @test adcontrolJacSP[row, col] ≈ controlJacWrap[row, col]
        @test ancontrolJacSP[row, col] ≈ controlJacWrap[row, col]
    end
end

# Static Jacobian
(rs,cs,_)       = findnz(staticSP)
adstaticJac     = zeros(6,1)
anstaticJac     = zeros(6,1)
adstaticJacSP   = sparse(rs,cs,ones(length(rs)))
anstaticJacSP   = sparse(rs,cs,ones(length(rs)))
staticJacWrap   = zeros(6,1)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Static(),
    adpf, adstaticJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Static(),
    anpf, anstaticJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Static(),
    adpf, adstaticJacSP, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Static(),
    anpf, anstaticJacSP, states, controls, static, time)
staticJac!(staticJacWrap, states, controls, static, time)
for row in 1:6
    @test adpf.staticSP[row, 1] == staticSP[row, 1]
    @test adstaticJac[row, 1] ≈ staticJacWrap[row, 1]
    @test anstaticJac[row, 1] ≈ staticJacWrap[row, 1]
    @test adstaticJacSP[row, 1] ≈ staticJacWrap[row, 1]
    @test anstaticJacSP[row, 1] ≈ staticJacWrap[row, 1]
end

# Time Jacobian
adtimeJac       = zeros(6,1)
antimeJac       = zeros(6,1)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Time(),
    adpf, adtimeJac, states, controls, static, time)
TrajectoryOptimization.EvaluateJacobian(TrajectoryOptimization.Time(),
    anpf, antimeJac, states, controls, static, time)
for row in 1:6
    @test adtimeJac[row, 1] ≈ 0.0
    @test antimeJac[row, 1] ≈ 0.0
end


text:
export Model, Flag

abstract type Model end
abstract type Flag end


text:
"""For use in a macro accepting arguments like
@m a,b,c = f(d,e)
returns a Vector{Symbol} containing [:a,:b,:c]"""
function lhs(expr::Expr)
    @assert expr.head == :(=)
    a = expr.args[1]
    isa(a,Symbol) ? Symbol[a] : convert(Vector{Symbol},a.args)
end
"""For use in a macro accepting arguments like
@m a,b,c = f(d,e)
returns a tuple (:f, [:d,:e])"""
function rhs(expr::Expr)
    @assert expr.head == :(=)
    @assert expr.args[2].head == :call
    expr.args[2].args[1], convert(Vector{Symbol},expr.args[2].args[2:end])
end
"""For use in a macro accepting arguments like
f(a)>b
returns a tuple(:f, :a, :b)"""
function conditionsyms(expr::Expr)
    @assert expr.head == :comparison
    @assert expr.args[2] == :>
    expr.args[1].args[1], expr.args[1].args[2], expr.args[3]
end
"""Accepts something like
@perpulse `filt_value = do_filter(filter, pulse_record)`
Converts it to `PerPulseStep(do_filter, [:filter, :pulse_record], [:filt_value])`."""
macro perpulse(expr)
f,args = rhs(expr)
outs = lhs(expr)
	quote
		PerPulseStep($f, $args, $outs)
	end
end
""" Accepts something like
@threshold calibration = calibrate_nofit(filt_value_dc_hist, known_energies) when counted(filt_value_dc_hist)>3000
ThresholdStep(calibrate_nofit, [:filt_value_dc_hist,:known_energies],[:calibration],:filt_value_dc_hist, counted, 3000, true))
"""
macro threshold(expr, when, condition)
    @assert when == :when "2nd argument to @threshold must be :when, this is just for readability"
    f,args = rhs(expr)
    outs = lhs(expr)
    cf, a, b = conditionsyms(condition)
    aa=string(a)
    quote
        ThresholdStep($f, $args, $outs, Symbol($aa), $cf, $b, true)
    end
end
""" Accepts something like
@histogram update_histogram!(filt_value_hist, selection_good, filt_value)
Converts it to HistogramStep(update_histogram!, [:selection_good, :filt_value], [:filt_value_hist])
"""
macro histogram(expr)
    f = expr.args[1]
    args = expr.args[2:end]
    quote
        HistogramStep($f, $(args[2:end]), $(args[1:1]))
    end
end


text:
"""
    svdsolve(A::AbstractMatrix, [howmany = 1, which = :LR, T = eltype(A)]; kwargs...)
    svdsolve(f, m::Int, n::Int, [howmany = 1, which = :LR, T = Float64]; kwargs...)
    svdsolve(f, x₀, y₀, [howmany = 1, which = :LM]; kwargs...)
    svdsolve(f, x₀, y₀, howmany, which, algorithm)

Compute `howmany` singular values from the linear map encoded in the matrix `A` or by the
function `f`. Return singular values, left and right singular vectors and a
`ConvergenceInfo` structure.

### Arguments:
The linear map can be an `AbstractMatrix` (dense or sparse) or a general function or callable
object. Since both the action of the linear map and its adjoint are required in order to compute
singular values, `f` can either be a tuple of two callable objects (each accepting a single argument),
representing the linear map and its adjoint respectively, or, `f` can be a single callable object
that accepts two input arguments, where the second argument is a flag that indicates whether
the adjoint or the normal action of the linear map needs to be computed. The latter form still
combines well with the `do` block syntax of Julia, as in
```julia
vals, lvecs, rvecs, info = svdsolve(x₀, y₀, howmany, which; kwargs...) do (x, flag)
    if flag
        # y = compute action of adjoint map on x
    else
        # y = compute action of linear map on x
    end
    return y
end
```

For a general linear map encoded using either the tuple or the two-argument form, the best approach
is to provide a start vector `x₀` (in the domain of the linear map). Alternatively, one can specify
the number `n` of columns of the linear map, in which case `x₀ = rand(T, n)` is used, where the
default value of `T` is `Float64`, unless specified differently. If an `AbstractMatrix` is used,
a starting vector `x₀` does not need to be provided; it is chosen as `rand(T, size(A,1))`.

The next arguments are optional, but should typically be specified. `howmany` specifies how many
singular values and vectors should be computed; `which` specifies which singular values should
be targetted. Valid specifications of `which` are
*   `LR`: largest singular values
*   `SR`: smallest singular values
However, the largest singular values tend to converge more rapidly.

### Return values:
The return value is always of the form `vals, lvecs, rvecs, info = svdsolve(...)` with
*   `vals`: a `Vector{<:Real}` containing the singular values, of length at least `howmany`,
    but could be longer if more singular values were converged at the same cost.
*   `lvecs`: a `Vector` of corresponding left singular vectors, of the same length as `vals`.
*   `rvecs`: a `Vector` of corresponding right singular vectors, of the same length as `vals`.
    Note that singular vectors are not returned as a matrix, as the linear map could act on any
    custom Julia type with vector like behavior, i.e. the elements of the lists `lvecs` (`rvecs`)
    are objects that are typically similar to the starting guess `y₀` (`x₀`), up to a possibly
    different `eltype`. When the linear map is a simple `AbstractMatrix`, `lvecs` and `rvecs`
    will be `Vector{Vector{<:Number}}`.
*   `info`: an object of type [`ConvergenceInfo`], which has the following fields
    -   `info.converged::Int`: indicates how many singular values and vectors were actually
        converged to the specified tolerance `tol` (see below under keyword arguments)
    -   `info.residual::Vector`: a list of the same length as `vals` containing the residuals
        `info.residual[i] = A * rvecs[i] - vals[i] * lvecs[i]`.
    -   `info.normres::Vector{<:Real}`: list of the same length as `vals` containing the norm
        of the residual `info.normres[i] = norm(info.residual[i])`
    -   `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A` or `A'`.
    -   `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)
!!! warning "Check for convergence"
    No warning is printed if not all requested eigenvalues were converged, so always check
    if `info.converged >= howmany`.

### Keyword arguments:
Keyword arguments and their default values are given by:
*   `krylovdim`: the maximum dimension of the Krylov subspace that will be constructed.
    Note that the dimension of the vector space is not known or checked, e.g. `x₀` should not
    necessarily support the `Base.length` function. If you know the actual problem dimension
    is smaller than the default value, it is useful to reduce the value of `krylovdim`, though
    in principle this should be detected.
*   `tol`: the requested accuracy according to `normres` as defined above. If you work in e.g.
    single precision (`Float32`), you should definitely change the default value.
*   `maxiter`: the number of times the Krylov subspace can be rebuilt; see below for further
    details on the algorithms.
*   `orth`: the orthogonalization method to be used, see [`Orthogonalizer`](@ref)

### Algorithm
The last method, without default values and keyword arguments, is the one that is finally called,
and can also be used directly. Here the algorithm is specified, though currently only [`GKL`](@ref)
is available. `GKL` refers to the the partial Golub-Kahan-Lanczos bidiagonalization which forms
the basis for computing the approximation to the singular values. This factorization is dynamically
shrunk and expanded (i.e. thick restart) similar to the Krylov-Schur factorization for eigenvalues.
"""
function svdsolve end

function svdsolve(A::AbstractMatrix, howmany::Int = 1, which::Selector = :LR, T::Type = eltype(A); kwargs...)
    svdsolve(A, rand(T, size(A,1)), howmany, which; kwargs...)
end
function svdsolve(f, n::Int, howmany::Int = 1, which::Selector = :LR, T::Type = Float64; kwargs...)
    svdsolve(f, rand(T, n), howmany, which; kwargs...)
end
function svdsolve(f, x₀, howmany::Int = 1, which::Symbol = :LR; kwargs...)
    which == :LR || which == :SR || error("invalid specification of which singular values to target: which = $which")
    alg = GKL(; kwargs...)
    svdsolve(f, x₀, howmany, which, alg)
end

function svdsolve(A, x₀, howmany::Int, which::Symbol, alg::GKL)
    krylovdim = alg.krylovdim
    maxiter = alg.maxiter
    howmany > krylovdim && error("krylov dimension $(krylovdim) too small to compute $howmany singular values")

    ## FIRST ITERATION: setting up
    numiter = 1
    # Compute Lanczos factorization
    iter = GKLIterator(svdfun(A), x₀, alg.orth)
    fact = initialize(iter)
    numops = 2
    sizehint!(fact, krylovdim)
    β = normres(fact)
    tol::eltype(β) = alg.tol
    while length(fact) < krylovdim
        fact = expand!(iter, fact)
        numops += 2
        normres(fact) < tol && length(fact) >= howmany && break
    end

    # Process
    # allocate storage
    HH = fill(zero(eltype(fact)), krylovdim+1, krylovdim)
    PP = fill(zero(eltype(fact)), krylovdim, krylovdim)
    QQ = fill(zero(eltype(fact)), krylovdim, krylovdim)

    # initialize
    β = normres(fact)
    m = length(fact)
    P = copyto!(view(PP, 1:m, 1:m), I)
    Q = copyto!(view(QQ, 1:m, 1:m), I)
    f = view(HH, m+1, 1:m)
    B = rayleighquotient(fact) # Bidiagional (lower)

    # compute singular value decomposition
    P, S, Q = bidiagsvd!(B, P, Q)
    if which == :SR
        reversecols!(P)
        reverserows!(S)
        reverserows!(Q)
    elseif which != :LR
        error("invalid specification of which singular values to target: which = $which")
    end
    mul!(f, view(Q', m, :), β)

    converged = 0
    while converged < length(fact) && abs(f[converged+1]) < tol
        converged += 1
    end

    ## OTHER ITERATIONS: recycle
    while numiter < maxiter && converged < howmany
        numiter += 1

        # Determine how many to keep
        keep = div(3*krylovdim + 2*converged, 5) # strictly smaller than krylovdim since converged < howmany <= krylovdim, at least equal to converged

        # Update basis by applying P and Q using Householder reflections
        U = basis(fact, :U)
        basistransform!(U, view(P, :, 1:keep))
        # for j = 1:m
        #     h, ν = householder(P, j:m, j)
        #     lmul!(h, view(P, :, j+1:krylovdim))
        #     rmul!(U, h')
        # end
        V = basis(fact, :V)
        basistransform!(V, view(Q', :, 1:keep))
        # for j = 1:m
        #     h, ν = householder(Q, j, j:m)
        #     rmul!(view(Q, j+1:krylovdim, :), h)
        #     rmul!(V, h)
        # end

        # Shrink GKL factorization (no longer strictly GKL)
        r = residual(fact)
        U[keep+1] = rmul!(r, 1/normres(fact))
        H = fill!(view(HH, 1:keep+1, 1:keep), 0)
        @inbounds for j = 1:keep
            H[j,j] = S[j]
            H[keep+1,j] = f[j]
        end

        # Restore bidiagonal form in the first keep columns
        @inbounds for j = keep:-1:1
            h, ν = householder(H, j+1, 1:j, j)
            H[j+1,j] = ν
            H[j+1,1:j-1] .= 0
            rmul!(view(H, 1:j, :), h')
            rmul!(V, h')
            h, ν = householder(H, 1:j, j, j)
            H[j,j] = ν
            @inbounds H[1:j-1,j] .= 0
            lmul!(h, view(H, :, 1:j-1))
            rmul!(U, h')
        end
        @inbounds for j = 1:keep
            fact.αs[j] = H[j,j]
            fact.βs[j] = H[j+1,j]
        end
        fact = shrink!(fact, keep)

        # GKL factorization: recylce fact
        while length(fact) < krylovdim
            fact = expand!(iter, fact)
            numops += 2
            normres(fact) < tol && length(fact) >= howmany && break
        end

        # post process
        β = normres(fact)
        m = length(fact)
        P = copyto!(view(PP, 1:m, 1:m), I)
        Q = copyto!(view(QQ, 1:m, 1:m), I)
        f = view(HH, m+1, 1:m)
        B = rayleighquotient(fact) # Bidiagional (lower)

        # compute singular value decomposition
        P, S, Q = bidiagsvd!(B, P, Q)
        if which == :SR
            reversecols!(P)
            reverserows!(S)
            reverserows!(Q)
        elseif which != :LR
            error("invalid specification of which singular values to target: which = $which")
        end
        mul!(f, view(Q', m, :), β)
        converged = 0
        while converged < length(fact) && abs(f[converged+1]) < tol
            converged += 1
        end
    end

    if converged > howmany
        howmany = converged
    end
    values = S[1:howmany]

    # Compute schur vectors
    Pv = view(P,:, 1:howmany)
    Qv = view(Q, 1:howmany, :)

    # Compute convergence information
    leftvectors = let U = basis(fact, :U)
        [U*v for v in cols(Pv)]
    end
    rightvectors = let V = basis(fact, :V)
        [V*v for v in cols(Qv')]
    end
    residuals = let r = residual(fact)
        [r*last(v) for v in cols(Qv')]
    end
    normresiduals = let f = f
        map(i->abs(f[i]), 1:howmany)
    end

    return values, leftvectors, rightvectors, ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
end

svdfun(A::AbstractMatrix) = (x,flag) -> flag ? A'*x : A*x
svdfun((f,fadjoint)::Tuple{Any,Any}) = (x,flag) -> flag ? fadjoint(x) : f(x)
svdfun(f) = f


text:
module DistributionsAD

using PDMats,
      LinearAlgebra,
      Distributions,
      Random,
      SpecialFunctions,
      StatsFuns,
      Compat,
      Requires,
      ZygoteRules,
      ChainRules,  # needed for `ChainRules.chol_blocked_rev`
      FillArrays

using SpecialFunctions: logabsgamma, digamma
using LinearAlgebra: copytri!, AbstractTriangular
using Distributions: AbstractMvLogNormal,
                     ContinuousMultivariateDistribution
using Base.Iterators: drop

import StatsFuns: logsumexp,
                  binomlogpdf,
                  nbinomlogpdf,
                  poislogpdf,
                  nbetalogpdf
import Distributions: MvNormal,
                      MvLogNormal,
                      poissonbinomial_pdf_fft,
                      logpdf,
                      quantile,
                      PoissonBinomial,
                      Binomial,
                      BetaBinomial,
                      Erlang

export TuringScalMvNormal,
       TuringDiagMvNormal,
       TuringDenseMvNormal,
       TuringMvLogNormal,
       TuringPoissonBinomial,
       TuringWishart,
       TuringInverseWishart,
       arraydist,
       filldist

include("common.jl")
include("univariate.jl")
include("multivariate.jl")
include("mvcategorical.jl")
include("matrixvariate.jl")
include("flatten.jl")
include("arraydist.jl")
include("filldist.jl")

include("zygote.jl")

@init begin
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
        using .ForwardDiff: @define_binary_dual_op # Needed for `eval`ing diffrules here
        include("forwarddiff.jl")

        # loads adjoint for `poissonbinomial_pdf_fft`
        include("zygote_forwarddiff.jl")
    end

    @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
        include("reversediff.jl")
    end

    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
        using DiffRules
        using SpecialFunctions
        using LinearAlgebra: AbstractTriangular
        using .Tracker: Tracker, TrackedReal, TrackedVector, TrackedMatrix,
                        TrackedArray, TrackedVecOrMat, track, @grad, data
        include("tracker.jl")
    end
end

end


text:
# Method for function with no Jacobian.
function fixedpoint(f,
    initial_x::AbstractArray{T};
    method::Symbol = :anderson,
    xtol::Real = zero(T),
    ftol::Real = convert(T,1e-8),
    iterations::Integer = 1_000,
    store_trace::Bool = false,
    show_trace::Bool = false,
    extended_trace::Bool = false,
    linesearch = LineSearches.Static(),
    factor::Real = one(T),
    autoscale::Bool = true,
    m::Integer = 5,
    beta::Real = 1,
    droptol::Real = 0,
    autodiff::Symbol = :central,
    inplace::Bool = !applicable(f, initial_x)) where T
    # Check for weird case. (Causes to hang for now)
    # typeof(f) <: Union{InplaceObjective, NotInplaceObjective} ? error("Union{InplaceObjective, NotInplaceObjective} Case") : true;
    # Wrapping
    if inplace
        function g!(out, x)
            f(out, x);
            out .-= x;
        end
        dg = OnceDifferentiable(g!, initial_x, copy(initial_x), autodiff)
    else
        g(x) = f(x) - x;
        dg = OnceDifferentiable(g, initial_x, copy(initial_x); autodiff = autodiff, inplace = inplace)
    end

    return nlsolve(dg,
            initial_x, method = method, xtol = xtol, ftol = ftol,
            iterations = iterations, store_trace = store_trace,
            show_trace = show_trace, extended_trace = extended_trace,
            linesearch = linesearch, factor = factor, autoscale = autoscale,
            m = m, beta = beta, droptol = droptol)
end