text:
@testset "Rosenbrock Tests" begin
## Breakout these since no other test of their adaptivity

using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra, LinearSolve
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear,
                              prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear

dts = (1/2) .^ (6:-1:3)
testTol = 0.2

### Rosenbrock23()

prob = prob_ode_linear

sim = test_convergence(dts,prob,Rosenbrock23())
@test sim.𝒪est[:final] ≈ 2 atol=testTol

sol = solve(prob,Rosenbrock23())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Rosenbrock23())
@test sim.𝒪est[:final] ≈ 2 atol=testTol

sol = solve(prob,Rosenbrock23())
@test length(sol) < 20

prob = prob_ode_bigfloat2Dlinear

sim = test_convergence(dts,prob,Rosenbrock23(linsolve=QRFactorization()))
@test sim.𝒪est[:final] ≈ 2 atol=testTol

sol = solve(prob,Rosenbrock23(linsolve=QRFactorization()))
@test length(sol) < 20

### Rosenbrock32()

prob = prob_ode_linear

sim = test_convergence(dts,prob,Rosenbrock32())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,Rosenbrock32())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Rosenbrock32())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,Rosenbrock32())
@test length(sol) < 20

### ROS3P()

prob = prob_ode_linear

sim = test_convergence(dts,prob,ROS3P())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS3P())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,ROS3P())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS3P())
@test length(sol) < 20

### Rodas3()

prob = prob_ode_linear

sim = test_convergence(dts,prob,Rodas3())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,Rodas3())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Rodas3())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,Rodas3())
@test length(sol) < 20

println("4th order Rosenbrocks")

### RosShamp4

prob = prob_ode_linear

sim = test_convergence(dts,prob,RosShamp4())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,RosShamp4())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,RosShamp4())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,RosShamp4())
@test length(sol) < 20

### Veldd4

prob = prob_ode_linear

sim = test_convergence(dts,prob,Veldd4())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,Veldd4())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Veldd4())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,Veldd4())
@test length(sol) < 20

### Velds4

prob = prob_ode_linear

sim = test_convergence(dts,prob,Velds4())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,Velds4())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Velds4())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,Velds4())
@test length(sol) < 20

### GRK4T

prob = prob_ode_linear

sim = test_convergence(dts,prob,GRK4T())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,GRK4T())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,GRK4T())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,GRK4T())
@test length(sol) < 20

### GRK4A
dts = (1/2) .^ (7:-1:4)

prob = prob_ode_linear

sim = test_convergence(dts,prob,GRK4A())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,GRK4A())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,GRK4A())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,GRK4A())
@test length(sol) < 20

### Ros4LStab

prob = prob_ode_linear

sim = test_convergence(dts,prob,Ros4LStab())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,Ros4LStab())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Ros4LStab())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,Ros4LStab())
@test length(sol) < 20

### Rosenbrock-W Algorithms

println("Rosenbrock-W")

### ROS34PW1a
prob = prob_ode_linear

sim = test_convergence(dts,prob,ROS34PW1a())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS34PW1a())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,ROS34PW1a())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS34PW1a())
@test length(sol) < 20

### ROS34PW1b
prob = prob_ode_linear

sim = test_convergence(dts,prob,ROS34PW1b())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS34PW1b())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,ROS34PW1b())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS34PW1b())
@test length(sol) < 20

### ROS34PW2
prob = prob_ode_linear

sim = test_convergence(dts,prob,ROS34PW2())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS34PW2())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,ROS34PW2())
@test sim.𝒪est[:final] ≈ 3 atol=testTol

sol = solve(prob,ROS34PW2())
@test length(sol) < 20

### ROS34PW3
prob = prob_ode_linear

sim = test_convergence(dts,prob,ROS34PW3())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,ROS34PW3())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,ROS34PW3())
@test sim.𝒪est[:final] ≈ 4 atol=testTol

sol = solve(prob,ROS34PW3())
@test length(sol) < 20

### RosenbrockW6S4OS
sim = test_convergence(dts,prob,RosenbrockW6S4OS())#test inplace
@test sim.𝒪est[:final] ≈ 4 atol=testTol

prob = prob_ode_linear

sim = test_convergence(dts,prob,RosenbrockW6S4OS())#test non-inplace
@test sim.𝒪est[:final] ≈ 4 atol=testTol

### Rodas4 Algorithms

println("RODAS")

dts = (1/2) .^ (7:-1:4)

prob = prob_ode_linear

sim = test_convergence(dts,prob,Rodas4(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4(autodiff=false),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4(autodiff=false))
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas42(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 5 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas42())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4P(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4P())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4P2(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4P2())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Rodas4(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4())
@test length(sol) < 20

println("Rodas4 with finite diff")

sim = test_convergence(dts,prob,Rodas4(autodiff=false),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4(autodiff=false))
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4(autodiff=false,
                       diff_type=Val{:forward}),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4(autodiff=false,diff_type=Val{:forward}))
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4(autodiff=false,
                       diff_type=Val{:complex}),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4(autodiff=false,diff_type=Val{:complex}))
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas42(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 5 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas42())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4P(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4P())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4P2(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4P2())
@test length(sol) < 20

println("Rodas4P2 with finite diff")

sim = test_convergence(dts,prob,Rodas4P2(autodiff=false),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas4P2(autodiff=false))
@test length(sol) < 20

### Rodas5
println("Rodas5")

prob = prob_ode_linear

