text:
using LinearAlgebra
using LoopVectorization

function RCCSDpT(Alg::ijk)
    val = Options.get("return_ints")
    Options.set("return_ints", true)
    ccsd, moints = RCCSD()
    Options.set("return_ints", val)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(ccsd::RCCSD, moints::IntegralHelper{T,E,O}, Alg::ijk) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    output("\n   • Perturbative Triples Started\n")

    T1 = permutedims(ccsd.T1, (2,1))
    T2 = permutedims(ccsd.T2, (4,3,2,1))

    # Invert order, but retains chemist's notation
    Vvvvo = permutedims(moints["OVVV"], (4,3,2,1))

    # Switch to physicist notation for better memory layout
    Vvooo = permutedims(moints["OOOV"], (4,1,2,3))
    Vvvoo = permutedims(moints["OVOV"], (2,4,1,3))

    v,o = size(T1)

    fo = moints["Fii"]
    fv = moints["Faa"]

    # Pre-allocate Intermediate arrays
    Ws  = [Array{T}(undef, v,v,v) for _ = 1:Threads.nthreads()]
    Vs  = [Array{T}(undef, v,v,v) for _ = 1:Threads.nthreads()]
    Evals = zeros(T, Threads.nthreads())

    output("Computing energy contribution from occupied orbitals:")
    BLAS_THREADS = BLAS.get_num_threads()
    BLAS.set_num_threads(1)

    t = @elapsed begin
    Threads.@threads for i in 1:o
    @inbounds begin
        id = Threads.threadid()
        W = Ws[id]
        V = Vs[id]

        # T1 views → V
        @views T1_i   = T1[:,i]

        # T2 views → W
        @views T2_i   = T2[:,:,:,i]

        @views Vvvvo_i = Vvvvo[:,:,:,i]

        for j in 1:i
            Dij = fo[i] + fo[j]

            # T1 views → V
            @views T1_j  = T1[:,j]

            # T2 views → W
            @views T2_j  = T2[:,:,:,j]
            ## Twice
            @views T2_ji = T2[:,:,j,i]
            @views T2_ij = T2[:,:,i,j]

            @views Vvvvo_j    = Vvvvo[:,:,:,j]

            @views Vvvoo_ij = Vvvoo[:,:,i,j]

            @views Vvooo_ij = Vvooo[:,:,i,j]
            @views Vvooo_ji = Vvooo[:,:,j,i]

            δij = i == j
            for k in 1:j
                Dijk = Dij + fo[k]

                # T1 views → V
                @views T1_k = T1[:,k]

                # T2 views → W
                @views T2_k = T2[:,:,:,k]
                ## Twice
                @views T2_jk = T2[:,:,j,k]
                @views T2_ik = T2[:,:,i,k]
                @views T2_ki = T2[:,:,k,i]
                @views T2_kj = T2[:,:,k,j]

                @views Vvvvo_k  = Vvvvo[:,:,:,k]

                @views Vvooo_jk = Vvooo[:,:,j,k]
                @views Vvooo_kj = Vvooo[:,:,k,j]
                @views Vvooo_ki = Vvooo[:,:,k,i]
                @views Vvooo_ik = Vvooo[:,:,i,k]

                @views Vvvoo_jk = Vvvoo[:,:,j,k]
                @views Vvvoo_ik = Vvvoo[:,:,i,k]

                # Build intermediates W and V for the set i j k
                @tensor begin
                    W[a,b,c] =  Vvvvo_i[d,b,a]*T2_jk[d,c] - Vvooo_ij[b,l]*T2_k[a,c,l]
                    W[a,b,c] += Vvvvo_j[d,a,b]*T2_ik[d,c] - Vvooo_ji[a,l]*T2_k[b,c,l]
                    W[a,b,c] += Vvvvo_k[d,a,c]*T2_ij[d,b] - Vvooo_ki[a,l]*T2_j[c,b,l]
                    W[a,b,c] += Vvvvo_k[d,b,c]*T2_ij[a,d] - Vvooo_kj[b,l]*T2_i[c,a,l]
                    W[a,b,c] += Vvvvo_i[d,c,a]*T2_jk[b,d] - Vvooo_ik[c,l]*T2_j[a,b,l]
                    W[a,b,c] += Vvvvo_j[d,c,b]*T2_ik[a,d] - Vvooo_jk[c,l]*T2_i[b,a,l]

                    V[a,b,c]  = W[a,b,c] + T1_i[a]*Vvvoo_jk[b,c] + Vvvoo_ik[a,c]*T1_j[b] + Vvvoo_ij[a,b]*T1_k[c]
                end

                # Compute Energy contribution
                δjk = j === k
                for a in 1:v
                    Dijka = Dijk - fv[a]
                    for b in 1:a
                        Dijkab = Dijka - fv[b]
                        δab = a === b
                        @fastmath for c in 1:b
                            Dd = Dijkab - fv[c]
                            δbc = b === c
                            X = (W[a,b,c]*V[a,b,c] + W[a,c,b]*V[a,c,b] + W[b,a,c]*V[b,a,c] + W[b,c,a]*V[b,c,a] + W[c,a,b]*V[c,a,b] + W[c,b,a]*V[c,b,a])
                            Y = (V[a,b,c] + V[b,c,a] + V[c,a,b])
                            Z = (V[a,c,b] + V[b,a,c] + V[c,b,a])
                            Ef = (Y - 2*Z)*(W[a,b,c] + W[b,c,a] + W[c,a,b]) + (Z - 2*Y)*(W[a,c,b]+W[b,a,c]+W[c,b,a]) + 3*X
                            Evals[id] += Ef*(2 - δij - δjk) / (Dd*(1 + δab + δbc))
                        end
                    end
                end
            end
        end
        output("  Orbital {} ✔", i)
    end # Threads (i loop)
    end # inbounds
    end # time

    BLAS.set_num_threads(BLAS_THREADS)
    Et = sum(Evals)
    output("Finished in {:5.5f} s", t)
    output("Final (T) contribution: {:15.10f}", Et)
    output("CCSD(T) energy:         {:15.10f}", Et+ccsd.energy)
    return RCCSDpT{T}(ccsd, Et+ccsd.energy, Et)
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# This file contains configuration-related functions and types

# default term file string
const tmpprefix = randstring(8)

# Default font, size and background for each supported terminal
const TerminalDefaults = Dict("wxt" => Dict(:font       => "Sans,10",
                                            :size       => "640,384",
                                            :background => "white"),
                              "qt" => Dict(:font       => "Sans,9",
                                           :size       => "640,480",
                                           :background => ""),
                              "x11" => Dict(:font       => "",
                                            :size       => "640,480",
                                            :background => ""),
                              "aqua" => Dict(:font       => "Times-Roman,14",
                                             :size       => "846,594",
                                             :background => ""),
                              "dumb" => Dict(:font       => "",
                                             :size       => "79,24",
                                             :background => ""),
                              "sixelgd" => Dict(:font       => "Sans,12",
                                                :size       => "640,480",
                                                :background => "white"),
                              "svg" => Dict(:font       => "Arial,12",
                                            :size       => "640,384",
                                            :background => "white"),
                              "gif" => Dict(:font       => "Sans,12",
                                            :size       => "640,480",
                                            :background => "white"),
                              "pngcairo" => Dict(:font       => "Sans,12",
                                                 :size       => "640,480",
                                                 :background => "white"),
                              "pdfcairo" => Dict(:font       => "Sans,12",
                                                 :size       => "5,3",
                                                 :background => "white"),
                              "epscairo" => Dict(:font       => "Sans,12",
                                                 :size       => "5,3",
                                                 :background => "white")
)

# Dicts to store user-specified configuration
default_config() = Dict(:mode => IsJupyterOrJuno ? "ijulia" : "normal",
                        :term => Dict(:terminal => IsJupyterOrJuno ? "svg" : "qt",
                                      :font => "",
                                      :size => "",
                                      :linewidth => "",
                                      :background => "",
                                      :termopts => ""),
                        :axes => Dict(:axis => "",
                                      :xrange => "",
                                      :yrange => "",
                                      :zrange => "",
                                      :fillstyle => "",
                                      :grid => "",
                                      :boxwidth => "",
                                      :xzeroaxis => "",
                                      :yzeroaxis => "",
                                      :zzeroaxis => "",
                                      :keyoptions => "",
                                      :palette => "",
                                      :onlyimpulses => false),
                        :curve => Dict(:plotstyle => "",
                                       :linecolor => "",
                                       :linestyle => "",
                                       :pointtype => "",
                                       :pointsize => "",
                                       :fillcolor => "",
                                       :fillstyle => ""),
                        :print => Dict(:print_term => "pdfcairo",
                                       :print_font => "",
                                       :print_size => "",
                                       :print_linewidth => "",
                                       :print_background => "",
                                       :print_outputfile => ""))

# Set any of Gaston's configuration variables
function set(;reset = false, terminal=config[:term][:terminal],
             mode = config[:mode], kw...)
    global config

    if reset
        config = default_config()
        return nothing
    end

    t = terminal
    mode == "ijulia" && (t = "svg")
    mode == "null" && (t = "dumb")
    if mode == "normal"
        terminal == "pdf" && (t = "pdfcairo")
        terminal == "pnf" && (t = "pnfcairo")
        terminal == "eps" && (t = "epscairo")
    end
    valid_terminal(t)
    config[:term][:terminal] = t
    config[:mode] = mode

    for k in keys(kw)
        k == :plotstyle && valid_plotstyle(kw[k])
        k == :linestyle && valid_linestyle(kw[k])
        k == :pointtype && valid_pointtype(kw[k])
        k == :axis && valid_axis(kw[k])
        k == :xrange && valid_range(kw[k])
        k == :yrange && valid_range(kw[k])
        k == :zrange && valid_range(kw[k])
        flag = true
        for i in [:term, :axes, :curve, :print]
            c = config[i]
            haskey(c, k) && (flag=false; c[k] = kw[k])
        end
        flag && throw(MethodError(set, "invalid setting"))
    end

    return nothing
end

### Encode terminal capabilities
# supports multiple windows
const term_window = ["qt", "wxt", "x11", "aqua"]
# outputs text
const term_text = ["dumb", "sixelgd"]
# outputs to a file
const term_file = ["svg", "gif", "pngcairo", "pdfcairo", "epscairo"]
# supports size
const term_sup_size = ["qt", "wxt", "x11", "sixelgd", "svg", "gif",
                       "dumb", "pngcairo", "pdfcairo", "epscairo"]
# supports font
const term_sup_font = ["qt", "wxt", "x11", "aqua", "sixelgd", "svg",
                       "gif","pngcairo", "pdfcairo", "epscairo"]
# supports linewidth
const term_sup_lw = ["qt", "wxt", "x11", "aqua", "sixelgd", "svg",
                     "gif", "pngcairo", "pdfcairo", "epscairo"]
# supports background color
const term_sup_bkgnd = ["sixelgd", "svg", "wxt", "gif", "pdfcairo", "pngcairo", "epscairo"]

# List of valid configuration values
const supported_terminals = ["", "qt", "wxt", "x11", "aqua", "dumb", "sixelgd",
                             "svg", "gif", "pngcairo", "pdfcairo", "epscairo"]
const supported_2Dplotstyles = ["", "lines", "linespoints", "points",
                                "impulses", "boxes", "errorlines", "errorbars",
                                "dots", "steps", "fsteps", "fillsteps",
                                "financebars"]
const supported_3Dplotstyles = ["", "lines", "linespoints", "points", "labels",
                                "impulses", "pm3d", "image", "rgbimage", "dots"]
const supported_plotstyles = vcat(supported_2Dplotstyles, supported_3Dplotstyles)
const supported_axis = ["", "normal", "semilogx", "semilogy", "semilogz",
                        "loglog"]
const supported_pointtypes = ["", "+", "x", "*", "esquare", "fsquare",
                          "ecircle", "fcircle", "etrianup", "ftrianup",
                          "etriandn", "ftriandn", "edmd", "fdmd"]
# List of plotstyles that support points
const ps_sup_points = ["linespoints", "points"]

#
# Validation functions
#

function valid_file_term(s)
    s ∈ term_file && return true
    throw(DomainError(s,"supported terminals are: $term_file"))
end

function valid_terminal(s)
    s ∈ supported_terminals && return true
    throw(DomainError(s,"supported terminals are: $supported_terminals"))
end
function valid_plotstyle(s)
    s ∈ supported_plotstyles && return true
    throw(DomainError(s,"supported plotstyles are: $supported_plotstyles"))
end
function valid_2Dplotstyle(s)
    s ∈ supported_2Dplotstyles && return true
    throw(DomainError(s,"supported 2-D plotstyles are: $supported_2Dplotstyles"))
end
function valid_3Dplotstyle(s)
    s ∈ supported_3Dplotstyles && return true
    throw(DomainError(s,"supported 3-D plotstyles are: $supported_3Dplotstyles"))
end
function valid_pointtype(s)
    s ∈ supported_pointtypes && return true
    length(s) == 1 && return true
    throw(DomainError(s,"supported point types are: $supported_pointtypes or single-character UTF-8 strings"))
end
function valid_axis(s)
    s ∈ supported_axis && return true
    throw(DomainError(s,"supported axis types are: $supported_axis"))
