text:
using Plots, Suppressor
includet("draftpolicy.jl")
ENV["GKSwstype"]=100
#N=3 and G=5 is feasible ~10 minutes

N = 3
G = 4
abilities = zeros(N).+1/N
s = Season(N,G,abilities)
pi_t = 10
pi_b = 1
nsims = pi_t - pi_b+1
obj_opt = zeros(nsims)
obj_ntd = zeros(nsims)

#a,b, opt_wts, ntd_wts = calc_objective(s,5)
for mult in pi_b:pi_t
    opt_wts = zeros(length(s.matches))
    ntd_wts= zeros(length(s.matches))
    i = mult-pi_b+1
    obj_opt[i], obj_ntd[i], opt_wts, ntd_wts = calc_objective(s, mult)
    if mult ==10
        println(obj_opt)
        println(obj_ntd)
        plot([ opt_wts ntd_wts],xlabel="Match",ylabel="Weight",ylims=[0,0.2],labels=["Optimal Weight" "T-IC Weight"])
        savefig("figures/weights10.pdf")
    end
end

obj_unif= zeros(nsims) .+ (1- 1/(N-1))^2 .+ (N-2)*(0-1/(N-1))^2
obj_ante= zeros(nsims) .+ (1- 1/(N))^2 .+ (N-1)*(0-1/(N))^2
plot(pi_b:pi_t,[ obj_ante obj_unif obj_opt obj_ntd], ylims=[0,1.1], xlabel="B", ylabel="Objective",label=["Ex-Ante Uniform" "Ex-Post Uniform" "Optimal" "T-IC"])
savefig("figures/simulated_gain.pdf")

obj_ante

max_stars_repo_name: evanmunro/draft-policy

max_issues_repo_name: evanmunro/draft-policy

max_forks_repo_name: evanmunro/draft-policy


text:
using Revise
using NNFEM
using PyPlot
using LinearAlgebra
using ADCME
using ADCMEKit
using MAT

NT = 100
Δt = 1/NT

node, elem = load_mesh("holes")

elements = SmallStrainContinuum[]
prop = Dict("name"=> "PlaneStress", "rho"=> 8000,  "E"=>200e9, "nu"=>0.3)
for i = 1:size(elem,1)
    nodes = node[elem[i,:], :]
    elnodes = elem[i,:]
    ngp = 3
    push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
end

Edge_Traction_Data = Array{Int64}[]
for i = 1:length(elements)
    elem = elements[i]
    for k = 1:4
        if elem.coords[k,1]>99.5 && elem.coords[k+1>4 ? 1 : k+1,1]>99.5
            push!(Edge_Traction_Data, [i, k, 1])
        end
    end
end
Edge_Traction_Data = hcat(Edge_Traction_Data...)'|>Array

EBC = zeros(Int64, size(node, 1), 2)
FBC = zeros(Int64, size(node, 1), 2)
g = zeros(size(node, 1), 2)
f = zeros(size(node, 1), 2)
for i = 1:size(node, 1)
    if node[i,1]<0.05
        EBC[i,:] .= -1
    end
end

domain = Domain(node, elements, 2, EBC, g, FBC, f, Edge_Traction_Data)

Dstate = zeros(domain.neqs)
state = zeros(domain.neqs)
velo = zeros(domain.neqs)
acce = zeros(domain.neqs)
EBC_func = nothing
FBC_func = nothing
Body_func = nothing
# Construct Edge_func
function Edge_func(x, y, t, idx)
  return [1e6*ones(length(x)) zeros(length(x))]
end
globaldata = GlobalData(state, Dstate, velo, acce, domain.neqs, EBC_func, FBC_func,Body_func, Edge_func)

d0 = zeros(2domain.nnodes)
Fext = getExternalForce(domain, globaldata)

H = constant(domain.elements[1].mat[1].H)

sess = Session(); init(sess)

# incremental load method
# d = ImplicitStaticSolver(globaldata, domain, d0, Δt, NT, H, Fext, method="incremental")
# d_ = run(sess, d)
# visualize_displacement(d_*25199.1, domain)

# # newton's method
# d = ImplicitStaticSolver(globaldata, domain, d0, Δt, NT, H, Fext, method="newton")
# sess = Session(); init(sess)
# d_ = run(sess, d)
# visualize_displacement(reshape(d_, 1, :)*25199.1, domain)

# newton's method (also work for nonlinear equations)

H = domain.elements[1].mat[1].H
function nn(ε, θ)
    H = reshape(θ, (3,3))
    Hs = repeat(reshape(H, (-1,)), getNGauss(domain))
    Hs = reshape(Hs, (getNGauss(domain), 3, 3) )
    ε*H, Hs
end
ADCME.options.newton_raphson.verbose = true
ADCME.options.newton_raphson.rtol = 1e-6
ADCME.options.newton_raphson.tol = 1e-6
d = ImplicitStaticSolver(globaldata, domain, d0,  nn, H[:], Fext)
d_ = run(sess, d)
visualize_displacement(reshape(d_, 1, :)*25199.1, domain)

# println("Press any key to learn the linear elastic matrix...")
# readline()

# H = spd(Variable(rand(3,3))) * 1e10
# d = ImplicitStaticSolver(globaldata, domain, d0, Δt, NT, H, Fext)
# loss = sum((d - d_)^2)
# init(sess)
# BFGS!(sess, loss)

# println("Estimated Stiffness Matrix:\n$(run(sess, H))")
# println("Exact Stiffness Matrix:\n$(domain.elements[1].mat[1].H)")

ext: jl


text:
function preprocess(PData::PDataSpectral, H, g, n1, n2)
    #global V, D, rho, ncomp

    # Δ, V = eig(H)
    X = eigen(H)
    Δ = X.values
    V = X.vectors

    l_m, = findmin(Δ)
    g̃ = V' * g
    n_g = norm(g)
    ɛ = 1.0e-10
    λ = max(-l_m, 0.0)

    PData.V = V
    PData.Δ = Δ
    PData.g̃ = g̃
    PData.λ = λ
    PData.success = true
    PData.OK = true

    return PData
end

function AInv(X::PDataSpectral, d̃::Array{Float64,1})
    return X.V * d̃
end


text:
function AngularIonosphericDelayModifier(arg0::IonosphericModel, arg1::jdouble)
    return AngularIonosphericDelayModifier((IonosphericModel, jdouble), arg0, arg1)
end