dts = (1/2) .^ (6:-1:3)
sim = test_convergence(dts,prob,Rodas5(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 5 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas5())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Rodas5(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 5 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol

sol = solve(prob,Rodas5())
@test length(sol) < 20

prob = ODEProblem((u,p,t)->0.9u, 0.1, (0., 1.0))
@test_nowarn solve(prob, Rosenbrock23(autodiff=false))
end


text:
using KendallTau
using Test

@testset "KendallTau.jl" begin
  include("rankcorr.jl")
  include("corkendall.jl")
  include("skipmissingpairs.jl")
end

ext: jl


text:
# # General example
# This function takes an SDP in standard form, reduces it, formulates it as (hermitian) SDP, and solves it with JuMP

using LinearAlgebra
using SDPSymmetryReduction
using JuMP
using CSDP
using SparseArrays

function reduceAndSolve(C, A, b;
    objSense = MathOptInterface.MAX_SENSE,
    verbose = false,
    complex = false,
    limitSize = 3000)

    tmd = @timed admPartSubspace(C, A, b, verbose)
    P = tmd.value
    jordanTime = tmd.time

    if P.n <= limitSize

        tmd = @timed blockDiagonalize(P, verbose; complex = complex)
        blkD = tmd.value
        blkDTime = tmd.time

        if blkD === nothing
            ## Either random/rounding error, or complex numbers needed
            return nothing
        end

        ## solve with solver of choice
        m = nothing
        if verbose
            m = Model(CSDP.Optimizer)
        else
            m = Model(optimizer_with_attributes(CSDP.Optimizer, "MSK_IPAR_LOG" => 0))
        end

        ## >= 0 because the SDP-matrices should be entry-wise nonnegative
        x = @variable(m, x[1:P.n] >= 0)

        PMat = hcat([sparse(vec(P.P .== i)) for i = 1:P.n]...)

        ## Reduce the number of constraints
        newConstraints = Float64.(hcat(A * PMat, b))
        newConstraints = sparse(svd(Matrix(newConstraints)').U[:, 1:rank(newConstraints)]')
        droptol!(newConstraints, 1e-8)

        newA = newConstraints[:, 1:end-1]
        newB = newConstraints[:, end]
        newC = C' * PMat

        @constraint(m, newA * x .== newB)
        @objective(m, objSense, newC * x)

        for i = 1:length(blkD[1])
            blkExpr =
                x[1] .* (
                    complex ?
                    [
                        real(blkD[2][1][i]) -imag(blkD[2][1][i])
                        imag(blkD[2][1][i]) real(blkD[2][1][i])
                    ] :
                    blkD[2][1][i]
                )
            for j = 2:P.n
                add_to_expression!.(
                    blkExpr,
                    x[j] .* (
                        complex ?
                        [
                            real(blkD[2][j][i]) -imag(blkD[2][j][i])
                            imag(blkD[2][j][i]) real(blkD[2][j][i])
                        ] :
                        blkD[2][j][i]
                    ),
                )
            end
            if size(blkExpr, 1) > 1
                @constraint(m, blkExpr in PSDCone())
            else
                @constraint(m, blkExpr .>= 0)
            end
        end

        tmd = @timed optimize!(m)
        optTime = tmd.time

        if Int64(termination_status(m)) != 1
            @show termination_status(m)
            @error("Solve error.")
        end
        return (
            jTime = jordanTime,
            blkTime = blkDTime,
            solveTime = optTime,
            optVal = newC * value.(x),
            blkSize = blkD[1],
            originalSize = size(P.P, 1),
            newSize = P.n
        )
    end
    return (
        jTime = jordanTime,
        blkTime = 0,
        solveTime = 0,
        optVal = 0,
        blkSize = 0,
        originalSize = size(P.P, 1),
        newSize = P.n
    )

end


text:
module Caesar

# import RoME: initfg # collision on RoME.initfg() since no parameters are given in both RoME and Caesar
import Distributions: Normal
import RoME: getRangeKDEMax2D
import IncrementalInference: getSample, initfg
import DistributedFactorGraphs: getManifold

using Reexport
using Requires
using Dates

@reexport using RoME
@reexport using IncrementalInference
@reexport using KernelDensityEstimate
@reexport using Distributions

using
  Pkg,
  DelimitedFiles,
  Distributed,
  Statistics,
  LinearAlgebra,
  IncrementalInference,
  TransformUtils,
  CoordinateTransformations,
  Rotations,
  JSON,
  JSON2,
  FileIO,
  DataStructures,
  ProgressMeter,
  ImageMagick, # figure out who else is using this and move to requires
  ImageCore,
  DocStringExtensions,
  Unmarshal,
  YAML,
  FFTW,
  TimeZones,
  TensorCast

export
  GenericInSituSystem,  # insitu components
  makeGenericInSituSys,
  InSituSystem,
  makeInSituSys,
  triggerPose,
  poseTrigAndAdd!,
  processTreeTrackersUpdates!,
  advOdoByRules,
  SLAMWrapper,

  # servers
  tcpStringSLAMServer,
  tcpStringBRTrackingServer,

  # user functions
  identitypose6fg,
  projectrbe,
  hasval,

  # Robot Utils
  getRangeKDEMax2D,

  # sas-slam
  CBFFilterConfig,
  CZTFilter,
  prepCZTFilter,
  getCBFFilter2Dsize,
  constructCBFFilter2D!,
  CBF2D_DelaySum!,
  MatchedFilter,
  SASBearing2D,
  PackedSASBearing2D,
  compare,
  SASDebug,
  reset!,
  prepMF,
  loadConfigFile,
  prepareSAS2DFactor,
  wrapRad,
  phaseShiftSingle!,
  liebf!,
  SASDebug

const NothingUnion{T} = Union{Nothing, T}

include("BearingRangeTrackingServer.jl")

include("SlamServer.jl")
include("DataUtils.jl")
include("UserFunctions.jl")

# Configuration
include("config/CaesarConfig.jl")

include("Deprecated.jl")

# Multisession operation
# include("attic/multisession/Multisession.jl")

# SAS-SLAM
include("beamforming/czt.jl")
include("beamforming/CBF.jl")
include("beamforming/MatchedFilter.jl")
include("beamforming/SASBearing2D.jl")
include("beamforming/SASUtils.jl")

# conditional loading for ROS
function __init__()
  # ZMQ server and endpoints
  @require ZMQ="c2297ded-f4af-51ae-bb23-16f91089e4e1" include("zmq/ZmqCaesar.jl")
  @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
    @info "Loading Caesar PyCall specific utilities (using PyCall)."
    @eval using .PyCall
    @require RobotOS="22415677-39a4-5241-a37a-00beabbbdae8" include("ros/Utils/RosbagSubscriber.jl")
  end
  @require AprilTags="f0fec3d5-a81e-5a6a-8c28-d2b34f3659de" begin
    include("images/apriltags.jl")
    @require ImageDraw="4381153b-2b60-58ae-a1ba-fd683676385f" include("images/AprilTagDrawingTools.jl")
  end
  @require ImageMagick="6218d12a-5da1-5696-b52f-db25d2ecc6d1" include("images/imagedata.jl")
  @require Images="916415d5-f1e6-5110-898d-aaa5f9f070e0" include("images/images.jl")
end

end


text:
struct LwplsrdaAgg
    X::Array{Float64}
    y::AbstractMatrix
    fm
    metric::String
    h::Real
    k::Int
    nlv::String
    tol::Real
    verbose::Bool
    lev::AbstractVector
    ni::AbstractVector
end

"""
    lwplsrda_agg(X, y; nlvdis, metric, h, k, nlv,
        tol = 1e-4, verbose = false)
Aggregation of kNN-LWPLSR-DA models with different numbers of LVs.
* `X` : X-data.
* `y` : y-data (class membership).
* `nlvdis` : Number of latent variables (LVs) to consider in the global PLS
    used for the dimension reduction before calculating the dissimilarities.
    If `nlvdis = 0`, there is no dimension reduction.
* `metric` : Type of dissimilarity used to select the neighbors.
    Possible values are "eucl" (default; Euclidean distance)
    and "mahal" (Mahalanobis distance).
* `h` : A scalar defining the shape of the weight function. Lower is h,
    sharper is the function. See function `wdist`.
* `k` : The number of nearest neighbors to select for each observation to predict.
* `nlv` : A character string such as "5:20" defining the range of the numbers of LVs
    to consider ("5:20": the predictions of models with nb LVS = 5, 6, ..., 20
    are averaged). Syntax such as "10" is also allowed ("10": correponds to
    the single model with 10 LVs).
* `tol` : For stabilization when very close neighbors.
* `verbose` : If true, fitting information are printed.

Ensemblist method where the predictions are calculated by "averaging"
the predictions of a set of models built with different numbers of
latent variables (LVs).

For instance, if argument `nlv` is set to `nlv = "5:10"`, the prediction for
a new observation is the most occurent class within the predictions
returned by the models with 5 LVS, 6 LVs, ... 10 LVs, respectively.
"""
function lwplsrda_agg(X, y; nlvdis, metric, h, k, nlv,
    tol = 1e-4, verbose = false)
    X = ensure_mat(X)
    y = ensure_mat(y)
    ztab = tab(y)
    if nlvdis == 0
        fm = nothing
    else
        fm = plskern(X, dummy(y).Y; nlv = nlvdis)
    end
    LwplsrdaAgg(X, y, fm, metric, h, k, nlv, tol, verbose,
        ztab.keys, ztab.vals)
end

"""
    predict(object::LwplsrdaAgg, X)
Compute the y-predictions from the fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::LwplsrdaAgg, X)
    X = ensure_mat(X)
    m = size(X, 1)
    ### Getknn
    if isnothing(object.fm)
        res = getknn(object.X, X; k = object.k, metric = object.metric)
    else
        Tnew = transform(object.fm, X)
        res = getknn(object.fm.T, Tnew; k = object.k, metric = object.metric)
    end
    listw = copy(res.d)
    for i = 1:m
        w = wdist(res.d[i]; h = object.h)
        w[w .< object.tol] .= object.tol
        listw[i] = w
    end
    ### End
    pred = locw(object.X, object.y, X;
        listnn = res.ind, listw = listw, fun = plsrda_agg, nlv = object.nlv,
        verbose = object.verbose).pred
    (pred = pred, listnn = res.ind, listd = res.d, listw = listw)
end


text:
"""
Golden-section Search (https://en.wikipedia.org/wiki/Golden-section_search)
"""
module GoldenSectionSearch

#=
I don't know how to prove the optimality of golden-section search. The explanation on Wikipedia is
not plausible. Referring to the diagram on Wikipedia, let $s = x_2 - x_1, t = x_3 - x_4$, and
without loss of generality, suppose that `x_3 - x_1 = 1`. Also suppose the minimum distributes
uniformly in $(a, b)$. The expected interval length after evaluating $f(x_2)$ and $f(x_4)$ is

    (s + 1/2 * (1 - s - t)) * (1 - t) + (t + 1/2 * (1 - s - t)) * (1 - s)
    s.t. 0 < s < 1, 0 < t < 1 and s + t < 1

The new interval length reaches its minimum when $s = t = 1/2 - ϵ$. However, this setting is greedy
as the search is carried out interatively so it doesn't necessarily lead to the minimum after N
steps. Also note that when $x_2$ is fixed, let $t = min(x_2 - x_1, x_3 - x_2)$, the greedy solution
of $s$ is

$$
s = t + 1/2,  if 0 < t <= 1/4
s = 1 - t,  if 1/4 < t <= 1/2
$$

and the corresponding minimum length of interval is

$$
7 / 8 - t,  if 0 < t <= 1/4
2 t^2 - 2 t + 1,  if 1/4 < t <= 1/2
$$

With this idea one can possibly derive the optimal search scheme but I get stuck here.
=#

const GOLDEN_RATIO = (sqrt(5.0) + 1.0) / 2.0  # ϕ
const GOLDEN_RATIO_INV = GOLDEN_RATIO - 1.0  # 1 / ϕ, or Φ
const GOLDEN_RATIO_INV2 = 1 - GOLDEN_RATIO_INV  # 1 / ϕ²

"""
Find the minimum of a unimodal function `f`.

# Pre-conditions

- The minimum lies within the interval `(a, b)`.
"""
function golden_section_search_minimum(f::Function, a::T, b::T)::T where T <: Number
    tolerance = eps(T)
    # tolerance = 1e-6
    @assert a < b
    l = b - a
    n = Int(ceil(log(l / tolerance) / log(GOLDEN_RATIO)))

    ll = l * GOLDEN_RATIO_INV2
    aa = a + ll
    bb = b - ll
    faa = f(aa)
    fbb = f(bb)

    for i in 1:n
        if faa > fbb
            a = aa
            aa = bb; faa = fbb

            # The following code preserves precision better than
            #
            #     bb = b - (aa - a); fbb = f(bb)
            #
            # Also, if computing the new interval using `+`/`-` you have to break the loop when
            # `a <= aa <= bb <= b` is violated.
            l *= GOLDEN_RATIO_INV
            bb = a + l * GOLDEN_RATIO_INV
            fbb = f(bb)
        else
            b = bb
            bb = aa; fbb = faa
            l *= GOLDEN_RATIO_INV

            # It seems better than
            #
            #     aa = b - l * GOLDEN_RATIO_INV
            #
            # but I don't know the reason.
            aa = a + l * GOLDEN_RATIO_INV2

            faa = f(aa)
        end
    end

    if faa > fbb
        (aa + b) / 2.0
    else
        (a + bb) / 2.0
    end
end

# Smoke tests
@assert golden_section_search_minimum(x -> (x - 2.0) ^ 2, 0.0, 10.0) ≈ 2.0f0
@assert golden_section_search_minimum(x -> x^2 - x - 1.0, -10.0, 10.0) ≈ 0.5f0

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
include("euler/euler.jl")
using .Primes: prime_sieve
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    sieve = prime_sieve(n)
    prime_numbers = findall(sieve)
    result = 0
    for prime ∈ prime_numbers
        str_prime = string(prime)
        result += 1
        for _ ∈ 1:length(str_prime)
            str_prime = string(str_prime[2:end], str_prime[1])
            if !sieve[parse(Int, str_prime)]
                result -= 1
                break
            end
        end
    end
    return result
end

compute(100)

compute(1_000_000)

@benchmark compute(1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
include("euler/euler.jl")
using .Math: fast_min
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(grid::Array{Int,2})::Int
    edge = size(grid)[1]
    for i ∈ edge - 1:-1:1
        grid[edge, i] += grid[edge, i + 1]
        grid[i, edge] += grid[i + 1, edge]
    end
    for i ∈ edge - 1:-1:1 for j ∈ edge - 1:-1:1 grid[i, j] += fast_min(grid[i + 1, j], grid[i, j + 1]) end end
    return grid[1, 1]
end

matrix = [131 673 234 103 18;
          201 96 342 965 150;
          630 803 746 422 111;
          537 699 497 121 956;
          805 732 524 37 331]

path = "problems/0081/p081_matrix.txt"

compute(matrix)

compute(parse.(Int, hcat(split.(split(read(path, String), "\n"), ",")...)))

@benchmark compute(parse.(Int, hcat(split.(split(read(path, String), "\n"), ",")...)))


text:
#=============
Filter methods
==============#

#=-
A. Entry points
-=#

function Base.filter!(ss::LinearStateSpace{KFP} where KFP<:KFParms{A}, Y::B) where {A<:AbstractArray, B<:AbstractArray}
    (T, p, m, r, offset) = Filthy.checkfilterdataconsistency(ss, Y)
    #rY = reshape(Y, (p, T))
    Filthy.resizecontainer!(ss, T+offset)
    for j in 1:size(Y, 2)
        onlinefilter_set!(ss, Y[:,j], offset)
    end
end

function Base.filter!(ss::LinearStateSpace{KFP} where KFP<:KFParms{A}, Y::Vector) where A<:AbstractFloat
    (T, p, m, r, offset) = Filthy.checkfilterdataconsistency(ss, Y)
    #rY = reshape(Y, (p, T))
    Filthy.resizecontainer!(ss, T+offset)
    for j in 1:length(Y)
        onlinefilter_set!(ss, Y[j], offset)
    end
    nothing
end

function onlinefilter_set!(ss::LinearStateSpace, y, offset::Int)
    hasdiffended = first(ss.f.hasdiffended)::Bool
    onlinefilterstep_set!(ss, y, offset, Val{!hasdiffended})
    ss.t[] += 1
    nothing
end

function onlinefilter!(ss::LinearStateSpace, y)
    hasdiffended = first(ss.f.hasdiffended)::Bool
    onlinefilterstep!(ss, y, Val{!hasdiffended})
    ss.t[] += 1
    nothing
end

function onlinefilterstep!(ss::LinearStateSpace, y, v::Type{Val{false}})
    a, P = currentstate(ss, Val{false})
    @unpack Z, H, T, R, Q = ss.p
    a′, P′, v, F, F⁻¹ = filterstep(a, P, Z, H, T, R, Q, y)
    storepush!(ss, a′, P′, v, F, F⁻¹, y)
end

function onlinefilterstep!(ss::LinearStateSpace, y, v::Type{Val{true}})
    a, P, Pinf, Pstar = currentstate(ss, Val{true})
    @unpack Z, H, T, R, Q = ss.p
    a′, P′, v, F, F⁻¹, Pinf′, Pstar′, flag = filterstep(a, P, Pinf, Pstar, Z, H, T, R, Q, y)
    ss.f.hasdiffended[1] = norm(Pinf′) <= 1e-08
    storepush!(ss, a′, P′, v, F, F⁻¹, Pinf′, Pstar′, y, flag)
end

function onlinefilterstep_set!(ss::LinearStateSpace, y, offset, v::Type{Val{false}})
    a, P = currentstate(ss, Val{false})
    @unpack Z, H, T, R, Q = ss.p
    a′, P′, v, F, F⁻¹ = filterstep(a, P, Z, H, T, R, Q, y)
    #@show a′, P′, v, F, F⁻¹
    storeset!(ss, a′, P′, v, F, F⁻¹, y, offset)
    updatelik!(ss, v, F, F⁻¹)
end

function onlinefilterstep_set!(ss::LinearStateSpace, y, offset, v::Type{Val{true}})
    a, P, Pinf, Pstar = currentstate(ss, Val{true})
    @unpack Z, H, T, R, Q = ss.p
    a′, P′, v, F, F⁻¹, Pinf′, Pstar′, flag = filterstep(a, P, Pinf, Pstar, Z, H, T, R, Q, y)
    ss.f.hasdiffended[1] = norm(Pinf′) <= 1e-06
    storeset!(ss, a′, P′, v, F, F⁻¹, Pinf′, Pstar′, y, offset, flag)
    updatelik!(ss, v, F, F⁻¹, flag)
end

#=-
B. Actual steps
-=#
function filterstep(a, P, Pinf, Pstar, Z, H, T, R, Q, y)
    Finf = Z*Pinf*Z'
    (Fflag, Finfinv) = Filthy.safeinverse(Finf)
    diffusefilterstep(a, P, Pinf, Pstar, Finf, Finfinv, Z, H, T, R, Q, y, Val{Fflag})
end

function filterstep(a::AA, P::PP, Z, H, T, R, Q, y) where {AA, PP}
    ## Non-diffuse step
    v = y .- Z*a
    F = Z*P*Z' .+ H
    Finv = inv(F)
    K = T*P*Z'*Finv
    L = T.-K*Z
    a′ = T*a .+ K*v
    P′ = T*P*L' .+ R*Q*R'
    (a′, P′, v, F, Finv)
end

function diffusefilterstep(a, P, Pinf, Pstar, Finf, Finfinv, Z, H, T, R, Q, y, ::Type{Val{1}})
    ## Finf is invertible
    v     = y .- Z*a
    Kinf  = T*Pinf*Z'*Finfinv
    Linf  = T - Kinf*Z
    Fstar = Z*Pstar*Z' + H
    Kstar = (T*Pstar*Z' + Kinf*Fstar)*Finfinv
    a′     = T*a + Kinf*v
    Pinf′  = T*Pinf*Linf'
    Pstar′ = T*Pstar*Linf' + Kinf*Finf*Kstar' + R*Q*R'
    P′     = Pstar′
    (a′, P′, v, Finf, Finfinv, Pinf′, Pstar′, 1)
end

function diffusefilterstep(a, P, Pinf, Pstar, Finf, Finfinv, Z, H, T, R, Q, y, ::Type{Val{0}})
    ## Finf is not > 0
    v     = y - Z*a
    Fstar = Z*Pstar*Z' + H
    Fstarinv = inv(Fstar)
    Kstar = T*Pstar*Z'*Fstarinv
    Lstar = T - Kstar*Z
    a′     = T*a + Kstar*v
    Pinf′  = T*Pinf*T'
    Pstar′ = T*Pstar*Lstar' + R*Q*R'
    P′     = Pstar′
    (a′, P′, v, Fstar, Fstarinv, Pinf′, Pstar′, 0)
end

# #=-
# C. In-place filter step
# -=#

# # Standard filtering equation
# function filterstep!(a′, P′, v, vp, C, M, ZM, TM, K, KZ, KZp, TPL, L, a::AA, P::PP, Z, H, T, R, Q, y) where {AA, PP}
#     ## v  is (p x 1)
#     ## v′ is (p x 1)
#     ## C is (p x p)
#     ## M  is (m x p)
#     ## ZM is (p x p)
#     ## TM is (m x p)
#     ## K  is (m x p)
#     ## KZ is (m x m)
#     ## KZ′ is (m x m)
#     ## TPL is (m x m)
#     ## L  is (m x m)

#     ## [v]
#     A_mul_B!(vp, Z, a)
#     @inbounds @simd for i in eachindex(v)
#         v[i] = y[i] - vp[i]
#     end
#     ## [M]
#     A_mul_Bt!(M, P, Z)
#     ## [ZPZ'] -> ZM
#     A_mul_B!(ZM, Z, M)
#     ## [ZPZ' + H] -> C
#     @inbounds @simd for i in eachindex(ZM)
#         C[i] = ZM[i] + H[i]
#     end
#     Cinv = inv(factorize(C))
#     ## [T*M] -> TM
#     A_mul_B!(TM, T, M)
#     ## [T*M*Cinv] -> K
#     A_mul_B!(K, TM, Cinv)
#     ## [KZ] -> KZ
#     A_mul_B!(KZ, K, Z)
#     ## [T-KZ] -> L
#     @inbounds @simd for i in eachindex(L)
#         L[i] = T[i] - KZ[i]
#     end
#     ## [T*a] -> Ta
#     A_mul_B!(a′, T, a)
#     A_mul_B!(vp, K, v)
#     @inbounds @simd for i in eachindex(a′)
#         a′[i] = a′[i] + vp[i]
#     end
#     ## [TPL']
#     A_mul_B!(KZp, T, P)
#     A_mul_Bt!(KZ, KZp, L)
#     ## [RQR]
#     A_mul_B!(K, R, Q)
#     A_mul_Bt!(TM, K, R)
#     ## [TPL'+RQR] -> P
#     @inbounds @simd for i in eachindex(P′)
#         P′[i] = KZ[i] + TM[i]
#     end
#     (a′, P′, v, C, Cinv)
# end

# ## diffuse filter version
# ## F>0 and invertible
# function filterstep!(a′, P′, v, vp, C, M, ZM, TM, K, KZ, KZp, L, a::AA, P::PP, Pinf, Pstar, Finf, Finfinv, Z, H, T, R, Q, y, ::Type{Val{1}}) where {AA, PP}
#     ## v  is (p x 1)
#     ## v′ is (p x 1)
#     ## C is (p x p)
#     ## M  is (m x p)
#     ## ZM is (p x p)
#     ## TM is (m x p)
#     ## K  is (m x p)
#     ## KZ is (m x m)
#     ## KZ′ is (m x m)
#     ## L  is (m x m)

#     ## [v]
#     A_mul_B!(vp, Z, a)
#     @inbounds @simd for i in eachindex(v)
#         v[i] = y[i] - vp[i]
#     end
#     ## [M]
#     A_mul_Bt!(M, Pstar, Z)
#     ## [ZPZ'] -> ZM
#     A_mul_B!(ZM, Z, M)
#     # [ZPZ' + H] -> C (Fstar)
#     @inbounds @simd for i in eachindex(ZM)
#         C[i] = ZM[i] + H[i]
#     end
#     ## [T*M] -> TM
#     A_mul_B!(TM, T, M)
#     ## [T*M*Cinv] -> K (Kinf)
#     A_mul_B!(K, TM, Cinv)
#     ## Kstar
#     A_mul_B!(KZ, K, C)
#     @inbounds @simd for i in eachindex(Cinv)
#         Cinv[i] = TM[i] + KZ[i]
#     end
#     A_mul_B!(ZM, Cinv, Finfinv)
#     ## [KZ] -> KZ
#     A_mul_B!(KZ, K, Z)
#     ## [T-KZ] -> L
#     @inbounds @simd for i in eachindex(L)
#         L[i] = T[i] - KZ[i]
#     end

#     ## [T*a] -> Ta
#     A_mul_B!(a′, T, a)
#     A_mul_B!(vp, K, v)
#     @inbounds @simd for i in eachindex(a′)
#         a′[i] = a′[i] + vp[i]
#     end
#     ## [TPL']
#     A_mul_B!(KZp, T, Pinf)
#     A_mul_Bt!(KZ, KZp, L)  ## Pinf'
#     ## [RQR]
#     A_mul_B!(K, R, Q)
#     A_mul_Bt!(TM, K, R)
#     ## [TPstarL'+RQR] ->
#     A_mul_B!(KZp, T, Pstar)
#     A_mul_Bt!(TPL, KZp, L)

#     @inbounds @simd for i in eachindex(P′)
#         P′[i] = KZ[i] + TPL[i] + TM[i] ## Pstar'
#     end

#     (a′, P′, v, Fstar, Fstarinv, Kz, P′, 0)
# end

# function filterstep!(a′, P′, v, vp, C, M, ZM, TM, K, KZ, KZp, L, a::AA, P::PP, Pinf, Pstar, Finf, Finfinv, Z, H, T, R, Q, y, ::Type{Val{0}}) where {AA, PP}
#     ## v  is (p x 1)
#     ## v′ is (p x 1)
#     ## C is (p x p)
#     ## M  is (m x p)
#     ## ZM is (p x p)
#     ## TM is (m x p)
#     ## K  is (m x p)
#     ## KZ is (m x m)
#     ## KZ′ is (m x m)
#     ## TPL is (m x m)
#     ## L  is (m x m)

#     ## [v]
#     A_mul_B!(vp, Z, a)
#     @inbounds @simd for i in eachindex(v)
#         v[i] = y[i] - vp[i]
#     end
#     ## [M]
#     A_mul_Bt!(M, Pstar, Z)
#     ## [ZPZ'] -> ZM
#     A_mul_B!(ZM, Z, M)
#     # [ZPZ' + H] -> C (Fstar)
#     @inbounds @simd for i in eachindex(ZM)
#         C[i] = ZM[i] + H[i]
#     end
#     Fstarinv = inv(factor(C))

#     ## Kstar
#     A_mul_B!(TM, T, M)
#     ## [T*M*Cinv] -> K
#     A_mul_B!(K, TM, Fstarinv)

#     ## Lstar

#     ## [T*a] -> Ta
#     A_mul_B!(a′, T, a)
#     A_mul_B!(vp, K, v)
#     @inbounds @simd for i in eachindex(a′)
#         a′[i] = a′[i] + vp[i]
#     end
#     ## [TPL']
#     A_mul_B!(KZp, T, Pinf)
#     A_mul_Bt!(KZ, KZp, L)  ## Pinf'
#     ## [RQR]
#     A_mul_B!(K, R, Q)
#     A_mul_Bt!(TM, K, R)
#     ## [TPstarL'+RQR] ->
#     A_mul_B!(KZp, T, Pstar)
#     A_mul_Bt!(TPL, KZp, L)

#     @inbounds @simd for i in eachindex(P′)
#         P′[i] = KZ[i] + TPL[i] + TM[i] ## Pstar'
#     end

#     (a′, P′, v, Fstar, Fstarinv, Kz, P′, 0)
# end

#=-
D. Storage
-=#

function storepush!(ss, a::T, P, v, F, F⁻¹, Pinf, Pstar, y, flag) where T<:Real
    push!(ss.f.a, a...)
    push!(ss.f.P, P...)
    push!(ss.c.F⁻¹, F⁻¹)
    push!(ss.c.F, F)
    push!(ss.c.Y, y)
    push!(ss.f.Pinf, Pinf)
    push!(ss.f.Pstar, Pstar)
    ss.f.d[1] += 1
    push!(ss.f.flagF, flag)
    nothing
end

function storepush!(ss, a::T, P, v, F, F⁻¹, Pinf, Pstar, y, flag) where T<:AbstractArray
    push!(ss.f.a, a)
    push!(ss.f.P, P)
    push!(ss.c.F⁻¹, F⁻¹)
    push!(ss.c.F, F)
    push!(ss.c.Y, y)
    push!(ss.f.Pinf, Pinf)
    push!(ss.f.Pstar, Pstar)
    ss.f.d[1] += 1
    push!(ss.f.flagF, flag)
    nothing
end

function storepush!(ss, a::T, P, v, F, F⁻¹, y) where T<:Real
    push!(ss.f.a, a...)
    push!(ss.f.P, P...)
    push!(ss.c.F⁻¹, F⁻¹)
    push!(ss.c.F, F)
    push!(ss.c.Y, y)
    nothing
end

function storepush!(ss, a::T, P, v, F, F⁻¹, y) where T<:AbstractArray
    push!(ss.f.a, a)
    push!(ss.f.P, P)
    push!(ss.c.F⁻¹, F⁻¹)
    push!(ss.c.F, F)
    push!(ss.c.Y, y)
    nothing
end

function storeset!(ss, a, P, v, F, F⁻¹, Pinf, Pstar, y, offset, flag)
    idx = ss.t[] + offset
    ss.f.a[idx]    = a
    ss.f.P[idx]    = P
    ss.c.F⁻¹[idx]  = F⁻¹
    ss.c.F[idx]    = F
    ss.c.Y[idx]    = y
    push!(ss.f.Pinf, Pinf)
    push!(ss.f.Pstar, Pstar)
    ss.f.d[1] += 1
    push!(ss.f.flagF, flag)
    nothing
end

function storeset!(ss, a, P, v, F, F⁻¹, y, offset)
    idx = ss.t[] + offset
    ss.f.a[idx] = a
    ss.f.P[idx] = P
    ss.c.F⁻¹[idx] = F⁻¹
    ss.c.F[idx] = F
    ss.c.Y[idx] = y
    # s = loglikpart(v, F, F⁻¹)
    # p, _ = size(ss.p)
    # ss.loglik[1] += - p/2*log(2π) - s
    nothing
end


text:
# NSWindow Styles
const NSBorderlessWindowMask             = 0
const NSTitledWindowMask                 = 1 << 0
const NSClosableWindowMask               = 1 << 1
const NSMiniaturizableWindowMask         = 1 << 2
const NSResizableWindowMask              = 1 << 3
const NSTexturedBackgroundWindowMask     = 1 << 8
const NSUnifiedTitleAndToolbarWindowMask = 1 << 12

# NSWindow Backing
const NSBackingStoreRetained      = 0
const NSBackingStoreNonretained   = 1
const NSBackingStoreBuffered      = 2

# NSWindow collection behaviours
const NSWindowCollectionBehaviorDefault  = 0
const NSWindowCollectionBehaviorCanJoinAllSpaces  = 1 << 0
const NSWindowCollectionBehaviorMoveToActiveSpace  = 1 << 1
const NSWindowCollectionBehaviorManaged  = 1 << 2
const NSWindowCollectionBehaviorTransient  = 1 << 3
const NSWindowCollectionBehaviorStationary  = 1 << 4
const NSWindowCollectionBehaviorParticipatesInCycle  = 1 << 5
const NSWindowCollectionBehaviorIgnoresCycle  = 1 << 6
const NSWindowCollectionBehaviorFullScreenPrimary  = 1 << 7
const NSWindowCollectionBehaviorFullScreenAuxiliary  = 1 << 8

# NSApplication Activation policies
const NSApplicationActivationPolicyRegular    = 0
const NSApplicationActivationPolicyAccessory  = 1
const NSApplicationActivationPolicyProhibited = 2

# NSApplication presentation styles
const NSApplicationPresentationDefault                    = 0
const NSApplicationPresentationAutoHideDock               = 1 <<  0
const NSApplicationPresentationHideDock                   = 1 <<  1
const NSApplicationPresentationAutoHideMenuBar            = 1 <<  2
const NSApplicationPresentationHideMenuBar                = 1 <<  3
const NSApplicationPresentationDisableAppleMenu           = 1 <<  4
const NSApplicationPresentationDisableProcessSwitching    = 1 <<  5
const NSApplicationPresentationDisableForceQuit           = 1 <<  6
const NSApplicationPresentationDisableSessionTermination  = 1 <<  7
const NSApplicationPresentationDisableHideApplication     = 1 <<  8
const NSApplicationPresentationDisableMenuBarTransparency = 1 <<  9
const NSApplicationPresentationFullScreen                 = 1 << 10
const NSApplicationPresentationAutoHideToolbar            = 1 << 11


text:
export model

"""
    model(which::Symbol; nbit::Int, V::Int, B::Int=4096, nlayer::Int=5)

predefined models, `which` should be one of :random, :u1, :su2.
* `nbit` is the system size (length of MPS),
* `V` is the number of virtual qubits,
* `B` is the batch size.
* `nlayer` is the number of layers in a block.
"""
model(which::Symbol, args...; kwargs...) = model(Val(which), args...; kwargs...)

include("general_circuit.jl")
include("u1_circuit.jl")
include("su2_circuit.jl")
include("twoqubit_circuit.jl")
include("cluster.jl")

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


text:
"""
    L(x, y, k) = -log(mean(k(xi, y) for xi in x))
"""
L(x, y, k) = -log(mean(k(xi, y) for xi in x))
"""
    KernelLoss(; loss = L, K = 50,
                 kerneltype = ModifiedGaussian,
                 gamma = .1,
                 heuristic = MedianHeuristic(),
                 prior = nothing,
                 distance = euclidean,
                 kernel = Kernel(type = kerneltype,
                                 heuristic = heuristic,
                                 gamma = gamma,
                                 distance = distance))

Constructs an structure, used to compare `K` samples generated from a model
to some data with a `kernel` and loss function `loss`. The loss function has the
signature `L(list_of_samples, data, kernel)`; by default it is the negative log
of the average distance between samples and data
(see [LikelihoodfreeInference.L](@ref)).
The `KernelLoss` can be used as a loss in a [`PointEstimator`](@ref).

The returned structure is callable.
# Example
```julia
model(x) = randn() .+ x
data = [2.]
k = KernelLoss()
k([1.], model, data)
```

`K` can be an integer or a function that returns integers, to implement a schedule.
# Example
```julia
using Statistics
model(x) = randn() .+ x
data = [2.]
k = KernelLoss(K = (n) -> n > 20 ? 10^3 : 10)
std(k([1.], model, data, 1) for _ in 1:100)   # large std, because K = 10
std(k([1.], model, data, 21) for _ in 1:100)  # small std, because K = 10^3
```

This `KernelLoss` is inspired by
[Bertl, J., Ewing, G., Kosiol, C., and Futschik, A., (2017)
Approximate maximum likelihood estimation for population genetic inference,
Statistical Applications in Genetics and Molecular Biology. 16:5-6
](http://dx.doi.org/10.1515/sagmb-2017-0016)
"""
struct KernelLoss{P,Tl,TK,Tk}
    loss::Tl
    K::TK
    kernel::Tk
    prior::P
end
function KernelLoss(; loss = L, K = 50,
                      kerneltype = ModifiedGaussian,
                      gamma = .1,
                      heuristic = MedianHeuristic(),
                      prior = nothing,
                      distance = euclidean,
                      kernel = Kernel(type = kerneltype,
                                      heuristic = heuristic,
                                      gamma = gamma,
                                      distance = distance))
    KernelLoss(loss, K, kernel, prior)
end
"""
    MMDLoss(; single_observation = false, kwargs...)

Returns a [`KernelLoss`](@ref) structure with the maximum mean discrepancy (MMD)
loss; to be used as a `loss` in a [`PointEstimator`](@ref). `kwargs` are passed
to the [`KernelLoss`](@ref) constructor.
If `single_observation = true` the data is expected to be a single vector of
numbers; otherwise a vector of vectors.

This `MMDLoss` is inspired by
[Briol, F.-X., Barp, A., Duncan, A. B. and Girolami, M., (2019),
Statistical Inference for Generative Models with Maximum Mean Discrepancy,
arXiv:1906.05944]
(https://arxiv.org/abs/1906.05944)
"""
function MMDLoss(; single_observation = false, kwargs...)
    loss = single_observation ? (S, data, k) -> mmd(S, [data], k) :
                                (S, data, k) -> mmd(vcat(S...), data, k)
    KernelLoss(; loss = loss, kwargs...)
end
computeK(K::Number, ::Any) = K
computeK(K::Function, n) = K(n)
computeK(K::Base.RefValue, ::Any) = K[]
function (k::KernelLoss)(θ, model, data, n = 1)
    S = [model(θ) for _ in 1:computeK(k.K, n)]
    n == 1 && update!(k.kernel, S, data)
    k.loss(S, data, k.kernel) - logpdf(k.prior, θ)
end
mutable struct QDLoss{D,EPS,Tk,P}
    distance::D
    ϵ::EPS
    K::Tk
    prior::P
    mins::Vector{Float64}
end
"""
    QDLoss(; K = 50, epsilon = 2/K, prior = nothing, distance = euclidean)

Constructs an structure, used to compare `K` samples generated from a model
to some data by computing the first epsilon-quantile of the `distance` between
samples and data.
The `QDLoss` can be used as a loss in a [`PointEstimator`](@ref).

The returned structure is callable.
# Example
```julia
model(x) = randn() .+ x
data = [2.]
q = QDLoss()
q([1.], model, data)
```

`K` can be an integer or a function that returns integers, to implement a schedule.
Similarly, `epsilon` can be a number or a function that returns numbers.
# Example
```julia
using Statistics
model(x) = randn() .+ x
data = [2.]
q = QDLoss(K = (n) -> n > 20 ? 10^3 : 10)
std(q([1.], model, data, 1) for _ in 1:100)   # large std, because K = 10
std(q([1.], model, data, 21) for _ in 1:100)  # small std, because K = 10^3
```
"""
function QDLoss(; K = 50, epsilon = 2/computeK(K, 1),
                  prior = nothing, distance = euclidean)
    QDLoss(distance, epsilon, K, prior, Float64[])
end
function lquantile(k::QDLoss{typeof(euclidean)}, S,
                   data::AbstractVector{<:AbstractVector}, n)
    length(k.mins) == 0 && (k.mins = zeros(length(data)))
    ds = [[k.distance(x[i], data[i]) for x in S] for i in eachindex(data)]
    sum(lquantile(ds[i], computeK(k.ϵ, n), Ref(k.mins, i)) for i in eachindex(data))
end
function lquantile(k, S, data, n)
    length(k.mins) == 0 && (k.mins = [0.])
    d = [k.distance(x, data) for x in S]
    lquantile(d, computeK(k.ϵ, n), Ref(k.mins, 1))
end
function lquantile(d, ϵ, min::Base.RefArray)
    q = quantile(d, ϵ)
    q == 0 && return -log(mean(d .== 0)/ϵ) + min[]
    lq = log(q)
    lq < min[] && (min[] = lq)
    lq
end
function (k::QDLoss)(θ, model, data, n = 1)
    S = [model(θ) for _ in 1:computeK(k.K, n)]
    lquantile(k, S, data, n) - 1/length(θ) * logpdf(k.prior, θ)
end
function objective(model, data, loss;
                   callback = (x, y) -> nothing)
    cb = hasmethod(callback, Tuple{Vector{<:Number}, Number}) ?
         callback : (x, y) -> callback()
    n = 0
    θ -> begin
        n += 1
        result = loss(θ, model, data, n)
        cb(θ, result)
        result
    end
end
struct EnergyLoss{TK}
    K::TK
end
"""
    EnergyLoss(K = 50)

Return an energy loss minization structure, that can be used as a `loss` in
[`PointEstimator`](@ref).

See e.g. section 6.1 of
[Székely, G. J. and Rizzo, M. L., (2013)
Energy statistics: A class of statistics based on distances,
Journal of Statistical Planning and Inference, 8:1249-1272]
(http://dx.doi.org/10.1016/j.jspi.2013.03.018).
"""
function EnergyLoss(; prior = nothing, K = 50)
    prior !== nothing && @warn "Priors are ignored in `EnergyLoss`."
    EnergyLoss(K)
end
function (e::EnergyLoss)(θ, model, data::T, n = 1) where T
    S = [model(θ) for _ in 1:computeK(e.K, n)]
    if T <: AbstractVector{<:AbstractVector}
        energydistance(vcat(S...), data)
    elseif T <: AbstractVector{<:Number}
        energydistance(S, [data])
    end
end
struct KLLoss{TK}
    K::TK
end
"""
    KLLoss(K = 50)

Return a KL-loss minization structure, that can be used as a `loss` in
[`PointEstimator`](@ref).

[Jiang, B., Wu, T.-Y. and Wong, W. H., (2018)
Approximate Bayesian Computation with Kullback-Leibler Divergence as Data Discrepancy,
Proceedings of Machine Learning Research,
(http://proceedings.mlr.press/v84/jiang18a.html).
"""
function KLLoss(; prior = nothing, K = 50)
    prior !== nothing && @warn "Priors are ignored in `KLLoss`."
    KLLoss(K)
end
function (e::KLLoss)(θ, model, data::T, n = 1) where T
    S = [model(θ) for _ in 1:computeK(e.K, n)]
    if T <: AbstractVector{<:AbstractVector}
        kldistance(vcat(S...), data)
    elseif T <: AbstractVector{<:Number}
        kldistance(S, data)
    end
end
function (e::KLLoss)(θ, model, data::AbstractVector{<:Number}, n = 1)
    S = [model(θ) for _ in 1:computeK(e.K, n)]
    energydistance(S, [data])
end
"""
    PointEstimator(optimizer, loss)
"""
struct PointEstimator{O,L} <: AbstractPointABC
    optimizer::O
    loss::L
end
"""
    PointEstimator(; prior = nothing, lower = nothing, upper = nothing,
                     optimizer = nlopt(prior),
                     losstype = QDLoss, kwargs...)

Creates a `PointEstimator` with given `prior`, `optimizer` and
`loss = losstype(; prior = prior, kwargs...)`. If the prior is `nothing` a
uniform prior is assumed. In this case the optimizer has to be initialized with
bounds, , e.g. `nlopt(lower, upper)` where `lower` and `upper` are arrays.

# Example
```julia
model(x) = randn() .+ x
data = [2.]
p = PointEstimator(lower = [-5], upper = [5], losstype = QDLoss, K = 10^3)
res = run!(p, model, data, maxfevals = 10^5)
```
"""
function PointEstimator(; prior = nothing, lower = nothing, upper = nothing,
                        optimizer = prior === nothing ? (lower === nothing || upper ===nothing) ? error("Please set `prior` or `lower` and `upper`") : nlopt(lower, upper) :
                                                          nlopt(prior),
                          losstype = QDLoss, kwargs...)
    PointEstimator(optimizer, losstype(; prior = prior, kwargs...))
end

function run!(p::PointEstimator, model, data;
              maxfevals = Inf, callback = () -> nothing, verbose = true)
    _optimize(p.optimizer, objective(model, data, p.loss, callback = callback),
              maxfevals = max(div(maxfevals, computeK(p.loss.K, 1)), 1),
              verbose = verbose)
end


text:
function testSMC(model::SMCModel, N::Int64, numTrials::Int64,
  fullOutput = false, essThreshold::Float64 = 2.0)
  smcio = SMCIO{model.particle, model.pScratch}(N, model.maxn, 1, fullOutput,
    essThreshold)
  println("Running SMC. N = $N, nthreads = 1, fullOutput = $fullOutput, essThreshold = $essThreshold")
  for i=1:numTrials
    @time smc!(model, smcio)
    left = max(length(smcio.logZhats)-10,1)
    right = length(smcio.logZhats)
    println(smcio.logZhats[left:right])
    println(sign.(smcio.Vhat1s[left:right]) .*
      sqrt.(abs.(smcio.Vhat1s[left:right])))
  end
end

function testSMCParallel(model::SMCModel, N::Int64, numTrials::Int64,
  fullOutput = false, essThreshold::Float64 = 2.0)
  nthreads = Threads.nthreads()
  smcio = SMCIO{model.particle, model.pScratch}(N, model.maxn, nthreads,
    fullOutput, essThreshold)
  println("Running SMC. N = $N, nthreads = $nthreads, fullOutput = $fullOutput, essThreshold = $essThreshold")
  for i = 1:numTrials
    @time smc!(model, smcio)
    left = max(length(smcio.logZhats)-10,1)
    right = length(smcio.logZhats)
    println(smcio.logZhats[left:right])
    println(sign.(smcio.Vhat1s[left:right]) .*
      sqrt.(abs.(smcio.Vhat1s[left:right])))
  end
end


text:
import MittagLefflerFunctions.Quadrature.QSumP
import MittagLefflerFunctions.Quadrature.QSumH
using PyPlot

figure(1)
marker = ("^", "o", "s", "v")
N = (5, 10, 15, 20)
for k = 1:4
    qs = QSumP(Float64, N[k])
    plot(real.(qs.w), imag.(qs.w), marker[k], label="N = $(N[k])")
end
legend()
grid(true)
axis((-50, 10, -5, 50))
xlabel(L"$\Re w$", fontsize=12)
ylabel(L"$\Im w$", fontsize=12)
savefig("figure04L.pdf")

figure(2)
for k = 1:4
    qs = QSumH(Float64, N[k])
    plot(real.(qs.w), imag.(qs.w), marker[k], label="N = $(N[k])")
end
legend()
grid(true)
axis((-50, 10,-5, 50))
xlabel(L"$\Re w$", fontsize=12)
ylabel(L"$\Im w$", fontsize=12)
savefig("figure04R.pdf")

figure(4)
for k = 1:4
    qs = QSumH(Float64, N[k])
    semilogy(collect(0:N[k]), abs.(qs.C), marker[k], label="N = $(N[k])")
end
legend()
grid(true)
xlabel(L"$n$", fontsize=12)
ylabel(L"$|C_n|$", fontsize=12)
xylims = axis()
xticks(collect(0:4:20))
savefig("figure05R.pdf")

figure(3)
for k = 1:4
    qs = QSumP(Float64, N[k])
    semilogy(collect(0:N[k]), abs.(qs.C), marker[k], label="N = $(N[k])")
end
legend()
grid(true)
#axis((-1.0, 21.0, 7.74879598668249e-22, 6034.9888883983585))
axis(xylims)
xlabel(L"$n$", fontsize=12)
ylabel(L"$|C_n|$", fontsize=12)
xticks(collect(0:4:20))
savefig("figure05L.pdf")


text:
# coordinates.jl
# Maps having to do with coordinate transforms.

"""
A Cartesion to Polar map. First dimension is interpreted as radial distance, second as an angle.
The unit circle is mapped to a square [-1,1]x[-1,1]

"""
struct CartToPolarMap{T} <: AbstractMap{SVector{2,T}, SVector{2,T}}
end

(m::CartToPolarMap)(x) = applymap(m, x)

applymap(map::CartToPolarMap, x::SVector{2,T}) where {T} = SVector{2,T}(sqrt(x[1]^2+x[2]^2)*2-1, atan(x[2],x[1])/pi)

inv(map::CartToPolarMap{T}) where {T} = PolarToCartMap{T}()

"""
A Polar to Cartesian map. The angle is mapped to the second dimension, radius to the first.
A square [-1,1]x[-1,1] is mapped to the unit circle
"""
struct PolarToCartMap{T} <: AbstractMap{SVector{2,T}, SVector{2,T}}
end

(m::PolarToCartMap)(x) = applymap(m, x)

applymap(map::PolarToCartMap, x::SVector{2,T}) where {T} = SVector{2,T}((x[1]+1)/2*cos(pi*x[2]), (x[1]+1)/2*sin(pi*x[2]))

inv(map::PolarToCartMap{T}) where {T} = CartToPolarMap{T}()


text:
######################################################################
# Pia Portal - Pia's standard discovery protocol.
# -----
# Licensed under Apache 2.0
module PortalProtocol
using ..Pia
using ExtraFun

const REGISTRY_STORE_VERSION = 0x1

struct PortalListing{S}
    address::Address{S}
    endpoint::Optional{Endpoint}
    last_seen::Float64
    scores::Dict{Symbol, Float32}
end
PortalListing(address::Address{S}) where S = PortalListing{S}(address, nothing, 0, Dict())

export PiaPortal
mutable struct PiaPortal{S} <: DiscoveryProtocol
    basedir::String
    lock::ReentrantLock

    function PiaPortal{S}(basedir, lock) where S
        @assert S isa Integer
        UInt32(S) # Ensures S is a *value* which can be re/stored
        new{S}(basedir, lock)
    end
end
PiaPortal(S::Integer, basedir::AbstractString = "./.Pia.portal") = PiaPortal{S}(basedir, ReentrantLock())

function resolve(portal::PiaPortal, transport::TransportProtocol, addr::Address)::Endpoint
    try
        return getpeer(portal, addr).endpoint::Endpoint
    catch; end

    error("not implemented")
end

function Pia.connect(transport::TransportProtocol, portal::PiaPortal, app::AbstractApplication)
    # TODO: Attempt to connect to rated peers

    return Pia.connect(transport, app)
end

end # module PortalProtocol


text:
using Documenter, NeuroCore

makedocs(;
    modules=[NeuroCore],
    format=Documenter.HTML(),
    pages = [
        "Introduction" => "index.md",
        "Sub-Modules" => [
            "Anatomical API" => "anatomical_api.md",
            "Spatial API" => "spatial_api.md",  # TODO
            "Color Channels" => "color_channels.md",  # TODO
            "Properties" => [
                "Introduction to Properties" => "introduction_properties.md",
                "Imaging" => "imaging_metadata.md",
                "Electrophysiology" => "electrophysiology.md",
            ],
        ],
        "Interfaceing With NeuroCore" => "working_with_interface.md",
        "Plots" => "plots.md"
    ],
    repo="https://github.com/JuliaNeuroscience/NeuroCore.jl/blob/{commit}{path}#L{line}",
    sitename="NeuroCore.jl",
    authors="Zachary P. Christensen",
)

deploydocs(repo = "github.com/JuliaNeuroscience/NeuroCore.jl.git")


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Test
using FiniteDifferences
using StaticArrays
using LinearAlgebra
using Suppressor
using Random
using Unitful
using Plots
using DataFrames

using OpticSim
# Geometry Module
using OpticSim.Geometry
# curve/surface imports
using OpticSim: findspan, makemesh, knotstoinsert, coefficients, inside, quadraticroots, tobeziersegments, evalcsg, makiemesh
# interval imports
using OpticSim: α, halfspaceintersection, positivehalfspace, lower, upper, EmptyInterval, rayorigininterval, intervalcomplement, intervalintersection, intervalunion, RayOrigin, Infinity, Intersection
include("TestData/TestData.jl")
# bounding box imports
using OpticSim: doesintersect
# RBT imports
using OpticSim: rotmat, rotmatd
# lens imports
using OpticSim: reflectedray, snell, refractedray, trace, intersection, pathlength, power
#Bounding volume hierarchy imports
# using OpticSim: partition!

const COMP_TOLERANCE = 25 * eps(Float64)
const RTOLERANCE = 1e-10
const ATOLERANCE = 10 * eps(Float64)
const SEED = 12312487
const ALL_TESTS = isempty(ARGS) || "all" in ARGS || "All" in ARGS || "ALL" in ARGS
const TESTSET_DIR = "testsets"

"""Evalute all functions not requiring arguments in a given module and test they don't throw anything"""
macro test_all_no_arg_functions(m)
    quote
        for n in names($m, all = true)
            # get all the valid functions
            if occursin("#", string(n)) || string(n) == string(nameof($m))
                continue
            end
            f = Core.eval($m, n)
            # get the methods of this function
            for meth in methods(f)
                # if this method has no args then try and evaluate it
                if meth.nargs == 1
                    # suppress STDOUT to prevent loads of stuff spamming the console
                    @test (@suppress_out begin
                        f()
                    end;
                    true)
                end
            end
        end
    end
end

"""
    @wrappedallocs(expr)
https://github.com/JuliaAlgebra/TypedPolynomials.jl/blob/master/test/runtests.jl

Given an expression, this macro wraps that expression inside a new function
which will evaluate that expression and measure the amount of memory allocated
by the expression. Wrapping the expression in a new function allows for more
accurate memory allocation detection when using global variables (e.g. when
at the REPL).
For example, `@wrappedallocs(x + y)` produces:
```julia
function g(x1, x2)
    @allocated x1 + x2
end
g(x, y)
```
You can use this macro in a unit test to verify that a function does not
allocate:
```
@test @wrappedallocs(x + y) == 0
```
"""
macro wrappedallocs(expr)
    argnames = [gensym() for a in expr.args]
    quote
        function g($(argnames...))
            @allocated $(Expr(expr.head, argnames...))
        end
        $(Expr(:call, :g, [esc(a) for a in expr.args]...))
    end
end

include("Benchmarks/Benchmarks.jl")

alltestsets = [
    "Repeat",
    "ParaxialAnalysis",
    "Paraxial",
    "JuliaLang",
    # "BVH",
    # "Examples", # slow
    "General",
    "SurfaceDefs",
    "Intersection",
    "Lenses",
    "OpticalSystem",
    "Emitters",
    "Comparison",
    "Visualization",
    "Allocations",
    "GlassCat",
    "Transform"
]

runtestsets = ALL_TESTS ? alltestsets : intersect(alltestsets, ARGS)
include.([joinpath(TESTSET_DIR, "$(testset).jl") for testset in runtestsets])


text:
"""
6, [1, 4] => [2, 3, 5, 6]
"""
function complement(N::Int, indices::Vector{Int})
    L = length(indices)
    x = Vector{Int}(undef, N - L)
    i_ = 1 # Position in the x vector
    j = 1 # Position in indices vector
    for i=1:N
        if j > L || indices[j]!=i
            x[i_] = i
            i_ += 1
        else
            j += 1
        end
    end
    x
end

"""
[1, 4, 5], [2, 4, 7] => [1, 5]
"""
function remove(ind1::Vector{Int}, ind2::Vector{Int})
    x = Int[]
    for i in ind1
        if i ∉ ind2
            push!(x, i)
        end
    end
    x
end

"""
[1, 4, 5], [2, 4, 7] => [1, 3]
"""
function shiftremove(ind1::Vector{Int}, ind2::Vector{Int})
    x = Int[]
    for i in ind1
        if i ∉ ind2
            counter = 0
            for i2 in ind2
                if i2 < i
                    counter += 1
                else
                    break
                end
            end
            push!(x, i-counter)
        end
    end
    x
end

function reducedindices(I_::Vector{Int}, I::Vector{Int})
    N = length(I_)
    x = Vector{Int}(undef, N)
    for n in 1:N
        x[n] = findfirst(isequal(I_[n]), I)
    end
    x
end

function reducedindices!(I_::Vector{Int}, I::Vector{Int})
    for n in 1:length(I_)
        I_[n] = findfirst(isequal(I_[n]), I)
    end
end

"""
Check if all indices are unique and smaller than or equal to imax.
"""
function check_indices(imax::Int, indices::Vector{Int})
    N = length(indices)
    for n=1:N
        i = indices[n]
        @assert 0 < i <= imax
        for m in n+1:N
            @assert i != indices[m]
        end
    end
end

"""
Check if the indices are sorted, unique and smaller than or equal to imax.
"""
function check_sortedindices(imax::Int, indices::Vector{Int})
    N = length(indices)
    if N == 0
        return nothing
    end
    i_ = indices[1]
    @assert 0 < i_ <= imax
    for i in indices[2:end]
        @assert 0 < i <= imax
        @assert i > i_
    end
end

"""
    check_embed_indices(indices::Array)

Determine whether a collection of indices, written as a list of (integers or lists of integers) is unique.
This assures that the embedded operators are in non-overlapping subspaces.
"""
function check_embed_indices(indices::Array)
    # short circuit return when `indices` is empty.
    length(indices) == 0 && return true

    err_str = "Variable `indices` comes in an unexpected form. Expecting `Array{Union{Int, Array{Int, 1}}, 1}`"
    @assert all(x isa Array || x isa Int for x in indices) err_str

    # flatten the indices and check for uniqueness
    # use a custom flatten because it's ≈ 4x  faster than Base.Iterators.flatten
    flatten(arr::Vector) = mapreduce(x -> x isa Vector ? flatten(x) : x, append!, arr, init=[])
    allunique(flatten(indices))
end


text:
# Julia wrapper for header: /usr/local/include/libavutil/audioconvert.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_channel_layout,
    av_get_channel_layout_string,
    av_get_channel_layout_nb_channels

function av_get_channel_layout(name)
    ccall((:av_get_channel_layout,libavutil),Uint64,(Ptr{Uint8},),name)
end

function av_get_channel_layout_string(buf,buf_size::Integer,nb_channels::Integer,channel_layout::Uint64)
    ccall((:av_get_channel_layout_string,libavutil),Void,(Ptr{Uint8},Cint,Cint,Uint64),buf,buf_size,nb_channels,channel_layout)
end

function av_get_channel_layout_nb_channels(channel_layout::Uint64)
    ccall((:av_get_channel_layout_nb_channels,libavutil),Cint,(Uint64,),channel_layout)
end


text:
# Linear layer, perform linear transformation to input array
# x₁ = softplus.(W) * x₀
struct LinearLayer{T,MT<:AbstractArray{T}}
    W::MT
end
@functor LinearLayer

LinearLayer(in_dim, out_dim) = LinearLayer(randn(out_dim, in_dim))

(lin::LinearLayer)(x) = softplus.(lin.W) * x

function Base.show(io::IO, layer::LinearLayer)
    return print(io, "LinearLayer(", size(layer.W, 2), ", ", size(layer.W, 1), ")")
end

# Product function, given an 2d array whose size is M×N, product layer will
# multiply every m neighboring rows of the array elementwisely to obtain
# an new array of size (M÷m)×N
function product(x, step=2)
    m, n = size(x)
    m % step == 0 || error("the first dimension of inputs must be multiple of step")
    new_x = reshape(x, step, m ÷ step, n)
    return .*([new_x[i, :, :] for i in 1:step]...)
end

# Primitive layer, mainly act as a container to hold basic kernels for the neural kernel network
struct Primitive{T}
    kernels::T
    Primitive(ks...) = new{typeof(ks)}(ks)
end
@functor Primitive

# flatten k kernel matrices of size Mk×Nk, and concatenate these 1d array into a k×(Mk*Nk) 2d array
_cat_kernel_array(x) = vcat([reshape(x[i], 1, :) for i in 1:length(x)]...)

# NOTE, though we implement `ew` & `pw` function for Primitive, it isn't a subtype of Kernel
# type, I do this because it will facilitate writing NeuralKernelNetwork
ew(p::Primitive, x) = _cat_kernel_array(map(k -> kernelmatrix_diag(k, x), p.kernels))
pw(p::Primitive, x) = _cat_kernel_array(map(k -> kernelmatrix(k, x), p.kernels))

function ew(p::Primitive, x, x′)
    return _cat_kernel_array(map(k -> kernelmatrix_diag(k, x, x′), p.kernels))
end
pw(p::Primitive, x, x′) = _cat_kernel_array(map(k -> kernelmatrix(k, x, x′), p.kernels))

function Base.show(io::IO, layer::Primitive)
    print(io, "Primitive(")
    join(io, layer.kernels, ", ")
    return print(io, ")")
end

"""
    NeuralKernelNetwork(primitives, nn)

Constructs a Neural Kernel Network (NKN) [1].

`primitives` are the based kernels, combined by `nn`.

```julia
k1 = 0.6 * (SEKernel() ∘ ScaleTransform(0.5))
k2 = 0.4 * (Matern32Kernel() ∘ ScaleTransform(0.1))
primitives = Primitive(k1, k2)
nkn = NeuralKernelNetwork(primitives, Chain(LinearLayer(2, 2), product))
```

[1] - Sun, Shengyang, et al. "Differentiable compositional kernel learning for Gaussian
    processes." International Conference on Machine Learning. PMLR, 2018.
"""
struct NeuralKernelNetwork{PT,NNT} <: Kernel
    primitives::PT
    nn::NNT
end
@functor NeuralKernelNetwork

# use this function to reshape the 1d array back to kernel matrix
_rebuild_kernel(x, n, m) = reshape(x, n, m)
_rebuild_diag(x) = reshape(x, :)

(κ::NeuralKernelNetwork)(x, y) = only(kernelmatrix(κ, [x], [y]))

function kernelmatrix_diag(nkn::NeuralKernelNetwork, x::AbstractVector)
    return _rebuild_diag(nkn.nn(ew(nkn.primitives, x)))
end

function kernelmatrix(nkn::NeuralKernelNetwork, x::AbstractVector)
    return _rebuild_kernel(nkn.nn(pw(nkn.primitives, x)), length(x), length(x))
end

function kernelmatrix_diag(nkn::NeuralKernelNetwork, x::AbstractVector, x′::AbstractVector)
    return _rebuild_diag(nkn.nn(ew(nkn.primitives, x, x′)))
end

function kernelmatrix(nkn::NeuralKernelNetwork, x::AbstractVector, x′::AbstractVector)
    return _rebuild_kernel(nkn.nn(pw(nkn.primitives, x, x′)), length(x), length(x′))
end

function kernelmatrix_diag!(K::AbstractVector, nkn::NeuralKernelNetwork, x::AbstractVector)
    K .= kernelmatrix_diag(nkn, x)
    return K
end

function kernelmatrix!(K::AbstractMatrix, nkn::NeuralKernelNetwork, x::AbstractVector)
    K .= kernelmatrix(nkn, x)
    return K
end

function kernelmatrix_diag!(
    K::AbstractVector, nkn::NeuralKernelNetwork, x::AbstractVector, x′::AbstractVector
)
    K .= kernelmatrix_diag(nkn, x, x′)
    return K
end

function kernelmatrix!(
    K::AbstractMatrix, nkn::NeuralKernelNetwork, x::AbstractVector, x′::AbstractVector
)
    K .= kernelmatrix(nkn, x, x′)
    return K
end

function Base.show(io::IO, kernel::NeuralKernelNetwork)
    print(io, "NeuralKernelNetwork(")
    join(io, [kernel.primitives, kernel.nn], ", ")
    return print(io, ")")
end


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

df1 = DataFrame(x=[1,2,3], y=[1, 2, 3])
df2 = DataFrame(x=[2, 4, 6], y=[2, 4, 6])

plot(df1, x=:x, y=:y, layer(df2, Geom.point))


text:
using LAR

#""" The `spacePartition` function takes as input a **non-valid** (with the meaning used in solid modeling field --- see~\cite{Requicha:1980:RRS:356827.356833}) `LAR` model of dimension d-1, i.e.~a triple `(V,FV,EV)`, The `spacePartition` function returns the **valid** `LAR` boundary model `(W,FW,EW)` of the space partition induced by `FV`. First an array `buckets` indexed on faces is computed. It contains the subsets of faces with greatest probability of intersecting each indexing face, respectively.
#"""

""" Half-line crossing test """
function crossingTest(new,old,count,status)
    if status == 0
        status = new
        count += 0.5
    else
        if status == old
        	count += 0.5
        else
        	count -= 0.5
        end
        status = 0
    end
end

function setTile(box)
	tiles = [[9,1,5],[8,0,4],[10,2,6]]
	b1,b2,b3,b4 = box
	function tileCode(point)
		x,y = point
		code = 0
		if y>b1 code=code|1 end
		if y<b2 code=code|2 end
		if x>b3 code=code|4 end
		if x<b4 code=code|8 end
		return code
	end
	return tileCode
end

# Point in polygon classification #
function pointInPolygonClassification(V,EV)
	V = V[1:2,:]
    function pointInPolygonClassification0(pnt)
        x,y = pnt
        xmin,xmax,ymin,ymax = x,x,y,y
        tilecode = setTile([ymax,ymin,xmax,xmin])
        count,status = 0,0

        for k in size(EV,2)
        	edge = EV[:,k]
            p1,p2 = V[:,edge[1]],V[:,edge[2]]
            (x1,y1),(x2,y2) = p1,p2
            c1,c2 = tilecode(p1),tilecode(p2)
            c_edge, c_un, c_int = c1$c2, c1|c2, c1&c2

            if (c_edge == 0) & (c_un == 0) return "p_on"
            elseif (c_edge == 12) & (c_un == c_edge) return "p_on"
            elseif c_edge == 3
                if c_int == 0 return "p_on"
                elseif c_int == 4 count += 1 end
            elseif c_edge == 15
                x_int = ((y-y2)*(x1-x2)/(y1-y2))+x2
                if x_int > x count += 1
                elseif x_int == x return "p_on" end
            elseif (c_edge == 13) & ((c1==4) | (c2==4))
                    crossingTest(1,2,status,count)
            elseif (c_edge == 14) & ((c1==4) | (c2==4))
                    crossingTest(2,1,status,count)
            elseif c_edge == 7 count += 1
            elseif c_edge == 11 count = count
            elseif c_edge == 1
                if c_int == 0 return "p_on"
                elseif c_int == 4 crossingTest(1,2,status,count) end
            elseif c_edge == 2
                if c_int == 0 return "p_on"
                elseif c_int == 4 crossingTest(2,1,status,count) end
            elseif (c_edge == 4) & (c_un == c_edge) return "p_on"
            elseif (c_edge == 8) & (c_un == c_edge) return "p_on"
            elseif c_edge == 5
                if (c1==0) | (c2==0) return "p_on"
                else crossingTest(1,2,status,count) end
            elseif c_edge == 6
                if (c1==0) | (c2==0) return "p_on"
                else crossingTest(2,1,status,count) end
            elseif (c_edge == 9) & ((c1==0) | (c2==0)) return "p_on"
            elseif (c_edge == 10) & ((c1==0) | (c2==0)) return "p_on"
            end
        end
        if (round(count)%2)==1
        	return "p_in"
        else
        	return "p_out"
        end
    end
    return pointInPolygonClassification0
end

function crossRelation(EV::Array{Int64,2},FV::Array{Array{Int64,1},1})
	sparseEV = cellComplex(EV)
	sparseFV = cellComplex(FV)
	sparseEF = (sparseFV' * sparseEV)'
	I,J,V = findnz(sparseEF)
	EF = [Int64[] for k=1:size(EV,2)]
	for (i,j,v) in zip(I,J,V)
		if v==2
			push!(EF[i],j)
		end
	end
	return EF
end

function crossRelation(FV::Array{Int64,2},EV::Array{Int64,2})
	FW = [FV[:,k] for k=1:size(FV,2)]
	out = crossRelation(FW,EV)
end

function crossRelation(FV::Array{Array{Int64,1},1},EV::Array{Int64,2})
	sparseFV = cellComplex(FV)
	sparseEV = cellComplex(EV)
	sparseEF = (sparseFV' * sparseEV)'
	I,J,V = findnz(sparseEF)
	FE = [Int64[] for k=1:length(FV)]
	for (i,j,v) in zip(I,J,V)
		if v==2
			push!(FE[j],i)
		end
	end
	return FE
end

function subModel(V::Array{Float64,2}, FV::Array{Int64,2}, EV::Array{Int64,2},
				bucket::Array{Int64,1},f::Int64,FE::Array{Array{Int64,1},1})
	FW = [FV[:,k] for k=1:size(FV,2)]
	out = subModel(V,FW,EV,bucket,f,FE)
end

function subModel(V::Array{Float64,2}, FV::Array{Array{Int64,1},1},
					EV::Array{Int64,2},bucket::Array{Int64,1},f::Int64,
					FE::Array{Array{Int64,1},1})
	FW = [FV[f] for f in bucket]
	EW = hcat(collect(Set([EV[:,e] for f in bucket for e in FE[f]]))...)

	oldVerts = Set(Int64[])
	for k=1:length(bucket)
		union!(oldVerts, Set(FW[k]))
	end
	olds = collect(oldVerts)
	vdict = Dict(zip(olds, 1:length(olds)))
	fz = [[vdict[v] for v in FW[k]] for k=1:length(FW)]
	ez = hcat([[vdict[v] for v in EW[:,k]] for k=1:size(EW,2)]...)
	z = hcat([V[:,v] for v in olds]...)
	pivot = find(bucket .== f)[1]
	return z,fz,ez,pivot
end

function visualize(Z,FZ::Array{Int64,2},EZ,f)
	FW = [FZ[:,k] for k=1:size(FZ,2)]
	out = visualize(Z,FW,EZ,f)
end

function visualize(Z,FZ::Array{Array{Int64,1},1},EZ,f)
	bucket = lar2hpc(Z,FZ)
	bucketEdges = lar2hpc(Z,EZ)
	params = PyObject(pyeval("list([1.,0.,0.,0.1,  0.,1.,0.,0.1,
				0.,0.,1.,0.1, 0.,0.,0.,0.1, 100.])"))
	glass = p.MATERIAL(params)
	pivot = p.COLOR(p.RED)(p.JOIN(lar2hpc(Z,[FZ[f]])))
	p.VIEW(p.STRUCT([glass(bucket),bucketEdges,pivot]))
end

function meetZero( Z, vpair )
	(u,v) = vpair
    testValue = Z[:,u][3] * Z[:,v][3]
    if testValue > abs(10^-6.)
        return false
    else
    	return true
    end
end

function crossZero( Z, vpair )
	(u,v) = vpair
    testValue = sign(vcode(Z[:,u])[3]) * sign(vcode(Z[:,v])[3])
    if testValue == -1
        return true
    else
    	return false
    end
end

function edgefilter(Z,EZ,e)
	cross = crossZero(Z,EZ[:,e])
	side1 = (abs(Z[3,EZ[1,e]])<10^-6.)
	side2 = (abs(Z[3,EZ[2,e]])<10^-6.)
	#@show cross,side1,side2
	out = cross | (side1 $ side2)
	#if out @show out end
	return out
end

function intersectSegmentWithZero(p1::Array{Float64,1}, p2::Array{Float64,1})
	x1,y1,z1 = p1
	x2,y2,z2 = p2
	# z = z1 + alpha(z2 - z1)
	alpha_0 = (0 - z1)/(z2 - z1)
	point = Float64[x1 + alpha_0*(x2 - x1) , y1 + alpha_0*(y2 - y1)]
	return point
end

function remap(W,FW,EW)
	W,close,clusters,vmap = p.pruneVertices(W')
	if typeof(W) == Array{Any,1}
		W = Array{Float64,2}(hcat(W...))
	else
		W = Array{Float64,2}(W')
	end
	vmap = vmap + 1
	EW = [sort([vmap[EW[1,k]],vmap[EW[2,k]]]) for k in 1:size(EW,2)]
	EW = collect(Set([EW[k] for k=1:length(EW)]))
	EW = hcat(EW...)
	FW = [sort([vmap[FW[k][h]] for h=1:length(FW[k])]) for k=1:length(FW)]
	FW = collect(Set([FW[k] for k=1:length(FW)]))
	return W,FW,EW
end

function checkLines4equalPoints(lines)
	out = []
	for line in lines
		if length(line) > 4
			b = reshape(line,2,Int(length(line)/2))
			line = vcat(p.pruneVertices(b')[1]...)
		end
		if (line[1,:]!=line[3,:]) | (line[2,:]!=line[4,:])
			push!(out,line)
		end
	end
	out
end

function onlyOne_INorON_vertex(Z,EZ,pivotEdges,lines)
	ev = hcat([EZ[:,e] for e in pivotEdges]...)
	classify = pointInPolygonClassification(Z,ev)
	out = Array{Array{Float64,1},1}()
	for line in lines
		p1,p2 = [line[1],line[2]],[line[3],line[4]]
		code1 = classify(p1)
		code2 = classify(p2)
		#@show code1,code2
		if !((code1=="p_out") & (code2=="p_out"))
			push!(out,line)
		end
	end
	return out
end

function spacePartition(V::Array{Float64,2}, FV::Array{Int64,2}, EV::Array{Int64,2},
						debug=false)
	FW = [FV[:,k] for k=1:size(FV,2)]
	out = spacePartition(V,FW,EV,debug)
end

function spacePartition(V::Array{Float64,2}, FV::Array{Array{Int64,1},1},
						EV::Array{Int64,2},debug=false)
	""" input: face index f; candidate incident faces F[f]; """
	boxes = lar2boxes(V,FV)
	buckets = boxBucketing(boxes)
	FE = crossRelation(FV,EV)

	Vertices = Array{Float64,2}()
	Edges = Array{Int64,2}()
	Faces = Array{Array{Int64,1},1}()
	nverts = 0

	for (f,F) in enumerate(buckets)
		@show (f,F)
		""" F[f] submodel extraction from (V,FV,EV) """
		Z,FZ,EZ,pivot = subModel(V,FV,EV,F,f,FE)
		if debug visualize(Z,FZ,EZ,pivot) end

		""" Computation of submanifold map M moving f to z=0 """
		M = submanifoldMapping(V,FV,f)

		""" Transformation of F(f) by M, giving (W,EW,FW) := M(F(f)) """
		Z1 = vcat(Z,ones((1,size(Z,2))))
		Y = M * Z1
		Z = Y[ 1:3, : ]
		if debug visualize(Z,FZ,EZ,pivot) end

		""" filtering of EW edges traversing z=0, """
		fe = crossRelation(FZ,EZ)
		edges4face = [[e for e in face] for (k,face) in enumerate(fe) if k!=pivot]

		pivotEdges = fe[pivot]
		crossingEdges = [[e for e in edges if edgefilter(Z,EZ,e)] for edges in edges4face]
		edges = [edgeList for edgeList in crossingEdges if edgeList!=[]]
		EW = vcat(edges...,pivotEdges...)
		ez = hcat([EZ[:,e] for e in EW]...)
		if debug visualize(Z,FZ,ez,pivot) end

		""" for each face in EW, computate the aligned set of points p(z=0) """
		lines = Array{Array{Float64,1},1}()
		for face in edges
			#assert(length(face)==2) # may fail with non convex faces TODO: make general
			line = Array{Float64,1}()
			for e in face
				v1,v2 = EZ[:,e]
				px,py = vcode(intersectSegmentWithZero(Z[:,v1],Z[:,v2]))[1:2]
				push!(line,px); push!(line,py)
			end
			push!(lines,line)
		end

		""" Further filtering of lines with one non-external vertex """
		lines = onlyOne_INorON_vertex(Z,EZ,pivotEdges,lines)

		for e in pivotEdges
			line = Array{Float64,1}()
			push!(line,vcode(Z[:,EZ[1,e]])[1]); push!(line,vcode(Z[:,EZ[1,e]])[2])
			push!(line,vcode(Z[:,EZ[2,e]])[1]); push!(line,vcode(Z[:,EZ[2,e]])[2])
			push!(lines,line)
		end
		lines = checkLines4equalPoints(lines)
		lineArray = hcat(lines...)
		W,FW,EW = larFromLines(lineArray')

		""" remove external cycle """
		chains = boundary(W,EW)
		if length(chains) == 2
			FW = [FW[1]]
		else
			areaPairing = chainAreas(W,EW,chains)
			maxAreas = maximum([abs(area) for area in areaPairing])
			FW = [face for (k,face) in enumerate(FW) if abs(areaPairing[k]) != maxAreas]
		end
		if debug viewLarIndices(W,EW,FW,3) end

		""" Apply the inverse submanifold transform """
		W1 = vcat(W,zeros((1,size(W,2))),ones((1,size(W,2))))
		Y = inv(M) * W1
		Z = Y[ 1:3, : ]
		if debug viewLarIndices(Z,EW,FW,3) end

		""" Accumulate the submodel parts """
		n = size(Z,2)
		if f==1
			Vertices = copy(Z)
			Edges = copy(EW)
			Faces = copy(FW)
			nverts = n
		else
			Vertices = hcat(Vertices,copy(Z))
			Edges = hcat(Edges,copy(EW)+nverts)
			Faces = vcat(Faces,copy(FW)+nverts)
			nverts += n
		end
	end
	""" return the **valid** `LAR` 2-skeleton model `(W,FW,EW)` """
	U,FU,EU = remap(Vertices,Faces,Edges)
	@show length(U),length(FU),length(EU)
end


text:
function count_increases(numbers::Vector{Int}, window::Int)
    count(numbers[1:end-window] .< numbers[1+window:end])
end

numbers = [parse(Int, number) for number in readlines()]
println(count_increases(numbers, 1))
println(count_increases(numbers, 3))


text:
using FIX
using Base.Test

@test FIX.checksum("123ioн") == 171

smpl = Dict(1 => "g", 2=> "h")
@test FIX.fixjoin(smpl, Char(1)) == "2=h" * Char(1) * "1=g" * Char(1)
smpl_back = FIX.fixparse(UInt8.([x for x in FIX.fixjoin(smpl, Char(1))]))
@test smpl_back == smpl
@test FIX.fixconvert(FIX.TAGS_INT_STRING, smpl_back) == Dict("AdvId" => "h", "Account" => "g")


text:
import ReverseDiff

# Gen versions of certain ReverseDiff functions
# this is because the generated code for the static IR is currently loaded in
# the Main module. to avoid requiring the Main module environment to import the
# name ReverseDiff, these aliases allow us to only require that the Main module
# environmment have Gen imported.
new_tape() = ReverseDiff.InstructionTape()
track(value, tape) = ReverseDiff.track(value, tape)
value(maybe_tracked) = ReverseDiff.value(maybe_tracked)
deriv!(maybe_tracked, deriv) = ReverseDiff.deriv!(maybe_tracked, deriv)
reverse_pass!(tape) = ReverseDiff.reverse_pass!(tape)
deriv(tracked) = ReverseDiff.deriv(tracked)
istracked(maybe_tracked) = ReverseDiff.istracked(maybe_tracked)
record!(tape, instruction_type, args...) = ReverseDiff.record!(tape, instruction_type, args...)
increment_deriv!(arg, deriv) = ReverseDiff.increment_deriv!(arg, deriv)
seed!(tracked) = ReverseDiff.seed!(tracked)
unseed!(tracked) = ReverseDiff.unseed!(tracked)

using ReverseDiff: InstructionTape, TrackedReal, TrackedArray

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
"""
    int(x) = trunc(Int, x)
"""
int(x) = trunc(Int, x)

"""
    dropmean(X, d) = dropdims(mean(X, dims=d), dims=d)
"""
function dropmean(X, d)
    if ndims(X) == 1
        mean(X; dims = d)
    else
        dropdims(mean(X; dims = d); dims = d)
    end
end

"""
    squeeze(X::AbstractArray)

removing singleton dimensions
"""
function squeeze(X::AbstractArray)
    keepd = Tuple(d for d in size(X) if d != 1)
    return reshape(X, keepd)
end

"""
    detrend!(data, n)

(in place) Linear detrend of signals along first axis removing the n-th order polynomial.
This detrend function is limited to linear orders (0th and 1st order).

### Arguments

  - `data::AbstractArray`: N-dim array where signal is in column-major order
  - `n::Integer`: `n = 0` subtracts mean from data, `n = 1` removes linear trend

**Note**: shape of data must be (signal length, ...)
"""
function detrend!(data::AbstractArray, n::Integer)
    original_shape = size(data)
    nsamp = size(data, 1)  # number of samples

    A = Array{Float64}([ones(nsamp) Array(1:nsamp)])

    data = reshape(data, (nsamp, :))  # reshaping data
    if n == 0
        data .-= mean(data; dims = 1)
    elseif n == 1
        data .-= A * (A \ data)
    end
    return reshape(data, original_shape)
end

"""
    window = hanning_fun(N)

Hanning window similar to MATLAB `hanning` implementation
"""
function hanning_fun(N::Integer)
    x = [range(0.0, 1.0; length = N + 2);]
    window = 0.5 .* (1 .- cospi.(2 .* x))
    window = (window + window[end:-1:1]) ./ 2  # forcing symmetry
    return window[2:(end-1)]  # excluding the zero values
end

"""
    review(X, dim)

Reshape + View
"""
review(X, dim) = reshape(view(X, :, dim), :, 1)

"""
    pink_noise(n, c)

Pink noise with `n` samples and pinkness of:
```math
P_s \\propto \\frac{1}{f^c} .
```

### Arguments

  - `n::Integer`: Number of samples
  - `c::Float64`: pinkness order (frequency power)

"""
function pink_noise(n::Integer, c::Float64)
    nf = int(n / 2)  # Nyquist frequency
    x = reshape(range(start = 0.0, stop = 1.0, length = n), 1, :)
    f = reshape(range(start = 1, stop = nf), :, 1)
    phi = rand(nf)
    y = f * x .+ phi
    z = sin.(2 * pi * y) ./ (f .^ c)
    return sum(z; dims = 1)[1, :]
end


text:
using Documenter, PrivateExample4

makedocs(modules = [PrivateExample4], sitename = "PrivateExample4.jl")

deploydocs(repo = "gitlab-test.eastus2.cloudapp.azure.com/vdayanand/PrivateExample4.jl.git")


text:
using FactCheck
include(Pkg.dir("PseudoGenomes", "src","pseudogenome-deletions.jl"))

ref_fasta = joinpath( Pkg.dir(),"PseudoGenomes","testdata","test.fasta")

# remember to change the coordinate strings in the strain fields if you want to play around with length of deletion masked.
SVs_filename = Pkg.dir("PseudoGenomes","testdata","28strains.REL-1410-SV.sdp.SAMPLE.MOD.tab")
outpath = joinpath( Pkg.dir(),"PseudoGenomes","testdata","out.maskdeletions.fasta")

# -- test
SV_df = read_SV_file(SVs_filename, gzip=false)

# -- test
#"PWK_PhJ"
strain_SVs_df = get_strain_SVs(SV_df, "129P2_OlaHsd",Regex("DEL"))

# -- test
mask_fasta_file_with_strain_SV_deletions(ref_fasta,outpath,strain_SVs_df)

# -- test
#mask_region!(sequence,SV_deletions_df)

#mask_fasta_file_with_strain_SV_deletions(mm10path,outpath,strain_SVs_df)
#mask_fasta_file_with_strain_SV_deletions(mm10path,outpath,SVs_filename,"CAST_EiJ")


text:
# using Test

using PaddedMatrices, DifferentiableObjects

function rosenbrock(x::AbstractArray{T}) where T
    out = zero(T)
    @inbounds for i ∈ 2:2:length(x)
        out += 100abs2(abs2(x[i-1]) - x[i]) + abs2(x[i-1] - 1)
    end
    out
end

P = 6

state = DifferentiableObjects.BFGSState2(Val(P));
initial_x = fill!(PaddedMatrices.MutableFixedSizeVector{P,Float64}(undef), 3.2);
ls2 = DifferentiableObjects.BackTracking2(Val(2));
ls3 = DifferentiableObjects.BackTracking2(Val(3));
obj = OnceDifferentiable(rosenbrock, initial_x);
DifferentiableObjects.optimize_light!(state, obj, initial_x, ls2), state.x_old
DifferentiableObjects.optimize_light!(state, obj, initial_x, ls3), state.x_old
DifferentiableObjects.optimize_scale!(state, obj, initial_x, ls2), state.x_old
DifferentiableObjects.optimize_scale!(state, obj, initial_x, ls3), state.x_old

using BenchmarkTools
@benchmark DifferentiableObjects.optimize_light!($state, $obj, $initial_x, $ls2)
@benchmark DifferentiableObjects.optimize_light!($state, $obj, $initial_x, $ls3)
@benchmark DifferentiableObjects.optimize_scale!($state, $obj, $initial_x, $ls2)
@benchmark DifferentiableObjects.optimize_scale!($state, $obj, $initial_x, $ls3)

using PaddedMatrices, DifferentiableObjects
Base.IndexStyle(::PaddedMatrices.AbstractPaddedVector) = IndexLinear()

function rosenbrock(x::AbstractArray{T}) where {T}
    out = zero(T)
    @inbounds for i ∈ 2:2:length(x)
        out += 100abs2(abs2(x[i-1]) - x[i]) + abs2(x[i-1] - 1)
    end
    out
end

P = 6
initial_x = fill!(PaddedMatrices.MutableFixedSizeVector{P,Float64}(undef), 3.2);
obj2 = TwiceDifferentiable(rosenbrock, initial_x);

using ForwardDiff, DiffResults, BenchmarkTools
x = initial_x;
c = obj2.config;

hess = DiffResults.hessian(c.result);
grad = DiffResults.gradient(c.result);
cfg = c.jacobian_config;

# hessian!(obj2.config, initial_x)

@which ForwardDiff.jacobian!(hess, c, grad, x, cfg, Val{false}())

result = hess;
f! = c;
y = grad;
@which ForwardDiff.vector_mode_jacobian!(result, f!, y, x, cfg)

# function vector_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
#     ydual = vector_mode_dual_eval(f!, y, x, cfg)
#     map!(d -> value(T,d), y, ydual)
#     extract_jacobian!(T, result, ydual, N)
#     extract_value!(T, result, y, ydual)
#     return result
# end

T, V, N = ForwardDiff.Tag{typeof(rosenbrock),Float64},Float64,6,Tuple{MutableFixedSizeArray{Tuple{6},ForwardDiff.Dual{ForwardDiff.Tag{typeof(rosenbrock),Float64},Float64,6},1,Tuple{1},6},MutableFixedSizeArray{Tuple{6},ForwardDiff.Dual{ForwardDiff.Tag{typeof(rosenbrock),Float64},Float64,6},1,Tuple{1},6}};

@which ForwardDiff.vector_mode_dual_eval(f!, y, x, cfg)

###
### Here, I can run random benchamrks without a problem.
###

ydual, xdual = cfg.duals;
typeof(ydual), typeof(xdual)

ForwardDiff.seed!(xdual, x, cfg.seeds)
ForwardDiff.seed!(ydual, y)
# seem to be in the clear, based on repeated
@benchmark ForwardDiff.seed!(xdual, x, cfg.seeds)
@benchmark ForwardDiff.seed!(ydual, y)

# Segfaults
# f!(ydual, xdual)

# segfaults
@which ForwardDiff.gradient!(c.inner_result, c.f, xdual, c.gradient_config, Val{false}())
@which ForwardDiff.vector_mode_gradient!(c.inner_result, c.f, xdual, c.gradient_config)
@which ForwardDiff.vector_mode_dual_eval(c.f, xdual, c.gradient_config)

xdualdual = c.gradient_config.duals;
# underlying tuple of length 0!!!
f, x = rosenbrock, initial_x;
result = DifferentiableObjects.HessianDiffResult(x);
chunk = DifferentiableObjects.Chunk(Val{P}())
# tag = ForwardDiff.Tag(f, T)

# here we construct a JacobianConfig so that yduals is a PtrVector to DiffResults.gradient(inner_result)
tag_instance = ForwardDiff.Tag(f, T)
tag_type = DifferentiableObjects.extract_tag(tag_instance)
# @show tag
seeds = ForwardDiff.construct_seeds(ForwardDiff.Partials{P,Float64})
# inner_result = GradientDiffResult(jacobian_config.duals[2])
inner_result = DifferentiableObjects.GradientDiffResult{ForwardDiff.Dual{tag_type,Float64,P},P}(undef)
duals = (DiffResults.gradient(inner_result), similar(x, ForwardDiff.Dual{tag_type,Float64,P}))
jacobian_config = ForwardDiff.JacobianConfig{tag_type,Float64,P,typeof(duals)}(seeds, duals)

# jacobian_config = ForwardDiff.JacobianConfig((f,ForwardDiff.gradient), DiffResults.gradient(result), x, chunk, tag)
# jacobian_config = ForwardDiff.JacobianConfig((f,ForwardDiff.gradient), DiffResults.gradient(result), x, chunk, tag)
gradient_config = ForwardDiff.GradientConfig(f, jacobian_config.duals[2], chunk, tag_instance)

# segfault
# @benchmark ForwardDiff.seed!(xdualdual, xdual, c.gradient_config.seeds)
@which ForwardDiff.seed!(xdualdual, xdual, c.gradient_config.seeds)
# xdualdual's tuple is of length 0!!!

c.f(xdual)

# segfaults
# @benchmark ForwardDiff.vector_mode_dual_eval(c.f, xdual, c.gradient_config)

ydualdual = ForwardDiff.vector_mode_dual_eval(c.f, xdual, c.gradient_config)

# Segfaults on call
@which ForwardDiff.extract_gradient!(ForwardDiff.Tag{typeof(rosenbrock),Float64}, c.inner_result, ydualdual)

# @which -> Segfaults!!! Yikes!?!?
v = ForwardDiff.value(ForwardDiff.Tag{typeof(rosenbrock),Float64}, ydualdual)
@benchmark ForwardDiff.value(ForwardDiff.Tag{typeof(rosenbrock),Float64}, ydualdual)

# body of extract_gradient! in:
# extract_gradient!(::Type{T}, result::DiffResults.DiffResult, dual::ForwardDiff.Dual) where T in ForwardDiff at /home/chriselrod/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:71
inner_result = DiffResults.value!(c.inner_result, ForwardDiff.value(ForwardDiff.Tag{typeof(rosenbrock),Float64}, ydualdual));

partials = ForwardDiff.partials(ForwardDiff.Tag{typeof(rosenbrock),Float64}, ydualdual);
grad_inner = DiffResults.gradient(c.inner_result);

typeof(grad_inner), typeof(partials)

function Base.copyto!()

end

@benchmark copyto!(grad_inner, partials)

@which DiffResults.gradient!(c.inner_result, partials)
# forwards to derivative!(r::Union{DifferentiableObjects.GradientDiffResult{V,P,R}, DifferentiableObjects.HessianDiffResult{V,P,R,L} where L} where R where P where V, x::AbstractArray) in DifferentiableObjects at /home/chriselrod/.julia/dev/DifferentiableObjects/src/forward_diff_differentiable.jl:61
# DifferentiableObjects/src/forward_diff_differentiable.jl:61: DiffResults.derivative!
#
grad_inner = DiffResults.gradient(c.inner_result);
# @code_warntype copyto!(grad_inner, partials)
@which copyto!(grad_inner, partials)

DiffResults.derivative!(c.inner_result, partials);
# benchmark segfaults

DiffResults.gradient!(c.inner_result, partials);
# segfault

@which copyto!(grad_inner, partials)
copyto!(grad_inner, partials);
# segfault

# segfaults
inner_result = DiffResults.gradient!(c.inner_result, ForwardDiff.partials(ForwardDiff.Tag{typeof(rosenbrock),Float64}, ydualdual))

result = ForwardDiff.extract_gradient!(ForwardDiff.Tag{typeof(rosenbrock),Float64}, c.inner_result, ydualdual)

# Segfaults
ForwardDiff.gradient!(c.inner_result, c.f, xdual, c.gradient_config, Val{false}())

DiffResults.value!(c.result, ForwardDiff.value(DiffResults.value(c.inner_result)))

ydual

ydual = ForwardDiff.vector_mode_dual_eval(f!, y, x, cfg)

map!(d -> ForwardDiff.value(T,d), y, ydual)
ForwardDiff.extract_jacobian!(T, result, ydual, N)
ForwardDiff.extract_value!(T, result, y, ydual)

ForwardDiff.jacobian!(hess, c, grad, x, jconfig, Val{false}())

hessian!(obj2, initial_x)
hessian(obj2)

@which hessian!(obj2, initial_x)
@which hessian(obj2)

using StaticArrays, StaticOptim

sx = @SVector fill(3.2, P);
soptimize(rosenbrock, sx)
@benchmark soptimize(rosenbrock, $sx)
@benchmark soptimize(rosenbrock, $sx, StaticOptim.Order3())

using Optim, LineSearches
xa = fill(3.2, 6);
hz = Optim.BFGS();
backtrack2 = Optim.BFGS(linesearch = LineSearches.BackTracking(order = 2));
backtrack3 = Optim.BFGS(linesearch = LineSearches.BackTracking(order = 3));
@benchmark optimize(rosenbrock, $xa, $hz)
@benchmark optimize(rosenbrock, $xa, $hz; autodiff = :forward)
@benchmark optimize(rosenbrock, $xa, $backtrack2)
@benchmark optimize(rosenbrock, $xa, $backtrack2; autodiff = :forward)
@benchmark optimize(rosenbrock, $xa, $backtrack3)
@benchmark optimize(rosenbrock, $xa, $backtrack3; autodiff = :forward)

# @inline rosenbrock(x) =  (@inbounds out = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2; out)

# initial_x = zero(SizedSIMDVector{2,Float64});
# d = OnceDifferentiable(rosenbrock, initial_x);
# # method = BFGS(;linesearch = DifferentiableObjects.BackTracking())
# method = BFGS();
# options = LightOptions();
# state = DifferentiableObjects.initial_state(method, options, d, initial_x);
#
# optimum, optimum_val = optimize_light(d, initial_x, method, options, state)
#
#
#
# using BenchmarkTools
# opt_bench = @benchmarkable optimize_light($d, $initial_x, $method, $options, init_state) setup=(init_state=DifferentiableObjects.initial_state($method, $options, $d, $initial_x); $initial_x .= 0)
# run(opt_bench)
#
# using LinearAlgebra, ForwardDiff
# dphi_0 = real(dot(gradient(d), state.s))
# # @show dphi_0
# # reset the direction if it becomes corrupted
# if dphi_0 >= zero(dphi_0) && DifferentiableObjects.reset_search_direction!(state, d, method)
#     dphi_0 = real(dot(gradient(d), state.s)) # update after direction reset
# end
# phi_0  = value(d)
# obj, x = d, state.x
#
# @benchmark DifferentiableObjects.update_g!($d, $state, $method)
# @benchmark DifferentiableObjects.value_gradient!($d, $state.x)
# @benchmark DifferentiableObjects.fdf($d, $state.x)
# @code_warntype DifferentiableObjects.fdf(d, state.x)
#
# using Profile
# @profile optimize_light(d, initial_x, method, options, state)
# Profile.print()
#
#
#
# @benchmark DifferentiableObjects.linesearch!($method.linesearch!, $d, $state.x, $state.s, $state.alpha, $state.x_ls, $phi_0, $dphi_0)
# @code_warntype DifferentiableObjects.linesearch!(method.linesearch!, d, state.x, state.s, state.alpha, state.x_ls, phi_0, dphi_0)
#
# @benchmark $method.linesearch!($d, $state.x, $state.s, $state.alpha, $state.x_ls, $phi_0, $dphi_0)
# @code_warntype method.linesearch!(d, state.x, state.s, state.alpha, state.x_ls, phi_0, dphi_0)
#
#
# @show optimum
# @test all(i -> optimum[i] ≈ 1, 1:length(optimum))
#
#
# using StaticArrays, StaticOptim
#
# sx = @SVector zeros(P)
# soptimize(rosenbrock, sx)
# @benchmark soptimize(rosenbrock, $sx)


text:
@testset "examples" begin

function find_sources(path::String, sources=String[])
    if isdir(path)
        for entry in readdir(path)
            find_sources(joinpath(path, entry), sources)
        end
    elseif endswith(path, ".jl")
        push!(sources, path)
    end
    sources
end

examples_dir = joinpath(@__DIR__, "..", "examples")
examples = find_sources(examples_dir)
filter!(file -> readline(file) != "# EXCLUDE FROM TESTING", examples)
filter!(file -> !occursin("Kaleidoscope", file), examples)

cd(examples_dir) do
    @testset for example in examples
        mod = @eval module $(gensym()) end
        @eval mod begin
            redirect_stdout(devnull) do
                include($example)
            end
        end
    end
end

end


text:
msg = "Hello world!!!"
println(msg)

max_stars_repo_name: rbpope/Site

max_issues_repo_name: rbpope/Site

max_forks_repo_name: rbpope/Site


text:
cuallocator()=true      # set to true to use the CUDA.jl allocator, false to use Knet allocator

# KnetPtr type holds a gpu allocated pointer.  We try to minimize the number of actual
# allocations, which are slow, by reusing preallocated but garbage collected pointers.

mutable struct KnetPtr
    ptr                         # actual pointer, removed type ::Cptr for serialization
    len::Int                    # size in bytes
    dev::Int                    # id of the device the pointer belongs to
    parent	                # used to implement shared memory pointers
end

# This is the low level KnetPtr constructor, it adds the finalizer and sets parent to
# nothing, which is only needed for shared pointers.

function KnetPtr(ptr::Cptr,len::Int,dev::Int)
    kp = KnetPtr(ptr,len,dev,nothing)
    finalizer(freeKnetPtr, kp)
end

# This constructor is used to create a shared pointer.  We need to keep the parent field to
# prevent premature gc of the parent.  The child does not need a special finalizer.

function KnetPtr(parent::KnetPtr, offs::Int, len::Int)
    if len < 0 || offs < 1 || offs+len-1 > parent.len; throw(BoundsError()); end
    KnetPtr(parent.ptr+offs-1, len, parent.dev, parent)
end

# When Julia gc reclaims a KnetPtr object, the following special finalizer does not actually
# release the memory, but inserts it back in the appropriate pool for reuse.

function freeKnetPtr(p::KnetPtr)
    if p.ptr == C_NULL
        # already freed, do nothing
    elseif p.parent isa Nothing
        @dbg (push!(arraysizes,-p.len); push!(blocksizes,-p.len))
        mem = KnetMems[p.dev+1]
        mem.bfree += p.len
        mem.kfree += 1
        push!(mem.pools[p.len].free, p.ptr)
        p.ptr = C_NULL # to avoid double free by gcnode then gc.
    elseif p.parent isa KnetPtr
        # subarray, do nothing
    else # p.parent isa CuArray
        freeKnetPtrCu(p)
    end
end

# We use the KnetPool type to keep track of allocated and garbage collected pointers: We
# keep one KnetPool struct per size per device.

mutable struct KnetPool
    free::Vector{Cptr}          # pointers available for reuse
    nptr::Int                   # number of allocated pointers
    KnetPool()=new(Cptr[],0)
end

const KNETMEMINIT = 1<<24       # initial gpu memory limit

# KnetMem type keeps memory information for one device.
mutable struct KnetMem
    pools::Dict{Int,KnetPool}   # pointers of a given size
    limit::Int                  # current memory limit
    bytes::Int                  # total bytes allocated (inuse + avail)
    bfree::Int                  # total bytes freed and available
    kptrs::Int                  # number of arrays allocated (inuse + avail)
    kfree::Int                  # number of arrays freed and available
    gc::Int                     # number of times GC.gc called
    knetgc::Int                 # number of times Knet.gc called
    gctime::UInt
end
KnetMem()=KnetMem(Dict{Int,KnetPool}(),KNETMEMINIT,0,0,0,0,0,0,0)

# KnetMems[dev+1] holds memory information for device dev.
KnetMems = nothing
initKnetMems() = (global KnetMems = [ KnetMem() for i in 1:gpuCount() ])
knetmem(dev=gpu()) = (if KnetMems == nothing; initKnetMems(); end; KnetMems[dev+1])

# Blocksize determines the actual allocation size given the array size in bytes, and can be
# larger than what the array needs for increased reuse.
# blocksize(n::Int)=n
# function blocksize(n::Int)
#     b = sizeof(n)<<3
#     z = leading_zeros(n-1)
#     1<<(b-z)
# end
blocksize(n::Int,b=cbrt(2))=floor(Int,b^ceil(log(b,n)))

# The following used for debugging and record every request
arraysizes = Int[]; allocs = Int[]; blocksizes = Int[]

# This the main KnetPtr constructor.  It tries to avoid actual allocation which is slow.
# Reusing a preallocated and garbage collected pointer is very fast.
# Allocating a new pointer is about 0.5ms.
# GC.gc() is about 100ms.
# Knet.gc() is about 250ms.

gc_interval() = 2*10^8  # gc interval in ns, optimized on seq2seq model, balancing costs of alloc, GC.gc, Knet.gc
putc(c)=nothing         # putc(c)=print(c) to observe GC.gc, Knet.gc and inclimit

function KnetPtr(arraybytes::Int)
    @assert gpu() >= 0 "Cannot use KnetArray without a GPU."
    cuallocator() && return KnetPtrCu(arraybytes)
    dev = gpu(); @assert dev >= 0 "KnetPtr: bad device id $dev."
    mem = knetmem(dev)
    blockbytes = blocksize(arraybytes)
    @dbg (push!(arraysizes,arraybytes); push!(blocksizes,blockbytes))
    pool = get!(KnetPool,mem.pools,blockbytes)

    ptr = reuse(mem, pool, blockbytes, 0) # 0. best case we have one available in pool
    ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    if pool.nptr > 0
        ptr = reuse(mem, pool, blockbytes, 1, trygc=false) # 1. try fast gc (~0.5 ms)
        ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    end
    if mem.bytes + blockbytes <= mem.limit # 2. allocate if within limit (~0.5 ms)
        ptr = alloc(mem, pool, blockbytes, 2)
        ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    end
    if pool.nptr > 0 && time_ns() - mem.gctime > gc_interval() # 3. try slow gc (~100 ms) if enough time passed
        ptr = reuse(mem, pool, blockbytes, 3, trygc=true)
        ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    end
    if maybe_inclimit!(mem, max(mem.bytes + blockbytes*2, mem.limit*6÷5)) # 4. try to increase limit
        ptr = alloc(mem, pool, blockbytes, 4)
        ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    end
    if pool.nptr > 0  # 5. try slow gc (~100ms)
        ptr = reuse(mem, pool, blockbytes, 5, trygc=true)
        ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    end
    Knet.gc()  # 6. last ditch effort: ~250 ms + future cost of lost pools
    if mem.bytes + blockbytes <= mem.limit
        ptr = alloc(mem, pool, blockbytes, 6)
        ptr != nothing && return KnetPtr(ptr,blockbytes,dev)
    end
    error("Out of gpu memory")
end

function alloc(mem, pool, blockbytes, dbg)
    ptr = knetMalloc(blockbytes)       # ~584μs
    if ptr != nothing
        @dbg push!(allocs, dbg)
        mem.bytes += blockbytes
        mem.kptrs += 1
        pool.nptr += 1
        return ptr
    else
        return nothing
    end
end

function reuse(mem, pool, blockbytes, dbg; trygc=nothing)
    if trygc !== nothing
        GC.gc(trygc)            # gc(true)=slow, gc(false)=fast
        if trygc
            mem.gctime = time_ns(); mem.gc += 1; putc('-')
        end
    end
    if !isempty(pool.free)
        @dbg push!(allocs, dbg)
        mem.bfree -= blockbytes
        mem.kfree -= 1
        return pop!(pool.free)
    else
        return nothing
    end
end

function maybe_inclimit!(m::KnetMem, minlimit=m.limit*6÷5)
    maxlimit = m.bytes + gpufree() - 500_000_000 # gpumem()[1] - 500_000_000 # m.bytes + gpufree()
    if m.limit < minlimit <= maxlimit
        putc('^')
        m.limit = minlimit
        @debug kmeminfo()
        return true
    else
        return false
    end
end

KnetPtr(n::Integer)=KnetPtr(Int(n))

# This does the actual allocation, returns `nothing` in case of error
function knetMalloc(nbytes::Int) # 584μs
    ptr = Cptr[0]
    # ret = @cudart1(cudaMalloc,(Ptr{Cptr},Csize_t),ptr,nbytes)
    ret = CUDA.unsafe_cuMemAlloc_v2(ptr, nbytes)
    ret == 0 ? ptr[1] : nothing
end

"""
    Knet.gc(dev=gpu())

cudaFree all pointers allocated on device `dev` that were previously allocated and garbage
collected. Normally Knet holds on to all garbage collected pointers for reuse. Try this if
you run out of GPU memory.
"""
function gc(dev=gpu())
    if KnetMems == nothing; return; end
    putc('+')
    mem = knetmem(dev)
    mem.knetgc += 1
    GC.gc(); GC.enable(false)
    for (n,v) in mem.pools
        for p in v.free
            # @cudart(cudaFree,(Cptr,),p)
            CUDA.cuMemFree_v2(CuPtr{Nothing}(UInt(p)))
        end
        v.nptr -= length(v.free)
        mem.kptrs -= length(v.free)
        mem.kfree -= length(v.free)
        mem.bfree -= n * length(v.free)
        mem.bytes -= n * length(v.free)
        empty!(v.free)
    end
    GC.enable(true); GC.gc()
end

function knetgc()
    @warn "knetgc is deprecated, use Knet.gc instead" maxlog=1
    Knet.gc()
end

# Some utilities

meminfo(i=gpu())=[(k,v.nptr,length(v.free)) for (k,v) in knetmem(i).pools]

using Printf
kmeminfo(i=gpu())=(m=knetmem(i); @sprintf("knetgc=%g gc=%g pools=%g kptrs=%g kfree=%g limit=%g bytes=%g bfree=%g", m.knetgc, m.gc, length(m.pools), m.kptrs, m.kfree, m.limit, m.bytes, m.bfree))

function gpuinfo(msg="",dev=gpu();n=10)
    msg != "" && print("$msg ")
    if nvmlfound # Libdl.find_library(["libnvidia-ml"],[]) != ""
        g = nvmlDeviceGetMemoryInfo(nvmlid(dev))
        println((:dev,dev,:total,g[1],:free,g[2],:used,g[3]))
    else
        dev0 = gpu(); gpu(dev)
        (free,total) = cudaMemGetInfo()
        gpu(dev0)
        println((:dev,dev,:total,total,:free,free,:used,total-free))
    end
    bytes = bfree = ptrs = k = 0
    for (s,u,f) in sort(meminfo(dev), by=(x->x[1]*x[2]), rev=true)
        bytes += s*u; bfree += s*f; ptrs += u; k += 1
        if k <= n; println((:bytes,s*u,:size,s,:alloc,u,:bfree,f)); end
    end
    if n < k; println('⋮'); end
    println((:totbytes,bytes,:bfree,bfree,:ptrs,ptrs))
    if KnetMems != nothing
        mem = KnetMems[dev+1]
        println((:membytes,mem.bytes,:bfree,mem.bfree,:limit,mem.limit,:pools,length(mem.pools)))
    end
end

function memdbg(msg="")
    if nvmlfound # Libdl.find_library(["libnvidia-ml"],[]) != ""
        m = nvmlDeviceGetMemoryInfo()
    else
        m = (0,0,0)
    end
    c = cudaMemGetInfo()
    x = [0,0]
    for (s,u,f) in meminfo()
        x[1] += s*u
        x[2] += s*f
    end
    println("""$msg:
cudaMemGetInfo: ctotal: $(c[2]) cfree: $(c[1]) ctotal-cfree: $(c[2]-c[1])
nvmlDeviceGetMemoryInfo: ntotal: $(m[1]) nfree: $(m[2]) nused: $(m[3]) nfree+used: $(m[2]+m[3])
KnetPtr: ktotal: $(x[1]) kavail: $(x[2]) ktotal-kavail: $(x[1]-x[2])
nused-ktotal: $(m[3]-x[1])
""")
end