end

function valid_linestyle(s)
    invalid = false
    s == "" && return true # allow empty string
    c = collect(s)
    # make sure only allowed characters are passed
    issubset(c, Set([' ', '-', '_', '.'])) || (invalid = true)
    # but do not allow spaces only
    unique(c) != [' '] || (invalid = true)
    invalid && throw(DomainError(s,"line style pattern accepts: space, dash, underscore and dot"))
    return true
end

# Validate that a given range follows gnuplot's syntax.
function valid_range(s::String)
    s == "" && return true # allow empty strings
    # floating point, starting with a dot
    f1 = "[-+]?\\.\\d+([eE][-+]?\\d+)?"
    # floating point, starting with a digit
    f2 = "[-+]?\\d+(\\.\\d*)?([eE][-+]?\\d+)?"
    # floating point
    f = "($f1|$f2)"
    # autoscale directive (i.e. `*` surrounded by
    # optional bounds lb < * < ub)
    as = "(($f\\s*<\\s*)?\\*(\\s*<\\s*$f)?)"
    # full range item: a floating point, or an
    # autoscale directive, or nothing
    it = "(\\s*($as|$f)?\\s*)"

    # empty range
    er = "\\[\\s*\\]"
    # full range: two colon-separated items
    fr = "\\[$it:$it\\]"

    # range regex
    rx = Regex("^\\s*($er|$fr)\\s*\$")

    if occursin(rx, s)
        return true
    end

    throw(DomainError(s,"range must have have the form of [x:y]"))
end

# Validate coordinates
function valid_coords(x,y;err=ErrorCoords(),fin=FinancialCoords())
    invalid = false
    length(x) != length(y) && (invalid = true)
    (err.valid && length(x) != length(err.ylow)) && (invalid = true)
    (err.valid && !isempty(err.yhigh) && length(x) !=
        length(err.yhigh)) && (invalid = true)
    (fin.valid && length(x) != length(fin.open)) && (invalid = true)
    (fin.valid && length(x) != length(fin.low)) && (invalid = true)
    (fin.valid && length(x) != length(fin.high)) && (invalid = true)
    (fin.valid && length(x) != length(fin.close)) && (invalid = true)

    invalid && throw(DimensionMismatch("input vectors must have the same nu mber of elements."))

    return true
end

ext: jl


text:
#
#      Gr(k,V) - The space of all k-dimensional linear subspaces
#           of the n-dimensional vector space C^n
#  A point on this manifold is a vector subspace of dimension k embedded in C^n
#  represented by a n×k matrix with full rank.
#  A point in an tangential space is a matrix.
#
import LinearAlgebra: norm, nullspace, det, tr, qr, triu, rank, svd, Diagonal, diag
import Base: exp, log, show, cat, rand, Matrix, real, atan
export Grassmannian, GrPoint, GrTVector, getValue
export distance, dot, exp, log, manifoldDimension, norm, retraction, inverseRetraction
export parallelTransport, randomTVector, randomMPoint, validateMPoint, validateTVector
export projection, zeroTVector
#
# Type definitions
#

@doc doc"""
    Grassmannian <: Manifold

The manifold $\mathcal M = \mathrm{Gr}(k,n)$ of the set of k-dimensional
subspaces in $\mathbb{K}^{n}$ represented by $n\times k$
orthonormal matrices, that are either real- or complex-valued.

# Abbreviation
`Gr`

# Constructor
    Grassmannian(k, n[, d=Float64])

generates the manifold $\mathrm{Gr}(k,n)$ where the integer `n` is the
number of rows and `k` is the number of columns of the matrices and the optional
parameter `d` sets the `DataType` of the matrix entries.
"""
struct Grassmannian{T<:Union{U, Complex{U}} where U<:AbstractFloat} <: Manifold
  name::String
  dimensionsubspace::Int
  dimensionvecspace::Int
  abbreviation::String
  function Grassmannian{T}(dimensionsubspace::Int, dimensionvecspace::Int) where T<:Union{U, Complex{U}} where U<:AbstractFloat
    if dimensionsubspace > dimensionvecspace
      throw(ErrorException("dimensionsubspace can't be bigger than dimensionvecspace: $dimensionsubspace > $dimensionvecspace"))
    else
     new("of the set of $dimensionsubspace-dimensional subspaces in K$dimensionvecspace in $T", dimensionsubspace, dimensionvecspace,"Gr($dimensionsubspace,$dimensionvecspace)")
    end
  end
end
Grassmannian(dimensionsubspace::Int, dimensionvecspace::Int, D::DataType=Float64) = Grassmannian{D}(dimensionsubspace, dimensionvecspace)

@doc doc"""
    GrPoint <: MPoint

A point $x$ on the manifold $\mathcal M = \mathrm{Gr}(k,n)$ is an orthonormal
matrix of size $ n×k$. The matrix is a numerical representation of the vector
subspace its columns span.

# Constructor

    GrPoint(p)

where `p::Matrix` is an orthonormal matrix of dimension $n×k$.
"""
struct GrPoint{T<:Union{U, Complex{U}} where U<:AbstractFloat} <: MPoint
  value::Matrix{T}
  GrPoint{T}(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = new(value)
end
GrPoint(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = GrPoint{T}(value)
getValue(x::GrPoint) = x.value;

@doc doc"""
    GrTVector <: TVector

A tangent vector $\xi \in T_x\mathcal M$ on the manifold
$\mathcal M = \mathrm{Gr}(k,n)$. The tangent space is given by as

```math
T_x\mathrm{Gr}(k,n) = \bigl\{
\xi \in \mathbb{K}^{n\times k} \ \big|
\ {\bar ξ}^\mathrm{T}x = {\bar x}^\mathrm{T}ξ = 0_{k×k} \bigr\}.
```

# Constructor

    GrTVector(ξ)

where `ξ` is an $n\times k$ `Matrix` that satisfies the above.
"""
struct GrTVector{T<:Union{U, Complex{U}} where U<:AbstractFloat} <: TVector
  value::Matrix{T}
  GrTVector{T}(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = new(value)
end
GrTVector(value::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat = GrTVector{T}(value)
getValue(ξ::GrTVector) = ξ.value;

@traitimpl IsMatrixM{Grassmannian}
@traitimpl IsMatrixP{GrPoint}
@traitimpl IsMatrixTV{GrTVector}

@doc doc"""
    dot(M,x,ξ,ν)

compute the Riemannian inner product for two [`GrTVector`](@ref)s `ξ` and `ν`
from $T_x\mathcal M$ of the [`Grassmannian`](@ref) manifold `M` given by

$\langle \xi, \nu \rangle_x = \operatorname{Re}(\operatorname{tr}(\xi^{\mathrm{T}} ν)),$
"""
function dot(M::Grassmannian{T}, x::GrPoint{T}, ξ::GrTVector{T}, ν::GrTVector{T})::Float64 where T<:Union{U, Complex{U}} where U<:AbstractFloat
  return real(tr(getValue(ξ)' * getValue(ν)))
end

@doc doc"""
    distance(M,x,y)

compute the Riemannian distance on [`Grassmannian`](@ref) manifold `M`$= \mathrm{Gr}(k,n)$ embedded in
$\mathbb R$. Let $USV = {\bar x}^\mathrm{T}y$ denote the SVD decomposition of
$x'y$. Then we compute

$d_{\mathrm{GR}(k,n)}(x,y) = \operatorname{norm}(\operatorname{Re}(b))$

where

$b_{i}=\begin{cases} 0 & \text{if} \; S_i≧1 \\ \operatorname{acos}(S_i) & \, \text{if} \; S_i<1 \end{cases}.$
"""
function distance(M::Grassmannian{T},x::GrPoint{T},y::GrPoint{T})::Float64 where T<:Union{U, Complex{U}} where U<:AbstractFloat
  if x==y
	return 0
  else
    z = getValue(x)'*getValue(y)
    a = svd(z).S
    b = zero(a)
    b[a.<1] = acos.(a[a.<1])
	return norm(real(b), 2)
  end
end

@doc doc"""
    exp(M,x,ξ,[t=1.0])

compute the exponential map on the [`Grassmannian`](@ref) manifold `M`$= \mathrm{Gr}(k,n)$ with
respect to the [`GrPoint`](@ref) `x` and the [`GrTVector`](@ref) `ξ`, which can
be shortened with `t` to `tξ`. Let $USV = t\xi$ denote the SVD decomposition of
$t\xi$. Then we compute

$A = x\cdot V\cdot \cos(S)\cdot {\bar V}^\mathrm{T} + U \cdot \sin(S) \cdot {\bar V}^\mathrm{T},$

where cosine and sine are applied element wise to the diagonal entries of $S$.
The resulting point $y$ of the exponential map is then the matrix $Q$ of the
QR decomposition $A=QR$ of $A$.
"""
function exp(M::Grassmannian{T},x::GrPoint{T},ξ::GrTVector{T},t::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  ξ = t*ξ
  d = svd(getValue(ξ))
  Y = getValue(x) * d.V * cos(Diagonal(d.S)) * (d.V)' + (d.U) * sin(Diagonal(d.S)) * (d.V)'
  GrPoint( Matrix( qr(Y).Q ) )
end

@doc doc"""
    inverseRetraction(M,x,y)

return a [`GrTVector`](@ref) `ξ` of the tagent space $T_x\mathrm{Gr}(k,n)$
with which the [`GrPoint`](@ref) `y` can be reached on the
[`Grassmannian`](@ref) manifold `M`by the
[`retraction`](@ref) from the [`GrPoint`](@ref) `x` after time 1.
The formula reads

$\xi = y \cdot (x^{\mathrm{T}}y)^{-1} - x$

This function is implemented only for the case $\mathbb{K}=\mathbb{R}$.
This is also the standard retraction.
"""
function inverseRetraction(M::Grassmannian{T},x::GrPoint{T},y::GrPoint{T}) where T<:AbstractFloat
  U = getValue(y)/(transpose(getValue(x))*getValue(y))
  A = U-getValue(x)
  GrPoint{T}(A)
end

@doc doc"""
    log(M,x,y)

compute the logarithmic map on the [`Grassmannian`](@ref) manifold
$\mathcal M=\mathrm{Gr}(k,n)$, i.e. the [`GrTVector`](@ref) whose corresponding
[`geodesic`](@ref) starting from [`GrPoint`](@ref) `x` reaches the
[`GrPoint`](@ref) `y` after time 1 on the [`Grassmannian`](@ref) manifold `M`.
The formula reads

$\log_x y = V\cdot \operatorname{atan}(S) \cdot {\bar U}^\mathrm{T}$

where $U$ and $V$ are the unitary matrices and $S$ is a diagonal matrix containing the
singular values of the SVD-decomposition of

$USV = ({\bar y}^\mathrm{T}x)^{-1} ( {\bar y}^\mathrm{T} - {\bar y}^\mathrm{T}x{\bar x}^\mathrm{T} )$

and the $\operatorname{atan}$ is meant elementwise.
"""
function log(M::Grassmannian{T},x::GrPoint{T},y::GrPoint{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  Z = getValue(y)'*getValue(x)
  if det(Z)≠0
   W = Z\(getValue(y)'-Z*getValue(x)')
   d = svd(W, full = false)
   A = (d.V)*atan(Diagonal(svd(W, full = false).S))*(d.U)'
   GrTVector{T}(A)
  else
   throw( ErrorException("The points $x and $y are antipodal, thus these input parameters are invalid.") )
  end
end

@doc doc"""
    manifoldDimension(x)

return the dimension of the [`Grassmannian`](@ref) manifold `M`$= \mathrm{Gr}(k,n)$, the
[`GrPoint`](@ref) `x`, itself embedded in $\mathbb{K}^{n\times k}$, belongs to.
The dimension for $\mathbb{K}=\mathbb{R}$ is defined by

$k(n-k)$

and for $\mathbb{K}=\mathbb{C}$

$2k \cdot (n-k).$
"""
function manifoldDimension(x::GrPoint{T}) where T<:AbstractFloat
	return Int(size(getValue(x),2) * (size(getValue(x),1) - size(getValue(x),2)))
end

function manifoldDimension(x::GrPoint{T}) where T<:Complex{U} where U<:AbstractFloat
	return Int(2 * size(getValue(x),2) * (size(getValue(x),1) - size(getValue(x),2)))
end

@doc doc"""
    manifoldDimension(M)

return the dimension of the [`Grassmannian`](@ref) manifold `M`.
The dimension for $\mathbb{K}=\mathbb{R}$ is defined by

$k(n-k)$

and for $\mathbb{K}=\mathbb{C}$

$2k \cdot (n-k).$
"""
function manifoldDimension(M::Grassmannian{T}) where T<:AbstractFloat
	return Int(M.dimensionsubspace * (M.dimensionvecspace - M.dimensionsubspace))
end

function manifoldDimension(M::Grassmannian{T}) where T<:Complex{U} where U<:AbstractFloat
	return Int(2 * M.dimensionsubspace * (M.dimensionvecspace - M.dimensionsubspace))
end

@doc doc"""
    norm(M,x,ξ)

compute the norm of the [`GrTVector`](@ref) `ξ` in the tangent space
$T_x\mathcal M$ at [`GrPoint`](@ref) `x` of the [`Grassmannian`](@ref) manifold `M`.

$\lVert \xi \rVert_x = \sqrt{\sum_{i,j=0}^n \xi_{ij}^2}$

where $\xi_{ij}$ are the entries of the matrix `ξ`, i.e. the norm is given by
the Frobenius norm of `ξ`.
"""
function norm(M::Grassmannian{T}, x::GrPoint{T}, ξ::GrTVector{T})::Float64 where T<:Union{U, Complex{U}} where U<:AbstractFloat
  norm(getValue(ξ))
end

@doc doc"""
    parallelTransport(M,x,y,ξ)

compute the paralllel transport of the [`GrTVector`](@ref) `ξ` from
the tangent space $T_x\mathcal M$ at [`GrPoint`](@ref) `x` to
$T_y\mathcal M$ at [`GrPoint`](@ref) `y` on the [`Grassmannian`](@ref) manifold `M` provided
that the corresponding [`geodesic`](@ref) $g(\cdot;x,y)$ is unique.
The formula reads

$P_{x\to y}(\xi) = \operatorname{proj}_{\mathcal M}(y,\xi).$

where $\operatorname{proj}_{\mathcal M}$ is the [`projection`](@ref) onto the
tangent space $T_y\mathcal M$.
"""
function parallelTransport(M::Grassmannian{T}, x::GrPoint{T}, y::GrPoint{T}, ξ::GrTVector{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
 GrTVector(projection(M,y,getValue(ξ)))
end

@doc doc"""
    projection(M,x,q)

project a matrix q orthogonally on the [`GrPoint`](@ref) `x` of the manifold
[`Grassmannian`](@ref) manifold `M`. The formula reads

$\operatorname{proj}_{\mathcal M}(x,q) = q-x({\bar x}^\mathrm{T}q),$

i.e. the difference matrix of the image and the output matrix lies in
the orthogonal complement of all [`GrTVector`](@ref)s from the tangent space
$T_x\mathcal M$ at [`GrPoint`](@ref) `x`.

# see also
[`parallelTransport`](@ref), [`randomTVector`](@ref)
"""
function projection(M::Grassmannian{T},x::GrPoint{T},q::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  A = getValue(x)'*q
  B = q - getValue(x)*A
  return B
end

@doc doc"""
    randomMPoint(M [,type=:Gaussian, σ=1.0])

return a random [`GrPoint`](@ref) `x` on [`Grassmannian`](@ref) manifold `M` by
generating a random (Gaussian) matrix with standard deviation `σ` in matching
size, which is orthonormal.
"""
function randomMPoint(M::Grassmannian{T}, ::Val{:Gaussian}, σ::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  V = σ * randn(T, (M.dimensionvecspace, M.dimensionsubspace))
  A = qr(V).Q[:,1:M.dimensionsubspace]
  GrPoint{T}( A )
end

@doc doc"""
    randomTVector(M,x [,type=:Gaussian, σ=1.0])

return a (Gaussian) random vector [`GrTVector`](@ref) in the tangential space
$T_x\mathrm{Gr}(k,n)$ with mean zero and standard deviation `σ` by projecting
a random Matrix onto the  [`GrPoint`](@ref) `x` with [`projection`](@ref).
"""
function randomTVector(M::Grassmannian{T},x::GrPoint{T}, ::Val{:Gaussian}, σ::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  y = σ * randn(T, (M.dimensionvecspace,M.dimensionsubspace))
  Y = projection(M, x, y)
  Y = Y/norm(Y)
  GrTVector(Y)
end

@doc doc"""
    retraction(M,x,ξ,[t=1.0])

move the [`GrPoint`](@ref) `x` in the direction of the [`GrTVector`](@ref) `ξ`
on the  [`Grassmannian`](@ref) manifold `M`. This SVD-based retraction is an
approximation of the exponential map [`exp`](@ref). Let

$USV = x + tξ$

be the singular value decomposition, then the formula reads

$\operatorname{retr}_x\xi = U{\bar V}^\mathrm{T}.$
"""
function retraction(M::Grassmannian{T},x::GrPoint{T},ξ::GrTVector{T},t::Float64=1.0) where T<:Union{U, Complex{U}} where U<:AbstractFloat
	y = getValue(x) + t * getValue(ξ)
    S = svd(y)
    A = S.U * S.V'
    GrPoint{T}(A)
end

@doc doc"""
    zeroTVector(M,x)

return a zero tangent vector in the tangent space of the [`GrPoint`](@ref) on
the [`Grassmannian`](@ref) manifold `M`.
"""
function zeroTVector(M::Grassmannian{T}, x::GrPoint{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  return GrTVector{T}(zeros(M.dimensionvecspace,M.dimensionsubspace))
end

function validateMPoint(M::Grassmannian{T}, x::GrPoint{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  if size(getValue(x), 1) ≠ M.dimensionvecspace
    throw( ErrorException("The dimension of $x must be $(M.dimensionvecspace) × $(M.dimensionsubspace) but it is $(size(getValue(x), 1)) × $(size(getValue(x), 2)).") )
  end
  if size(getValue(x), 2) ≠ M.dimensionsubspace
    throw( ErrorException("The dimension of $x must be $(M.dimensionvecspace) × $(M.dimensionsubspace) but it is $(size(getValue(x), 1)) × $(size(getValue(x), 2)).") )
  end
  if rank(getValue(x)) ≠ M.dimensionsubspace
    throw( ErrorException("$x must have  full rank.") )
  end

  z=getValue(x)'*getValue(x)
  if z ≉ one(z)
    throw( ErrorException("$x has to be orthonormal but it's not.") )
  end
end

function validateTVector(M::Grassmannian{T}, x::GrPoint{T}, ξ::GrTVector{T},) where T<:Union{U, Complex{U}} where U<:AbstractFloat
  if size(getValue(ξ), 1) ≠ M.dimensionvecspace
    throw( ErrorException("The dimension of $ξ must be $(M.dimensionvecspace) × $(M.dimensionsubspace) but it is $(size(getValue(ξ), 1)) × $(size(getValue(ξ), 2))") )
  end
  if size(getValue(ξ), 2) ≠ M.dimensionsubspace
    throw( ErrorException("The dimension of $ξ must be $(M.dimensionvecspace) × $(M.dimensionsubspace) but it is $(size(getValue(ξ), 1)) × $(size(getValue(ξ), 2))") )
  end
  if norm(getValue(x)'*getValue(ξ) + getValue(ξ)'*getValue(x)) > 10^(-15)
    throw( ErrorException("The matrix $x'$ξ must be skew-symmetric!") )
  end
end


text:
using Tensors
using BenchmarkTools
using ForwardDiff

const SUITE = BenchmarkGroup()
const ALL_DIMENSIONS = true
const MIXED_SYM_NONSYM = true
const MIXED_ELTYPES = true

const dT = ForwardDiff.Dual{Nothing,Float64,4}

function create_tensors()
    tensor_dict = Dict{Tuple{Int, Int, DataType}, AbstractTensor}()
    symtensor_dict = Dict{Tuple{Int, Int, DataType}, AbstractTensor}()
    for dim in 1:3
        for order in (1,2,4)
            for T in (Float32, Float64)
                tensor_dict[(dim, order, T)] = rand(Tensor{order, dim, T})
                if order != 1
                    symtensor_dict[(dim, order, T)] = rand(SymmetricTensor{order, dim, T})
                else
                    symtensor_dict[(dim, order, T)] = rand(Tensor{order, dim, T})
                end
            end
            tensor_dict[(dim, order, dT)] = Tensor{order, dim, dT}(([ForwardDiff.Dual(rand(5)...,) for i in 1:length(rand(Tensor{order, dim}))]...,))
            if order != 1
                symtensor_dict[(dim, order, dT)] = SymmetricTensor{order, dim, dT}(([ForwardDiff.Dual(rand(5)...,) for i in 1:length(rand(SymmetricTensor{order, dim}).data)]...,))
            else
                symtensor_dict[(dim, order, dT)] = Tensor{order, dim, dT}(([ForwardDiff.Dual(rand(5)...,) for i in 1:length(rand(Tensor{order, dim}).data)]...,))
            end
        end
    end
    return tensor_dict, symtensor_dict
end

tensor_dict, symtensor_dict = create_tensors()

include("benchmark_functions.jl")
include("benchmark_ad.jl")


text:
export GoToDoor

using Random

mutable struct GoToDoor{W<:GridWorldBase, R} <: AbstractGridWorld
    world::W
    agent_pos::CartesianIndex{2}
    agent::Agent
    target::Door
    target_reward::Float64
    penalty::Float64
    reward::Float64
    rng::R
end

function GoToDoor(;n=8, agent_start_pos=CartesianIndex(2,2), rng=Random.GLOBAL_RNG)
    doors = [Door(c) for c in COLORS[1:4]]
    objects = (EMPTY, WALL, doors...)

    world = GridWorldBase(objects, n, n)

    world[EMPTY, 2:n-1, 2:n-1] .= true

    target_reward = 1.0
    penalty = -1.0
    reward = 0.0

    env = GoToDoor(world, agent_start_pos, Agent(dir=RIGHT), doors[begin], target_reward, penalty, reward, rng)

    reset!(env)

    return env
end

function (env::GoToDoor)(::MoveForward)
    env.reward = 0.0
    dir = get_dir(env.agent)
    dest = dir(env.agent_pos)
    if dest ∈ CartesianIndices((size(env.world, 2), size(env.world, 3))) && !env.world[WALL,dest]
        env.agent_pos = dest
        if get_terminal(env)
            if env.world[env.target, env.agent_pos]
                env.reward = env.target_reward
            else
                env.reward = env.penalty
            end
        end
    end
    env
end

RLBase.get_state(env::GoToDoor) = (get_agent_view(env), env.target)

RLBase.get_terminal(env::GoToDoor) = any([env.world[x, env.agent_pos] for x in env.world.objects[end-3:end]])

function RLBase.reset!(env::GoToDoor)
    n = size(env.world)[end]

    env.world[WALL, [1,n], 1:n] .= true
    env.world[WALL, 1:n, [1,n]] .= true

    doors = env.world.objects[end-3:end]
    for door in doors
        env.world[door, :, :] .= false
    end

    env.target = rand(env.rng, doors)
    env.reward = 0.0

    door_pos = [CartesianIndex(rand(env.rng, 2:n-1),1),
                CartesianIndex(rand(env.rng, 2:n-1),n),
                CartesianIndex(1,rand(env.rng, 2:n-1)),
                CartesianIndex(n,rand(env.rng, 2:n-1))]

    rp = randperm(env.rng, length(door_pos))

    for (door, pos) in zip(doors, door_pos[rp])
        env.world[door, pos] = true
        env.world[WALL, pos] = false
    end

    env
end


text:
# This file was generated, do not modify it. # hide
k_range = range(krb, :(knn_model.K), lower=2, upper=100, scale=:log10)
l_range = range(krb, :(ridge_model.lambda), lower=1e-4, upper=10, scale=:log10)
w_range = range(krb, :(knn_weight), lower=0.1, upper=0.9)

ranges = [k_range, l_range, w_range]


text:
using ForwardDiff, Roots, Plots, LinearAlgebra, Optim

function strong_backtracking(f, ∇, x, d; α=1, β=0.1, σ=0.4)
    y0, g0, y_prev, α_prev = f(x), ∇(x)⋅d, NaN, 0
    αlo, αhi = NaN, NaN
    # bracket phase
    while true
        y = f(x + α*d)
        if y > y0 + β*α*g0 || (!isnan(y_prev) && y ≥ y_prev)
            αlo, αhi = α_prev, α
            break
        end
        g = ∇(x + α*d)⋅d
        if abs(g) ≤ -σ*g0
            return α
        elseif g ≥ 0
            αlo, αhi = α, α_prev
            break
        end
        y_prev, α_prev, α = y, α, 2α
    end

    # zoom phase
    ylo = f(x + αlo*d)
    while true
        α = (αlo + αhi)/2
        y = f(x + α*d)
        if y > y0 + β*α*g0 || y ≥ ylo
            αhi = α
        else
            g = ∇(x + α*d)⋅d
            if abs(g) ≤ -σ*g0
                return α
            elseif g*(αhi - αlo) ≥ 0
                αhi = αlo
            end
            αlo = α
        end
    end
end

f(x) = ((x[1]-1.0)^2) + (0.1*(x[2]-4)^2) + sin(5*pi*x[1])+sin(5*pi*x[2])
grad_f(x) = ForwardDiff.gradient(f,x)

# Ponto inicial
xk = [-1.0,1.0]

#direcao de busca
dk = -grad_f(xk)/norm(grad_f(xk))

beta = 0.1
alpha_f = 1.0

phi(alpha) = f(xk+alpha*dk)
d_phi(alpha) = grad_f(xk+alpha*dk)'dk

l_alpha(alpha) = phi(0)+alpha*beta*d_phi(0)
sigma = 0.4

wolf = []

for alpha = 0:0.01:alpha_f
    if (phi(alpha) <= l_alpha(alpha)) && (abs(grad_f(xk+alpha*dk)'dk) <= -sigma*grad_f(xk)'dk)
        push!(wolf,alpha)
    end
end

wolfbt = strong_backtracking(f,grad_f,xk,dk)

plot(phi,0,alpha_f,lw=2,lab="phi")
plot!(l_alpha,0,alpha_f,lw=2,lab="L_phi",c=:black)
scatter!(wolf,phi.(wolf),lab="StrongWolfe",ms=3,c=:red)
scatter!([wolfbt],[phi(wolfbt)],lab="StrongWolf - BT",c=:yellow)
plot!(legend=:outerbottomright)
plot!(title="beta = 0.1 e sigma = 0.4",xlabel = "alpha", ylabel = "phi(alpha), l_alpha(alpha")


text:
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(
    joinpath(pkgpath, "test", "3_PowerSimulations_examples", "01_operations_problems.jl"),
);

using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

ed_template = template_economic_dispatch(network = StandardPTDFModel)

ed_template.devices[:Hydro] = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)

PTDF_matrix = PTDF(sys)

problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 4,
    optimizer = solver,
    balance_slack_variables = true,
    constraint_duals = [:CopperPlateBalance, :network_flow],
    PTDF = PTDF_matrix,
)

res = solve!(problem);

λ = convert(Array, res.dual_values[:CopperPlateBalance])
μ = convert(Array, res.dual_values[:network_flow])

buses = get_components(Bus, sys)
congestion_lmp = Dict()
for bus in buses
    congestion_lmp[get_name(bus)] = μ * PTDF_matrix[:, get_number(bus)]
end
congestion_lmp = DataFrame(congestion_lmp)

LMP = λ .- congestion_lmp

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
#
# Advent of Code
# Day 7: Internet Protocol Version 7
#
println("Advent of Code")
println("Day 7: Internet Protocol Version 7")
println("")

# Data file
# Input data has been saved to this file
inputDataFileName = "day07-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

# part 1 test data
#inputData = """abba[mnop]qrst
#abcd[bddb]xyyx
#aaaa[qwer]tyui
#ioxxoj[asdfgh]zxcvbn"""

# part 2 test data
#inputData = """aba[bab]xyz
#xyx[xyx]xyx
#aaa[kek]eke
#zazbz[bzb]cdb"""

inputDataLines = split(strip(inputData), '\n')

function isABBA(window::SubString{String})
  # If string matches ABBA, return true
  # NB: AAAA is not allowed

  if window[1] == window[4]
    if window[2] == window[3]
      if window[1] != window[2]
        return true
      end
    end
  end

  return false
end

function isABA(window::SubString{String})
  # If string matches ABA, return true
  # NB: AAA is not allowed

  if window[1] == window[3]
    if window[1] != window[2]
      return true
    end
  end
  return false
end

function isBAB(ABA::SubString{String}, window::SubString{String})
  # If string is inverse of ABA: BAB, return true
  #
  if ABA[1] == window[2]
    if ABA[2] == window[1]
      if ABA[2] == window[3]
        if ABA[3] == window[2]
          return true
        end
      end
    end
  end
  return false
end

nFoundPart1 = 0
nFoundPart2 = 0

for line in inputDataLines
  inHypernet = false
  supportsTLS = false
  supportsSSL = false
  ABA = ""
  hasABA = false
  rejected = false

  for counter in 1:length(line)-2
    # step through all the characters,
    # setting state if we're in the Hypernet card
    # or not
    if line[counter] == '['
      inHypernet = true
    elseif line[counter] == ']'
      inHypernet = false
    end

    # check part 1 ABBA pattern first
    if counter <= length(line) - 3
      if isABBA(line[counter:counter+3])
        # ABBA patterns are only "counted" outside the Hypernet cards
        if !inHypernet
          # don't break out, there may be a pattern in a
          # subsequent Hypernet card that rejects this IP
          supportsTLS = true
        else
          # rejected. Can't break out, due to part 2 processing
          supportsTLS = false
          rejected=true
        end
      end
    end

    # only evaluate "SSL" ABA:BAB patterns if we haven't found
    # one yet. This avoids double-counting.

    if !supportsSSL
      # Find each ABA pattern, one at a time
      if isABA(line[counter:counter+2])
        if !inHypernet
          # now we need to search for matching BAB patterns
          # in the Hypernet cards
          hasABA = true
          ABA = line[counter:counter+2]
          inSubHypernet = false # need to keep track of this separately

          # It seems awful to loop within the loop but we need
          # to capture cases where the ABA occurs after the BAB-containing
          # hypernet cards
          for subCounter in 1:length(line)-2
            if line[subCounter] == '['
              inSubHypernet = true
            elseif line[subCounter] == ']'
              inSubHypernet = false
            end

            # BABs only occur in the Hypernet cards
            if inSubHypernet
              if isBAB(ABA, line[subCounter:subCounter+2])
                # it's okay to break out here, there's nothing that
                # will invalidate this if we find a BAB match
                supportsSSL = true
                nFoundPart2 += 1
                break
              end
            end
          end
        end
      end
    end
  end

  # part 1 summing
  if supportsTLS && !rejected
    nFoundPart1 += 1
  end
end

println("Found $nFoundPart1 TLS-supporting IPs.")
println("Found $nFoundPart2 SSL-supporting IPs.")


text:
# Defines the Curve and UncertainCurve types for handling (uncertain) x-y data

"""
    AbstractCurve

Abstract supertype for all curves (xy-data).
Subtypes need to implement fields `x` and `y`.
"""
abstract type AbstractCurve end

Base.eltype(s::AbstractCurve) = eltype(s.x)
Base.length(s::AbstractCurve) = length(s.x)
Base.lastindex(s::AbstractCurve) = length(s)
Base.getindex(s::AbstractCurve, i::Integer) = (s.x[i], s.y[i])
Base.getindex(s::AbstractCurve, r::UnitRange) = typeof(s)(s.x[r], s.y[r])
Base.iterate(s::AbstractCurve, i::Int64) = i > length(s) ? nothing : (s[i], i + 1)
Base.iterate(s::AbstractCurve) = iterate(s, 1)
Base.IteratorSize(itertype::Type{AbstractCurve}) = Base.HasLength()
Base.minimum(s::AbstractCurve) = s[argmin(s.y)]
Base.maximum(s::AbstractCurve) = s[argmax(s.y)]

Base.:(==)(c0::AbstractCurve, c1::AbstractCurve) = c0.x == c1.x && c0.y == c1.y
Base.:(+)(c::AbstractCurve, y) = typeof(c)(c.x, c.y .+ y)
Base.:(-)(c::AbstractCurve, y) = typeof(c)(c.x, c.y .- y)
Base.:(*)(c::AbstractCurve, y) = typeof(c)(c.x, c.y .* y)
Base.:(/)(c::AbstractCurve, y) = typeof(c)(c.x, c.y ./ y)
Base.:(+)(y, c::AbstractCurve) = typeof(c)(c.x, y .+ c.y)
Base.:(-)(y, c::AbstractCurve) = typeof(c)(c.x, y .- c.y)
Base.:(*)(y, c::AbstractCurve) = typeof(c)(c.x, y .* c.y)
Base.:(/)(y, c::AbstractCurve) = typeof(c)(c.x, y ./ c.y)
Base.vcat(c0::AbstractCurve, c1::AbstractCurve) = typeof(c0)(vcat(c0.x, c1.x), vcat(c0.y, c1.y))

"""
    stitch(curves...)

Append curves one after the other.

Useful when cutting out signals from `Curves` for baseline correction purposes.
"""
stitch(curves...) = reduce(vcat, curves)

function Base.sort(c::T) where {T <: AbstractCurve}
    x, y = c.x, c.y
    ind = sortperm(x)
    return T(x[ind], y[ind])
end

function Base.show(io::IO, c::AbstractCurve)
    println(io, "$(typeof(c)), $(length(c)) datapoints")
    if length(c) > 10
        ps = [c[1:5]; c[end-4:end]]
        indicator = "    ⋮\n"
    else
        ps = c
        indicator = ""
    end
    for (i, p) in enumerate(ps)
        println(io, p)
        i == 5 ? print(io, indicator) : nothing
    end
end

"""
    Curve{T} <: AbstractCurve

Datatype holding x-y data.
x and y vectors have to have the same length.

# Fields
- `x :: Vector{T}` : x-grid
- `y :: Vector{T}` : y-values

# Constructors

    Curve(x, y)
    Curve(y)

# Note

Both `x` and `y` must be convertabile to vectors with the same element type.

# Examples

```jldoctest
julia> Curve([2.0, 3.0], [4.0, 5.0])
Curve{Float64}, 2 datapoints
(2.0, 4.0)
(3.0, 5.0)

julia> Curve([2.0, 3.0], [5.0])
ERROR: ArgumentError: x and y need to have the same length.
[...]

julia> Curve([4, 5, 6])
Curve{Int64}, 3 datapoints
(1, 4)
(2, 5)
(3, 6)

julia> Curve(1:3, [1., 2., 3.])
Curve{Float64}, 3 datapoints
(1.0, 1.0)
(2.0, 2.0)
(3.0, 3.0)
```
"""
struct Curve{T} <: AbstractCurve
    x::Vector{T}
    y::Vector{T}
    function Curve{T}(x, y) where {T <: Number}
        verify_same_length(x, y)
        return new{T}(x, y)
    end
end

function Curve(y)
    T = eltype(y)
    return Curve{T}(collect(T, 1:length(y)), y)
end

Curve(x::Vector{T}, y::Vector{T}) where T = Curve{T}(x, y)

function Curve(x, y)
    verify_same_length(x, y)
    xypairs = promote.(x, y)
    x = map(first, xypairs)
    y = map(last, xypairs)
    T = eltype(y)
    Curve{T}(x, y)
end

"""
    UncertainCurve{T} <: AbstractCurve

Datatype holding x-y data where y-data holds a vector of particles to model uncertainty.
`x` and `y` vectors have to have the same length.

# Fields
- `x::Vector{T}` : x-grid
- `y::Vector{Particles{T, N}}` : uncertain y-values

# Constructors

    UncertainCurve(x::Vector{T}, y::Vector{Particles{N, T}})
    UnertainCurve(y::Vector{T})

# Examples

```jldoctest
julia> using MonteCarloMeasurements

julia> UncertainCurve([2.0, 3.0], [4.0 ± 1.0 , 5.0 ± 1.0])
UncertainCurve{Float64, 2000}, 2 datapoints
(2.0, 4.0 ± 1.0)
(3.0, 5.0 ± 1.0)

julia> UncertainCurve([2.0], [4.0 ± 1.0 , 5.0 ± 1.0])
ERROR: ArgumentError: x and y need to have the same length.
[...]

julia> UncertainCurve([3.0 ± 1.0, 4.0 ± 1.0 , 5.0 ± 1.0])
UncertainCurve{Float64, 2000}, 3 datapoints
(1.0, 3.0 ± 1.0)
(2.0, 4.0 ± 1.0)
(3.0, 5.0 ± 1.0)
```
"""
struct UncertainCurve{T, N} <: AbstractCurve
    x::Vector{T}
    y::Vector{Particles{T, N}}
    function UncertainCurve{T, N}(x::Vector{T}, y::Vector{Particles{T, N}}) where {T, N}
        verify_same_length(x, y)
        return new{T, N}(x, y)
    end
end
UncertainCurve(x::Vector{T}, y::Vector{Particles{T, N}}) where {T, N} = UncertainCurve{T, N}(x, y)
UncertainCurve(y::Vector{Particles{T, N}}) where {T, N} = UncertainCurve(collect(T, 1:length(y)), y)

"""
    get_draw(n, p::Particles)

Retrieve the `n`th particle of the particles stored in `p`.
"""
get_draw(n, p::Particles) = p.particles[n]

"""
    get_draw(n, uc::UncertainCurve)

Retrieve the `n`th sample of the samples stored in `UncertainCurve` `uc`.
"""
get_draw(n, uc::UncertainCurve) = Curve(uc.x, get_draw.(n, uc.y))

"""
    mean(uc::UncertainCurve)

Retrieve the mean of the `UncertainCurve` `uc`.
"""
Statistics.mean(uc::UncertainCurve) = Curve(uc.x, Statistics.mean(uc.y))

MonteCarloMeasurements.:(±)(c::Curve, v) = UncertainCurve(c.x, c.y ± v)

# utility functions

"""
    crop(s::AbstractCurve, left, right)

Crop a slice from a curve between `left` and `right` and return a new curve.

# Examples

```jldoctest
julia> c = Curve(1:0.5:10, 11:0.5:20)
Curve{Float64}, 19 datapoints
(1.0, 11.0)
(1.5, 11.5)
(2.0, 12.0)
(2.5, 12.5)
(3.0, 13.0)
    ⋮
(8.0, 18.0)
(8.5, 18.5)
(9.0, 19.0)
(9.5, 19.5)
(10.0, 20.0)

julia> crop(c, 2, 3)
Curve{Float64}, 3 datapoints
(2.0, 12.0)
(2.5, 12.5)
(3.0, 13.0)

julia> crop(c, 2.1, 3.5)
Curve{Float64}, 3 datapoints
(2.5, 12.5)
(3.0, 13.0)
(3.5, 13.5)

julia> crop(c, 2.1, 3.9)
Curve{Float64}, 3 datapoints
(2.5, 12.5)
(3.0, 13.0)
(3.5, 13.5)

julia> crop(c, 2.1, 4.0)
Curve{Float64}, 4 datapoints
(2.5, 12.5)
(3.0, 13.0)
(3.5, 13.5)
(4.0, 14.0)
```
"""
function crop(s::AbstractCurve, left, right)
    T = eltype(s.x)
    S = typeof(s)
    i = searchsorted(s.x, T(left))  |> first
    j = searchsorted(s.x, T(right)) |> last
    return S(s.x[i:j], s.y[i:j])
end


text:
module cash # Cash

Quantity :: Number #number of a bankNote
Quality ::  Number     # how much of a bankNote

#TODO: add bankNotes (functionality)

# Dr cash

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using StochasticPrograms
using JuMP
using GLPKMathProgInterface
using LinearAlgebra
using Test

import StochasticPrograms: probability, expected

struct SPResult
    x̄::Vector{Float64}
    VRP::Float64
    EWS::Float64
    EVPI::Float64
    VSS::Float64
    EV::Float64
    EEV::Float64
end

problems = Vector{Tuple{StochasticProgram,SPResult,String}}()
@info "Loading test problems..."
@info "Loading simple..."
include("simple.jl")
@info "Loading instant simple..."
include("instant_simple.jl")
@info "Loading deferred simple..."
include("deferred.jl")
@info "Loading farmer..."
include("farmer.jl")
@info "Loading sampler..."
include("sampler.jl")
@info "Test problems loaded. Starting test sequence."

@testset "Stochastic Programs" begin
    @testset "SP Constructs: $name" for (sp,res,name) in problems
        @test optimize!(sp) == :Optimal
        @test norm(optimal_decision(sp)-res.x̄) <= 1e-2
        @test abs(optimal_value(sp)-res.VRP) <= 1e-2
        @test abs(EWS(sp)-res.EWS) <= 1e-2
        @test abs(EVPI(sp)-res.EVPI) <= 1e-2
        @test abs(VSS(sp)-res.VSS) <= 1e-2
        @test abs(EV(sp)-res.EV) <= 1e-2
        @test abs(EEV(sp)-res.EEV) <= 1e-2
    end
    @testset "Inequalities: $name" for (sp,res,name) in problems
        @test EWS(sp) <= VRP(sp)
        @test VRP(sp) <= EEV(sp)
        @test VSS(sp) >= 0
        @test EVPI(sp) >= 0
        @test VSS(sp) <= EEV(sp)-EV(sp)
        @test EVPI(sp) <= EEV(sp)-EV(sp)
    end
    @testset "Deferred model creation" begin
        @test decision_length(deferred) == 0
        @test nscenarios(deferred) == 2
        @test nsubproblems(deferred) == 0
        @test optimize!(deferred) == :Optimal
        @test decision_length(deferred) == 2
        @test nscenarios(deferred) == 2
        @test nsubproblems(deferred) == 2
        @test abs(optimal_value(deferred)+855.83) <= 1e-2
    end
    @testset "Copying: $name" for (sp,res,name) in problems
        sp_copy = copy(sp)
        add_scenarios!(sp_copy, scenarios(sp))
        @test nscenarios(sp_copy) == nscenarios(sp)
        generate!(sp_copy)
        @test nsubproblems(sp_copy) == nsubproblems(sp)
        @test optimize!(sp_copy) == :Optimal
        optimize!(sp)
        @test norm(optimal_decision(sp_copy)-optimal_decision(sp)) <= 1e-2
        @test abs(optimal_value(sp_copy)-optimal_value(sp)) <= 1e-2
        @test abs(EWS(sp_copy)-EWS(sp)) <= 1e-2
        @test abs(EVPI(sp_copy)-EVPI(sp)) <= 1e-2
        @test abs(VSS(sp_copy)-VSS(sp)) <= 1e-2
        @test abs(EV(sp_copy)-EV(sp)) <= 1e-2
        @test abs(EEV(sp_copy)-EEV(sp)) <= 1e-2
    end
    @testset "Sampling" begin
        sampled_sp = SAA(simple_model, sampler, 100, solver=GLPKSolverLP())
        @test nscenarios(sampled_sp) == 100
        @test nsubproblems(sampled_sp) == 100
        @test abs(stage_probability(sampled_sp)-1.0) <= 1e-6
        sample!(sampled_sp, sampler, 100)
        @test nscenarios(sampled_sp) == 200
        @test nsubproblems(sampled_sp) == 200
        @test abs(stage_probability(sampled_sp)-1.0) <= 1e-6
    end
    @testset "Confidence intervals" begin
        glpk = GLPKSolverLP()
        CI = confidence_interval(simple_model, sampler, N = 100, M = 10, confidence = 0.95, solver = glpk)
        @test lower(CI) <= upper(CI)
        sol = optimize(simple_model, sampler, solver = glpk, confidence = 0.95)
        @test lower(confidence_interval(sol)) <= upper(confidence_interval(sol))
    end
end

@info "Starting distributed tests..."

include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
disttestfile = joinpath(@__DIR__, "run_dtests.jl")
push!(test_exeflags.exec,"--color=yes")
cmd = `$test_exename $test_exeflags $disttestfile`

if !success(pipeline(cmd; stdout=stdout, stderr=stderr)) && ccall(:jl_running_on_valgrind, Cint, ()) == 0
    @error "Distributed test failed, cmd : $cmd"
end


text:
using DashBootstrapComponents, DashHtmlComponents

function make_button(placement)
    return dbc_button(
        "Tooltip on $placement",
        id="tooltip-target-$placement",
        className="mx-2",
        n_clicks=0,
    )
end;

function make_tooltip(placement)
    return dbc_tooltip(
        "Tooltip on $placement",
        target="tooltip-target-$placement",
        placement=placement,
    )
end;

tooltips = html_div(
    [[make_button(p) for p in ["top", "left", "right", "bottom"]] ; [make_tooltip(p) for p in ["top", "left", "right", "bottom"]]],
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using Flux: testmode!
using Flux.Tracker: data

@testset "Dropout" begin
  x = [1.,2.,3.]
  @test x == testmode!(Dropout(0.1))(x)
  @test x == Dropout(0)(x)
  @test zero(x) == Dropout(1)(x)

  x = rand(100)
  m = Dropout(0.9)
  y = m(x)
  @test count(a->a==0, y) > 50
  testmode!(m)
  y = m(x)
  @test count(a->a==0, y) == 0
  testmode!(m, false)
  y = m(x)
  @test count(a->a==0, y) > 50

  x = rand(100)
  m = Chain(Dense(100,100),
            Dropout(0.9))
  y = m(x)
  @test count(a->a == 0, y) > 50
  testmode!(m)
  y = m(x)
  @test count(a->a == 0, y) == 0
end

@testset "BatchNorm" begin
  let m = BatchNorm(2), x = param([1 3 5;
                                   2 4 6])

    @test m.β.data == [0, 0]  # initβ(2)
    @test m.γ.data == [1, 1]  # initγ(2)
    # initial m.σ is 1
    # initial m.μ is 0
    @test m.active

    # @test m(x).data ≈ [-1 -1; 0 0; 1 1]'
    m(x)

    # julia> x
    #  2×3 Array{Float64,2}:
    #  1.0  3.0  5.0
    #  2.0  4.0  6.0
    #
    # μ of batch will be
    #  (1. + 3. + 5.) / 3 = 3
    #  (2. + 4. + 6.) / 3 = 4
    #
    # ∴ update rule with momentum:
    #  .1 * 3 + 0 = .3
    #  .1 * 4 + 0 = .4
    @test m.μ ≈ reshape([0.3, 0.4], 2, 1)

    # julia> .1 .* var(x, dims = 2, corrected=false) .* (3 / 2).+ .9 .* [1., 1.]
    # 2×1 Array{Float64,2}:
    #  1.3
    #  1.3
    @test m.σ² ≈ .1 .* var(x.data, dims = 2, corrected=false) .* (3 / 2).+ .9 .* [1., 1.]

    testmode!(m)
    @test !m.active

    x′ = m(x).data
    @test isapprox(x′[1], (1 .- 0.3) / sqrt(1.3), atol = 1.0e-5)
  end

  # with activation function
  let m = BatchNorm(2, sigmoid), x = param([1 3 5;
                                            2 4 6])
    @test m.active
    m(x)

    testmode!(m)
    @test !m.active

    y = m(x).data
    @test isapprox(y, data(sigmoid.((x .- m.μ) ./ sqrt.(m.σ² .+ m.ϵ))), atol = 1.0e-7)
  end

  let m = BatchNorm(2), x = param(reshape(1:6, 3, 2, 1))
    y = reshape(permutedims(x, [2, 1, 3]), 2, :)
    y = permutedims(reshape(m(y), 2, 3, 1), [2, 1, 3])
    @test m(x) == y
  end

  let m = BatchNorm(2), x = param(reshape(1:12, 2, 3, 2, 1))
    y = reshape(permutedims(x, [3, 1, 2, 4]), 2, :)
    y = permutedims(reshape(m(y), 2, 2, 3, 1), [2, 3, 1, 4])
    @test m(x) == y
  end

  let m = BatchNorm(2), x = param(reshape(1:24, 2, 2, 3, 2, 1))
    y = reshape(permutedims(x, [4, 1, 2, 3, 5]), 2, :)
    y = permutedims(reshape(m(y), 2, 2, 2, 3, 1), [2, 3, 4, 1, 5])
    @test m(x) == y
  end

  let m = BatchNorm(32), x = randn(Float32, 416, 416, 32, 1);
    m(x)
    @test (@allocated m(x)) <  100_000_000
  end
end

@testset "InstanceNorm" begin
  # helper functions
  expand_inst = (x, as) -> reshape(repeat(x, outer=[1, as[length(as)]]), as...)
  # begin tests
  let m = InstanceNorm(2), sizes = (3, 2, 2),
      x = param(reshape(collect(1:prod(sizes)), sizes))

      @test m.β.data == [0, 0]  # initβ(2)
      @test m.γ.data == [1, 1]  # initγ(2)

      @test m.active

      m(x)

      #julia> x
      #[:, :, 1] =
      # 1.0  4.0
      # 2.0  5.0
      # 3.0  6.0
      #
      #[:, :, 2] =
      # 7.0  10.0
      # 8.0  11.0
      # 9.0  12.0
      #
      # μ will be
      # (1. + 2. + 3.) / 3 = 2.
      # (4. + 5. + 6.) / 3 = 5.
      #
      # (7. + 8. + 9.) / 3 = 8.
      # (10. + 11. + 12.) / 3 = 11.
      #
      # ∴ update rule with momentum:
      # (1. - .1) * 0 + .1 * (2. + 8.) / 2 = .5
      # (1. - .1) * 0 + .1 * (5. + 11.) / 2 = .8
      @test m.μ ≈ [0.5, 0.8]
      # momentum * var * num_items / (num_items - 1) + (1 - momentum) * sigma_sq
      # julia> reshape(mean(.1 .* var(x.data, dims = 1, corrected=false) .* (3 / 2), dims=3), :) .+ .9 .* 1.
      # 2-element Array{Float64,1}:
      #  1.
      #  1.
      @test m.σ² ≈ reshape(mean(.1 .* var(x.data, dims = 1, corrected=false) .* (3 / 2), dims=3), :) .+ .9 .* 1.

      testmode!(m)
      @test !m.active

      x′ = m(x).data
      @test isapprox(x′[1], (1 - 0.5) / sqrt(1. + 1f-5), atol = 1.0e-5)
  end
  # with activation function
  let m = InstanceNorm(2, sigmoid), sizes = (3, 2, 2),
      x = param(reshape(collect(1:prod(sizes)), sizes))

    affine_shape = collect(sizes)
    affine_shape[1] = 1

    @test m.active
    m(x)

    testmode!(m)
    @test !m.active

    y = m(x).data
    @test isapprox(y, data(sigmoid.((x .- expand_inst(m.μ, affine_shape)) ./ sqrt.(expand_inst(m.σ², affine_shape) .+ m.ϵ))), atol = 1.0e-7)
  end

  let m = InstanceNorm(2), sizes = (2, 4, 1, 2, 3),
      x = param(reshape(collect(1:prod(sizes)), sizes))
    y = reshape(permutedims(x, [3, 1, 2, 4, 5]), :, 2, 3)
    y = reshape(m(y), sizes...)
    @test m(x) == y
  end

  # check that μ, σ², and the output are the correct size for higher rank tensors
  let m = InstanceNorm(2), sizes = (5, 5, 3, 4, 2, 6),
      x = param(reshape(collect(1:prod(sizes)), sizes))
    y = m(x)
    @test size(m.μ) == (sizes[end - 1], )
    @test size(m.σ²) == (sizes[end - 1], )
    @test size(y) == sizes
  end

  # show that instance norm is equal to batch norm when channel and batch dims are squashed
  let m_inorm = InstanceNorm(2), m_bnorm = BatchNorm(12), sizes = (5, 5, 3, 4, 2, 6),
      x = param(reshape(collect(1:prod(sizes)), sizes))
    @test m_inorm(x) == reshape(m_bnorm(reshape(x, (sizes[1:end - 2]..., :, 1))), sizes)
  end

  let m = InstanceNorm(32), x = randn(Float32, 416, 416, 32, 1);
    m(x)
    @test (@allocated m(x)) <  100_000_000
  end

end


text:
#coast.msh geometry, created by JIGSAW
#adapted by @miguelraz 28/03/2020 from
#https://github.com/dengwirda/jigsaw-matlab
#mshid=1
#ndims=2
#point=303
const lakes = [(-8.915414699999999,1.661592)
(-8.884706400000001,1.5538653)
(-8.7270571,1.4396306)
(-8.6780329,1.3310253)
(-8.2763869,1.3642085)
(-7.7589437,1.5521208)
(-7.3399811,1.5866412)
(-7.2064877,1.7400111)
(-6.987161,1.7838988)
(-6.9190874,1.6751018)
(-6.6594362,1.823788)
(-6.4705078,2.0285109)
(-6.3982459,2.0257761)
(-6.1776033,2.1236043)
(-5.9683636,1.9041636)
(-6.0500631,1.6420917)
(-6.2751037,1.4382319)
(-6.2083449,1.2768143)
(-6.3593873,1.1233594)
(-6.4379782,0.9672913)
(-6.4389687,0.940829)
(-6.3308569,0.9103191)
(-5.921668,1.1606956)
(-5.885353,1.1594465)
(-5.5306989,0.88275)
(-5.4571118,0.9068911)
(-5.2776995,0.8218321)
(-4.6882658,1.0696109)
(-4.4325551,1.115784)
(-3.9867787,1.5284986)
(-3.6232783,1.573356)
(-3.5156848,1.5181578)
(-3.1151141,1.6162773)
(-2.6031886,1.9253444)
(-2.4201615,2.0814581)
(-2.0576784,2.182551)
(-1.9469493,2.3929843)
(-1.6564457,2.6016249)
(-1.4385855,2.8113568)
(-1.185718,2.9681802)
(-0.896576,3.2839699)
(-0.4656316,3.4939341)
(-0.0715676,3.6520868)
(0.3221563,3.5994536)
(0.6089023,3.4943945)
(0.6092879,3.3885375)
(0.3942456,3.3878929)
(0.2866788,3.4141384)
(0.0717038,3.3345153)
(0.0896864,3.2286668)
(0.1076576,3.1757488)
(-0.1436793,3.01699)
(-0.3235837,2.8584538)
(-0.4676923,2.7529346)
(-0.9391942,2.0668111)
(-1.0495193,1.7498742)
(-1.0882363,1.3531376)
(-1.0158452,1.3262275)
(-0.9424102,1.4845979)
(-0.615238,1.7477533)
(-0.5061946000000001,1.9061816)
(-0.2890742,2.0115316)
(-0.506667,1.7473959)
(-0.6160022000000001,1.5360392)
(-0.397972,1.8000363)
(-0.1808403,1.8525961)
(0.2894343,1.7998173)
(0.5431103,1.6681146)
(0.9058832,1.5373304)
(0.9432855999999999,1.3258124)
(1.0534256,1.1147328)
(1.1641909,0.8507940000000001)
(1.3109203,0.6930545)
(1.3117251,0.5871977)
(1.421469,0.5351345)
(1.9317215,0.5400872)
(2.2226054,0.5965222)
(2.3340522,0.4391798)
(2.4073552,0.4137156)
(2.4449497,0.334835)
(2.5909063,0.3369578)
(2.6249835,0.4962917)
(2.7343488,0.4979874)
(2.8454551,0.3938967)
(2.936647,0.3954205)
(2.9201944,0.2892562)
(2.9931921,0.2905022)
(3.3135431,0.7199108)
(3.3863589,0.7213257)
(3.9960329,1.1315466)
(4.0705567,1.0538628)
(4.6490632,1.174476)
(4.903998,1.1550887)
(5.4495532,1.1450567)
(5.9850466,1.4278414)
(6.1653805,1.4607402)
(6.3476152,1.4409063)
(6.6355959,1.5048353)
(6.5309084,1.3948027)
(6.4623697,1.2862027)
(6.4803085,0.8098819)
(6.4094979,0.7542354)
(6.410009,0.7012564)
(6.4882671,0.5981826)
(6.7068647,0.6065329)
(6.8546905,0.559329)
(6.9296957,0.5093111)
(7.3625916,0.6332439)
(7.4195041,0.5826584)
(7.4810256,0.4262426)
(7.6622,0.4606821)
(7.7339418,0.4903887)
(7.6927729,0.5946095)
(7.6199031,0.5913803)
(7.5788199,0.6956243)
(7.4649955,0.7967085)
(7.4604098,0.9025555)
(7.6668095,1.176764)
(7.770946,1.2875089)
(7.7298813,1.3917132)
(7.6707366,1.4951115)
(7.6356734,1.4670289)
(7.6246141,1.3074531)
(7.5532118,1.2777792)
(7.5280861,1.4357493)
(7.5345391,1.7011603)
(7.7131429,1.762169)
(7.8603274,1.7158544)
(7.9228745,1.9308907)
(7.7348899,2.0813402)
(7.3713821,2.1182087)
(7.2698318,1.9547751)
(7.1885795,2.1633926)
(7.1152657,2.1868122)
(7.0441691,2.15734)
(6.9336956,2.2058041)
(6.9676141,2.2602054)
(7.0648074,2.529299)
(7.0897106,2.7954185)
(7.2917362,3.1221593)
(7.3230152,3.2295634)
(7.2074219,3.4101256)
(7.0987913,3.4319941)
(7.056259,3.5892427)
(6.8349672,3.7390777)
(6.5426013,3.8863685)
(6.3538119,4.1440578)
(6.6493349,4.7917061)
(6.5385166,4.8932842)
(6.7052993,5.1650313)
(6.7009764,5.2708721)
(5.7441174,5.2345817)
(5.5686945,5.1755508)
(4.788607,5.1509481)
(4.4636817,5.3536511)
(4.2075761,5.6646354)
(3.607578,7.1858876)
(3.3873644,7.6576363)
(3.2804384,7.7613007)
(3.1762847,7.7327191)
(3.0381962,7.6506117)
(2.8286881,7.6467634)
(2.8598071,7.8590921)
(2.6845531,7.9089847)
(2.5468001,7.8008409)
(2.4437671,7.6933492)
(2.3390419,7.6917801)
(2.1983122,7.7691861)
(2.1616242,7.9010276)
(2.0918985,7.9000964)
(2.0235167,7.7933405)
(1.9540682,7.7660064)
(1.570508,7.7353016)
(1.5002116,7.7875578)
(1.2547448,7.9442315)
(0.730965,8.1526824)
(0.4522013,8.2575018)
(0.1738657,8.3098811)
(-0.0694533,8.521515300000001)
(-0.1737495,8.4157382)
(-0.2084645,8.4422444)
(-0.347615,8.363095700000001)
(-0.4518238,8.389823099999999)
(-0.5558115,8.4695476)
(-0.5215963,8.310643900000001)
(-0.5221188,8.1518584)
(-0.4879606,7.9400337)
(-0.2092654,7.833566)
(-0.2097517,7.463066)
(-0.5257651,7.0403597)
(-0.5604465,7.1463352)
(-0.7716285000000001,6.9354925)
(-0.9473082,6.8835005)
(-0.9833619,6.7249256)
(-1.1593398,6.6731632)
(-1.2287946,6.7795396)
(-1.1574372,6.9378046)
(-0.6994015,7.4115277)
(-0.5933112,7.7286747)
(-0.6973173,7.8878837)
(-0.9759166,7.9422776)
(-1.1175557,7.6256239)
(-1.1194069,7.3609825)
(-1.2951661,7.2564438)
(-1.5483098,6.411761)
(-1.9366363,6.31007)
(-1.8637764,6.5209552)
(-1.7582839,6.5197737)
(-1.7202965,6.7840348)
(-1.5427314,6.9939696)
(-1.5417155,7.0998256)
(-1.5762336,7.1530937)
(-1.9985176,6.9989869)
(-2.3506579,6.8978706)
(-2.5645144,6.6893724)
(-2.5351899,6.3183276)
(-2.7562921,5.7395724)
(-3.7980607,5.0191439)
(-4.3373509,4.7676755)
(-4.9834948,4.5740477)
(-5.9256532,4.1282314)
(-7.0382285,3.1643631)
(-7.4517045,2.7580044)
(-8.203390300000001,2.1245258)
(-8.7235361,1.8639192)
(-6.0222359,1.402651)
(-5.9152018,1.5578835)
(-5.8726343,1.5299217)
(-5.7920161,1.5536639)
(-5.8345898,1.5816055)
(-5.8318773,1.6609928)
(-5.663557,1.7083083)
(-5.5966615,1.7590827)
(-5.6808641,1.8413432)
(-5.7929717,1.7391469)
(-6.0818599,1.5637339)
(-6.0865686,1.4314222)
(-5.7324455,1.9225485)
(-5.5683675,1.9594866)
(-5.5335991,1.9159672)
(-5.5318777,1.9688923)
(-5.449816,2.0457126)
(-5.4914569,2.1000333)
(-5.7297704,2.0019357)
(-5.7677013,1.9502324)
(-5.2706487,2.2254789)
(-5.2257103,2.277064)
(-5.2196558,2.3563407)
(-5.191165,2.4614151)
(-5.2924029,2.4539566)
(-5.3166868,2.3699543)
(-5.3193455,2.2852741)
(-5.344493,2.174812)
(-2.628509,4.8904087)
(-2.5231726,4.8093509)
(-2.4524965,4.7818143)
(-2.166425,4.9101487)
(-1.8814086,4.9860466)
(-2.0588911,4.9881727)
(-2.093713,5.0415487)
(-1.9150532,5.1452403)
(-1.7011793,5.2487523)
(-1.5239785,5.2470098)
(-1.1676824,5.5087407)
(-1.2022868,5.6148538)
(-0.9185046,5.7718016)
(-0.882891,5.8245348)
(-0.9532116,5.8778612)
(-1.3075201,5.7215257)
(-1.5563954,5.5649116)
(-1.6276655,5.5126727)
(-1.7340945,5.4872998)
(-1.840557,5.4619959)
(-2.0900008,5.3326514)
(-2.3756722,5.1777288)
(-2.590489,5.0486276)
(-0.1745905,7.648274)
(-0.1048759,7.4629629)
(-0.0174793,7.4629295)
(-0.0349355,7.5687895)
(0.0349007,7.7275752)
(0.2267794,7.7806613)
(0.3839078,7.7280333)
(0.558318,7.7550129)
(0.6632222,7.7024853)
(0.8201638,7.7296796)
(0.732548,7.8086476)
(0.6273773,7.9405222)
(0.4880418,7.9135694)
(0.1393711,7.9922753)
(0.0348428,7.9922181)
(-0.1743299,7.8864526)
(4.5964504,4.4039327)
(4.8814214,4.4120893)
(4.9866866,4.468199)
(5.1274771,4.5254771)
(5.1614199,4.5795088)
(5.0886198,4.6302255)
(4.9818982,4.6269721)
(4.8364965,4.7285905)
(4.5535316,4.6675729)
(4.3785074,4.5569131)
(4.3813116,4.4510625)]


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -32.35079088597315 <= q <= 33.04707617725232)

                     add_NL_constraint(m, :(sigmoid(0.7516814527569986 + 0.07211062597860263*sigmoid(0.800956367833014 + 0.1575422515944389*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.28099609060437203*sigmoid(-0.8241959750660866 + 0.48752067922981546*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.19401527858092393*sigmoid(-0.014178711859067938 + -0.6496362145686132*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.43991506890783594*sigmoid(-0.38559964579408224 + -0.6882459841770032*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + sigmoid(0.5442086655639207 + -0.859033717055464*sigmoid(0.800956367833014 + 0.1575422515944389*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.6994107397450926*sigmoid(-0.8241959750660866 + 0.48752067922981546*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.21447765318835765*sigmoid(-0.014178711859067938 + -0.6496362145686132*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9247182813449877*sigmoid(-0.38559964579408224 + -0.6882459841770032*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + sigmoid(-0.33457802196328057 + 0.8511359109817849*sigmoid(0.800956367833014 + 0.1575422515944389*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + 0.9056512987027685*sigmoid(-0.8241959750660866 + 0.48752067922981546*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.901983268277478*sigmoid(-0.014178711859067938 + -0.6496362145686132*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.02854306019224717*sigmoid(-0.38559964579408224 + -0.6882459841770032*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) + sigmoid(0.85223750045835 + -0.7468516096711739*sigmoid(0.800956367833014 + 0.1575422515944389*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.30181017969165413*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.5024447348308962*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + -0.7463877606567912*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.1269236474608202*sigmoid(-0.8241959750660866 + 0.48752067922981546*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + -0.8942226763869914*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + -0.1882956347370066*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.0613423792438379*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.5407492362268407*sigmoid(-0.014178711859067938 + -0.6496362145686132*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5744456236018198*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.06990002879133872*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.11738131892224501*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))) + -0.8670343803826994*sigmoid(-0.38559964579408224 + -0.6882459841770032*sigmoid(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4])) + 0.5869261088566509*sigmoid(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4])) + 0.7286086900917508*sigmoid(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4])) + 0.9385542676650069*sigmoid(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function load_a1at_aa_rel_ab()
    _A1AT_AA_REL_AB = Dict()
    datfile = _a1at_aa_rel_abundance_proc_file()
    !isfile(datfile) && error("data file missing ", datfile)
    df = CSV.read(_a1at_aa_rel_abundance_proc_file(), DataFrame)
    for (m1, m2) in zip(df[!,1], df[!,2])
        _A1AT_AA_REL_AB[m1] = m2
    end
    return _A1AT_AA_REL_AB
end


text:
function initial_velocity(data, np)
    constante = -2*data.R*data.T/data.m
    p=[rand() for i in 1:np]
    vmod =  [sqrt(constante * log(1-p[i])) for i in 1:np]
    angle = [2*pi*rand() for i in 1:np]
return [ [vmod[i]*cos(angle[i]), vmod[i]*sin(angle[i])] for i in 1:np ], vmod
end


text:
"""
    plate_uniaxial_local_buckling(A::Array{<:Real,2}, D::Array{<:Real,2}, b::Real)
Calculates local buckling load and strains. `b` is the panel width.

assumes: balanced, symmetric, large aspect ratio, simply supported,
uniaxial compression, flat rectangular plate. buckling strain is output as a
positive value

see chapter on Structural Component Design Techniques from Alastair Johnson
section 6.2: Design of composite panels

probably assumes wave numbers are continuous in order to eliminate them
"""
function plate_uniaxial_local_buckling(A::AbstractArray{<:Real,2}, D::AbstractArray{<:Real,2},
    b::Real)
    bucklingload = 2.0*(pi/b)^2.0*(sqrt(D[1,1]*D[2,2])+D[1,2]+2*D[3,3])
    # bucklingstrain = bucklingload/(A[1,1]-A[1,2]^2.0/A[2,2])
    return bucklingload
end #buckling

# Lekhnitskii (see LOCAL BUCKLING OF COMPOSITE BEAMS G. Tarján, A. Sapkás and L.P. Kollár)
function plate_bending_local_buckling(D::AbstractArray{<:Real,2}, b::Real)
    bucklingload = (pi/b)^2.0*(13.4*sqrt(D[1,1]*D[2,2])+10.4*(D[1,2]+2*D[3,3]))
    return bucklingload
end #buckling

"""
    plate_shear_local_buckling(A::Array{<:Real,2}, D::Array{<:Real,2}, b::Real)
Calculates local buckling load and strains. `b` is the panel width.

assumes: balanced, symmetric, large aspect ratio, simply supported,
uniaxial compression, flat rectangular plate

see chapter on Structural Component Design Techniques from Alastair Johnson
section 6.2: Design of composite panels

probably assumes wave numbers are continuous in order to eliminate them
"""
function plate_shear_local_buckling(A::AbstractArray{<:Real,2}, D::AbstractArray{<:Real,2},
    b::Real)
    D1 = D[1,1]
    D2 = D[2,2]
    D3 = D[1,2] + 2*D[3,3]
    c = (D1*D2)^(1/2)/D3

    # # non-smooth version
    # if c >= 1.0
    #     fc = c^(1/2)*(0.62+0.38/c)
    # else
    #     fc = 0.89 + 0.04*c + 0.07*c^2
    # end

    # smooth version
    if (1 + tanh(20*(c-1))) == 0.0
        fc = 0.89 + 0.04*c + 0.07*c^2
    else
        fc = (0.89 + 0.04*c + 0.07*c^2) + (1 + tanh(20*(c-1)))/2*(c^(1/2)*(0.62+0.38/c)-(0.89 + 0.04*c + 0.07*c^2))
    end

    K2 = (D2*D3)^(1/2)*fc/D1
    bucklingload = 52*K2*D1/b^2
    # bucklingstrain = -bucklingload/A[3,3]
    return bucklingload
end #buckling

"""
    cylinder_uniaxial_local_buckling(A::Array{<:Real,2}, D::Array{<:Real,2}, R::Real, h::Real, correlation::Bool=true)
Calculates local buckling load. `R` is the cylinder radius. `h` is the total laminate thickness.

assumes: infinitely long cylinder, simply supported, uniaxial compression,
wave numbers are real-valued rather than integers (this allows them to be
eliminated as variables)

see An Approximate Solution and Master Curves for Buckling of Symmetrically
Laminated Composite Cylinders by Michael P. Nemeth

also see Simple Formulas and Results for Buckling Resistance and Stiffness Design of
Compression-Loaded Laminated-Composite Cylinders by Michael P. Nemeth

The correlation factor from SP-8001 is also applied, but can be disabled by
setting `correlation=false`.

"""
function cylinder_uniaxial_local_buckling(A, D, R, h, correlation=true)
    nu_m = A[1,2]/sqrt(A[1,1]*A[2,2])
    rho = R/h*sqrt(1-nu_m^2)*(A[1,1]*A[2,2]*h^4/(144*D[1,1]*D[2,2]))^(1/4)
    beta = (D[1,2] + 2*D[3,3])/sqrt(D[1,1]*D[2,2])
    I = (A[1,1]*D[2,2]/(D[1,1]*A[2,2]))^(1/2)
    pcr_tilde = 4*sqrt(3*I)*rho/(pi^2)
    bucklingload = pcr_tilde*pi^2*sqrt(D[1,1]*D[2,2])/(R^2)

    if correlation
        psi = 1/29.8*(R/(D[1,1]*D[2,2]/(A[1,1]*A[2,2]))^(1/4))^(1/2)
        gamma = 1-0.901*(1-exp(-psi))
        bucklingload = gamma*bucklingload
    end

    return bucklingload
end

"""
    cylinder_bending_local_buckling(A::Array{<:Real,2}, D::Array{<:Real,2}, R::Real, h::Real)
Calculates local buckling load. `R` is the cylinder radius. `h` is the total laminate thickness.

Follows approach suggested in SP-8001, which uses the theoretical axial compression cylinder
buckling equation with a correction factor specific for bending cases, except
assumes an infinitely long cylinder and that wave numbers are real-valued rather
than integers (in order to elimate them as variables)

see An Approximate Solution and Master Curves for Buckling of Symmetrically
Laminated Composite Cylinders by Michael P. Nemeth for details on the equation
for the buckling of an infinitely long orthotropic cylinder

"""
function cylinder_bending_local_buckling(A, D, R, h)

    bucklingload = cylinder_bending_local_buckling(A, D, R, h, correlation=false)

    # apply correlation
    psi = 1/29.8*(R/(D[1,1]*D[2,2]/(A[1,1]*A[2,2]))^(1/4))^(1/2)
    gamma = 1-0.731*(1-exp(-psi))
    bucklingload = gamma*bucklingload

    return bucklingload
end


text:
module Uchiyama

using LinearAlgebra
using Random

include("particles.jl")
include("squares.jl")
include("collisions.jl")
include("periodic_collisions.jl")
include("periodic_uchiyama.jl")
include("hard_spheres.jl")
include("event_driven.jl")
include("periodic_hard_spheres.jl")

end


text:
using LinearAlgebra
using Flux
import Flux.gpu, Flux.cpu

# Output layer
mutable struct output{T}
    C2::Union{Matrix{T},CuMatrix{T}}
    D21::Union{Matrix{T},CuMatrix{T}}
    D22::Union{Matrix{T},CuMatrix{T}}
    by::Union{Vector{T},CuVector{T}}
end

function output{T}(nu::Int, nx::Int, nv::Int, ny::Int) where T
    C2 = randn(T, ny, nx) / convert(T, sqrt(ny + nx))
    D21 = randn(T, ny, nv) / convert(T, sqrt(ny + nv))
    D22 = zeros(T, ny, nu)
    by = randn(T, ny) / convert(T, sqrt(ny))
    return output{T}(C2, D21, D22, by)
end

(layer::output)(x, w, u) = layer.C2 * x + layer.D21 * w + layer.D22 * u .+ layer.by
Flux.trainable(layer::output) = (layer.C2, layer.D21, layer.by)

function Flux.gpu(M::output{T}) where T
    if T != Float32
        println("Moving type: ", T, " to gpu may not be supported. Try Float32!")
    end
    return output{T}(gpu(M.C2), gpu(M.D21), gpu(M.D22), gpu(M.by))
end

function Flux.cpu(M::output{T}) where T
    return output{T}(cpu(M.C2), cpu(M.D21), cpu(M.D22), cpu(M.by))
end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
function constrain_sdp(model, s::ConstrainedContinuousIdentitySystem, Qi, Ui, Zi, nQ, nZ)
    Ei = s.X[1]; Di = s.X[2];
    @constraint(model, Qi - Ei'*Ui*Ei - Di'*Zi*Di - Array(I(nQ)) in PSDCone())
    @constraint(model, Zi - Array(I(nZ)) in PSDCone())
end

function constrain_sdp(model, s::ConstrainedLinearMap, Qi, Qj, Uij, Zij, nZ, γ)
    Ei = s.X[1]; Di = s.X[2]; Ai = s.A
    @constraint(model, γ^2*Qi - Ai'*Qj*Ai - Ei'*Uij*Ei - Di'*Zij*Di in PSDCone())
    @constraint(model, Zij - Array(I(nZ)) in PSDCone())
end

function sdplyap_feasible_prog(γ, s::StateDepDiscreteSwitchedLinearSystem, optimizer; verbose=0)
    modes=1:nstates(s)
    model = Model(optimizer)
    nA = size(s.resetmaps[1].A, 1)
    nE = size(s.resetmaps[1].X[1], 1)
    nD = size(s.resetmaps[1].X[2], 1)
    if verbose ≤ 0
        set_silent(model)
    end
    Q_vrefs = Dict(state => @variable(model, [i=1:nA, j=1:nA]) for state in modes)
    Ui_vrefs = Dict(state => @variable(model, [i=1:nE, j=1:nE], lower_bound=0) for state in modes)
    Zi_vrefs = Dict(state => @variable(model, [i=1:nD, j=1:nD]) for state in modes)

    for i in modes
        Qi = Q_vrefs[i]
        constrain_sdp(model, HybridSystems.mode(s, i), Qi, Ui_vrefs[i], Zi_vrefs[i], nA, nD)
        for t in out_transitions(s, i)
            j = target(s, t)
            Qj = Q_vrefs[j]
            Uij_ref = @variable(model, [i=1:nE, j=1:nE], lower_bound=0)
            Zij_ref = @variable(model, [i=1:nD, j=1:nD])
            constrain_sdp(model, resetmap(s, t),  Qi, Qj, Uij_ref, Zij_ref, nD, γ)
        end
    end

    if verbose ≥ 2
        print(model)
    end

    optimize!(model)

    if verbose ≥ 3
        print(model)
    end
    if verbose ≥ 2
        try
            @show JuMP.solution_summary(model, verbose=true)
        catch err
            println("solution_summary requires JuMP version 0.21.7 or later")
        end
    elseif verbose ≥ 1
        try
            @show JuMP.solution_summary(model)
        catch err
            println("solution_summary requires JuMP version 0.21.7 or later")
        end
    end

    status = (JuMP.termination_status(model), JuMP.primal_status(model), JuMP.dual_status(model), JuMP.raw_status(model))
end


text:
export BatchedSampler

mutable struct BatchedGradSampler{BN, P, S, Sc, Sv, Pv, Gv, Gvv, Gva, LC, Lv, Lgv, Pa, Pc, Pd} <: AbstractIterativeSampler
    bnet::BN
    Ĉ::P

    sampler::S
    sampler_cache::Sc
    samples::Sv

    ψvals::Pv
    ∇vals::Gv
    ∇vals_vec::Gvv
    ∇vec_avg::Gva

    accum::LC

    local_vals::Lv
    ∇local_vals::Lgv

    precond_alg::Pa
    precond_cache::Pc

    observables_sampler

    parallel_cache::Pd
end

function BatchedGradSampler(net,
                        sampl,
                        prob,
                        algo=Gradient();
                        batch_sz=2^4,
                        local_batch_sz=batch_sz,
                        obs_sampler=sampl,
                        par_type=automatic_parallel_type())
    if net isa CachedNet
        throw("Only takes standard network.")
    end

    par_cache      = parallel_execution_cache(par_type)

    bnet           = cached(net, batch_sz)
    v              = state(prob, bnet)
    sampler_cache  = cache(sampl, basis(prob), bnet, par_cache)

    ch_len         = chain_length(sampl, sampler_cache)

    samples        = state(prob, bnet, ch_len)
    ψvals          = similar(trainable_first(bnet), out_type(bnet), 1, batch_sz, ch_len)
    ∇vals, ∇vec    = grad_cache(bnet, batch_sz, ch_len)
    ∇vec_avg       = tuple([similar(∇vec_i, size(∇vec_i, 1)) for ∇vec_i=∇vec]...)

    local_acc      = AccumulatorObsGrad(net, basis(prob), v, local_batch_sz)
    Llocal_vals    = collect(similar(ψvals, size(ψvals)[2:end]...)) #ensure it's on cpu #tocheck
    ∇Llocal_vals   = tuple([similar(ψvals, size(∇vec_i, 1), size(ψvals)[2:end]...) for ∇vec_i=∇vec]...)

    precond        = algorithm_cache(algo, prob, net, par_cache)

    if prob isa KLocalLiouvillian
        obs        = BatchedObsDMSampler(bnet, obs_sampler, basis(prob), batch_sz=batch_sz)
    else
        obs        = BatchedObsKetSampler(samples, ψvals, local_acc)
    end

    nq = BatchedGradSampler(bnet, prob,
            sampl, sampler_cache, samples,
            ψvals, ∇vals, ∇vec, ∇vec_avg,
            local_acc, Llocal_vals, ∇Llocal_vals,
            algo, precond,
            obs,
            par_cache)

    return nq
end

function sample!(is::BatchedGradSampler; sample=true)
    ch_len         = chain_length(is.sampler, is.sampler_cache)
    batch_sz       = size(is.local_vals, 1)

    # Sample the new configurations
    sample && _sample_state!(is)

    # Compute logψ and ∇logψ
    logψ_and_∇logψ!(is.∇vals, is.ψvals, is.bnet, is.samples)

    # Compute LdagL
    for i=1:ch_len
        for j = 1:batch_sz
            σv = state_uview(is.samples, j, i)
            init!(is.accum, σv, is.ψvals[1,j,i], vec_data(is.∇vals[i])[1][:,j] )
            accumulate_connections!(is.accum, is.Ĉ, σv)
            L_loc, ∇L_loc = NeuralQuantum.finalize!(is.accum)
            is.local_vals[j, i]   = L_loc
            uview(is.∇local_vals[1], :, j, i) .= vec_data(∇L_loc)[1]
        end
    end

    Ĉ2 = abs2.(is.local_vals)

    L_stat = stat_analysis(Ĉ2, is.parallel_cache)

    # Center the gradient so that it has zero-average
    _center_gradient!(is)

    # Compute the gradient
    # Flatten the batches and the iterations
    ∇vr = reshape(is.∇vals_vec[1], :, ch_len*batch_sz)
    ∇Ĉr = reshape(is.∇local_vals[1], :, ch_len*batch_sz)
    Ĉr  = reshape(is.local_vals, 1, :)
    #Ĉ2r = reshape(Ĉ2, 1, :)

    # Compute the gradient
    ∇C  = Ĉr*∇Ĉr'
    ∇C ./= (ch_len*batch_sz)
    ∇C .-= L_stat.mean .* is.∇vec_avg[1]'
    ∇C .= conj.(∇C)
    workers_mean!(∇C, is.parallel_cache)

    setup_algorithm!(is.precond_cache, reshape(∇C,:), ∇vr, is.parallel_cache)

    return L_stat, is.precond_cache
end

function compute_observables(is::BatchedGradSampler)
    return compute_observables(is.observables_sampler)
end

function add_observable!(is::BatchedGradSampler, name::String, obs)
    return add_observable!(is.observables_sampler, name, obs)
end

Base.show(io::IO, is::BatchedGradSampler) = print(io,
    "BatchedGradSampler for :"*
    "\n\tnet\t\t: $(is.bnet)"*
    "\n\tproblem\t: $(is.Ĉ)"*
    "\n\tsampler\t: $(is.sampler)")


text:
a_data_props = [:pos]
m_data_props = []


text:
# Inflación de Percentiles Equiponderados

# InflationPercentileEq.jl - Definición de percentiles equiponderados

"""
    InflationPercentileEq <: InflationFunction
Función de inflación para computar el percentil equiponderado `k`.

## Utilización

    function (inflfn::InflationPercentileEq)(base::VarCPIBase{T}) where T

Define cómo opera `InflationPercentileEq` sobre un objeto de tipo VarCPIBase.

    function (inflfn::InflationPercentileEq)(cs::CountryStructure)

Define cómo opera `InflationPercentileEq` sobre un objeto de tipo CountryStructure.

### Ejemplo
Cálculo del percentil 70 de la distribución de variaciones intermensuales ponderadas

```julia-repl
julia> percEqfn = InflationPercentileEq(70)
(::InflationPercentileEq) (generic function with 5 methods)
julia>percEqfn(gtdata) #gtdata es de tipo UniformCountryStructure
231-element Vector{Float32}:
9.812069
10.271847
10.382021
10.323441
10.427952
⋮
1.9042492
1.9850135
1.9871473
2.0293117
2.1856546
```
"""
Base.@kwdef struct InflationPercentileEq <: InflationFunction
    k::Float32
end

# Métodos para crear funciones de inflación a partir de enteros y flotantes
"""
    InflationPercentileEq(k::Int) = InflationPercentileEq(k = Float32(k) / 100)
    InflationPercentileEq(q::T) where {T <: AbstractFloat}
Permite utilizar valores de `k` que no necesariamente son de tipo de punto flotante, como por ejemplo: enteros, fracciones, Float32 y Float64.

# Ejemplo
Se usa 0.70 en lugar de 70

```julia-repl
julia> percEqfn = InflationPercentileEq(0.70)
(::InflationPercentileEq) (generic function with 5 methods)

julia>percEqfn(gtdata)
231-element Vector{Float32}:
9.812069
10.271847
10.382021
10.323441
10.427952
⋮
1.9042492
1.9850135
1.9871473
2.0293117
2.1856546
```
Se obtienen los mismos resultados.
"""
InflationPercentileEq(k::Int) = InflationPercentileEq(k = Float32(k) / 100)
function InflationPercentileEq(q::T) where {T <: AbstractFloat}
    q < 1.0 && return InflationPercentileEq(convert(Float32, q))
    InflationPercentileEq(convert(Float32, q) / 100)
end

"""
    measure_name(inflfn::InflationPercentileEq)

Indica qué medida se utiliza para una instancia de una función de inflación.

# Ejemplo
```julia-repl
julia> percEqfn = InflationPercentileEq(0.70)
julia> measure_name(percEqfn)
"Percentil equiponderado 70.0"
```
"""
measure_name(inflfn::InflationPercentileEq) = "Percentil equiponderado " * string(round(100inflfn.k, digits=2))

# Parámetro de la función de inflación
params(inflfn::InflationPercentileEq) = (inflfn.k, )

# Las funciones sobre VarCPIBase resumen en variaciones intermensuales
function (inflfn::InflationPercentileEq)(base::VarCPIBase{T}) where T
    # InflationPercentileEq k de la distribución de variaciones intermensuales
    k = inflfn.k

    # Obtener el percentil k de la distribución intermensual
    rows = size(base.v, 1)

    k_interm = Vector{T}(undef, rows)
    Threads.@threads for r in 1:rows
        row = @view base.v[r, :]
        k_interm[r] = quantile(row, k)
    end

    k_interm
end

function InflationPercentileEq(vec::Vector{<:Real})
    x = vec[1]
    return InflationPercentileEq(x)
end


text:
def arrcpy(a) {
    cpy = []
    for i=0, i<|a|, i=i+1
        cpy :<= a[i]
    cpy
}

def partition(a, lo, hi) {
	{p, i} = (a[hi], lo)
	for j=lo, j < hi, j=j+1
		if a[j] <= p {
			a[i] <-> a[j]
			i = i+1
		}

	a[i] <-> a[hi]
	i
}

def quicksort(a, lo, hi)
	if lo < hi {
		p = partition(a, lo, hi)
		quicksort(a, lo, p-1)
		quicksort(a, p+1, hi)
	}

def qsort(a) quicksort(a, 0, |a|-1)

def qsorted(a) {
    sorted = arrcpy(a)
    qsort(sorted)
    sorted
}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
# cat_data
"""
    cat_data(vec_data::Data...)

Do this will concatinate all model in `vec_data`, returning a variable in `Data_Sum` type
"""
function cat_data(vec_data::Data...)
    data = Data_Sum(vec_data[1])
    num_data = length(vec_data)
    if num_data == 1
        return
    else
        for id = 2 : num_data

            flag_str = 0
            # Str Info
            str_now = [typeof(data.vec_str[n]) for n = 1:length(data.vec_str)]
            #if !in(typeof(vec_data[id].data_str), str_now)
            if !in_str(vec_data[id].data_str, data.vec_str)
                flag_str = 1  # Differnet Structure
                typ_tilt = max(data.vec_atom, "typ")
                vec_data[id].data_str.vec_type_id .+= typ_tilt
                data.vec_str = vcat(data.vec_str, vec_data[id].data_str)
                add!(vec_data[id].vec_atom, typ_tilt, "typ")
            else
                flag_str = 0 # Same Structure
                typ_tilt = data.vec_str[findall(x->x==typeof(vec_data[id].data_str), str_now)[1]].vec_type_id[1]
                typ_tilt = typ_tilt - vec_data[id].vec_atom[1].typ
                add!(vec_data[id].vec_atom, typ_tilt, "typ")
            end

            # Basic Info
            data.data_basic.num_atoms += vec_data[id].data_basic.num_atoms
            data.data_basic.num_bonds += vec_data[id].data_basic.num_bonds
            data.data_basic.num_angles += vec_data[id].data_basic.num_angles
            data.data_basic.num_dihedrals += vec_data[id].data_basic.num_dihedrals
            data.data_basic.num_impropers += vec_data[id].data_basic.num_impropers
            if flag_str == 1
                data.data_basic.num_atom_types += vec_data[id].data_basic.num_atom_types
                data.data_basic.num_bond_types += vec_data[id].data_basic.num_bond_types
                data.data_basic.num_angle_types += vec_data[id].data_basic.num_angle_types
                data.data_basic.num_dihedral_types += vec_data[id].data_basic.num_dihedral_types
                data.data_basic.num_improper_types += vec_data[id].data_basic.num_improper_types
            end
            data.data_basic.box_size[:, 2] = max(data.data_basic.box_size[:, 2], vec_data[id].data_basic.box_size[:, 2])
            data.data_basic.box_size[:, 1] = min(data.data_basic.box_size[:, 1], vec_data[id].data_basic.box_size[:, 1])
            data.data_basic.box_tilt = conv(max(data.data_basic.box_tilt, vec_data[id].data_basic.box_tilt), 1)

            atom_tilt = max(data.vec_atom, "atom")
            # Atom Info
            add!(vec_data[id].vec_atom, atom_tilt, "atom")
            data.vec_atom = vcat(data.vec_atom, vec_data[id].vec_atom)

            # Bond Info
            if typeof(data.vec_bond) == Int64
                if typeof(vec_data[id].vec_bond) != Int64
                    add!(vec_data[id].vec_bond, atom_tilt, "atom")
                    data.vec_bond = vec_data[id].vec_bond
                end
            else
                if typeof(vec_data[id].vec_bond) != Int64
                    add!(vec_data[id].vec_bond, atom_tilt, "atom")
                    data.vec_bond = vcat(data.vec_bond, vec_data[id].vec_bond)
                end
            end

            # Angle Info
            if typeof(data.vec_angle) == Int64
                if typeof(vec_data[id].vec_angle) != Int64
                    add!(vec_data[id].vec_angle, atom_tilt, "atom")
                    data.vec_angle = vec_data[id].vec_angle
                end
            else
                if typeof(vec_data[id].vec_angle) != Int64
                    add!(vec_data[id].vec_angle, atom_tilt, "atom")
                    data.vec_angle = vcat(data.vec_angle, vec_data[id].vec_angle)
                end
            end

            # Dihedral Info
            if typeof(data.vec_dihedral) == Int64
                if typeof(vec_data[id].vec_dihedral) != Int64
                    add!(vec_data[id].vec_dihedral, atom_tilt, "atom")
                    data.vec_dihedral = vec_data[id].vec_dihedral
                end
            else
                if typeof(vec_data[id].vec_dihedral) != Int64
                    add!(vec_data[id].vec_dihedral, atom_tilt, "atom")
                    data.vec_dihedral = vcat(data.vec_dihedral, vec_data[id].vec_dihedral)
                end
            end

            # Improper Info
            if typeof(data.vec_improper) == Int64
                if typeof(vec_data[id].vec_improper) != Int64
                    add!(vec_data[id].vec_improper, atom_tilt, "atom")
                    data.vec_improper = vec_data[id].vec_improper
                end
            else
                if typeof(vec_data[id].vec_improper) != Int64
                    add!(vec_data[id].vec_improper, atom_tilt, "atom")
                    data.vec_improper = vcat(data.vec_improper, vec_data[id].vec_improper)
                end
            end
        end
    end
    data
end

# sort_data
"""
    sort_data!(data::Data_Unit, list_atom::Array)

Do this will rearrange the atom id to consecutive one after calling `remove!`
"""
function sort_data!(data::Data_Unit, list_atom::Array)
    fields = fieldnames(typeof(data))
    name_fields = [string(fields[n]) for n = 1:length(fields)]
    num_fields = length(fields)
    list_fields = findall(x->occursin("vec", x), name_fields)
    list_atom .-= 1    # Find atom next to the deleted one
    num_atoms = length(list_atom)

    for field in list_fields
        vec_now = getfield(data, fields[field])
        if typeof(vec_now) != Int64
            len = length(vec_now)
            dims = length(vec_now[1].atom)
            vec_id = [vec_now[n].atom[i] for n = 1 : len, i = 1 : dims]
            atom_now = 1
            for id = 1 : len-1
                for dim = 1 : dims
                    diff = vec_id[id+1, dim] - vec_id[id, dim]
                    if diff >= 1
                        vec_id[1:id, dim] .+= diff - 1
                    end
                end
            end
            vec_id .-= vec_id[1, 1] - 1
            set!(vec_now, vec_id, "atom")
        end
    end
    data
end