function get_parameters_drivers(obj::AngularIonosphericDelayModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::AngularIonosphericDelayModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
using Test
using SafeTestsets

@safetestset "Parametric entity" begin include("parametricentity_test.jl") end

@safetestset "Meshgen" begin include("meshgen_test.jl") end

max_stars_repo_name: WaveProp/ParametricSurfaces

max_issues_repo_name: WaveProp/ParametricSurfaces

max_forks_repo_name: WaveProp/ParametricSurfaces


text:
#__precompile__()
module SeuratRDS

using Pkg
using Conda
using Dates
using DelimitedFiles
using DataFrames

ENV["R_HOME"]="*";Pkg.build("RCall");using RCall # make sure we have a good r environment

Conda.add("r-matrix",joinpath("/",split(RCall.libR,"/")[2:end-4]...)) # make sure Matrix is in the new conda env

export loadSeur

# return features x barcodes data as tuple with data matrix, metadata, colnames, rownames, dataframe representation
function loadSeur(rdsPath::String,modality::String,assay::String,metadata::String)
    R"""
      rdsPath = $rdsPath;
      seur = readRDS(rdsPath);"""; # read in annotations

    # export counts and embedded labels
    R"""
      modality = $modality;
      assay = $assay;
      metadata = $metadata;
      modl = get($modality,slot(seur,"assays"));
      dat = slot(modl,$assay);
      met = get(metadata,slot(seur,"meta.data"));
      cnm = colnames(dat);
      rnm = rownames(dat);"""

    dat = rcopy(R"as.matrix(dat)")
    met = rcopy(R"as.matrix(met)")
    cnm = rcopy(R"as.matrix(cnm)")
    rnm = rcopy(R"as.matrix(rnm)")

    df = DataFrame(dat)
    rename!(df,Symbol.(reduce(vcat,cnm)))
    insertcols!(df,1,(:gene=>reduce(vcat,rnm)))

    (dat=dat,met=met,col=cnm,row=rnm,df=df)
end

# convert to loadSeur output to barcodes x features dataframe and add labels column
function bcFeatLabels(seurData::NamedTuple,labels::Vector)
  df = DataFrame(seurData.dat')
  rename!(df,Symbol.(reduce(vcat,seurData.row)))
  insertcols!(df,1,(:barcode=>reduce(vcat,seurData.col)))
  insertcols!(df,1,(:label=>reduce(vcat,labels)))
  df
end

# return barcodes x features data where the metadata::Dict includes extra features like:
# :featurename => metadata, where metadata corresponds to get(metadata,slot(seur,"meta.data"))
function loadSeur(rdsPath::String,
                  modality::String,assay::String,
                  metadata::Dict)
    R"""
      library(Matrix)
      rdsPath = $rdsPath;
      seur = readRDS(rdsPath);
      modality = $modality;
      assay = $assay;
      modl = get($modality,slot(seur,"assays"));
      dat = slot(modl,$assay);
      cnm = colnames(dat);
      rnm = rownames(dat);"""; # read in annotations

    dat = rcopy(R"as.matrix(dat)")
    cnm = rcopy(R"as.matrix(cnm)")
    rnm = rcopy(R"as.matrix(rnm)")

    df = DataFrame(dat')
    rename!(df,Symbol.(reduce(vcat,rnm)))
    insertcols!(df,1,(:barcode=>reduce(vcat,cnm)))

    # export counts and embedded labels
    for k in keys(metadata)
      seurkey = get(metadata,k,"")
      R"""
        met = get($seurkey,slot(seur,"meta.data"));"""
      met = rcopy(R"as.matrix(met)")
      insertcols!(df,1,(k=>reduce(vcat,met)))
    end
    df
end

end


text:
using COMTRADE
using Base.Test

z1 = read_comtrade("data/1991-ascii-sel")
z2 = read_comtrade("data/1999-ascii-pq")
z3 = read_comtrade("data/1999-binary-hif-sel")

@test z1.filename == "data/1991-ascii-sel"
@test z1.cfg.rev_year == 1991
@test z1.cfg.nA == 24
@test z1.cfg.A[2, :ch_id] == "IB"
@test z1.cfg.D[2, :ch_id] == "TRP"
@test z1.cfg.triggertime > z1.cfg.time
@test z1.dat[1, :IA] ≈ -270.999876
@test z1.dat[1, S"VA(kV)"] ≈ -33.3998801
@test !z1.dat[1, :TRP]

@test z2.filename == "data/1999-ascii-pq"
@test z2.cfg.rev_year == 1999
@test z2.cfg.nA == 6
@test z2.cfg.A[2, :ch_id] == "Ib"
@test z2.cfg.triggertime > z1.cfg.time
@test z2.dat[1, :Ia] ≈ 101.06138883816476

@test z3.filename == "data/1999-binary-hif-sel"
@test z3.cfg.rev_year == 1999
@test z3.cfg.nA == 18
@test z3.cfg.A[2, :ch_id] == "IBRMS"
@test z3.cfg.D[2, :ch_id] == "NTUNE_A"
@test z3.cfg.triggertime > z1.cfg.time
@test z3.dat[1, :IARMS] ≈ 0.0
@test z3.dat[1, :EN]


text:
using Base.Test

@testset "Construction and Initialization" begin
# ============================================= #
@testset "reshape" begin
    a::Array{Float64,3} = Array{Float64,3}(3,3,3)

    # - reshape - #
    a = reshape(1:27,3,3,3)
    # - same as the code below
    # value::Int64 = 1
    # for k in 1:3
    #     for j in 1:3
    #         for i in 1:3
    #             a[i,j,k] = value
    #             value += 1
    #         end
    #     end
    # end
    @test a[1,1,1] ≈ 1.
    @test isapprox(a[2,1,1], 2.)
    @test a[3,3,3] ≈ 27.
end
# ============================================= #
# - Comprehension - #
@testset "Comprehensions" begin
    # a =
    # [1.0 4.0 7.0;
    #  2.0 5.0 8.0;
    #  3.0 6.0 9.0]
    # [10.0 13.0 16.0;
    #  11.0 14.0 17.0;
    #  12.0 15.0 18.0]
    # [19.0 22.0 25.0;
    #  20.0 23.0 26.0;
    #  21.0 24.0 27.0]
    a::Array{Float64,3} = Float64[i+3(j-1)+9(k-1) for i=1:3, j=1:3, k=1:3]
    @test a[1,1,1] ≈ 1.
    @test isapprox(a[2,1,1], 2.)
    @test a[3,3,3] ≈ 27.
end
# ============================================= #
end


text:
# Classical Multidimensional Scaling

## convert Gram matrix to Distance matrix

function gram2dmat!{DT}(D::AbstractMatrix{DT}, G::AbstractMatrix)
    # argument checking
    m = size(G, 1)
    n = size(G, 2)
    m == n || error("D should be a square matrix.")
    size(D) == (m, n) ||
        throw(DimensionMismatch("Sizes of D and G do not match."))

    # implementation
    for j = 1:n
        for i = 1:j-1
            @inbounds D[i,j] = D[j,i]
        end
        D[j,j] = zero(DT)
        for i = j+1:n
            @inbounds D[i,j] = sqrt(G[i,i] + G[j,j] - 2 * G[i,j])
        end
    end
    return D
end

gram2dmat{T<:Real}(G::AbstractMatrix{T}) = gram2dmat!(similar(G, Base.momenttype(T)), G)

## convert Distance matrix to Gram matrix

function dmat2gram!{GT}(G::AbstractMatrix{GT}, D::AbstractMatrix)
    # argument checking
    m = size(D, 1)
    n = size(D, 2)
    m == n || error("D should be a square matrix.")
    size(G) == (m, n) ||
        throw(DimensionMismatch("Sizes of G and D do not match."))

    # implementation
    u = zeros(GT, n)
    s = 0.0
    for j = 1:n
        s += (u[j] = Base.sumabs2(view(D,:,j)) / n)
    end
    s /= n

    for j = 1:n
        for i = 1:j-1
            @inbounds G[i,j] = G[j,i]
        end
        for i = j:n
            @inbounds G[i,j] = (u[i] + u[j] - abs2(D[i,j]) - s) / 2
        end
    end
    return G
end

dmat2gram{T<:Real}(D::AbstractMatrix{T}) = dmat2gram!(similar(D, Base.momenttype(T)), D)

## classical MDS

function classical_mds{T<:Real}(D::AbstractMatrix{T}, p::Integer)
    n = size(D, 1)
    p < n || error("p must be less than n.")

    G = dmat2gram(D)
    E = eigfact!(Symmetric(G))
    ord = sortperm(E.values; rev=true)
    (v, U) = extract_kv(E, ord, p)
    for i = 1:p
        @inbounds v[i] = sqrt(v[i])
    end
    scale!(U, v)'
end


text:
export optimize!

function optimize!(loss::IDAPBCLoss{P}, paramvec::AbstractVector, data;
    η=0.001, batchsize=64, maxiters=1e4, tol=1e-4) where
    {J2,M,MD,V,VD<:Function,P<:IDAPBCProblem{J2,M,MD,V,VD}}

    batchgs = Vector{typeof(paramvec)}(undef, batchsize)
    dataloader = Flux.Data.DataLoader(data; batchsize=batchsize, shuffle=true)
    max_batch = round(Int, dataloader.imax / dataloader.batchsize, RoundUp)
    optimizer = Flux.ADAM(η)
    nepoch = 0;
    train_loss = 1/length(data) * pmap(x->loss(x,paramvec), data)
    losstypestr = typeof(loss).name.name |> string
    while nepoch < maxiters && train_loss > tol
        estatus(losstypestr, nepoch, train_loss, maxiters)
        nbatch = 1
        for batch in dataloader
            npoints = length(batch)
            batchloss = 1/npoints * pmap(x->loss(x,paramvec), batch)
            bstatus(nbatch, max_batch, batchloss)
            pmap!(batchgs, x->gradient(loss,x,paramvec), batch)
            grads = 1/npoints * sum(batchgs[1:npoints])
            if !any(isnan.(grads))
                Flux.Optimise.update!(optimizer, paramvec, grads)
                nbatch += 1
            else
                @warn "∂ℓ/∂θ computations yielded NaNs."
            end
        end
        nepoch += 1
        train_loss = 1/length(data) * pmap(x->loss(x,paramvec), data)
    end
    estatus(losstypestr, nepoch, 1/length(data) * pmap(x->loss(x,paramvec), data), maxiters)
end

function optimize!(loss::IDAPBCLoss{P}, paramvec::Flux.Params, data;
    η=0.001, batchsize=64, maxiters=1e4, tol=1e-4) where
    {J2,M,MD,V,VD<:Chain,P<:IDAPBCProblem{J2,M,MD,V,VD}}

    batchgs = Vector{Zygote.Grads}(undef, batchsize)
    dataloader = Flux.Data.DataLoader(data; batchsize=batchsize, shuffle=true)
    max_batch = round(Int, dataloader.imax / dataloader.batchsize, RoundUp)
    optimizer = Flux.ADAM(η)
    nepoch = 0;
    train_loss = 1/length(data) * pmap(x->loss(x), data)
    losstypestr = typeof(loss).name.name |> string
    while nepoch < maxiters && train_loss > tol
        estatus(losstypestr, nepoch, train_loss, maxiters)
        nbatch = 1
        for batch in dataloader
            npoints = length(batch)
            batchloss = 1/npoints * pmap(loss, batch)
            bstatus(nbatch, max_batch, batchloss)
            pmap!(batchgs, x->gradient(loss,x,paramvec), batch)
            grads = @. +(batchgs[1:npoints]...)
            foreach(grads) do x
                x = x/npoints
            end
            if !any(mapreduce(x->vec(isnan.(x)), vcat, grads))
                Flux.Optimise.update!(optimizer, paramvec, grads)
                nbatch += 1
            else
                @warn "∂ℓ/∂θ computations yielded NaNs."
            end
        end
        nepoch += 1
        train_loss = 1/length(data) * pmap(loss, data)
    end
    estatus(losstypestr, nepoch, 1/length(data) * pmap(loss, data), maxiters)
end

bstatus(b,bs,l) = @printf("Batch %05d/%05d | %10.4e\r", b, bs, l)
estatus(t,e,l,maxiters) = @printf("EPOCH %05d/%05d | TRAIN LOSS (%s) = %10.4e\n", e, maxiters, t, l)

function pmap(f, batch::Vector{Vector{T}}) where {T<:Real}
    l = Threads.Atomic{eltype(first(batch))}(0)
    n = length(batch)
    Threads.@threads for id in 1:n
        Threads.atomic_add!(l, f(batch[id]))
    end
    return l.value
end

function pmap!(out, f, batch::Vector{Vector{T}}) where {T<:Real}
    n = length(batch)
    Threads.@threads for id in 1:n
        out[id] = f(batch[id])
    end
end


text:
using Documenter, RatingCurves, Cairo, Fontconfig

CI = get(ENV, "CI", nothing) == "true"

makedocs(sitename = "RatingCurves.jl",
    format = Documenter.HTML(
    prettyurls = CI,
    ),
    pages = [
       "index.md",
       "Tutorial" =>["Getting started" => "tutorial/index.md",
            "Data" => "tutorial/gauging.md",
            "Rating curve fitting" => "tutorial/rcfit.md",
            "Compound rating curve fitting" => "tutorial/crcfit.md"],
       "contributing.md",
       "functions.md"]
)

# if CI
#     deploydocs(
#     repo   = "github.com/JuliaExtremes/RatingCurves.jl.git",
#     devbranch = "dev",
#     versions = ["stable" => "v^", "v#.#"],
#     push_preview = false,
#     target = "build"
#     )
# end

deploydocs(
    repo = "github.com/JuliaExtremes/RatingCurves.jl.git",
)


text:
render(ast) = ast |> to_math |> to_latex |> latexstring |> display

"""
Renders the value of `v`
$(SIGNATURES)
"""
renderval(v) = v |> repr |> parse |> render

"Render an expression. Eg `@render x^(2/(y+1))`"
macro render(ast)
	:(render($(Meta.quot(ast))))
end


text:
using ChainRulesCore

ChainRulesCore.debug_mode() = true

#####
##### `eigen`
#####

@Zygote.adjoint function LinearAlgebra.eigen(A::LinearAlgebra.RealHermSymComplexHerm)
    dU = eigen(A)
    return dU, function (Δ)
      d, U = dU
      d̄, Ū = Δ
      if Ū === nothing
        P = Diagonal(d̄)
      else
        F = inv.(d' .- d)
        P = F .* (U' * Ū)
        if d̄ === nothing
          P[diagind(P)] .= 0
        else
          P[diagind(P)] = d̄
        end
      end
      return (U * P * U',)
    end
  end

Zygote.refresh()

# function rrule(::typeof(eigen), X::AbstractMatrix{<:Real})
#     F = eigen(X)
#     function eigen_pullback(Ȳ::Composite)
#         # `getproperty` on `Composite`s ensures we have no thunks.
#         ∂X = eigen_rev(F, Ȳ.values, Ȳ.vectors)
#         return (NO_FIELDS, ∂X)
#     end
#     return F, eigen_pullback
# end

# function rrule(::typeof(getproperty), F::T, x::Symbol) where T <: Eigen
#     function getproperty_eigen_pullback(Ȳ)
#         C = Composite{T}
#         ∂F = if x === :values
#             C(values=Ȳ,)
#         elseif x === :vectors
#             C(vectors=Ȳ,)
#         # elseif x === :V
#         #     C(V=Ȳ,)
#         # elseif x === :Vt
#         #     # TODO: https://github.com/JuliaDiff/ChainRules.jl/issues/106
#         #     throw(ArgumentError("Vt is unsupported; use V and transpose the result"))
#         end
#         return NO_FIELDS, ∂F, DoesNotExist()
#     end
#     return getproperty(F, x), getproperty_eigen_pullback
# end

# # When not `Zero`s expect `λ̄ ::AbstractVector, Ū::AbstractMatrix`
# function eigen_rev(res::Eigen, λ̄ , Ū)
#     λ = res.values
#     U = res.vectors

#     k = length(λ)
#     T = eltype(λ)
#     F = T[i == j ? 0 : inv(@inbounds λ[j] - λ[i]) for i = 1:k, j = 1:k]
#     D = Diagonal(λ)
#     D̄ = λ̄  isa AbstractZero ? λ̄  : Diagonal(λ̄ )
#     Ut = U'
#     Ā = conj.(U) * ( D̄ + ( F .* ( Ut * Ū ) ) ) * Ut
#     return Ā
# end

######################################################################################

# Manifolds.jl Forward mode eigen AD by mateuszbaran
# https://github.com/JuliaManifolds/Manifolds.jl/pull/27

##
using LinearAlgebra
using ForwardDiff

import LinearAlgebra.eigen
function eigen(A::StridedMatrix{<:ForwardDiff.Dual})
    A_values = map(d -> d.value, A)
    A_values_eig = eigen(A_values)
    UinvAU = A_values_eig.vectors \ A * A_values_eig.vectors
    vals_diff = diag(UinvAU)
    F = similar(A_values)
    for i ∈ axes(A_values, 1), j ∈ axes(A_values, 2)
        if i == j
            F[i, j] = 0
        else
            F[i, j] = inv(A_values_eig.values[j] - A_values_eig.values[i])
        end
    end
    vectors_diff = A_values_eig.vectors * (F .* UinvAU)
    for i ∈ eachindex(vectors_diff)
        vectors_diff[i] = ForwardDiff.Dual{ForwardDiff.tagtype(vectors_diff[i])}(A_values_eig.vectors[i], vectors_diff[i].partials)
    end
    return Eigen(vals_diff, vectors_diff)
end

function f(A)
    A_eig = eigen(A)
    return A_eig.vectors * Diagonal(A_eig.values) / A_eig.vectors
end
##
ForwardDiff.derivative(t -> f([t+1 t; t t+1]), 1.0)
##
ForwardDiff.jacobian(x -> f(x), [1.0 2.0; 4.0 1.0])

##

######################################################################################

### Legume code

# def vjp_maker_eigsh(ans, x, **kwargs):
#     """Gradient for eigenvalues and vectors of a hermitian matrix."""
#     numeig = kwargs['k']
#     N = x.shape[-1]
#     w, v = ans              # Eigenvalues, eigenvectors.
#     vc = np.conj(v)

#     def vjp(g):
#         wg, vg = g          # Gradient w.r.t. eigenvalues, eigenvectors.
#         w_repeated = np.repeat(w[..., np.newaxis], numeig, axis=-1)

#         # Eigenvalue part
#         vjp_temp = np.dot(vc * wg[..., np.newaxis, :], T(v))

#         # Add eigenvector part only if non-zero backward signal is present.
#         # This can avoid NaN results for degenerate cases if the function
#         # depends on the eigenvalues only.
#         if np.any(vg):
#             off_diag = np.ones((numeig, numeig)) - np.eye(numeig)
#             F = off_diag / (T(w_repeated) - w_repeated + np.eye(numeig))
#             vjp_temp += np.dot(np.dot(vc, F * np.dot(T(v), vg)), T(v))

#         return vjp_temp

#     return vjp

# def vjp_maker_eigsh(ans, mat, **kwargs):
#     """Steven Johnson method extended to a Hermitian matrix
#     https://math.mit.edu/~stevenj/18.336/adjoint.pdf
#     """
#     numeig = kwargs['k']
#     N = mat.shape[0]

#     def vjp(g):
#         vjp_temp = np.zeros_like(mat)
#         for iv in range(numeig):
#             a = ans[0][iv]
#             v = ans[1][:, iv]
#             vc = np.conj(v)
#             ag = g[0][iv]
#             vg = g[1][:, iv]

#             # Eigenvalue part
#             vjp_temp += ag*np.outer(vc, v)

#             # Add eigenvector part only if non-zero backward signal is present.
#             # This can avoid NaN results for degenerate cases if the function
#             # depends on the eigenvalues only.
#             if np.any(vg):
#                 # Projection operator on space orthogonal to v
#                 P = np.eye(N, N) - np.outer(vc, v)
#                 Amat = T(mat - a*np.eye(N, N))
#                 b = P.dot(vg)

#                 # Initial guess orthogonal to v
#                 v0 = P.dot(np.random.randn(N))

#                 # Find a solution lambda_0 using conjugate gradient
#                 (l0, _) = sp.linalg.cg(Amat, b, x0=v0, atol=0)
#                 # Project to correct for round-off errors
#                 l0 = P.dot(l0)

#                 vjp_temp -= np.outer(l0, v)

#         return vjp_temp


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: auto_scaling_plans
using AWS.Compat
using AWS.UUIDs

"""
    CreateScalingPlan()

Creates a scaling plan.

# Required Parameters
- `ApplicationSource`: A CloudFormation stack or set of tags. You can create one scaling plan per application source.
- `ScalingInstructions`: The scaling instructions.
- `ScalingPlanName`: The name of the scaling plan. Names cannot contain vertical bars, colons, or forward slashes.

"""
create_scaling_plan(ApplicationSource, ScalingInstructions, ScalingPlanName; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("CreateScalingPlan", Dict{String, Any}("ApplicationSource"=>ApplicationSource, "ScalingInstructions"=>ScalingInstructions, "ScalingPlanName"=>ScalingPlanName); aws_config=aws_config)
create_scaling_plan(ApplicationSource, ScalingInstructions, ScalingPlanName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("CreateScalingPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationSource"=>ApplicationSource, "ScalingInstructions"=>ScalingInstructions, "ScalingPlanName"=>ScalingPlanName), args)); aws_config=aws_config)

"""
    DeleteScalingPlan()

Deletes the specified scaling plan. Deleting a scaling plan deletes the underlying ScalingInstruction for all of the scalable resources that are covered by the plan. If the plan has launched resources or has scaling activities in progress, you must delete those resources separately.

# Required Parameters
- `ScalingPlanName`: The name of the scaling plan.
- `ScalingPlanVersion`: The version number of the scaling plan.

"""
delete_scaling_plan(ScalingPlanName, ScalingPlanVersion; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("DeleteScalingPlan", Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion); aws_config=aws_config)
delete_scaling_plan(ScalingPlanName, ScalingPlanVersion, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("DeleteScalingPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion), args)); aws_config=aws_config)

"""
    DescribeScalingPlanResources()

Describes the scalable resources in the specified scaling plan.

# Required Parameters
- `ScalingPlanName`: The name of the scaling plan.
- `ScalingPlanVersion`: The version number of the scaling plan.

# Optional Parameters
- `MaxResults`: The maximum number of scalable resources to return. The value must be between 1 and 50. The default value is 50.
- `NextToken`: The token for the next set of results.
"""
describe_scaling_plan_resources(ScalingPlanName, ScalingPlanVersion; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlanResources", Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion); aws_config=aws_config)
describe_scaling_plan_resources(ScalingPlanName, ScalingPlanVersion, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlanResources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion), args)); aws_config=aws_config)

"""
    DescribeScalingPlans()

Describes one or more of your scaling plans.

# Optional Parameters
- `ApplicationSources`: The sources for the applications (up to 10). If you specify scaling plan names, you cannot specify application sources.
- `MaxResults`: The maximum number of scalable resources to return. This value can be between 1 and 50. The default value is 50.
- `NextToken`: The token for the next set of results.
- `ScalingPlanNames`: The names of the scaling plans (up to 10). If you specify application sources, you cannot specify scaling plan names.
- `ScalingPlanVersion`: The version number of the scaling plan. If you specify a scaling plan version, you must also specify a scaling plan name.
"""
describe_scaling_plans(; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlans"; aws_config=aws_config)
describe_scaling_plans(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlans", args; aws_config=aws_config)

"""
    GetScalingPlanResourceForecastData()

Retrieves the forecast data for a scalable resource. Capacity forecasts are represented as predicted values, or data points, that are calculated using historical data points from a specified CloudWatch load metric. Data points are available for up to 56 days.

# Required Parameters
- `EndTime`: The exclusive end time of the time range for the forecast data to get. The maximum time duration between the start and end time is seven days.  Although this parameter can accept a date and time that is more than two days in the future, the availability of forecast data has limits. AWS Auto Scaling only issues forecasts for periods of two days in advance.
- `ForecastDataType`: The type of forecast data to get.    LoadForecast: The load metric forecast.     CapacityForecast: The capacity forecast.     ScheduledActionMinCapacity: The minimum capacity for each scheduled scaling action. This data is calculated as the larger of two values: the capacity forecast or the minimum capacity in the scaling instruction.    ScheduledActionMaxCapacity: The maximum capacity for each scheduled scaling action. The calculation used is determined by the predictive scaling maximum capacity behavior setting in the scaling instruction.
- `ResourceId`: The ID of the resource. This string consists of the resource type and unique identifier.    Auto Scaling group - The resource type is autoScalingGroup and the unique identifier is the name of the Auto Scaling group. Example: autoScalingGroup/my-asg.   ECS service - The resource type is service and the unique identifier is the cluster name and service name. Example: service/default/sample-webapp.   Spot Fleet request - The resource type is spot-fleet-request and the unique identifier is the Spot Fleet request ID. Example: spot-fleet-request/sfr-73fbd2ce-aa30-494c-8788-1cee4EXAMPLE.   DynamoDB table - The resource type is table and the unique identifier is the resource ID. Example: table/my-table.   DynamoDB global secondary index - The resource type is index and the unique identifier is the resource ID. Example: table/my-table/index/my-table-index.   Aurora DB cluster - The resource type is cluster and the unique identifier is the cluster name. Example: cluster:my-db-cluster.
- `ScalableDimension`: The scalable dimension for the resource.
- `ScalingPlanName`: The name of the scaling plan.
- `ScalingPlanVersion`: The version number of the scaling plan.
- `ServiceNamespace`: The namespace of the AWS service.
- `StartTime`: The inclusive start time of the time range for the forecast data to get. The date and time can be at most 56 days before the current date and time.

"""
get_scaling_plan_resource_forecast_data(EndTime, ForecastDataType, ResourceId, ScalableDimension, ScalingPlanName, ScalingPlanVersion, ServiceNamespace, StartTime; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("GetScalingPlanResourceForecastData", Dict{String, Any}("EndTime"=>EndTime, "ForecastDataType"=>ForecastDataType, "ResourceId"=>ResourceId, "ScalableDimension"=>ScalableDimension, "ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion, "ServiceNamespace"=>ServiceNamespace, "StartTime"=>StartTime); aws_config=aws_config)
get_scaling_plan_resource_forecast_data(EndTime, ForecastDataType, ResourceId, ScalableDimension, ScalingPlanName, ScalingPlanVersion, ServiceNamespace, StartTime, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("GetScalingPlanResourceForecastData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EndTime"=>EndTime, "ForecastDataType"=>ForecastDataType, "ResourceId"=>ResourceId, "ScalableDimension"=>ScalableDimension, "ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion, "ServiceNamespace"=>ServiceNamespace, "StartTime"=>StartTime), args)); aws_config=aws_config)

"""
    UpdateScalingPlan()

Updates the specified scaling plan. You cannot update a scaling plan if it is in the process of being created, updated, or deleted.

# Required Parameters
- `ScalingPlanName`: The name of the scaling plan.
- `ScalingPlanVersion`: The version number of the scaling plan.

# Optional Parameters
- `ApplicationSource`: A CloudFormation stack or set of tags.
- `ScalingInstructions`: The scaling instructions.
"""
update_scaling_plan(ScalingPlanName, ScalingPlanVersion; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("UpdateScalingPlan", Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion); aws_config=aws_config)
update_scaling_plan(ScalingPlanName, ScalingPlanVersion, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = auto_scaling_plans("UpdateScalingPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion), args)); aws_config=aws_config)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule opencv_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("opencv")
JLLWrappers.@generate_main_file("opencv", UUID("296abffa-b662-51bb-8aac-8360b8502390"))
end  # module opencv_jll


text:
export get_state_timeseries, get_state_cumulative,
       get_vaccine_doses_administered_statewise, get_statewise_daily,
       get_active_cases_state

# 2. Statewise timeseries of Confirmed, Recovered and Deceased numbers.
function get_state_timeseries(state::String)
    df = read_CSV_url(states_csv_url)
    df.Active = df.Confirmed - df.Recovered - df.Deceased
    return df[df[!,:State].==state,:]
end

function get_state_timeseries()
    df = read_CSV_url(states_csv_url)
    df.Active = df.Confirmed - df.Recovered - df.Deceased
    return df
end

# 4. Statewise per day delta of Confirmed, Recovered and Deceased numbers.
function get_statewise_daily()
    df = read_CSV_url(state_wise_daily_csv_url)
    return df
end

# 5. Statewise cumulative numbers till date.
function get_state_cumulative(brief::Bool=true)
    df = read_CSV_url(state_wise_csv_url)
    if brief == true
        colindex = columnindex(df,:State_Notes)
        select!(df,Not(colindex))
    end
    return df
end

# Function to create state_code to state and vice versa dictionary objects.
function generate_dict_statecodes()
    df = get_state_cumulative()
    codeToStateDict = Dict{String,String}()
    stateToCodeDict = Dict{String,String}()
    nrows = size(df)[1]
    for i in 1:nrows
        state = df[i,:State]
        code = df[i,:State_code]
        codeToStateDict[code]=state
        stateToCodeDict[state]=code
    end

    return stateToCodeDict, codeToStateDict
end

state_to_code, code_to_state = generate_dict_statecodes()

# Function to calculate active cases on a particular date in all states.
function get_active_cases_state(date::Date)
    df = get_state_timeseries()
    df_today = df[df.Date.==date,:]
    sort!(df_today, :State)

    dateY = date-Dates.Day(1)# date yesterday
    df_yesterday = df[df.Date.==dateY,:]
    sort!(df_yesterday, :State)

    df_today.DeltaDeceased = df_today.Deceased - df_yesterday.Deceased
    df_today.DeltaActive = df_today.Active - df_yesterday.Active

    return df_today
end
# function calculate_active_cases_states(date::Date)
#     df = Covid19India.get_statewise_daily()
#     dfdate = df[df[!,:Date_YMD].==date,:]
#     dfout = DataFrame(State=String[], State_code=String[], Confirmed=Int[], Recovered=Int[], Deceased=Int[], Active=Int[], DeltaActive=Int[])
#     nstates = size(dfdate)[2]-3
#     colnames = names(dfdate)
#     for col_id in 1:nstates
#         state_code = colnames[col_id+3]
#         state=try
#             code_to_state[state_code]
#         catch e
#             continue
#         end
#         println(state)
#         confirmed = dfdate[dfdate.Status.=="Confirmed", state_code][1]
#         recovered = dfdate[dfdate.Status.=="Recovered", state_code][1]
#         deceased = dfdate[dfdate.Status.=="Deceased", state_code][1]

#         df_timeseries = Covid19India.get_state_timeseries(state)
#         #df_timeseries_today = df_timeseries[df_timeseries.Date.==date,:]
#         #df_timeseries_yesterday = df_timeseries[df_timeseries.Date=(date-Dates.Day(1))]
#         # active_today = df_timeseries_today[1,"Confirmed"]-df_timeseries_today[1,"Recovered"]-df_timeseries_today[1,"Deceased"]
#         # active_yesterday = df_timeseries_yesterday[1,"Confirmed"]-df_timeseries_yesterday[1,"Recovered"]-df_timeseries_yesterday[1,"Deceased"]
#         # active_delta = active_today - active_yesterday

#         # push!(dfout.State, code_to_state[state_code])
#         # push!(dfout.State_code, state_code)
#         # push!(dfout.Confirmed, confirmed)
#         # push!(dfout.Recovered, recovered)
#         # push!(dfout.Deceased, deceased)
#         # push!(dfout.Active, active_today)
#         # push!(dfout.DeltaActive, active_delta)
#     end

#     return df_timeseries
# end

# VACCINE DOSES ADMINISTERED.
# 12. Number of vaccine doses administered statewise [Reading]
function get_vaccine_doses_administered_timeseries_statewise()
    df = read_CSV_url(vaccine_doses_administered_statewise_csv_url)
    return df
end

function get_vaccine_doses_administered_timeseries_statewise(state::String)
    df = get_vaccine_doses_administered_timeseries_statewise()
    return df[df[!,:State].==state,:]
end

function get_vaccine_doses_administered_yesterday()
    df = get_vaccine_doses_administered_timeseries_statewise()
    #dttoday = Dates.format(Dates.now(), "dd/mm/yyyy")
    dtyesterday = Dates.format(Dates.now() - Dates.Day(1), "dd/mm/yyyy")
    dt_dby = Dates.format(Dates.now() - Dates.Day(2), "dd/mm/yyyy")
    dt_dbdby = Dates.format(Dates.now() - Dates.Day(3), "dd/mm/yyyy")
    df = df[!,["State",dt_dbdby, dt_dby,dtyesterday]]
    df.delta_today = df[!,4]-df[!,3]
    df.delta_yesterday = df[!,3]-df[!,2]
    df.delta_diff = df.delta_today - df.delta_yesterday
    sort!(df,:delta_today, rev=true)
    return df
end

# 13. Key data points from CoWin database at a state level
function get_cowin_vaccine_data_statewise(brief::Bool=false)
    df = read_CSV_url(cowin_vaccine_data_statewise_csv_url)
    if brief == true
        df = df[!,[]]
    end
    return df
end


text:
using CombinedParsers.Regexp
import CombinedParsers.Regexp: word, non_word
import CombinedParsers.Regexp: newline, inline, whitespace_maybe, whitespace_horizontal

@testset "CombinedParsers" begin
@testset "CharIn" begin
    @test parse(CharIn(isuppercase),"A") =='A'
    @test parse(map(v->length(v),re"a*"),"aaaa") == 4
end

@testset "map" begin
    @test parse(re"abc"[2:3],"abc")==('b','c')
    @test parse(map(v->length(v),re"a*"),"aaaa") == 4
    @test parse(re"^abc$"[2],"abc") == 'a'
end

@testset "NamedTuple" begin
    @test parse(Sequence(first = CharIn(isuppercase), second = CharIn(islowercase)),"Ab") == (first='A', second='b')
    @test parse(Sequence(CharIn(isuppercase), :second => CharIn(islowercase)),"Ab") == (second='b',)
end

@testset "Repeat" begin
    @test parse(join(Repeat('a'),","),"a,a") == ['a','a']
end

@testset "parse_all" begin
    @test collect(parse_all(Repeat("a"|"ab"),"aabab")) ==
        [ ["a","a"], ["a","ab","a"], ["a","ab","ab"], ["a","ab"], ["a"],[] ]
    @test collect(parse_all(re"^(a|ab|b)+$"[2],"abab")) ==
        [ ['a','b','a','b'], ['a','b', ('a','b')], [('a','b'),'a','b'], [('a','b'),('a','b')]]
end

@testset "Either" begin
    @test parse(Either("a","ab","ac")*AtEnd(),"ab") == ("ab", AtEnd())
    @test_throws ArgumentError parse(Atomic(Either("a","ab","ac"))*AtEnd(),"ab")
end

@testset "Repeat_until" begin
    @test parse(Sequence("(", Repeat_until(!!Either(word,non_word), ")")),
                "(balanced parenthesis)") ==
        ("(",["balanced", " ", "parenthesis"])
end
end

@testset "FlatMap" begin
    @test parse(
        after(with_log("left",Either("a","ab"))) do v
        parser(v)
        end,
        "abab")=="ab"
    # parse(pattern,with_options(Base.PCRE.MULTILINE,"^"))
end

attribute_parser =
    map(
        (v) -> (lowercase(v[1]) => "$(v[5])")::Pair{String,String},
        Sequence(
            !Repeat1(word), whitespace_maybe,"=", whitespace_maybe,
            Either(Sequence(2,"\"", Repeat_until(AnyChar(),"\"",wrap=JoinSubstring)),
                   Sequence(2,"'", Repeat_until(AnyChar(),"'",wrap=JoinSubstring)),
                   !re"[0-9]+%",
                   !re"[-+]?[0-9]+",
                   !Repeat1(word),
                   !re"#[0-9A-Fa-f]{6}")
        ));

attributes = join(attribute_parser, whitespace_maybe);

parse(attributes,"a = 1 b=6% font=\"+1asd\"")

@testset "continue options of last Either" begin
    @test parse(attributes, "size=10% class=1") == [ "size"=>"10%", "class"=>"1" ]
end

function html(tags::CombinedParser, inner::CombinedParser, attrs=attributes)
    html(result_type(inner), tags, attrs) do until
        Repeat_until(inner, until)
    end
end

function html(inner::Function, T::Type, tags::CombinedParser, attrs_parser=attributes)
    A = eltype(result_type(attrs_parser))
    function nested_html(x,)
        (tag,attrs) = x
        Either{Any}(map(parser("/>")) do v
                    (tag=tag, attrs=attrs, children=T[])
                    end,
                    Sequence(">",
                             inner(Sequence("</",tag,">"))) do v
                    (tag=tag, attrs=attrs, children=v[2])
                    end)
    end
    FlatMap{Any}(
        nested_html,
        Sequence(
            2,
            "<",
            Sequence(
                tags,
                Optional(Sequence(
                    2,
                    whitespace_horizontal,
                    attrs_parser,
                    whitespace_maybe)))))
end

@testset "html" begin
    inner = Either{Any}(Any[!Repeat(CharNotIn("<>"))]);
    pushfirst!(inner,html(!re"[[:alpha:]]+",inner,attributes));
    parse(inner,"ba")
    @test parse(inner,"iboldj") ==
        (tag="a", attrs=["font"=>"+1"], children=[
            "i",
            (tag="b", attrs=[], children=["bold"]),
            "j"])

    @test parse(inner,"") ==
        (tag="a", attrs=["font"=>"+1"], children=[])
end


text:
include("sample.jl")

using PyPlot
using Distances
###
### Box Constrained Minimization
###
using Optim, NLSolversBase #hide
import NLSolversBase: clear! #hide

fun(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;

x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)
lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)
res = optimize(df, dfc, x0, IPNewton())
if false
ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)
clear!(df)
res = optimize(df, dfc, x0, IPNewton())
con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    J
end
function con_h!(h, x, λ)
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
end;
lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
x0 = [0.2, 0.2]
res = optimize(df, dfc, x0, IPNewton())

lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res1 = optimize(df, dfc, x0, IPNewton())

function con2_c!(c, x)
    c[1] = x[1]^2 + x[2]^2     ## First constraint
    c[2] = x[2]*sin(x[1])-x[1] ## Second constraint
    c
end

function con2_jacobian!(J, x)
    # First constraint
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    # Second constraint
    J[2,1] = x[2]*cos(x[1])-1.0
    J[2,2] = sin(x[1])
    J
end

function con2_h!(h, x, λ)
    # First constraint
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
    # Second constraint
    h[1,1] += λ[2]*x[2]*-sin(x[1])
    h[1,2] += λ[2]*cos(x[1])
    # Symmetrize h
    h[2,1]  = h[1,2]
    h
end;

x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
                                     lx, ux, lc, uc)
res2 = optimize(df, dfc, x0, IPNewton())

if true
    f = x -> (x[1]^2 + x[2] -11)^2 + (x[1]+x[2]^2-7)^2
    figure(1)
    x = collect(Float16, range(-5,length=200,stop=5));
    y = collect(Float16, range(-5,length=200, stop=5));
    pairs = collect(Iterators.product(x, y));
    z = f.(pairs);
    #z = hcat(x, y); z = mapslices(fun, z, dims = 2)
    surf(x, y, z);
    (res, loc) = sample(f, [-5, -5], [5, 5]; num = 1000)
    scatter3D(loc[:, 1], loc[:, 2], res, color = "red")
    #contour(x, y, z)
    #scatter3D([0], [0], [100], color = "red")
end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
export ELFSymbols, ELFSymtabEntry, ELFSymtabEntry64, ELFSymtabEntry32, ELFSymbolRef

"""
    ELFSymbols

ELF symbol table, contains the list of symbols defined within the object file.
"""
struct ELFSymbols{H<:ELFHandle} <: Symbols{H}
    section_ref::SectionRef{H}
end

function Symbols(section::SectionRef{H}) where {H <: ELFHandle}
    return ELFSymbols(section)
end

"""
    Symbols(oh::ELFHandle)

Return the `ELFSymbols` object that is contained within the first `.dynsym` or
`.symtab` section within the ELF object.
"""
function Symbols(oh::H) where {H <: ELFHandle}
    sections = Sections(oh)
    # First, try to load `.symtab`.  If we can't load that guy, we'll fall back
    # to `.dynsym` which has less information, but is more likely to exist.
    dyn_sections = findall(sections, ".symtab")
    if !isempty(dyn_sections)
        return ELFSymbols(first(dyn_sections))
    end

    dyn_sections = findall(sections, ".dynsym")
    if !isempty(dyn_sections)
        return ELFSymbols(first(dyn_sections))
    end
    error("Could not find a .dynsym or .symtab section in object file")
end

"""
    Section(syms::ELFSymbols)

Return the `Section` that contains this ELF symbol table.
"""
Section(syms::ELFSymbols) = syms.section_ref
handle(syms::ELFSymbols) = handle(Section(syms))
symtab_entry_type(syms::ELFSymbols) = symtab_entry_type(handle(syms))

function lastindex(syms::ELFSymbols{H}) where {H <: ELFHandle}
    sect_size = section_size(Section(syms))
    sym_size = sizeof(symtab_entry_type(syms))
    return div(sect_size, sym_size)
end

# Add an StrTab() override to be able to load the symbol string table
"""
    StrTab(syms::ELFSymbols)

Given an ELF `Symbols` container, construct an `ELFStrTab` that refers to the
associated symbol string table.
"""
function StrTab(syms::ELFSymbols{H}) where {H <: ELFHandle}
    # The `sh_link` field in a symbol section always points to the symbol table
    section_idx = deref(Section(syms)).sh_link+1
    strtab_section = Sections(handle(syms))[section_idx]
    return ELFStrTab{H}(strtab_section)
end

"""
    ELFSymtabEntry

ELF Symtab Entry type, contains the data relevant to a symbol defined within
this ELF file, garnered from the symbol table (symtab)
"""
abstract type ELFSymtabEntry{H <: ELFHandle} <: SymtabEntry{H} end

@io struct ELFSymtabEntry32{H <: ELFHandle} <: ELFSymtabEntry{H}
    st_name::UInt32
    st_value::UInt32
    st_size::UInt32
    st_info::UInt8
    st_other::UInt8
    st_shndx::UInt16
end

@io struct ELFSymtabEntry64{H <: ELFHandle} <: ELFSymtabEntry{H}
    st_name::UInt32
    st_info::UInt8
    st_other::UInt8
    st_shndx::UInt16
    st_value::UInt64
    st_size::UInt64
end

function read(oh::ELFHandle, ::Type{ST}) where {ST <: ELFSymtabEntry}
    type_func = symtab_entry_type
    size_func = symtab_entry_size
    return read_struct(oh, type_func, size_func, "Symbol Entry")
end

# Define two helper functions that return ELF-specific symbol infos
st_bind(st_info::UInt8) = st_info >> 4
st_type(st_info::UInt8) = st_info & 0xf

# Define helper functions to detect symbols marked as "reserved"
function isreserved(sym::ELFSymtabEntry)
    return SHN_LORESERVE <= sym.st_shndx <= SHN_HIRESERVE
end

symbol_name(sym::ELFSymtabEntry) = string("strtab@", sym.st_name)
symbol_value(sym::ELFSymtabEntry) = sym.st_value
isundef(sym::ELFSymtabEntry) = (sym.st_shndx == SHN_UNDEF)
isglobal(sym::ELFSymtabEntry) = (st_bind(sym.st_info) & STB_GLOBAL) != 0
islocal(sym::ELFSymtabEntry) = !isglobal(sym)
isweak(sym::ELFSymtabEntry) = (st_bind(sym.st_info) & STB_WEAK) != 0

"""
    ELFSymbolRef

Contains a reference to an `ELFSymtabEntry`, as well as an `ELFSymbols`, etc...
"""
struct ELFSymbolRef{H<:ELFHandle} <: SymbolRef{H}
    syms::ELFSymbols{H}
    entry::ELFSymtabEntry{H}
    idx::UInt32
end
function SymbolRef(syms::ELFSymbols, entry::ELFSymtabEntry, idx)
    return ELFSymbolRef(syms, entry, UInt32(idx))
end

deref(sym::ELFSymbolRef) = sym.entry
Symbols(sym::ELFSymbolRef) = sym.syms

@derefmethod isreserved(sym::ELFSymbolRef)
symbol_number(sym::ELFSymbolRef) = sym.idx
function symbol_name(sym::ELFSymbolRef)
    # We look up the symbol name in the linked string table for this particular section
    strtab = StrTab(Symbols(sym))
    return strtab_lookup(strtab, deref(sym).st_name)
end
# Helper function to get the section a symbol refers to
symbol_section(sym::ELFSymbolRef) = Sections(handle(sym))[deref(sym).st_shndx + 1]

function symbol_value(sym::ELFSymbolRef)
    # First, grab the actual `st_value` entry from the SymtabEntry
    value = symbol_value(deref(sym))

    # If the given symbol is defined, and is not reserved
    if !isundef(sym) && !isreserved(sym)
        # Grab the address of the section this symbol refers to
        sh_addr = section_address(symbol_section(sym))

        # What to do next depends on the object kind. Shared libraries and
        # executable's st_value's are virtual addresses
        e_type = header(handle(sym)).e_type

        # If this is a plain relocatable file (`.o` file) just increment by
        # section address
        if sh_addr != 0 && e_type == ET_REL
            value += sh_addr
        end

        # If this is an executable or dynamic library
        if e_type in (ET_EXEC, ET_DYN)
            # Subtract off the virtual address of the first PT_LOAD segment
            P = [s for s in Segments(handle(sym)) if segment_type(s) == PT_LOAD]
            if !isempty(P)
                value -= segment_virtual_address(first(P))
            end
        end
    end

    # Return our ill-gotten goods
    return value
end


text:
module TestDiagrammaticPrograms
using Test

using Catlab, Catlab.Graphs, Catlab.CategoricalAlgebra
using Catlab.Programs.DiagrammaticPrograms
using Catlab.Programs.DiagrammaticPrograms: NamedGraph
using Catlab.Graphs.BasicGraphs: TheoryGraph, TheoryReflexiveGraph
using Catlab.Graphs.BipartiteGraphs: TheoryBipartiteGraph
using Catlab.WiringDiagrams.CPortGraphs: ThCPortGraph

@present TheoryDDS(FreeSchema) begin
  X::Ob
  Φ::Hom(X,X)
end

# Graphs
########

g = @graph begin
  s
  t
  s → t
  s → t
end
@test g == parallel_arrows(NamedGraph{Symbol,Union{Symbol,Nothing}}, 2,
                           V=(vname=[:s,:t],), E=(ename=[nothing,nothing],))

g = @graph NamedGraph{Symbol,Symbol} begin
  x, y
  (f, g): x → y
end
@test g == parallel_arrows(NamedGraph{Symbol,Symbol}, 2,
                           V=(vname=[:x,:y],), E=(ename=[:f,:g],))

tri_parsed = @graph NamedGraph{Symbol,Symbol} begin
  v0, v1, v2
  fst: v0 → v1
  snd: v1 → v2
  comp: v0 → v2
end
tri = @acset NamedGraph{Symbol,Symbol} begin
  V = 3
  E = 3
  src = [1,2,1]
  tgt = [2,3,3]
  vname = [:v0, :v1, :v2]
  ename = [:fst, :snd, :comp]
end
@test tri_parsed == tri

# Categories
############

Δ¹_parsed = @fincat begin
  V, E
  (δ₀, δ₁): V → E
  σ₀: E → V

  σ₀ ∘ δ₀ == id(V)
  σ₀ ∘ δ₁ == id(V)
end
Δ¹_graph = @acset NamedGraph{Symbol,Symbol} begin
  V = 2
  E = 3
  src = [1,1,2]
  tgt = [2,2,1]
  vname = [:V, :E]
  ename = [:δ₀, :δ₁, :σ₀]
end
Δ¹ = FinCat(Δ¹_graph, [ [1,3] => empty(Path, Δ¹_graph, 1),
                        [2,3] => empty(Path, Δ¹_graph, 1) ])
@test Δ¹_parsed == Δ¹

# Functors
##########

# Underlying graph of circular port graph.
F = @finfunctor TheoryGraph ThCPortGraph begin
  V => Box
  E => Wire
  src => src ⨟ box
  tgt => tgt ⨟ box
end
@test F == FinFunctor(Dict(:V => :Box, :E => :Wire),
                      Dict(:src => [:src, :box], :tgt => [:tgt, :box]),
                      TheoryGraph, ThCPortGraph)

# Incomplete definition.
@test_throws ErrorException begin
  @finfunctor TheoryGraph ThCPortGraph begin
    V => Box
    src => src ⨟ box
    tgt => tgt ⨟ box
  end
end

# Failure of functorality.
@test_throws ErrorException begin
  @finfunctor TheoryGraph ThCPortGraph begin
    V => Box
    E => Wire
    src => src
    tgt => tgt
  end
end

# Extra definitions.
@test_throws ErrorException begin
  @finfunctor TheoryGraph TheoryReflexiveGraph begin
    V => Box
    E => Wire
    src => src
    tgt => tgt
    refl => refl
  end
end

# GAT expressions.
F = @finfunctor TheoryDDS TheoryDDS begin
  X => X; Φ => id(X)
end
F′ = @finfunctor TheoryDDS TheoryDDS begin
  X => X; Φ => id{X}
end
@test F == F′

# Diagrams
##########

C = FinCat(TheoryGraph)
F_parsed = @diagram C begin
  v::V
  (e1, e2)::E
  (t: e1 → v)::tgt
  (s: e2 → v)::src
end
J = FinCat(@acset NamedGraph{Symbol,Union{Symbol,Nothing}} begin
  V = 3
  E = 2
  src = [2,3]
  tgt = [1,1]
  vname = [:v, :e1, :e2]
  ename = [:t, :s]
end)
@test dom(F_parsed) == J
F = FinFunctor([:V,:E,:E], [:tgt, :src], J, C)
@test F_parsed == F

F_parsed = @diagram TheoryGraph begin
  v => V
  (e1, e2) => E
  t: e1 → v => tgt
  s: e2 → v => src
end
@test F_parsed == F

F_parsed = @diagram TheoryGraph begin
  v::V
  (e1, e2)::E
  (e1 → v)::tgt
  (e2 → v)::src
end
J_parsed = dom(F_parsed)
@test src(graph(J_parsed)) == src(graph(J))
@test tgt(graph(J_parsed)) == tgt(graph(J))

F_parsed′ = @free_diagram TheoryGraph begin
  v::V
  (e1, e2)::E
  tgt(e1) == v
  v == src(e2)
end
@test F_parsed′ == F_parsed

F = @free_diagram TheoryGraph begin
  (e1, e2)::E
  tgt(e1) == src(e2)
end
@test is_functorial(F)
@test collect_ob(F) == [TheoryGraph[:E], TheoryGraph[:E], TheoryGraph[:V]]
@test collect_hom(F) == [TheoryGraph[:tgt], TheoryGraph[:src]]

F = @diagram TheoryDDS begin
  x::X
  (f: x → x)::(Φ⋅Φ)
end
@test only(collect_ob(F)) == TheoryDDS[:X]
@test only(collect_hom(F)) == compose(TheoryDDS[:Φ], TheoryDDS[:Φ])

# Migrations
############

# Pullback migration
#-------------------

# Graph with reversed edges.
F = @migration TheoryGraph TheoryGraph begin
  V => V
  E => E
  src => tgt
  tgt => src
end
@test F isa DataMigrations.DeltaSchemaMigration
@test F == FinFunctor(Dict(:V => :V, :E => :E),
                      Dict(:src => :tgt, :tgt => :src),
                      TheoryGraph, TheoryGraph)

# Variant where target schema is not given.
F = @migration TheoryGraph begin
  E => E
  V => V
  (src: E → V) => tgt
  (tgt: E → V) => src
end
J = FinCat(parallel_arrows(NamedGraph{Symbol,Union{Symbol,Nothing}}, 2,
                           V=(vname=[:E,:V],), E=(ename=[:src,:tgt],)))
@test F == FinDomFunctor([:E,:V], [:tgt,:src], J, FinCat(TheoryGraph))

# Conjunctive migration
#----------------------

# Graph with edges that are paths of length 2.
F = @migration TheoryGraph TheoryGraph begin
  V => V
  E => @join begin
    v::V
    (e₁, e₂)::E
    (e₁ → v)::tgt
    (e₂ → v)::src
  end
  src => e₁ ⋅ src
  tgt => e₂ ⋅ tgt
end
@test F isa DataMigrations.ConjSchemaMigration
F_E = diagram(ob_map(F, :E))
@test nameof.(collect_ob(F_E)) == [:V, :E, :E]
@test nameof.(collect_hom(F_E)) == [:tgt, :src]
F_tgt = hom_map(F, :tgt)
@test collect_ob(F_tgt) == [(3, TheoryGraph[:tgt])]

# Syntactic variant of above.
F′ = @migration TheoryGraph TheoryGraph begin
  V => V
  E => @join begin
    v::V
    (e₁, e₂)::E
    tgt(e₁) == v
    src(e₂) == v
  end
  src => src(e₁)
  tgt => tgt(e₂)
end
@test F′ == F

# Cartesian product of graph with itself.
F = @migration TheoryGraph TheoryGraph begin
  V => @product (v₁::V; v₂::V)
  E => @product (e₁::E; e₂::E)
  src => (v₁ => e₁⋅src; v₂ => e₂⋅src)
  tgt => (v₁ => e₁⋅tgt; v₂ => e₂⋅tgt)
end
F_V = diagram(ob_map(F, :V))
@test collect_ob(F_V) == fill(TheoryGraph[:V], 2)
F_src = hom_map(F, :src)
@test collect_ob(F_src) == [(1, TheoryGraph[:src]), (2, TheoryGraph[:src])]

# Reflexive graph from graph.
F = @migration TheoryReflexiveGraph TheoryGraph begin
  V => @join begin
    v::V
    ℓ::E
    (s: ℓ → v)::src
    (t: ℓ → v)::tgt
  end
  E => @join begin
    (v₁, v₂)::V
    (ℓ₁, ℓ₂, e)::E
    (s₁: ℓ₁ → v₁)::src
    (t₁: ℓ₁ → v₁)::tgt
    (s₂: ℓ₂ → v₂)::src
    (t₂: ℓ₂ → v₂)::tgt
    (s: e → v₁)::src
    (t: e → v₂)::tgt
  end
  refl => begin
    (v₁, v₂) => v
    (ℓ₁, ℓ₂, e) => ℓ
    (s₁, s₂, s) => s
    (t₁, t₂, t) => t
  end
  src => begin
    v => v₁; ℓ => ℓ₁; s => s₁; t => t₁
  end
  tgt => begin
    v => v₂; ℓ => ℓ₂; s => s₂; t => t₂
  end
end
F_tgt = hom_map(F, :tgt)
@test ob_map(F_tgt, :v) == (2, id(TheoryGraph[:V]))
@test hom_map(F_tgt, :t) |> edges |> only == 4

# Free/initial port graph on a graph.
# This is the left adjoint to the underlying graph functor.
F = @migration TheoryGraph begin
  Box => V
  Wire => E
  InPort => @join begin
    v::V
    e::E
    (t: e → v)::tgt
  end
  OutPort => @join begin
    v::V
    e::E
    (s: e → v)::src
  end
  (in_port_box: InPort → Box) => v
  (out_port_box: OutPort → Box) => v
  (src: Wire → OutPort) => begin
    v => src
  end
  (tgt: Wire → InPort) => begin
    v => tgt
  end
end
F_src = hom_map(F, :src)
@test collect_ob(F_src) == [(1, TheoryGraph[:src]), (1, id(TheoryGraph[:E]))]
@test collect_hom(F_src) == [id(shape(codom(F_src)), 1)]

# Gluing migration
#-----------------

# Coproduct of graph with itself.
F = @migration TheoryGraph TheoryGraph begin
  V => @cases (v₁::V; v₂::V)
  E => @cases (e₁::E; e₂::E)
  src => begin
    e₁ => v₁ ∘ src
    e₂ => v₂ ∘ src
  end
  tgt => begin
    e₁ => v₁ ∘ tgt
    e₂ => v₂ ∘ tgt
  end
end
@test F isa DataMigrations.GlueSchemaMigration
F_V = diagram(ob_map(F, :V))
@test collect_ob(F_V) == fill(TheoryGraph[:V], 2)
F_src = hom_map(F, :src)
@test collect_ob(F_src) == [(1, TheoryGraph[:src]), (2, TheoryGraph[:src])]

# Free reflexive graph on a graph.
F = @migration TheoryReflexiveGraph TheoryGraph begin
  V => V
  E => @cases (e::E; v::V)
  src => (e => src)
  tgt => (e => tgt)
  refl => v
end
F_tgt = hom_map(F, :tgt)
@test collect_ob(F_tgt) == [(1, TheoryGraph[:tgt]), (1, id(TheoryGraph[:V]))]

# Vertices in a graph and their connected components.
F = @migration TheoryGraph begin
  V => V
  Component => @glue begin
    e::E; v::V
    (e → v)::src
    (e → v)::tgt
  end
  (component: V → Component) => v
end
F_C = diagram(ob_map(F, :Component))
@test nameof.(collect_ob(F_C)) == [:E, :V]
@test nameof.(collect_hom(F_C)) == [:src, :tgt]

# Gluc migration
#---------------

# Graph with edges that are paths of length <= 2.
F = @migration TheoryGraph TheoryGraph begin
  V => V
  E => @cases begin
    v => V
    e => E
    path => @join begin
      v::V
      (e₁, e₂)::E
      (e₁ → v)::tgt
      (e₂ → v)::src
    end
  end
  src => begin
    e => src
    path => e₁⋅src
  end
  tgt => begin
    e => tgt
    path => e₂⋅tgt
  end
end
@test ob_map(F, :V) isa DataMigrations.GlucQuery
@test F isa DataMigrations.GlucSchemaMigration
F_src = hom_map(F, :src)
@test collect_ob(shape_map(F_src)) == [1,1,1]
F_src_v, F_src_e, F_src_path = components(diagram_map(F_src))
@test collect_ob(F_src_v) == [(1, id(TheoryGraph[:V]))]
@test collect_ob(F_src_e) == [(1, TheoryGraph[:src])]
@test collect_ob(F_src_path) == [(2, TheoryGraph[:src])]

# Graph with edges that are minimal paths b/w like vertices in bipartite graph.
F = @migration TheoryGraph TheoryBipartiteGraph begin
  V => @cases (v₁::V₁; v₂::V₂)
  E => @cases begin
    e₁ => @join begin
      v₂::V₂; e₁₂::E₁₂; e₂₁::E₂₁
      (e₁₂ → v₂)::tgt₂
      (e₂₁ → v₂)::src₂
    end
    e₂ => @join begin
      v₁::V₁; e₂₁::E₂₁; e₁₂::E₁₂
      (e₂₁ → v₁)::tgt₁
      (e₁₂ → v₁)::src₁
    end
  end
  src => begin
    e₁ => v₁ ∘ (e₁₂ ⋅ src₁)
    e₂ => v₂ ∘ (e₂₁ ⋅ src₂)
  end
  tgt => begin
    e₁ => v₁ ∘ (e₂₁ ⋅ tgt₁)
    e₂ => v₂ ∘ (e₁₂ ⋅ tgt₂)
  end
end
@test ob_map(F, :V) isa DataMigrations.GlucQuery
@test F isa DataMigrations.GlucSchemaMigration
F_src = hom_map(F, :src)
@test collect_ob(shape_map(F_src)) == [1,2]
F_src1, F_src2 = components(diagram_map(F_src))
@test collect_ob(F_src1) == [(2, TheoryBipartiteGraph[:src₁])]
@test collect_ob(F_src2) == [(2, TheoryBipartiteGraph[:src₂])]

# Box product of reflexive graph with itself.
F = @migration TheoryReflexiveGraph TheoryReflexiveGraph begin
  V => @product (v₁::V; v₂::V)
  E => @glue begin
    vv => @product (v₁::V; v₂::V)
    ev => @product (e₁::E; v₂::V)
    ve => @product (v₁::V; e₂::E)
    (refl₁: vv → ev) => (e₁ => v₁⋅refl; v₂ => v₂)
    (refl₂: vv → ve) => (v₁ => v₁; e₂ => v₂⋅refl)
  end
  src => begin
    ev => (v₁ => e₁⋅src; v₂ => v₂)
    ve => (v₁ => v₁; v₂ => e₂⋅src)
  end
  tgt => begin
    ev => (v₁ => e₁⋅tgt; v₂ => v₂)
    ve => (v₁ => v₁; v₂ => e₂⋅tgt)
  end
  refl => vv
end
@test ob_map(F, :V) isa DataMigrations.GlucQuery
@test F isa DataMigrations.GlucSchemaMigration
F_src = hom_map(F, :src)
@test collect_ob(shape_map(F_src)) == [1,1,1]
F_src_vv, F_src_ev, F_src_ve = components(diagram_map(F_src))
@test collect_ob(F_src_ev) == [(1, TheoryReflexiveGraph[:src]),
                               (2, id(TheoryReflexiveGraph[:V]))]

end

max_forks_repo_name: epatters/CompCat


text:
# assumption: already ran the sample code -

# what output are we looking at?
output_index = 4
#xlabel_text = "Clot time CT (min)"
#xlabel_text = "MCF (mm)"
#xlabel_text = "α (deg)"
xlabel_text = "AUC (mm*s)"
#xlabel_text = "Area under FIIa curve (nmol-min/L)"

#μ_bias = 1.2
#μ_bias = 38.32

dm = Normal(mean(S[:, output_index]), std(S[:, output_index]))
dd = Normal(mean(validation_data_array[:, output_index]), std(validation_data_array[:, output_index]))
μ_bias = mean(dd) - mean(dm)
#μ_bias = 0.0 # no mean correction

# sample -
number_of_samples = 10000
dms = rand(dm, number_of_samples) .+ μ_bias
dds = rand(dd, number_of_samples)

# plot -
number_of_bins = round(Int64, 0.02 * number_of_samples)
stephist(dms, bins = number_of_bins, label = "Model", lw = 2, normed = true)
stephist!(dds, bins = number_of_bins, label = "Measured", lw = 2, normed = true)

# labels -
xlabel!(xlabel_text, fontsize = 18)
ylabel!("Frequency of occurance (dimensionless)", fontsize = 18)

# dump fig -
fig_file_name = "Fig-distribution-OI-$(output_index)-N2-TPA-0.pdf"
savefig(joinpath(_PATH_TO_FIGS, fig_file_name))

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
using Borg

function dtlz2(vars, objs, consts)
    nvars = length(vars)
    nobjs = length(objs)

    k = nvars - nobjs + 1
    g = 0.0

    for i=nvars-k:nvars-1
        g += (vars[i+1] - 0.5)^2.0
    end

    for i=0:nobjs-1
        objs[i+1] = 1.0 + g

        for j=0:nobjs-i-1-1
            objs[i+1] *= cos(0.5*pi*vars[j+1])
        end

        if i != 0
            objs[i+1] *= sin(0.5*pi*vars[nobjs-i-1+1])
        end
    end
end

nvars = 11
nobjs = 2

lowerBounds = ones(nvars) * 0.0
upperBounds = ones(nvars) * 1.0
epsilons = ones(nobjs) * 0.01

problem = createProblem(lowerBounds, upperBounds, epsilons, 0, dtlz2)

result = solve(problem, 1000000)


text:
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Architectures: device
using Oceananigans.Operators: ∂xᶠᵃᵃ, ∂yᵃᶠᵃ, Δzᵃᵃᶠ, Δzᵃᵃᶜ
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions
using Oceananigans.Solvers: solve!
using Oceananigans.Fields

using Adapt
using KernelAbstractions: NoneEvent

struct ImplicitFreeSurface{E, G, B, V, R, I, S}
    η :: E
    gravitational_acceleration :: G
    barotropic_volume_flux :: B
    vertically_integrated_lateral_face_areas :: V
    implicit_step_right_hand_side :: R
    implicit_step_solver :: I
    solver_settings :: S
end

# User interface to ImplicitFreeSurface
ImplicitFreeSurface(; gravitational_acceleration=g_Earth, solver_settings...) =
    ImplicitFreeSurface(nothing, gravitational_acceleration, nothing, nothing, nothing, nothing, solver_settings)

Adapt.adapt_structure(to, free_surface::ImplicitFreeSurface) =
    ImplicitFreeSurface(Adapt.adapt(to, free_surface.η), free_surface.gravitational_acceleration,
                        nothing, nothing, nothing, nothing, nothing)

# Internal function for HydrostaticFreeSurfaceModel
function FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, arch, grid)
    η = FreeSurfaceDisplacementField(velocities, arch, grid)
    g = convert(eltype(grid), free_surface.gravitational_acceleration)

    # Initialize barotropic volume fluxes
    barotropic_x_volume_flux = ReducedField(Face, Center, Nothing, arch, grid; dims=3)
    barotropic_y_volume_flux = ReducedField(Center, Face, Nothing, arch, grid; dims=3)
    barotropic_volume_flux = (u=barotropic_x_volume_flux, v=barotropic_y_volume_flux)

    # Initialize vertically integrated lateral face areas
    ∫ᶻ_Axᶠᶜᶜ = ReducedField(Face, Center, Nothing, arch, grid; dims=3)
    ∫ᶻ_Ayᶜᶠᶜ = ReducedField(Center, Face, Nothing, arch, grid; dims=3)

    vertically_integrated_lateral_face_areas = (xᶠᶜᶜ = ∫ᶻ_Axᶠᶜᶜ, yᶜᶠᶜ = ∫ᶻ_Ayᶜᶠᶜ)

    compute_vertically_integrated_lateral_face_areas!(vertically_integrated_lateral_face_areas, grid, arch)

    implicit_step_solver = PreconditionedConjugateGradientSolver(implicit_free_surface_linear_operation!,
                                                                 template_field = η,
                                                                 maximum_iterations = grid.Nx * grid.Ny;
                                                                 free_surface.solver_settings...)

    implicit_step_right_hand_side = ReducedField(Center, Center, Nothing, arch, grid; dims=3)

    return ImplicitFreeSurface(η,
                               g,
                               barotropic_volume_flux,
                               vertically_integrated_lateral_face_areas,
                               implicit_step_right_hand_side,
                               implicit_step_solver,
                               free_surface.solver_settings)
end

@inline explicit_barotropic_pressure_x_gradient(i, j, k, grid, ::ImplicitFreeSurface) = 0
@inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, ::ImplicitFreeSurface) = 0

"""
Implicitly step forward η.
"""
ab2_step_free_surface!(free_surface::ImplicitFreeSurface, model, Δt, χ, velocities_update) =
    implicit_free_surface_step!(free_surface::ImplicitFreeSurface, model, Δt, χ, velocities_update)

function implicit_free_surface_step!(free_surface::ImplicitFreeSurface, model, Δt, χ, velocities_update)

    η = free_surface.η
    g = free_surface.gravitational_acceleration
    rhs = free_surface.implicit_step_right_hand_side
    ∫ᶻ_Q = free_surface.barotropic_volume_flux
    ∫ᶻ_A = free_surface.vertically_integrated_lateral_face_areas

    #=
    # Take an explicit step first to produce an improved initial guess for η for the iterative solver.
    event = explicit_ab2_step_free_surface!(free_surface, model, Δt, χ)
    wait(device(model.architecture), event)
    =#

    fill_halo_regions!(η, model.architecture)

    compute_vertically_integrated_volume_flux!(∫ᶻ_Q, model, velocities_update)

    compute_implicit_free_surface_right_hand_side!(rhs, model, g, Δt, ∫ᶻ_Q, η)

    fill_halo_regions!(rhs, model.architecture)

    # solve!(x, solver, b, args...) solves A*x = b for x.
    solve!(η, free_surface.implicit_step_solver, rhs, ∫ᶻ_A.xᶠᶜᶜ, ∫ᶻ_A.yᶜᶠᶜ, g, Δt)

    return NoneEvent()
end

function compute_implicit_free_surface_right_hand_side!(rhs, model, g, Δt, ∫ᶻ_Q, η)

    event = launch!(model.architecture,
                    model.grid,
                    :xy,
                    implicit_free_surface_right_hand_side!,
                    rhs,
                    model.grid,
                    g,
                    Δt,
                    ∫ᶻ_Q,
                    η,
                    dependencies=Event(device(model.architecture)))

    wait(device(model.architecture), event)

    fill_halo_regions!(rhs, model.architecture)

    return nothing
end

""" Compute the divergence of fluxes Qu and Qv. """
@inline flux_div_xyᶜᶜᵃ(i, j, k, grid, Qu, Qv) = δxᶜᵃᵃ(i, j, k, grid, Qu) + δyᵃᶜᵃ(i, j, k, grid, Qv)

@kernel function implicit_free_surface_right_hand_side!(rhs, grid, g, Δt, ∫ᶻ_Q, η)
    i, j = @index(Global, NTuple)
    @inbounds rhs[i, j, 1] = - Azᶜᶜᵃ(i, j, 1, grid) * η[i, j, 1] / (g * Δt^2) +
                               flux_div_xyᶜᶜᵃ(i, j, 1, grid, ∫ᶻ_Q.u, ∫ᶻ_Q.v) / (g * Δt)
end


text:
function RangeTroposphericDelayModifier(arg0::DiscreteTroposphericModel)
    return RangeTroposphericDelayModifier((DiscreteTroposphericModel,), arg0)
end

function get_parameters_drivers(obj::RangeTroposphericDelayModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::RangeTroposphericDelayModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
using GeometryTypes
function f(G, x)
    fx = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

    if !(G == nothing)
        G1 = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
        G2 = 200.0 * (x[2] - x[1]^2)
        gx = Point(G1,G2)

        return fx, gx
    else
        return fx
    end
end

minimize(f, Point(3.0,3.0), NearlyNewton.GradientDescent(InverseApprox()), I, NearlyNewton.OptOptions())
minimize(f, Point(3.0,3.0), NearlyNewton.BFGS(InverseApprox()), I, NearlyNewton.OptOptions())
minimize(f, Point(3.0,3.0), NearlyNewton.DFP(InverseApprox()), I, NearlyNewton.OptOptions())
minimize(f, Point(3.0,3.0), NearlyNewton.SR1(InverseApprox()), I, NearlyNewton.OptOptions())

minimize(f, Point(3.0,3.0), NearlyNewton.GradientDescent(DirectApprox()), I, NearlyNewton.OptOptions())
minimize(f, Point(3.0,3.0), NearlyNewton.BFGS(DirectApprox()), I, NearlyNewton.OptOptions())
minimize(f, Point(3.0,3.0), NearlyNewton.DFP(DirectApprox()), I, NearlyNewton.OptOptions())
minimize(f, Point(3.0,3.0), NearlyNewton.SR1(DirectApprox()), I, NearlyNewton.OptOptions())


text:
@testset "BernsteinPolynomial" begin
    @testset "evaluation" begin
        b = BernsteinPolynomial(1, 2, 3, 4)
        @test b(0) === 1
        @test b(1) === 4
    end

    @testset "utility" begin
        b = BernsteinPolynomial(1, 2, 3, 4)
        @test zero(b) === zero(typeof(b)) === BernsteinPolynomial(0, 0, 0, 0)
        @test constant(b) == 1
    end

    @testset "derivative" begin
        b = BernsteinPolynomial(1, 2, 3, 4)
        b′ = derivative(b)
        for t in range(0., 1., length=10)
            @test ForwardDiff.derivative(b, t) ≈ b′(t) atol=1e-10
        end
    end

    @testset "exponential integral" begin
        b = BernsteinPolynomial(1, 2, 3, 4)
        c = 5
        for t in range(0., 1., length=10)
            @test ForwardDiff.derivative(t -> exponential_integral(b, c, t), t) ≈ b(t) * exp(c * t) atol=1e-10
        end
        @test exponential_integral(b, c) ≈ exponential_integral(b, c, 1.0) atol=1e-10
    end

    @testset "arithmetic" begin
        b = BernsteinPolynomial(1, 2, 3, 4)
        c = 5
        for t in range(0., 1., length=10)
            @test (b + c)(t) ≈ b(t) + c
            @test (b - c)(t) ≈ b(t) - c
            @test (c + b)(t) ≈ c + b(t)
            @test (c - b)(t) ≈ c - b(t)
        end
        b1 = BernsteinPolynomial(1, 2, 3, 4)
        b2 = BernsteinPolynomial(2, 3, 4, 5)
        for t in range(0., 1., length=10)
            @test (b1 + b2)(t) ≈ b1(t) + b2(t)
            @test (b1 - b2)(t) ≈ b1(t) - b2(t)
        end
    end

    @testset "scaling" begin
        b1 = BernsteinPolynomial(2, 3, 4)
        @test b1 * 4 === 4 * b1 === BernsteinPolynomial(8, 12, 16)
        @test b1 / 2 === BernsteinPolynomial(2 / 2, 3 / 2, 4 / 2)
        for t in range(0., 1., length=10)
            @test (b1 * 4)(t) ≈ b1(t) * 4
            @test (b1 / 2)(t) ≈ b1(t) / 2
        end
        b2 = BernsteinPolynomial(ntuple(i -> 21 - i, Val(20)))
        allocs = let b2 = b2
            @allocated b2 * 4
        end
        @test allocs == 0
    end

    @testset "BernsteinPolynomial to Polynomial" begin
        b = BernsteinPolynomial(5, 4, 3, 2, 1, 10)
        p = Polynomial(b)
        for t in range(0., 1., length=10)
            @test b(t) ≈ p(t)
        end
        allocs = let b = b
            @allocated Polynomial(b)
        end
        @test allocs == 0
    end
end


text:
# This is a Cooley-Tukey FFT algorithm for any number type.
function fft_pow2(x::Vector{F}) where F
    n = length(x)
    T = mapreduce(eltype,promote_type,x)
    @assert ispow2(n)
    if n==1
        return x
    elseif n==2
        return F[x[1]+x[2];x[1]-x[2]]
    end
    even,odd = fft_pow2(x[1:2:end-1]),fft_pow2(x[2:2:end])
    twiddle = exp(-2im*convert(T,π)/n*collect(0:n-1))
    half1 = even + odd.*twiddle[1:div(n,2)]
    half2 = even + odd.*twiddle[div(n,2)+1:n]
    return vcat(half1,half2)
end
ifft_pow2(x::Vector) = conj(fft_pow2(conj(x)))/length(x)


text:
"""
    (fx,ft)=fft3d(z)

Compute spatial (`fx`) and temporal (`ft`) Fourier spectra of field `z`.
"""
        function fft3d(z)

@static        if VERSION < v"0.7.0-DEV.2005"
           function sum(a;dims=1)
              return Base.sum(a,dims);
           end
        end

#       function [fx,fy,ft]=fft3d(ss);
        nss=size(z); ns=nss[1];
        if(length(nss)>=3)
            nt=nss[3]
	else
	    nt=1
        end

        ns2=div(ns,2);
        nt2=div(nt,2);
        if((nt>1)&&(mod(nt,2)!=0))
            z=z[:,:,1:(nt-1)];
            nt=nt-1;
        end

        zf=abs.(fft(z)/(ns*ns*nt)).^2;
        zf0=zf;
        zf[ns2+1,:,:]=zf0[ns2+1,:,:]/2;
        zf[:,ns2+1,:]=zf0[:,ns2+1,:]/2;
        if(nt>1)
          zf[:,:,nt2+1]=zf0[:,:,nt2+1]/2;
        end

@compat        fs=reshape(sum(zf,dims=3),ns,ns);
@compat        ft=reshape(sum(sum(zf,dims=1),dims=2),nt,1);

        fs=fftshift(fs) ;
        fs=fs/nt;

        nn=zeros(ns,1);
        fx=zeros(ns,1);
        for j=-ns2:ns2-1
           for i=-ns2:ns2-1
              k2=sqrt(i*i+j*j);
              ik2=floor(Integer,k2+1.5);
              if(ik2>1)
                fx[ik2]=fx[ik2]+fs[i+ns2+1,j+ns2+1];
                nn[ik2]=nn[ik2]+1;
              end
           end
        end
    if(nt>1)
        ft=ft[2:nt2+1];
    end
   fx=fx[2:ns2+1]./nn[2:ns2+1];
   return (fx,ft,fs)
   end


text:
mutable struct OptimizationProblem{N, Ny, Nx, Nθ, M}
    ms::MultipleShooting{Vector{Float64}, N, Ny, Nx, Nθ, M}
    # Parameter vector
    θ_ext::Vector{Float64}
    θ_aux::Vector{Float64}
    x0_list_aux::Vector{Vector{Float64}}
    # Gradient
    grad::Vector{Float64}
    gradθ_aux::Vector{Float64}
    gradx0_aux::Vector{Vector{Float64}}
    # Multiplicative vector p
    p_θ::Vector{Float64}
    p_x0::Vector{Vector{Float64}}
    # Hessian
    hessp::Vector{Float64}
    hesspθ_aux::Vector{Float64}
    hesspx0_aux::Vector{Vector{Float64}}
    # Constraint
    constr::Vector{Float64}
    constr_aux::Vector{Vector{Float64}}
    # Jacobian
    jacobian::SparseMatrixCOO{Float64,Int64}
    jacobian_θ::Vector{Matrix{Float64}}
    jacobian_x0::Vector{Matrix{Float64}}
end

function OptimizationProblem(f::Function, g::Function,
                             x0_list::Vector{Vector{Float64}},
                             y::Vector{Vector{Float64}},
                             k0_list::Vector{Int}, θ;
                             loss=L2DistLoss())
   ms = MultipleShooting(f, g, x0_list, y,
                         k0_list, θ, loss)
   # Problem dimensions
   M = length(x0_list)
   N = length(y)
   Ny = length(y[1])
   Nx = length(x0_list[1])
   Nθ = length(θ)
   # Parameter Vector
   θ_ext = zeros(Nθ + M*Nx)
   θ_aux = zeros(Nθ)
   x0_list_aux = Vector{Vector{Float64}}(M)
   for i = 1:M
       x0_list_aux[i] = zeros(Nx)
   end
   build_extended_vector!(θ_ext, θ, x0_list, Nθ, M, Nx)
   # Gradient
   grad = zeros(Nθ + M*Nx)
   gradθ_aux = zeros(Nθ)
   gradx0_aux = Vector{Vector{Float64}}(M)
   for i = 1:M
       gradx0_aux[i] = zeros(Nx)
   end
   # Multiplicative_vector_p
   p_θ = zeros(Nθ)
   p_x0 = Vector{Vector{Float64}}(M)
   for i = 1:M
       p_x0[i] = zeros(Nx)
   end
   # Hessian
   hessp = zeros(Nθ + M*Nx)
   hesspθ_aux = zeros(Nθ)
   hesspx0_aux = Vector{Vector{Float64}}(M)
   for i = 1:M
       hesspx0_aux[i] = zeros(Nx)
   end
   # Constraint
   constr = Vector{Float64}((M-1)*Nx)
   constr_aux = Vector{Vector{Float64}}(M-1)
   for i = 1:M-1
       constr_aux[i] = zeros(Nx)
   end
   # Jacobian
   jacobian_θ = Vector{Matrix{Float64}}(M-1)
   for i = 1:M-1
       jacobian_θ[i] = ones(Nx, Nθ)
   end
   jacobian_x0 = Vector{Matrix{Float64}}(M-1)
   for i = 1:M-1
       jacobian_x0[i] = ones(Nx, Nx)
   end
   i, j = jacobian_indices(Nθ, M, Nx)
   nnz = length(i)
   v = Vector{Float64}(nnz)
   jacobian = SparseMatrixCOO{Float64,Int64}(Nx*(M-1), Nθ + M*Nx, i, j, v)
   jacobian_from_matrices!(jacobian, jacobian_θ, jacobian_x0,
                           Nθ, M, Nx; reset_constants=true)
   # Call constructor
   OptimizationProblem{N, Ny, Nx, Nθ, M}(
      ms, θ_ext, θ_aux, x0_list_aux, grad, gradθ_aux, gradx0_aux, p_θ, p_x0,
      hessp, hesspθ_aux, hesspx0_aux, constr, constr_aux, jacobian,
      jacobian_θ, jacobian_x0)
end

function cost_function{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M},
        θ_ext::Vector{Float64})
    update_simulation!(opt, θ_ext)
    return cost_function(opt.ms)
end

function gradient{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M},
        θ_ext::Vector{Float64})
    update_simulation!(opt, θ_ext)
    gradient!(opt.gradθ_aux, opt.ms, "θ")
    gradient!(opt.gradx0_aux, opt.ms, "x0")
    return build_extended_vector!(opt.grad, opt.gradθ_aux, opt.gradx0_aux,
                                  Nθ, M, Nx)
end

function hessian{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M},
        θ_ext::Vector{Float64})
    update_simulation!(opt, θ_ext)
    function hessp(p)
        read_extended_vector!(opt.p_θ, opt.p_x0, p, Nθ, M, Nx)
        hessian_approx!(opt.hesspθ_aux, opt.ms, opt.p_θ, "θ")
        hessian_approx!(opt.hesspx0_aux, opt.ms, opt.p_x0, "x0")
        return build_extended_vector!(opt.hessp, opt.hesspθ_aux, opt.hesspx0_aux,
                                      Nθ, M, Nx)
    end
    return LinearOperator(Nθ+M*Nx, Nθ+M*Nx, hessp)
end

function constraint{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M},
        θ_ext::Vector{Float64})
    update_simulation!(opt, θ_ext)
    constr!(opt.constr_aux, opt.ms)
    k = 1
    for i = 1:M-1
        for j = 1:Nx
            opt.constr[k] = opt.constr_aux[i][j]
            k +=1
        end
    end
    return opt.constr
end

function jacobian{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M},
        θ_ext::Vector{Float64})
    update_simulation!(opt, θ_ext)
    constr_jac!(opt.jacobian_θ, opt.ms, "θ")
    constr_jac!(opt.jacobian_x0, opt.ms, "x0")
    return jacobian_from_matrices!(opt.jacobian, opt.jacobian_θ, opt.jacobian_x0,
                                   Nθ, M, Nx; reset_constants=false)
end

function update_simulation!{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M}, θ_ext::Vector{Float64})
    if θ_ext != opt.θ_ext
        copy!(opt.θ_ext, θ_ext)
        read_extended_vector!(opt.θ_aux, opt.x0_list_aux, opt.θ_ext, Nθ, M, Nx)
        copy!(opt.θ_ext, θ_ext)
        new_simulation!(opt.ms, opt.x0_list_aux, opt.θ_aux)
        return true
    else
        return false
    end
end

function build_extended_vector!(θ_extended, θ, x0_list, Nθ, M, Nx)
    k = 1
    for i = 1:Nθ
        θ_extended[k] = θ[i]
        k += 1
    end
    for i = 1:M
        for j = 1:Nx
            θ_extended[k] = x0_list[i][j]
            k += 1
        end
    end
    return θ_extended
end

function read_extended_vector!(θ, x0_list, θ_extended, Nθ, M, Nx)
    k = 1
    for i = 1:Nθ
        θ[i] = θ_extended[k]
        k += 1
    end
    for i = 1:M
        for j = 1:Nx
            x0_list[i][j] = θ_extended[k]
            k += 1
        end
    end
    return θ, x0_list
end

function jacobian_indices(Nθ, M, Nx)
    # Number of nonzero elements
    nnz = (M-1)*Nx*(Nθ+Nx+1)
    # Initialize indeces vectors
    i = Vector{Int64}(nnz)
    j = Vector{Int64}(nnz)
    # Counter
    k = 1
    # Include elements from ``jacobian_θ`` matrices
    for rows_offset = 0:Nx:(M-2)*Nx
        for cols = 1:Nθ
            for rows = 1:Nx
                i[k] = rows + rows_offset
                j[k] = cols
                k += 1
            end
        end
    end
    # Include elements from ``jacobian_x0`` matrices
    cols_offset = Nθ
    for offset = 0:Nx:(M-2)*Nx
        for cols = 1:Nx
            for rows = 1:Nx
                i[k] = rows + offset
                j[k] = cols + offset + cols_offset
                k += 1
            end
        end
    end
    # Include indices relative to -eye(Nx) matrices
    cols_offset = Nθ+Nx
    for offset = 0:Nx:(M-2)*Nx
        for rows = 1:Nx
            i[k] = rows + offset
            j[k] = rows + offset + cols_offset
            k += 1
        end
    end
    return i, j
end

function jacobian_from_matrices!(jacobian::SparseMatrixCOO,
                                 jacobian_θ, jacobian_x0, Nθ, M, Nx;
                                 reset_constants=false)
    # Counter
    k = 1
    # Include elements from ``jacobian_θ`` matrices
    for matrix_number = 1:(M-1)
        for cols = 1:Nθ
            for rows = 1:Nx
                jacobian.v[k] = jacobian_θ[matrix_number][rows, cols]
                k += 1
            end
        end
    end
    # Include elements from ``jacobian_x0`` matrices
    for matrix_number = 1:(M-1)
        for cols = 1:Nx
            for rows = 1:Nx
                jacobian.v[k] = jacobian_x0[matrix_number][rows, cols]
                k += 1
            end
        end
    end
    # Include indices relative to -eye(Nx) matrices
    if reset_constants
        for matrix_number = 1:(M-1)
            for rows = 1:Nx
                jacobian.v[k] = -1
                k += 1
            end
        end
    end
    return jacobian
end

function solve{N, Ny, Nx, Nθ, M}(
        opt::OptimizationProblem{N, Ny, Nx, Nθ, M};
        callback=(x, state)->false, options=Dict(),
        kwargs...)
    fun(θ_ext) = copy(cost_function(opt, θ_ext))
    grad(θ_ext) = copy(gradient(opt, θ_ext))
    hess(θ_ext) = to_python(hessian(opt, θ_ext))
    constr(θ_ext) = copy(constraint(opt, θ_ext))
    jac(θ_ext) = to_python(jacobian(opt, θ_ext))
    nlconstr = scipy_opt["NonlinearConstraint"](constr, 0, 0,
                                                jac, "2-point")
    θ_ext0 = copy(opt.θ_ext)

    if M >= 2
        scipy_opt["minimize"](fun, θ_ext0, jac=grad, hess=hess,
                              options=options, callback=callback,
                              method="trust-constr",
                              constraints=nlconstr;
                              kwargs...)
    else
        scipy_opt["minimize"](fun, θ_ext0, jac=grad, hess=hess,
                              options=options, callback=callback,
                              method="trust-constr";
                              kwargs...)
    end
end


text:
## TODO: Documentation

include("prebuilt_arrays.jl")
include("SNAPParams.jl")
include("runtime_arrays.jl")


text:
using TimeSync
using Test

@testset "TimeSync.jl" begin
    # Write your tests here.
end


text:
"""
C = α*A + β*C

The bias tensor A must match the corresponding dimension of the destination tensor
C or must be equal to 1.
"""
function add!(α, A::CuArray{T}, β, C::CuArray{T}) where T
    h = gethandle()
    adesc = TensorDesc(A, 4)
    cdesc = TensorDesc(C, 4)
    @cudnn(:cudnnAddTensor,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr),
        h, T[α], adesc, A, T[β], cdesc, C)
    C
end


text:
include("./Legendre_polynomials.jl")

function Coupling_coefficients(K::Array{Tuple{Float64,Float64,Float64},1}, lmax::Int64; G::Float64 = 1.0)
    #### Returns an array of size n_pop * n_pop * lmax_half. Each cell holds one (even) coupling coefficient.    ####
    #### The difference with Coupling_coefficients_v2 is that here we take K the list of all orbital populations ####
    #### And return the list of all coupling coefficients. _v2 treats only 2 populations.                        ####

    # getting all the constants we will need
    n_pop = length(K)
    lmax_half = div(lmax, 2)
    n_bin = 10^5 # number of Riemann nodes we will use

    # preallocating all the arrays we will need
    J = zeros(n_pop, n_pop, lmax_half) # will store all the coupling coefficients
    P_l = zeros(lmax_half) # these will be the legendre polynomials evaluted in 0.0
    Coef_l = zeros(lmax_half) # these will be the multiplicative coefficients in front of the integrals (modified at each iteration)
    I_l = zeros(lmax_half) # these will store the integrals
    R₁ = zeros(n_bin) # will store the list of radii at Riemann nodes for the first annuli
    R₂ = zeros(n_bin) # will store the list of radii at Riemann nodes for the second annuli

    # Computing the prefactors
    Legendre_polynomials!(0.0, lmax, Pₗ = P_l) # P_l now holds the list of even Legendre polynomials at point 0.0

    # Computing the Riemann nodes
    cos_E = zeros(n_bin)
    for i = 1:n_bin
        cos_E[i] = cos(π * (2 * i - 1) / (2 * n_bin)) # These are the nodes of the Riemann sums (mid-point rule for faster convergence)
    end

    # Starting the integrals
    for star_1 = 1:n_pop
        m₁, a₁, e₁ = K[star_1]

        # getting the radii
        for i = 1:n_bin
            R₁[i] = a₁ * (1.0 - e₁ * cos_E[i])
        end

        for star_2 = 1:star_1 # We use the symmetry Jₗ[K₁, K₂] = Jₗ[K₂, K₁]
            m₂, a₂, e₂ = K[star_2]

            # getting the radii
            for i = 1:n_bin
                R₂[i] = a₂ * (1.0 - e₂ * cos_E[i])
            end

            # resetting the array of integrals
            for l_half = 1:lmax_half
                I_l[l_half] = 0.0
            end

            # Computing the prefactors
            for l_half = 1:lmax_half
                l = 2 * l_half
                Coef_l[l_half] = G * 4 * π * m₁ * m₂ / (2 * l + 1)
            end

            # starting the integrals -- We use the multipole algorithm
            for l_half = 1:lmax_half

                # let's start with the pⱼ
                wⱼ = 0
                pⱼ = 0.0
                for j = 1:n_bin
                    # first part of the recurrence : multiplication
                    if j > 1
                        pⱼ = pⱼ * (R₂[j - 1] / R₂[j])^((2 * l_half) + 1)
                    end

                    # second part of the recurrence : addition of δp
                    i = wⱼ
                    while ((i + 1 <= n_bin) && (R₁[i + 1] <= R₂[j]))
                        i += 1
                        μ₁ᵢ = R₁[i] / a₁
                        pⱼ += μ₁ᵢ * (R₁[i] / R₂[j])^(2 * l_half) / R₂[j]
                    end
                    wⱼ = i

                    μ₂ⱼ = R₂[j] / a₂
                    I_l[l_half] += μ₂ⱼ * pⱼ
                end

                # now continue with the qⱼ
                wⱼ = n_bin + 1
                qⱼ = 0.0
                for j = n_bin:-1:1
                    # first part of the recurrence : multiplication
                    if j < n_bin
                        qⱼ = qⱼ * (R₂[j] / R₂[j + 1])^(2 * l_half)
                    end

                    # second part of the recurrence : addition of δq
                    i = wⱼ
                    while ((i - 1 >= 1) && (R₁[i - 1] > R₂[j]))
                        i -= 1
                        μ₁ᵢ = R₁[i] / a₁
                        qⱼ += μ₁ᵢ * (R₂[j] / R₁[i])^(2 * l_half) / R₁[i]
                    end
                    wⱼ = i

                    μ₂ⱼ = R₂[j] / a₂
                    I_l[l_half] += μ₂ⱼ * qⱼ
                end
            end

            # Let's store the results of the integral in J
            for l_half = 1:lmax_half
                J_K₁_K₂_l = Coef_l[l_half] * I_l[l_half] * (P_l[l_half]^(2)) / (n_bin^2) # Jₗ[K₁, K₂] = G m₁ m₂ 4 π Pₗ(0)² / (2 * l + 1) * int et int = double Riemann sum / n_bin^2
                J[star_1, star_2, l_half] = J_K₁_K₂_l # We use the symmetry once again there
                J[star_2, star_1, l_half] = J_K₁_K₂_l
            end
        end
    end

    return(J)
end

function Coupling_coefficients_v2!(K₁::Tuple{Float64,Float64,Float64}, K₂::Tuple{Float64,Float64,Float64}, lmax::Int64;
                    G::Float64 = 1.0,
                    Pₗ = nothing, Cₗ = nothing, R₁ = nothing, R₂ = nothing, Iₗ = nothing, Jₗ = nothing)
    #### returns an array of size lmax_half. Each cell contains one (even) coupling coupling coefficient between K₁ and K₂. ####
    #### The difference with Coupling_coefficients is that here we take 2 populations K₁ and K₂ and return the coupling     ####
    #### coefficients between them, whereas Coupling_coefficients works with any number n_pop of populations.               ####

    # getting the constants we will need
    lmax_half = div(lmax, 2)
    n_bin = 10^3 # number of Riemann nodes we will use

    # let's preallocate the arrays that have not been given as argument
    if isnothing(Pₗ)
        Pₗ = zeros(lmax_half) # Will store the (even) Legendre polynomials evaluated in 0.0
    end
    if isnothing(Cₗ)
        Cₗ = zeros(lmax_half) # Will store the prefactors G m₁ m₂ 4 π / (2l + 1)
    end
    if isnothing(R₁)
        R₁ = zeros(n_bin) # Will store the radii at Riemann nodes for the first annuli
    end
    if isnothing(R₂)
        R₂ = zeros(n_bin) # Will store the radii at Riemann nodes for the second annuli
    end
    if isnothing(Iₗ)
        Iₗ = zeros(lmax_half) # Will store the Riemann integrals
    end
    if isnothing(Jₗ)
        Jₗ = zeros(lmax_half) # Will store the coupling coefficients
    end

    # let's compute the prefactors
    m₁, a₁, e₁ = K₁
    m₂, a₂, e₂ = K₂
    for l_half = 1:lmax_half
        Cₗ[l_half] = G * 4 * π * m₁ * m₂ / (2 * (2 * l_half) + 1)
    end

    Legendre_polynomials!(0.0, lmax, Pₗ = Pₗ) # Pₗ now holds the even Legendre polynomials at point 0.0

    # let's compute the radii at the nodes of the Riemann sum
    for j = 1:n_bin
        Ecc = π * (2 * j - 1) / (2 * n_bin)
        R₁[j] = a₁ * (1.0 - e₁ * cos(Ecc))
        R₂[j] = a₂ * (1.0 - e₂ * cos(Ecc))
    end

    # let's reset the array containing the integrals
    for l_half = 1:lmax_half
        Iₗ[l_half] = 0.0
    end

    # Now let's compute the integrals -- We use the multipole algorithm
    for l_half = 1:lmax_half
        # let's start with the pⱼ
        wⱼ = 0
        pⱼ = 0.0
        for j = 1:n_bin
            # first part of the recurrence : multiplication
            if j > 1
                pⱼ = pⱼ * (R₂[j - 1] / R₂[j])^((2 * l_half) + 1)
            end

            # second part of the recurrence : addition of δp
            i = wⱼ
            while ((i + 1 <= n_bin) && (R₁[i + 1] <= R₂[j]))
                i += 1
                μ₁ᵢ = R₁[i] / a₁
                pⱼ += μ₁ᵢ * (R₁[i] / R₂[j])^(2 * l_half) / R₂[j]
            end
            wⱼ = i

            μ₂ⱼ = R₂[j] / a₂
            Iₗ[l_half] += μ₂ⱼ * pⱼ / (n_bin^2)
        end

        # now continue with the qⱼ
        wⱼ = n_bin + 1
        qⱼ = 0.0
        for j = n_bin:-1:1
            # first part of the recurrence : multiplication
            if j < n_bin
                qⱼ = qⱼ * (R₂[j] / R₂[j + 1])^(2 * l_half)
            end

            # second part of the recurrence : addition of δq
            i = wⱼ
            while ((i - 1 >= 1) && (R₁[i - 1] > R₂[j]))
                i -= 1
                μ₁ᵢ = R₁[i] / a₁
                qⱼ += μ₁ᵢ * (R₂[j] / R₁[i])^(2 * l_half) / R₁[i]
            end
            wⱼ = i

            μ₂ⱼ = R₂[j] / a₂
            Iₗ[l_half] += μ₂ⱼ * qⱼ / (n_bin^2)
        end
    end

    # Let's save the integrals in J
    for l_half = 1:lmax_half
        Jₗ[l_half] = Cₗ[l_half] * Iₗ[l_half] * (Pₗ[l_half]^(2)) # Jₗ[K₁, K₂] = G m₁ m₂ 4 π Pₗ(0)² / (2 l + 1) * int
    end
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe