text:
using Unitful

type GaussBeamParams
    waist_radius
    waist_pos
    wavelength
end

"""
Returns the additive unit (zero) of the same type as the given value,
to transparently support unitful and symbolic use cases.
"""
zero_like(x) = x - x

"""
Returns the intensity at the given point of a Gaussian beam with the
specified parameters and total power.

Parameters
----------
power
    The optical power integrated over the whole beam profile.
waist
    The :math:`1 / e^2`-radius of the beam at its waist.
wavelength
    The wavelength of the radiation.
r_offset
    If given, the distance of the point from the focus in the radial direction.
z_offset
    If given, the distance of the point from the focus along the direction of
    propagation.
"""
function gauss_intensity(power, waist, wavelength; r_offset=zero_like(waist),
        z_offset=zero_like(waist))
    gauss_intensity(power, GaussBeamParams(waist, zero_like(waist), wavelength),
        r_offset=r_offset, z_pos=z_offset)
end

function gauss_intensity(power, beam::GaussBeamParams; r_offset=zero_like(waist),
        z_pos=zero_like(waist))
    raleigh_range = π * beam.waist_radius^2 / beam.wavelength
    radius = beam.waist_radius * sqrt(1 + ((z_pos - beam.waist_pos) / raleigh_range)^2)
    2 * power / (π * radius^2) * exp(-2 *(r_offset / radius)^2)
end

export GaussBeamParams, gauss_intensity


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "FluidSynth"
version = v"2.0.3"

# Collection of sources required to build FluidSynth
sources = [
    "https://github.com/FluidSynth/fluidsynth.git" =>
    "1bae9b2fe1a958f54f4910c802a79673e0df9850",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd fluidsynth/
mkdir build
cd build/
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DLIB_SUFFIX="" -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libfluidsynth", :libfluidsynth),
    ExecutableProduct(prefix, "fluidsynth", :fluidsynth)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/staticfloat/GlibBuilder/releases/download/v2.54.2-2/build.jl",
    "https://github.com/staticfloat/libffiBuilder/releases/download/v3.2.1-0/build.jl",
    "https://github.com/yukota/libsndfileBuilder/releases/download/v1.0.28-0/build_libsndfile.v1.0.28.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: yukota/FluidSynthBuilder

max_issues_repo_name: yukota/FluidSynthBuilder

max_forks_repo_name: yukota/FluidSynthBuilder


text:
# CodeTable32
# ===========

const CodeTable32 = CodeTable{32}

const BASE32_CODEPAD = 0x20  # PADding
const BASE32_CODEIGN = 0x21  # IGNore
const BASE32_CODEEND = 0x22  # END
const BASE32_CODEERR = 0xff  # ERRor

ignorecode(::Type{CodeTable32}) = BASE32_CODEIGN

"""
    CodeTable32(asciicode::String, pad::Char; casesensitive::Bool=false)

Create a code table for base32.
"""
function CodeTable32(asciicode::String, pad::Char; casesensitive::Bool=false)
    if !isascii(asciicode) || !isascii(pad)
        throw(ArgumentError("the code table must be ASCII"))
    elseif length(asciicode) != 32
        throw(ArgumentError("the code size must be 32"))
    end
    encodeword = Vector{UInt8}(undef, 32)
    decodeword = Vector{UInt8}(undef, 128)
    fill!(decodeword, BASE32_CODEERR)
    for (i, char) in enumerate(asciicode)
        bits = UInt8(i-1)
        code = UInt8(char)
        encodeword[bits+1] = code
        decodeword[code+1] = bits
        if !casesensitive
            if isuppercase(char)
                code = UInt8(lowercase(char))
                decodeword[code+1] = bits
            end
            if islowercase(char)
                code = UInt8(uppercase(char))
                decodeword[code+1] = bits
            end
        end
    end
    padcode = UInt8(pad)
    decodeword[padcode+1] = BASE32_CODEPAD
    return CodeTable32(encodeword, decodeword, padcode)
end

@inline function encode(table::CodeTable32, byte::UInt8)
    return table.encodeword[Int(byte & 0x1f) + 1]
end

@inline function decode(table::CodeTable32, byte::UInt8)
    return table.decodeword[Int(byte)+1]
end

"""
The standard base32 code table (cf. Table 3 of RFC4648).
"""
const BASE32_STD = CodeTable32("ABCDEFGHIJKLMNOPQRSTUVWXYZ234567", '=')

"""
The extended hex code table (cf. Table 4 of RFC4648).
"""
const BASE32_HEX = CodeTable32("0123456789ABCDEFGHIJKLMNOPQRSTUV", '=')


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license
# Portions derived from LAPACK, see below.

# invariant subspace computations using complex Schur decompositions

import LinearAlgebra: _ordschur!, _ordschur

function _ordschur(T::StridedMatrix{Ty},
                    Z::StridedMatrix{Ty},
                   select::Union{Vector{Bool},BitVector}) where Ty <: Complex
    _ordschur!(copy(T), copy(Z), select)
end

function _ordschur!(T::StridedMatrix{Ty},
                    Z::StridedMatrix{Ty},
                    select::Union{Vector{Bool},BitVector}) where Ty <: Complex
    # suppress most checks since this is an internal function expecting
    # components of a Schur object
    n = size(T,1)
    ks = 0
    for k=1:n
        if select[k]
            ks += 1
            if k != ks
                _trexchange!(T,Z,k,ks)
            end
        end
    end
    triu!(T)
    T,Z,diag(T)
end

"""
reorder `T` by a unitary similarity transformation so that `T[iold,iold]`
is moved to `T[inew,inew]`. Also right-multiply `Z` by the same transformation.
"""
function _trexchange!(T,Z,iold,inew)
    # this is Algorithm 7.6.1 in Golub & van Loan 4th ed.
    n = size(T,1)
    if (n < 1) || (iold == inew)
        return
    end
    if iold < inew
        krange = iold:inew-1
    else
        krange = iold-1:-1:inew
    end
    for k in krange
        Tkk = T[k,k]
        Tpp = T[k+1,k+1]
        G,_ = givens(T[k,k+1], Tpp - Tkk, k, k+1)
        if k+1 <= n
            lmul!(G,T)
        end
        rmul!(T,adjoint(G))
        Z === nothing || rmul!(Z,adjoint(G))
    end
end

# eigvalscond and subspacesep are derived from LAPACK's ztrsen.
# LAPACK is released under a BSD license, and is
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.

"""
    eigvalscond(S::Schur,nsub::Integer) => Real

Estimate the reciprocal of the condition number of the `nsub` leading
eigenvalues of `S`. (Use `ordschur` to move a subspace of interest
to the front of `S`.)

See the LAPACK User's Guide for details of interpretation.
"""
function eigvalscond(S::Schur{Ty},nsub) where Ty
    n = size(S.T,1)
    # solve T₁₁ R - R T₂₂ = σ T₁₂
    R = S.T[1:nsub,nsub+1:n] # copy
    R, scale = trsylvester!(view(S.T,1:nsub,1:nsub),view(S.T,nsub+1:n,nsub+1:n),
                    R)
    rnorm = norm(R) # Frobenius norm, as desired
    s = (rnorm > 0) ?
        scale / (sqrt(scale*scale / rnorm + rnorm ) * sqrt(rnorm)) :
        one(real(Ty))
    s
end

"""
    subspacesep(S::Schur,nsub::Integer) => Real

Estimate the reciprocal condition of the separation angle for the
invariant subspace corresponding
to the leading block of size `nsub` of a Schur decomposition.
(Use `ordschur` to move a subspace of interest
to the front of `S`.)

See the LAPACK User's Guide for details of interpretation.
"""
function subspacesep(S::Schur{Ty}, nsub) where Ty
    n = size(S.T,1)
    scale = one(real(Ty))
    function f(X)
        # solve T₁₁ R - R T₂₂ = σ X
        R, s = trsylvester!(view(S.T,1:nsub,1:nsub),
                            view(S.T,nsub+1:n,nsub+1:n),
                            reshape(X,nsub,n-nsub))
        scale = s
        R
    end
    function fH(X)
        # solve T₁₁ᴴ R - R T₂₂ᴴ = σ X
        R, s = adjtrsylvester!(view(S.T,1:nsub,1:nsub),
                               view(S.T,nsub+1:n,nsub+1:n),
                               reshape(X,nsub,n-nsub))
        scale = s
        R
    end
    est = norm1est!(f,fH,zeros(Ty,nsub*(n-nsub)))
    return scale / est
end


text:
const RTR = ReverseDiff.TrackedReal
const RTV = ReverseDiff.TrackedVector
const RTM = ReverseDiff.TrackedMatrix
const RTA = ReverseDiff.TrackedArray
using ReverseDiff: SpecialInstruction
import NaNMath
using ForwardDiff: Dual
import SpecialFunctions: logbeta
import Distributions: Gamma

Base.:*(A::Adjoint{<:Real, <:RTV{<:Real}}, B::AbstractVector{<:Real}) = dot(A, B)
Base.:*(A::Adjoint{<:Real, <:RTV{<:Real}}, B::RTV{<:Real}) = dot(A, B)
Base.:*(A::AbstractVector{<:Real}, B::Adjoint{<:Real, <:RTV{<:Real}}) = dot(A, B)
Base.:*(A::RTV{<:Real}, B::Adjoint{<:Real, <:RTV{<:Real}}) = dot(A, B)

Gamma(α::RTR, θ::Real; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::Real, θ::RTR; check_args=true) = pgamma(α, θ, check_args = check_args)
Gamma(α::RTR, θ::RTR; check_args=true) = pgamma(α, θ, check_args = check_args)
pgamma(α, θ; check_args=true) = Gamma(promote(α, θ)..., check_args = check_args)
Gamma(α::T; check_args=true) where {T <: RTR} = Gamma(α, one(T), check_args = check_args)
function Gamma(α::T, θ::T; check_args=true) where {T <: RTR}
    check_args && Distributions.@check_args(Gamma, α > zero(α) && θ > zero(θ))
    return Gamma{T}(α, θ)
end

# Work around to stop TrackedReal of Inf and -Inf from producing NaN in the derivative
function Base.minimum(d::LocationScale{T}) where {T <: RTR}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * minimum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(minimum)(d.ρ))
    end
end
function Base.maximum(d::LocationScale{T}) where {T <: RTR}
    if isfinite(minimum(d.ρ))
        return d.μ + d.σ * maximum(d.ρ)
    else
        return convert(T, ReverseDiff.@skip(maximum)(d.ρ))
    end
end

for T in (:RTV, :RTM)
    @eval begin
        function Distributions.logpdf(d::MvNormal{<:Any, <:PDMats.ScalMat}, x::$T)
            logpdf(TuringScalMvNormal(d.μ, d.Σ.value), x)
        end
        function Distributions.logpdf(d::MvNormal{<:Any, <:PDMats.PDiagMat}, x::$T)
            logpdf(TuringDiagMvNormal(d.μ, d.Σ.diag), x)
        end
        function Distributions.logpdf(d::MvNormal{<:Any, <:PDMats.PDMat}, x::$T)
            logpdf(TuringDenseMvNormal(d.μ, d.Σ.chol), x)
        end

        function Distributions.logpdf(d::MvLogNormal{<:Any, <:PDMats.ScalMat}, x::$T)
            logpdf(TuringMvLogNormal(TuringScalMvNormal(d.normal.μ, d.normal.Σ.value)), x)
        end
        function Distributions.logpdf(d::MvLogNormal{<:Any, <:PDMats.PDiagMat}, x::$T)
            logpdf(TuringMvLogNormal(TuringDiagMvNormal(d.normal.μ, d.normal.Σ.diag)), x)
        end
        function Distributions.logpdf(d::MvLogNormal{<:Any, <:PDMats.PDMat}, x::$T)
            logpdf(TuringMvLogNormal(TuringDenseMvNormal(d.normal.μ, d.normal.Σ.chol)), x)
        end
    end
end

# zero mean, dense covariance
MvNormal(A::RTM) = TuringMvNormal(A)

# zero mean, diagonal covariance
MvNormal(σ::RTV) = TuringMvNormal(σ)

# dense mean, dense covariance
MvNormal(m::AbstractVector{<:Real}, A::RTM{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::RTV{<:Real}, A::Matrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::RTV{<:Real}, A::RTM{<:Real}) = TuringMvNormal(m, A)

# dense mean, diagonal covariance
function MvNormal(
    m::RTV{<:Real},
    D::Diagonal{<:RTR, <:RTV{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:RTR, <:RTV{<:Real}},
)
    return TuringMvNormal(m, D)
end
function MvNormal(
    m::RTV{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvNormal(m, D)
end

# dense mean, diagonal covariance
MvNormal(m::RTV{<:Real}, σ::RTV{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::RTV{<:Real}, σ::AbstractVector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::RTV{<:Real}, σ::Vector{<:Real}) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::RTV{<:Real}) = TuringMvNormal(m, σ)

# dense mean, constant variance
MvNormal(m::RTV{<:Real}, σ::RTR) = TuringMvNormal(m, σ)
MvNormal(m::RTV{<:Real}, σ::Real) = TuringMvNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::RTR) = TuringMvNormal(m, σ)

# dense mean, constant variance
function MvNormal(m::RTV{<:Real}, A::UniformScaling{<:RTR})
    return TuringMvNormal(m, A)
end
function MvNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:RTR})
    return TuringMvNormal(m, A)
end
function MvNormal(m::RTV{<:Real}, A::UniformScaling{<:Real})
    return TuringMvNormal(m, A)
end

# zero mean,, constant variance
MvNormal(d::Int, σ::RTR) = TuringMvNormal(d, σ)

# zero mean, dense covariance
MvLogNormal(A::RTM) = TuringMvLogNormal(TuringMvNormal(A))

# zero mean, diagonal covariance
MvLogNormal(σ::RTV) = TuringMvLogNormal(TuringMvNormal(σ))

# dense mean, dense covariance
MvLogNormal(m::RTV{<:Real}, A::RTM{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::RTV{<:Real}, A::Matrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::AbstractVector{<:Real}, A::RTM{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))

# dense mean, diagonal covariance
function MvLogNormal(
    m::RTV{<:Real},
    D::Diagonal{<:RTR, <:RTV{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{<:RTR, <:RTV{<:Real}},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end
function MvLogNormal(
    m::RTV{<:Real},
    D::Diagonal{T, <:AbstractVector{T}} where {T<:Real},
)
    return TuringMvLogNormal(TuringMvNormal(m, D))
end

# dense mean, diagonal covariance
MvLogNormal(m::RTV{<:Real}, σ::RTV{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::RTV{<:Real}, σ::AbstractVector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::RTV{<:Real}, σ::Vector{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))
MvLogNormal(m::AbstractVector{<:Real}, σ::RTV{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, σ))

# dense mean, constant variance
function MvLogNormal(m::RTV{<:Real}, σ::RTR)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::RTV{<:Real}, σ::Real)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end
function MvLogNormal(m::AbstractVector{<:Real}, σ::RTR)
    return TuringMvLogNormal(TuringMvNormal(m, σ))
end

# dense mean, constant variance
function MvLogNormal(m::RTV{<:Real}, A::UniformScaling{<:RTR})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:RTR})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end
function MvLogNormal(m::RTV{<:Real}, A::UniformScaling{<:Real})
    return TuringMvLogNormal(TuringMvNormal(m, A))
end

# zero mean,, constant variance
MvLogNormal(d::Int, σ::RTR) = TuringMvLogNormal(TuringMvNormal(d, σ))

function LinearAlgebra.cholesky(A::RTM; check=true)
    factors, info = turing_chol(A, check)
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end

function turing_chol(x::ReverseDiff.TrackedArray{V,D}, check) where {V,D}
    tp = ReverseDiff.tape(x)
    x_value = ReverseDiff.value(x)
    check_value = ReverseDiff.value(check)
    C, back = pullback(_turing_chol, x_value, check_value)
    out = ReverseDiff.track(C.factors, D, tp)
    ReverseDiff.record!(tp, SpecialInstruction, turing_chol, (x, check), out, (back, issuccess(C)))
    return out, C.info
end

@noinline function ReverseDiff.special_reverse_exec!(instruction::SpecialInstruction{typeof(turing_chol)})
    output = instruction.output
    instruction.cache[2] || throw(PosDefException(C.info))
    input = instruction.input
    input_deriv = ReverseDiff.deriv(input[1])
    P = instruction.cache[1]
    input_deriv .+= P((factors = ReverseDiff.deriv(output),))[1]
    ReverseDiff.unseed!(output)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(instruction::SpecialInstruction{typeof(turing_chol)})
    output, input = instruction.output, instruction.input
    C = cholesky(ReverseDiff.value(input[1]), check = ReverseDiff.value(input[2]))
    ReverseDiff.value!(output, C.factors)
    return nothing
end

Distributions.Dirichlet(alpha::RTV) = TuringDirichlet(alpha)
Distributions.Dirichlet(d::Integer, alpha::RTR) = TuringDirichlet(d, alpha)

function Distributions.logpdf(d::MatrixBeta, X::AbstractArray{<:RTM{<:Real}})
    return mapvcat(x -> logpdf(d, x), X)
end

Distributions.Wishart(df::RTR, S::Matrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::RTR, S::AbstractMatrix{<:Real}) = TuringWishart(df, S)
Distributions.Wishart(df::Real, S::RTM) = TuringWishart(df, S)
Distributions.Wishart(df::RTR, S::RTM) = TuringWishart(df, S)

Distributions.InverseWishart(df::RTR, S::Matrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::RTR, S::AbstractMatrix{<:Real}) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::Real, S::RTM) = TuringInverseWishart(df, S)
Distributions.InverseWishart(df::RTR, S::RTM) = TuringInverseWishart(df, S)

function Distributions.logpdf(d::Wishart, X::RTM)
    return logpdf(TuringWishart(d), X)
end
function Distributions.logpdf(d::Wishart, X::AbstractArray{<:RTM})
    return logpdf(TuringWishart(d), X)
end

function Distributions.logpdf(d::InverseWishart, X::RTM)
    return logpdf(TuringInverseWishart(d), X)
end
function Distributions.logpdf(d::InverseWishart, X::AbstractArray{<:RTM})
    return logpdf(TuringInverseWishart(d), X)
end

# Modified from Tracker.jl

Base.vcat(xs::RTM...) = _vcat(xs...)
Base.vcat(xs::RTV...) = _vcat(xs...)
function _vcat(xs::Union{RTV{<:Any, D}, RTM{<:Any, D}}...) where {D}
    tp = ReverseDiff.tape(xs...)
    xs_value = ReverseDiff.value.(xs)
    out_value = vcat(xs_value...)
    function back(Δ)
        start = 0
        Δs = [begin
          x = map(_ -> :, size(xsi))
          i = isempty(x) ? x : Base.tail(x)
          d = Δ[start+1:start+size(xsi,1), i...]
          start += size(xsi, 1)
          d
        end for xsi in xs]
        return (Δs...,)
    end
    out = ReverseDiff.track(out_value, D, tp)
    ReverseDiff.record!(tp, SpecialInstruction, vcat, xs, out, (back,))
    return out
end

@noinline function ReverseDiff.special_reverse_exec!(instruction::SpecialInstruction{typeof(vcat)})
    output = instruction.output
    input = instruction.input
    input_derivs = ReverseDiff.deriv.(input)
    P = instruction.cache[1]
    jtvs = P(ReverseDiff.deriv(output))
    for i in 1:length(input_derivs)
        input_derivs[i] .+= jtvs[i]
    end
    ReverseDiff.unseed!(output)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(instruction::SpecialInstruction{typeof(vcat)})
    output, input = instruction.output, instruction.input
    out_value = vcat(ReverseDiff.value.(input)...)
    ReverseDiff.value!(output, out_value)
    return nothing
end

Base.hcat(xs::RTM...) = _hcat(xs...)
Base.hcat(xs::RTV...) = _hcat(xs...)
function _hcat(xs::Union{RTV{<:Any, D}, RTM{<:Any, D}}...) where {D}
    tp = ReverseDiff.tape(xs...)
    xs_value = ReverseDiff.value.(xs)
    out_value = hcat(xs_value...)
    function back(Δ)
        start = 0
        Δs = [begin
          d = if ndims(xsi) == 1
            Δ[:, start+1]
          else
            i = map(_ -> :, size(xsi)) |> Base.tail |> Base.tail
            Δ[:, start+1:start+size(xsi,2), i...]
          end
          start += size(xsi, 2)
          d
        end for xsi in xs]
        return (Δs...,)
    end
    out = ReverseDiff.track(out_value, D, tp)
    ReverseDiff.record!(tp, SpecialInstruction, hcat, xs, out, (back,))
    return out
end

@noinline function ReverseDiff.special_reverse_exec!(instruction::SpecialInstruction{typeof(hcat)})
    output = instruction.output
    input = instruction.input
    input_derivs = ReverseDiff.deriv.(input)
    P = instruction.cache[1]
    jtvs = P(ReverseDiff.deriv(output))
    for i in 1:length(input_derivs)
        input_derivs[i] .+= jtvs[i]
    end
    ReverseDiff.unseed!(output)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(instruction::SpecialInstruction{typeof(hcat)})
    output, input = instruction.output, instruction.input
    out_value = hcat(ReverseDiff.value.(input)...)
    ReverseDiff.value!(output, out_value)
    return nothing
end

Base.cat(Xs::RTA...; dims) = _cat(dims, Xs...)
Base.cat(Xs::RTV...; dims) = _cat(dims, Xs...)
function _cat(dims, Xs::Union{RTV{<:Any, D}, RTM{<:Any, D}}...) where {D}
    tp = ReverseDiff.tape(dims, Xs...)
    Xs_value = ReverseDiff.value.(Xs)
    out_value = cat(Xs_value...; dims = dims)
    function back(Δ)
        start = ntuple(i -> 0, Val(ndims(Δ)))
        Δs = [begin
          dim_xs = 1:ndims(xs)
          till_xs = ntuple((i -> i in dims ? (i in dim_xs ? size(xs,i) : 1) : 0), Val(ndims(Δ)))
          xs_in_Δ = ntuple(i -> till_xs[i] > 0 ? (start[i]+1:start[i]+till_xs[i]) : Colon(), Val(ndims(Δ)))
          d = reshape(Δ[xs_in_Δ...],size(xs))
          start = start .+ till_xs
          d
        end for xs in Xs]
        return (Δs...,)
    end
    out = ReverseDiff.track(out_value, D, tp)
    ReverseDiff.record!(tp, SpecialInstruction, cat, (dims, Xs...), out, (back,))
    return out
end

@noinline function ReverseDiff.special_reverse_exec!(instruction::SpecialInstruction{typeof(cat)})
    output = instruction.output
    input = instruction.input
    input_derivs = ReverseDiff.deriv.(Base.tail(input))
    P = instruction.cache[1]
    jtvs = P(ReverseDiff.deriv(output))
    for i in 1:length(jtvs)
        input_derivs[i] .+= jtvs[i]
    end
    ReverseDiff.unseed!(output)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(instruction::SpecialInstruction{typeof(cat)})
    output, input = instruction.output, instruction.input
    dims = ReverseDiff.value(input[1])
    Xs = ReverseDiff.value.(Base.tail(input))
    out_value = cat(Xs..., dims = dims)
    ReverseDiff.value!(output, out_value)
    return nothing
end

###########

# Broadcasting

using ReverseDiff: ForwardOptimize
using Base.Broadcast: Broadcasted
import Base.Broadcast: materialize
const RDBroadcasted{F, T} = Broadcasted{<:Any, <:Any, F, T}

_materialize(f, args) = broadcast(ForwardOptimize(f), args...)

for (M, f, arity) in ReverseDiff.DiffRules.diffrules()
    if arity == 1
        @eval @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTA}}) = _materialize(bc.f, bc.args)
    elseif arity == 2
        @eval begin
            @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTA, RTA}}) = _materialize(bc.f, bc.args)
            @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTA, RTR}}) = _materialize(bc.f, bc.args)
            @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTR, RTA}}) = _materialize(bc.f, bc.args)
        end
        for A in ReverseDiff.ARRAY_TYPES
            @eval begin
                @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{$A, RTA}}) = _materialize(bc.f, bc.args)
                @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTA, $A}}) = _materialize(bc.f, bc.args)
                @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{$A, RTR}}) = _materialize(bc.f, bc.args)
                @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTR, $A}}) = _materialize(bc.f, bc.args)
            end
        end
        for R in ReverseDiff.REAL_TYPES
            @eval begin
                @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{$R, RTA}}) = _materialize(bc.f, bc.args)
                @inline materialize(bc::RDBroadcasted{typeof($M.$f), <:Tuple{RTA, $R}}) = _materialize(bc.f, bc.args)
            end
        end
    end
end


text:
include("presets/slg.jl")
include("presets/tblg.jl")
include("presets/geom2D_utils.jl")


text:
@plugin struct SystemImage <: Plugin
    path::String="deps/lib"
    incremental::Bool=false
    filter_stdlibs::Bool=true
    cpu_target::String="x86-64"
end

@plugin struct Comonicon <: Plugin
    name::Union{Nothing, String} = nothing
    # install
    completion::Bool=true
    quiet::Bool=false
    compile::Union{Nothing, String} = nothing
    optimize::Int=2
end

customizable(::Type{<:Comonicon}) = [
    :name=>String,
    :completion=>Bool,
    :quiet=>Bool,
    :compile=>String,
    :optimize=>Int
]

customizable(::Type{<:SystemImage}) = [
    :path=>String,
    :incremental=>Bool,
    :filter_stdlibs=>Bool,
    :cpu_target=>String,
]

# set up deps
function prehook(::Comonicon, ::Template, pkg_dir::AbstractString)
    if !ispath(joinpath(pkg_dir, "deps"))
        mkpath(joinpath(pkg_dir, "deps"))
    end
end

# create Comonicon.toml and build.jl
function hook(p::Comonicon, t::Template, pkg_dir::AbstractString)
    suffix = ".jl"
    sysimg = nothing
    for each in t.plugins
        if (each isa Git) && (each.jl == false)
            suffix = ""
        end

        if each isa SystemImage
            sysimg = each
        end
    end

    pkg = basename(pkg_dir)
    toml = OrderedDict{String, Any}(
        "name" => isnothing(p.name) ? default_name(pkg) : p.name
    )

    toml["install"] = OrderedDict{String, Any}(
        "completion" => p.completion,
        "quiet" => p.quiet,
        "optimize" => p.optimize,
    )

    if !isnothing(p.compile)
        toml["install"]["compile"] = p.compile
    end

    # system image related files/configs
    if !isnothing(sysimg)
        toml["sysimg"] = OrderedDict(
            "path" => sysimg.path,
            "incremental" => sysimg.incremental,
            "filter_stdlibs" => sysimg.filter_stdlibs,
            "cpu_target" => sysimg.cpu_target,
        )

        toml["download"] = OrderedDict(
            "host" => t.host,
            "user" => t.user,
            "repo" => pkg * suffix,
        )
    end # system image related files/configs

    open(joinpath(pkg_dir, "Comonicon.toml"), "w") do f
        TOML.print(f, toml)
    end

    open(joinpath(pkg_dir, "deps", "build.jl"); append=true) do f
        println(f, "using $pkg; $pkg.comonicon_install()")
    end
end

function hook(p::SystemImage, t::Template, pkg_dir::AbstractString)
    any(x->(x isa Comonicon), t.plugins) || error("SystemImage plugin must be used with Comonicon")
    workflow_dir = joinpath(pkg_dir, ".github", "workflows")
    mkpath(workflow_dir)
    cp(
        joinpath(
            templates("package", "github", "workflows", "sysimg.yml")
        ),
        joinpath(workflow_dir, "sysimg.yml")
    )
end

gitignore(::Comonicon) = ["/deps/build.log"]
gitignore(::SystemImage) = ["/deps/lib", "/deps/precompile.jl"]
needs_username(::Comonicon) = true


text:
@testset "basic" begin
    @testset "GNNChain" begin
        n, din, d, dout = 10, 3, 4, 2

        g = GNNGraph(random_regular_graph(n, 4),
                    graph_type=GRAPH_T,
                    ndata= randn(Float32, din, n))

        gnn = GNNChain(GCNConv(din => d),
                       BatchNorm(d),
                       x -> tanh.(x),
                       GraphConv(d => d, tanh),
                       Dropout(0.5),
                       Dense(d, dout))

        testmode!(gnn)

        test_layer(gnn, g, rtol=1e-5)

        @testset "Parallel" begin
            AddResidual(l) = Parallel(+, identity, l)

            gnn = GNNChain(ResGatedGraphConv(din => d, tanh),
                           BatchNorm(d),
                           AddResidual(ResGatedGraphConv(d => d, tanh)),
                           BatchNorm(d),
                           Dense(d, dout))

            testmode!(gnn)

            test_layer(gnn, g, rtol=1e-5)
        end
    end
end


text:
"""
```julia
h = EllipsoidAltitude(Body::T, PosRsph::Vector{Float64}) where {T <: abstractCelestialBody}
```
Computes the altitude above the reference ellipsoid with equatorial radius `Re`,
and a flattening `flat`.
"""

function EllipsoidAltitude(Body::T, PosRsph::Vector{Float64}) where {T <: abstractCelestialBody}
  return PosRsph[1] - equatorial_radius(Body)*(1-flattening(Body)*sin(PosRsph[3])^2)
end


text:
# CORK data types

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test,Random
using LinearAlgebra

@bench @testset "CORK types" begin

    # Construct a low-rank NEP
    A0=[1.0 3.0; -1.0 2.0];
    v=reshape([-1.0 ; 1]/sqrt(2),2,1);
    dep=DEP([A0,v*v'],[0,1]);

    # Form a CORKPencil
    cp_org=CORKPencil(dep,IarCorkLinearization(d=9))
    # .. and compress it
    cplr=lowRankCompress(cp_org,1,1);
    (AA_LR,BB_LR)=buildPencil(cplr);

    # Manually construct a pencil
    Av=[-A0-v*v']
    Bv=[-one(A0)-v*v']
    BvLR=[v/2, -v/3, v/4, -v/5, v/6, -v/7,  v/8, -v/9]
    AvLR=zero.(BvLR);
    Z=v;

    d=9;
    M=diagm( 0 =>  ones(d) )[2:end,:]
    N=diagm( -1 =>  1 ./ (1:d-1) )[2:end,:]
    cplr2=CORKPencilLR(M,N,Av,AvLR,Bv,BvLR,Z);

    (AA2,BB2)=buildPencil(cplr2);

    # ... and check that it's the same as we had before
    @test AA2 ≈ AA_LR
    @test BB2 ≈ BB_LR

end


text:
using Propagators, Test

@testset "Propagators.jl" begin
  include("basics.jl")
end


text:
function lintdict(ex::Expr, ctx::LintContext)
    typed = isexpr(ex.args[1], :curly)
    ks = Set{Any}()
    ktypes = Set{Any}()
    vtypes = Set{Any}()
    for a in ex.args[2:end]
        if ispairexpr(a)
            keyexpr = lexicalfirst(a)
            lit = lexicalvalue(keyexpr)
            if lit !== nothing
                if keyexpr in ks
                    msg(ctx, :E334, keyexpr, "duplicate key in Dict")
                end
                push!(ks, keyexpr)
            end
            for (j,s) in [(lexicalfirst,ktypes), (lexicallast,vtypes)]
                kvexpr = j(a)
                typeguess = lexicaltypeof(kvexpr)
                if isconcretetype(typeguess)
                    push!(s, typeguess)
                elseif isexpr(kvexpr, :call) && in(kvexpr.args[1], [:Date, :DateTime])
                    # TODO: use the existing guesstype infrastructure
                    # we want to add more immutable types such as Date, DateTime, etc.
                    push!(s, kvexpr.args[1])
                else
                    if typed
                        push!(s, Any)
                    end
                end
            end
        end

        lintexpr(a, ctx)
    end

    if typed
        if length(ktypes) > 1 && ex.args[1].args[2] != :Any && !isexpr(ex.args[1].args[2], :call)
            msg(ctx, :E531, "multiple key types detected. Use Dict{Any,V}() for mixed type dict")
        end
        if length(vtypes) > 1 && ex.args[1].args[3] != :Any && !isexpr(ex.args[1].args[3], :call)
            msg(ctx, :E532, "multiple value types detected. Use Dict{K,Any}() for mixed type dict")
        end
    end
end


text:
"""
    JModels.fit(t, x; settings=NamedTuple())

Return a fitted model of type `t` on `data`. Implementing this function is encouraged but
optional. Without implementing this function, things such as model evaluation where your
model is instantiated, `fit`ted and used for prediction multiple times are not possible.

It is advised to assign default values to all keyword arguments. This makes it easier for
people to compare different models.
"""
function fit(t, data; settings=NamedTuple())
    verify_model(t)
    error("JModels.fit is not implemented for $t")
end

"""
    JModels.fit!(model, data; settings=NamedTuple())

Fit an existing `model` on `data` by mutating `model`. In contrast to [`fit`](@ref), this
method is more flexible in configuring the model since a predefined model can be passed to
be `fit`ted. Also, this method can offer more performance if the model is trained in
multiple steps.
"""
function fit!(model, data; settings=NamedTuple())
    verify_model(model)
    error("JModels.fit! is not implemented for $typeof(model)")
end

"""
    JModels.predict(fmodel, data; settings=NamedTuple())

Predict with fitted model `fmodel` on `data`. For example, a k-means clustering model can predict target
labels.
"""
function predict(fmodel, data; settings=NamedTuple())
    verify_model(fmodel)
    error("JModels.predict is not implemented for $typeof(fmodel)")
end

"""
    JModels.transform(fmodel, data; settings=NamedTuple())

Transform `data` via fitted model `fmodel`. For example, a k-means clustering model reduce
dimensionality.
"""
function transform(fmodel, data; settings=NamedTuple())
    verify_model(fmodel)
    error("JModels.transform is not implemented for $typeof(fmodel)")
end

"""
    JModels.inverse_transform(fmodel, data; settings=NamedTuple())

Inversely transform `data` via fitted model `fmodel`.
"""
function inverse_transform(fmodel, data; settings=NamedTuple())
    verify_model(fmodel)
    error("JModels.inverse_transform is not implemented for $typeof(fmodel)")
end

"""
    JModels.ismodel(x) -> Bool

Check if an object `x` has defined that it is a statistical model and has implemented the
JModels interface.

# Example
```julia
JModels.ismodel(::ExampleModel) = true
```
"""
function ismodel end

ismodel(x) = false

"""
    JModels.verify_model(x)

Throw an error if `!ismodel(x)`.
"""
function verify_model(m)
    if !ismodel(m)
        T = typeof(m)
        error("JModels.ismodel(::$T) = true not implemented for $T which is either a bug or means that $T doesn't satisfy the JModels interface")
    end
end


text:
@testset "Solves" begin
# Test max iterations
solver = ALTROSolver(Problems.Pendulum()...)
set_options!(solver,
    iterations=100,
    cost_tolerance_intermediate=1e-10,
    gradient_tolerance_intermediate=1e-10,
    constraint_tolerance=1e-10,
    show_summary=false
)
solve!(solver)
status(solver)
@test iterations(solver) == solver.opts.iterations
@test status(solver) == Altro.MAX_ITERATIONS

## test max cost
model = RobotZoo.Quadrotor()
N = 101
tf = 5.0
x0,u0 = zeros(model)
xf = MArray(x0); xf[1] = 100
Q = Diagonal(@SVector ones(13))*1e6
R = Diagonal(@SVector ones(4))
obj = LQRObjective(Q,R,Q*N,xf,N,uf=u0)
prob = Problem(model, obj, x0, tf, xf=xf)
U0 = [u0*1e1 for k = 1:N-1]
initial_controls!(prob, U0)

solver = ALTROSolver(prob, show_summary=false)
Altro.is_constrained(solver)
solve!(solver)
@test iterations(solver) == 1
@test status(solver) == Altro.MAXIMUM_COST

## test max state value
model = RobotZoo.Quadrotor()
dmodel = RD.DiscretizedDynamics{RD.RK4}(model)
N = 101
tf = 5.0
dt = tf / (N-1)
x0,u0 = zeros(model)
xf = MArray(x0); xf[3] = 200
Q = Diagonal(SA[10,10,10, 1,1,1,1, 10,10,10, 10,10,10.])*1e-1 * dt
R = Diagonal(@SVector ones(4)) * dt
obj = LQRObjective(Q,R,Q*N,xf,N,uf=u0)
prob = Problem(dmodel, obj, x0, tf, xf=xf)
U0 = [u0*2 for k = 1:N-1]
initial_controls!(prob, U0)
rollout!(prob)
cost(prob)

solver = ALTROSolver(prob, show_summary=false)
cost(solver)
set_options!(solver, max_state_value=100)
@test Altro.is_constrained(solver) == false
solve!(solver)
@test status(solver) == Altro.MAX_ITERATIONS

## Max outer iterations
prob, opts = Problems.DubinsCar(:parallel_park, N=N)
opts.projected_newton = false
solver = ALTROSolver(prob, opts, show_summary=false)
solve!(solver)
@test max_violation(solver) > opts.constraint_tolerance
status(solver) == Altro.MAX_ITERATIONS_OUTER
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export allknn

"""
    allknn(index::Index, dist::Function; k::Int=1)

Finds k-nearest-neighbors for all items in the index's dataset; removes the object itself from its knn
"""
function allknn(index::Index, dist::Function; k::Int=1)
    n = length(index.db)
    A = [KnnResult(k+1) for i in 1:n]

    for i in 1:n
        x = index.db[i]
        res = A[i]
        search(index, dist, x, res)
        if first(res).objID == i
            popfirst!(res)
        end

        if (i % 1000) == 1
            println(stderr, "computing all-knn for index $(typeof(index)); k=$(k); advance $i of n=$n")
        end
    end

    A
end

"""
    allknn(index::Index, dist::Function, X::AbstractVector; k::Int=1)

Finds k-nearest-neighbors for all items in the given X
"""
function allknn(index::Index, dist::Function, X::AbstractVector; k::Int=1)
    n = length(X)
    A = [KnnResult(k) for i in 1:n]

    for i in 1:n
        x = X[i]
        res = A[i]
        search(index, dist, x, res)

        if (i % 1000) == 1
            println(stderr, "computing all-knn for index $(typeof(index)); k=$(k); advance $i of n=$n")
        end
    end

    A
end

"""
    allknn(index::SearchGraph{T}, dist::Function; k::Int=1)

Finds k-nearest-neighbors for all items in the index's dataset; removes the object itself from its knn.
This function is optimized for the SearchGraph similarity-index
"""
function allknn(index::SearchGraph, dist::Function; k::Int=1)
    n = length(index.db)
    A = [KnnResult(k+1) for i in 1:n]
    # TODO: Test caching the distance function
    for i in 1:n
        q = index.db[i]
        res = A[i]
        hints = [p.objID for p in res]
        append!(hints, index.links[i])
        empty!(res)
        search(index, dist, q, res, hints=hints)
        if first(res).objID == i
            popfirst!(res)
        end

        for p in res
            if p.objID > i
                push!(A[p.objID], i, p.dist)
            end
        end

        if (i % 1000) == 1
            println(stderr, "computing all-knn for $(index.search_algo); k=$(k); advance $i of n=$n")
        end
    end

    A
end

## """
## Finds k-nearest-neighbors for all items in the given dataset (represented as a SearchGraph object).
##
## This function is optimized for the SearchGraph similarity-index
## """
## function allknn(index::SearchGraph{T}, dist::Function, X::SearchGraph{T}; k::Int=1, bsize::Int=4) where T
##     n = length(X.db)
##     A = Vector{KnnResult}(undef, n)
##     # TODO: Test caching the distance function
##
##     L = [(id=n, hints=rand(1:n, bsize))]
##     while length(L) > 0
##         t = pop!(L)
##         search(index, dist, X.db[i], )
##     end
##     hints = Int32[]
##     for i in 1:n
##         q = X.db[i]
##
##
##         if A[i] == undef
##             A[i] = KnnResult(k)
##         else
##             continue
##         end
##
##         # hints = rand(1:n, bsize)
##         res = search(index, dist, q, A[i], hints=hints)
##         push!(L, (id=i, res))
##         clear!(hints)
##         append!(index.links[first(res).objID])
##         if (i % 10000) == 1
##             println(stderr, "computing all-knn for $(index.search_algo); k=$(k); advance $i of n=$n")
##         end
##     end
##
##     A
## end


text:
module Hydrogen_Vh
    include("hydrogen_vh_module.jl")
    include("hydrogen_fem_julia.jl")
    using FastGaussQuadrature
    using LinearAlgebra
    using Match
    using MKL
    using Printf
    using .Hydrogen_FEM
    using .Hydrogen_Vh_module

    function construct(hfem_param)
        vh_param = Hydrogen_Vh_module.Hydrogen_Vh_param(100, "result.csv")
        vh_val = Hydrogen_Vh_module.Hydrogen_Vh_variables(
            Array{Float64}(undef, hfem_param.ELE_TOTAL, 2, 2),
            SymTridiagonal(Array{Float64}(undef, hfem_param.NODE_TOTAL), Array{Float64}(undef, hfem_param.NODE_TOTAL - 1)),
            Array{Float64}(undef, hfem_param.NODE_TOTAL),
            Array{Float64}(undef, hfem_param.ELE_TOTAL, 2),
            Array{Float64}(undef, hfem_param.NODE_TOTAL),
            Array{Float64}(undef, vh_param.INTEGTABLENUM),
            Array{Float64}(undef, vh_param.INTEGTABLENUM))

        vh_val.x, vh_val.w = gausslegendre(vh_param.INTEGTABLENUM)

        return vh_param, vh_val
    end

    function do_run(hfem_param, hfem_val, vh_val)
        # 要素行列とLocal節点ベクトルを生成
        make_element_matrix_and_vector(hfem_param, hfem_val, vh_val)

        # 全体行列と全体ベクトルを生成
        tmp_dv, tmp_ev = make_global_matrix_and_vector(hfem_param, hfem_val, vh_val)

        # 境界条件処理
        boundary_conditions(vh_val, hfem_param, tmp_dv, tmp_ev)

        # 連立方程式を解く
        vh_val.ug = vh_val.mat_A_glo \ vh_val.vec_b_glo
    end

    save_result(hfem_val, vh_param, vh_val) = let
        open(vh_param.RESULT_FILENAME, "w" ) do fp
            for i = 2:length(hfem_val.node_r_glo)
                r = hfem_val.node_r_glo[i]
                println(fp, @sprintf "%.14f, %.14f, %.14f" (r) (vh_val.ug[i] / r) (- (1.0 + 1.0 / r) * exp(-2.0 * r) + 1.0 / r))
            end
        end
    end

    function boundary_conditions(vh_val, hfem_param, tmp_dv, tmp_ev)
        a = 0.0
        tmp_dv[1] = 1.0
        vh_val.vec_b_glo[1] = a
        vh_val.vec_b_glo[2] -= a * tmp_ev[1]
        tmp_ev[1] = 0.0

        b = 1.0
        tmp_dv[hfem_param.NODE_TOTAL] = 1.0
        vh_val.vec_b_glo[hfem_param.NODE_TOTAL] = b
        vh_val.vec_b_glo[hfem_param.NODE_TOTAL - 1] -= b * tmp_ev[hfem_param.NODE_TOTAL - 1]
        tmp_ev[hfem_param.NODE_TOTAL - 1] = 0.0

        vh_val.mat_A_glo = SymTridiagonal(tmp_dv, tmp_ev)
    end

    function gl_integ(f, x1, x2, vh_val)
        xm = 0.5 * (x1 + x2);
        xr = 0.5 * (x2 - x1);

        s = sum(i -> vh_val.w[i] * f(xm + xr * vh_val.x[i]), eachindex(vh_val.x))

        return s * xr;
    end

    function make_element_matrix_and_vector(hfem_param, hfem_val, vh_val)
        # 要素行列とLocal節点ベクトルの各成分を計算
        @inbounds for e = 1:hfem_param.ELE_TOTAL
            for i = 1:2
                for j = 1:2
                    vh_val.mat_A_ele[e, i, j] = (-1) ^ i * (-1) ^ j / hfem_val.length[e]
                end

                vh_val.vec_b_ele[e, i] =
                    @match i begin
                        1 => gl_integ(r -> r * Hydrogen_FEM.rho(hfem_param, hfem_val, r) * (hfem_val.node_r_ele[e, 2] - r) / hfem_val.length[e],
                                hfem_val.node_r_ele[e, 1],
                                hfem_val.node_r_ele[e, 2],
                                vh_val)

                        2 => gl_integ(r -> r * Hydrogen_FEM.rho(hfem_param, hfem_val, r) * (r - hfem_val.node_r_ele[e, 1]) / hfem_val.length[e],
                                hfem_val.node_r_ele[e, 1],
                                hfem_val.node_r_ele[e, 2],
                                vh_val)

                        _ => 0.0
                    end
            end
        end
    end

    function make_global_matrix_and_vector(hfem_param, hfem_val, vh_val)
        tmp_dv = zeros(hfem_param.NODE_TOTAL)
        tmp_ev = zeros(hfem_param.NODE_TOTAL - 1)

        # 全体行列と全体ベクトルを生成
        @inbounds for e = 1:hfem_param.ELE_TOTAL
            for i = 1:2
                for j = 1:2
                    if hfem_val.node_num_seg[e, i] == hfem_val.node_num_seg[e, j]
                        tmp_dv[hfem_val.node_num_seg[e, i]] += vh_val.mat_A_ele[e, i, j]
                    elseif hfem_val.node_num_seg[e, i] + 1 == hfem_val.node_num_seg[e, j]
                        tmp_ev[hfem_val.node_num_seg[e, i]] += vh_val.mat_A_ele[e, i, j]
                    end
                end

                vh_val.vec_b_glo[hfem_val.node_num_seg[e, i]] += vh_val.vec_b_ele[e, i]
            end
        end

        return tmp_dv, tmp_ev
    end
end


text:
get_type(v::NapiValue) = @napi_call napi_typeof(v::NapiValue)::NapiValueType

instanceof(a, b) = @with_scope try
    @napi_call napi_instanceof(a::NapiValue, b::NapiValue)::Bool
catch _
    false
end

is_undefined(nv::NapiValue) = get_type(nv) == NapiTypes.napi_undefined
is_function(nv::NapiValue) = get_type(nv) == NapiTypes.napi_function

macro wrap_is_type(name)
    func_name = Symbol(:is_, name)
    napi_name = Symbol(:napi_, func_name)
    @eval $func_name(v::NapiValue) = @napi_call $napi_name(v::NapiValue)::Bool
end

@wrap_is_type array
@wrap_is_type arraybuffer
@wrap_is_type buffer
@wrap_is_type typedarray
@wrap_is_type dataview
@wrap_is_type date
@wrap_is_type error
@wrap_is_type promise

@global_js_const _JS_MAP = "Map"
is_map(v::NapiValue) = instanceof(v, _JS_MAP)
@global_js_const _JS_SET = "Set"
is_set(v::NapiValue) = instanceof(v, _JS_SET)
@global_js_const _JS_ITERATOR_SYMBOL = "Symbol.iterator" false
is_iterator(v::NapiValue) = is_function(get(v, _JS_ITERATOR_SYMBOL, nothing; result=RESULT_RAW))

for func in (
    :get_type,
    :is_array, :is_arraybuffer, :is_buffer,
    :is_typedarray, :is_dataview,
    :is_date, :is_error,
    :is_promise,
    :is_map, :is_set, :is_iterator
)
    @eval $func(v) = @with_scope $func(napi_value(v))
end

macro wrap_coerce_convert(name)
    func_name = Symbol(:to_, name)
    napi_name = Symbol(:napi_coerce_, func_name)
    @eval $func_name(v::ValueTypes) = @napi_call $napi_name(v::NapiValue)::NapiValue
end

@wrap_coerce_convert bool
@wrap_coerce_convert number
@wrap_coerce_convert object
@wrap_coerce_convert string

Base.show(io::IO, v::JsObjectType) = print(io, string(
    typeof(v), ": ",
    string_pointer(getfield(v, :ref))
    # value(String, getfield(v, :ref))
))

function Base.get(
    o::ValueTypes, key, default=nothing;
    result=RESULT_VALUE
)
    with_result(result, this=o) do
        nv = @napi_call napi_get_property(o::NapiValue, key::NapiValue)::NapiValue
        is_undefined(nv) && return default
        nv
    end
end
set!(o::ValueTypes, key, value) = @with_scope begin
    @napi_call napi_set_property(o::NapiValue, key::NapiValue, value::NapiValue)
end
Base.haskey(o::ValueTypes, key) = @with_scope begin
    @napi_call napi_has_property(o::NapiValue, key::NapiValue)::Bool
end
Base.delete!(o::ValueTypes, key) = @with_scope begin
    @napi_call napi_delete_property(o::NapiValue, key::NapiValue)::Bool
end
Base.in(key, o::ValueTypes) = haskey(o, key)

Base.get(
    o::ValueTypes, key::AbstractString, default=nothing;
    result=RESULT_VALUE
) = with_result(result, this=o) do
    nv = @napi_call napi_get_named_property(o::NapiValue, key::Cstring)::NapiValue
    is_undefined(nv) && return default
    nv
end
set!(o::ValueTypes, key::AbstractString, value) = @with_scope begin
    @napi_call napi_set_named_property(o::NapiValue, key::Cstring, value::NapiValue)
end
Base.haskey(o::ValueTypes, key::AbstractString) = @with_scope begin
    @napi_call napi_has_named_property(o::NapiValue, key::Cstring)::Bool
end
Base.keys(o::ValueTypes) = @with_scope begin
    value(Array, @napi_call napi_get_property_names(o::NapiValue)::NapiValue)
end

Base.get(
    o::ValueTypes, key::Integer, default=nothing;
    result=RESULT_VALUE
) = with_result(result, this=o) do
    nv = @napi_call napi_get_element(o::NapiValue, key::UInt32)::NapiValue
    is_undefined(nv) && return default
    nv
end
set!(o::ValueTypes, key::Integer, value) = @with_scope begin
    @napi_call napi_set_element(o::NapiValue, key::UInt32, value::NapiValue)
end
Base.haskey(o::ValueTypes, key::Integer) = @with_scope begin
    @napi_call napi_has_element(o::NapiValue, key::UInt32)::Bool
end
Base.delete!(o::ValueTypes, key::Integer) = @with_scope begin
    @napi_call napi_delete_element(o::NapiValue, key::UInt32)::Bool
end

Base.getproperty(o::ValueTypes, key::Symbol) = get(o, string(key))
Base.setproperty!(o::ValueTypes, key::Symbol, value) = set!(o, string(key), value)
Base.hasproperty(o::ValueTypes, key::Symbol) = haskey(o, string(key))
Base.propertynames(o::ValueTypes) = Tuple(Symbol.(keys(o)))
Base.getindex(o::ValueTypes, key) = get(o, key)
Base.setindex!(o::ValueTypes, value, key) = set!(o, key, value)

_convert_key(::AbstractDict{T}, k::T) where T = k
_convert_key(::AbstractDict{String}, k) = string(k)
_convert_key(::AbstractDict{String}, k::String) = k
_convert_key(::AbstractDict{T}, k::AbstractString) where T <: Number = T(
    if isconcretetype(T)
        parse(T, k)
    else
        parse(Float64, k)
    end
)
_convert_key(::AbstractDict{T}, k) where T = convert(T, k)

@noinline jlnode_getproperty(o, k) = getproperty(o, Symbol(k))
@noinline jlnode_getproperty(o::AbstractDict, k) = getindex(o, _convert_key(o, k))
@noinline jlnode_setproperty!(o, k, v) = setproperty!(o, Symbol(k), v)
@noinline jlnode_setproperty!(o::AbstractDict, k, v) = setindex!(o, v, _convert_key(o, k))
@noinline jlnode_propertynames(o) = string.(propertynames(o))
@noinline jlnode_propertynames(o::AbstractDict) = keys(o)
@noinline jlnode_propertynames(o::Module) = names(o)
@noinline jlnode_hasproperty(o, k) = hasproperty(o, Symbol(k))
@noinline jlnode_hasproperty(o::AbstractDict, k) = try
    haskey(o, _convert_key(o, k))
catch
    false
end


text:
using Test
using MagneticReadHead
using Revise: includet
# Note: The patching used to be done with Mocking.jl
# However, that had some weird interactions between modules
# and includes
# So now the patching is done directly.
# Luckily we only have to patch some basic stuff intended only for this

"""
    make_readline_patch(text_queue)

Patches `MagneticReadHead.better_readline(io)`, such that each time it is
called it will return succecutive elements from the `text_queue`.
And will error if it is called more than it has canned responses for.
"""
function make_readline_patch(text_queue)
    @eval MagneticReadHead text_state=nothing
    @eval MagneticReadHead function better_readline(io)
        global text_state
        ret = (text_state==nothing ?
            iterate($text_queue) :
            iterate($text_queue, text_state)
        )
        ret == nothing && error("Out of programmed responses")
        text, text_state = ret
        println(text)
        sleep(0.1)
        return text
    end
end

"""
    make_recording_breakpoint_hit_patch()

Patchs the `breakpoint_hit(meth, statement_ind, variables)` method in MagneticReadHead
so that it will record what breakpoints were hit.
this function returns a reference to that record `Vector`.
"""
function make_recording_breakpoint_hit_patch()
    record = []
    @eval MagneticReadHead function breakpoint_hit(meth, statement_ind, variables)
        push!($record, (
            f=MagneticReadHead.functiontypeof(meth).instance,
            method=meth,
            statement_ind=statement_ind,
            variables=variables
        ))
    end
    return record
end

"""
    reset_patched_functions()

Patchs `breakpoint_hit` and `better_readline` back to their original state.
"""
function reset_patched_functions!()
    @eval MagneticReadHead function breakpoint_hit(meth, statement_ind)
        return nothing
    end

    @eval MagneticReadHead function better_readline(io)
        # This is not technically the same as the better_readline defined in utils
        # but we will only be doing automated tests so that doesn't matter.
        return readline(io)
    end
end

includet("demo.jl")


text:
struct Item
    url::String
    data
    assets
end

# https://github.com/radiantearth/stac-spec/blob/master/item-spec/item-spec.md

for (prop,name) in ((:id, "identifier"),
                    (:bbox, "bounding box"),
                    (:links, "links"),
                    (:properties, "properties"))
    @eval begin
        @doc """
     data = $($prop)(item)

Get the $($name) of STAC `item`.
        """
        $prop(item::Item) = item.data[$prop]
        export $prop
    end
end

"""
    dt = DateTime(item)

Get the date time of STAC `item` as a `Dates.DateTime` (or `nothing`
if this properties is not specified).
"""
function DateTime(item::Item)
    if haskey(item.data.properties,:datetime)
        dt = item.data.properties[:datetime]
        return CFTime.parseDT(Dates.DateTime,dt)
    else
        return nothing
    end
end
export DateTime

"""
    data = geometry(item)

Get the geometry of STAC `item` as a GeoJSON object
"""
geometry(item::Item) = GeoJSON.dict2geo(item.data[:geometry])
export geometry

function Item(url)
    data = cached_resolve(url)
    assets = _assets(data)
    return Item(url,data,assets)
end

function Base.show(io::IO,item::Item)
    fmt(x) = @sprintf("%0.6f",x)
    west, south, east, north = fmt.(bbox(item))
    printstyled(io, id(item), "\n", bold=true, color=title_color[])
    printstyled(io, "bounding box:\n")
    print(io,"""
     ┌──────$(north )───────┐
     │                      │
$(west)                $(east)
     │                      │
     └──────$(south )───────┘

""")

    _printstyled(io, "date time: ",DateTime(item),"\n")
    _show_assets(io,item)
end


text:
module Calcium

using DiffEqBase: ODEProblem, ODEFunction
using ForwardDiff
using Parameters: @with_kw, @unpack
using StaticArrays: SVector
using Setfield: @lens
import Setfield

using ...Bifurcations: BifurcationProblem
using ...Codim2: DiffEqCodim2Problem

@with_kw struct CalciumParam{
        vlType, vcaType, iType, glType, gcaType, cType, v1Type, v2Type}
    vl::vlType = -60
    vca::vcaType = 120
    i::iType = -220.0
    gl::glType = 2
    gca::gcaType = 4
    c::cType = 20
    v1::v1Type = -1.2
    v2::v2Type = 18
end

# http://www2.gsu.edu/~matrhc/Tutorial.html
# https://github.com/robclewley/pydstool/blob/master/examples/Tutorial_Calcium.py
# http://docs.juliadiffeq.org/latest/analysis/bifurcation.html
function f(u::SVector, p::CalciumParam, t)
    @unpack vl, vca, i, gl, gca, c, v1, v2 = p
    v = u[1]
    w = u[2]
    dv = (i + gl * (vl - v) - gca * 0.5 * (1 + tanh((v-v1)/v2)) * (v-vca)) / c
    dw = v-w
    return SVector(dv, dw)
end

function f(du, u, p, t)
    du .= f(SVector{2}(u), p, t)
    nothing
end

function jac(J, u, p, t)
    J .= ForwardDiff.jacobian(x -> f(x, p, t), SVector{2}(u))
    nothing
end

const _f = f
const _jac = jac

make_prob(
        p = CalciumParam();
        f = _f,
        u0 = SVector(-170.0, -170.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis = (@lens _.i),
        t_domain = (-300.0, 100.0),
        kwargs...) =
    BifurcationProblem(ode, param_axis, t_domain;
                       kwargs...)

prob = make_prob()
ode = prob.p.de_prob
u0 = ode.u0
param_axis = prob.p.param_axis

function make_codim2_prob(
        p = CalciumParam();
        u0 = SVector(-170.0, -170.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis1 = (@lens _.i),
        param_axis2 = (@lens _.gca),
        t_domain = ([-300.0, 0.0], [100.0, 8.0]),
        kwargs...)
    return DiffEqCodim2Problem(
        ode,
        param_axis1,
        param_axis2,
        t_domain;
        kwargs...)
end

make_ode_jac(
    p = CalciumParam();
    u0 = [-170.0, -170.0],
    tspan = (0.0, 30.0),
    f = _f,
    jac = _jac,
) =
    ODEProblem(ODEFunction{true}(f; jac = jac), u0, tspan, p)

end  # module


text:
using Libfreenect2
using Cxx
using CVCore
using CVHighGUI
using CVImgProc

isesc(key) = key == 27

f = Freenect2()
num_devices = enumerateDevices(f)
if num_devices <= 0
    error("No device!")
end
device = openDefaultDevice(f, OpenCLPacketPipeline())
listener = SyncMultiFrameListenerPtr()
setIrAndDepthFrameListener(device, listener)
setColorFrameListener(device, listener)

start(device)

# NOTE: must be called after start(device)
registration = Registration(getIrCameraParams(device),
    getColorCameraParams(device))
undistorted = FramePtr(512, 424, 4, key=Libfreenect2.FRAME_DEPTH)
registered = FramePtr(512, 424, 4, key=Libfreenect2.FRAME_COLOR)

try
    while true
        frames = waitForNewFrame(listener)
        color = frames[FrameType.COLOR]
        ir = frames[FrameType.IR]
        depth = frames[FrameType.DEPTH]

        apply(registration, color, depth, undistorted, registered)

        # Convert to Julia Array
        colorarr = convert(Array{UInt8,3}, color)
        irarr = convert(Array{Float32,2}, ir)
        deptharr = convert(Array{Float32,2}, depth)
        registeredarr = convert(Array{UInt8,3}, registered)
        undistortedarr = convert(Array{Float32,2}, undistorted)

        # Scale array to range [0,1]
        scale!(1/65535, irarr)
        scale!(1/4500, deptharr)

        # resize color image since it's a bit large to draw
        colormat = Mat(colorarr)
        colormat = resize(colormat,(Libfreenect2.height(color)/3, Libfreenect2.width(color)/3))

        # imshow("color", colormat)
        # imshow("ir", irarr)
        imshow("depth", deptharr)

        # imshow("registered", registeredarr)
        # imshow("unistored", undistortedarr)

        release(listener, frames)

        key = waitKey(delay=1)
        isesc(key) && break

        # TODO: remove this
        rand() > 0.98 && gc(false)
    end
finally
    stop(device)
    close(device)
    destroyAllWindows()
end


text:
# A work in progress - beware, dragons!
# Feynman-Kleinert-1986-PRA
#
# All codes lead to MAPI...
# Processing Feynman-Kleinert double-well method to apply to
#     https://github.com/jarvist/Julia-SoftModeTISH-DeformationPotential

push!(LOAD_PATH,"../src/") # load module from local directory

using FeynmanKleinert

# Variables
# Potential energy g; either anharmonic strength or double-well setup
g=0.1976
# Thermodynamic Beta, natch
β=1000

# Potential energy curves - V(x) functions
# Double well potential as given in figure 2 + 3 captions
# Hold my beer - I'm going to try metaprogramming
macro doublewell(g)
    return :(x -> -0.5*x^2 + 0.25*g*x^4 + 0.25/g)
end

# Anharmonic potential, following Kleinhert's book (p. 471); Also Figure 1 Feynman-Kleinert
macro anharmonicwell(g)
    return :(x^2/2+x^4*g/4)
end

# Harmonic oscillator - simple soln
macro harmonicwell(g)
    return :(0.5*x^2)
end

#V=@doublewell(g)
#println("@doublewell test. g=$g β=$β ")
#println("V(0.0) = ",V(0.0))
#println("Va2(0.0,1.0) = ",Va2(0.0,1.0,V))
#println("Wtilde(0.0,4,4,V,β=$β) = ",Wtilde(0.0,1.0,4,V,β))

using Plots

### OK - let's go for the big one!
# Try and reproduce Fig 3; Double-well for varying g.

xrange=-2:0.05:2 # For the plots; this to agree with FeynmanKleinertPRA
#Ws=[ W(x,g,β) for x in xrange ]
#println(Ws)
#plot(Ws,xrange,label="W(x)")

# This is same data as Fig 3 in FeynmanKleinertPRA

plot(size=(1200,600),fmt=:png) # force aspect ratio;
    # use PNG as changing size seems to break whatever the default is (SVG?)
for g in [.2, .3,.35,.4,.6] # range of g used in Fig 3. (look carefully at unbroken lines)
    plot!(x->W(x/sqrt(g),@doublewell(g),β), xrange, label="g=$g")
        # Note scaling of x-range within fn. argument
    plot!(x->@doublewell(g)(x/sqrt(g)), xrange, label="", linestyle = :dash)
        # Bare potential = infinite Temperature
end
plot!(ylim=(0,1.25)) # same as Fig 3 scale

savefig("MAPI.png")

g=.3 # detailed analysis + partition fn. calc on this g=
β=50

import QuadGK.quadgk

# This is horribly small as minimum in energy is offset ^^^ upwards.
Zq,Zerr=quadgk(x->exp(-β*W(x/sqrt(g),@doublewell(g),β))/sqrt(2*π*β), -2, 2)
println("Quantum Partition Function Z=$Zq")

Zqfudge,Zerr=quadgk(x->exp(-β*(W(x/sqrt(g),@doublewell(g),β)-0.7))/sqrt(2*π*β), -2, 2)
# Did I actually just do that?
println("Quantum Partition W-0.7 Function Z=$Zqfudge")

Zc,Zerr=quadgk(x->exp(-β*@doublewell(g)(x/sqrt(g)))/sqrt(2*π*β), -2, 2)
println("Classical Partition Function Z=$Zc")

# Effective classical potential, eval'd for x
U=x->W(x/sqrt(g),@doublewell(g),β)

classical=x->@doublewell(g)(x/sqrt(g))

plot(size=(1200,600),fmt=:png) # force aspect ratio;
plot!(x->W(x/sqrt(g),@doublewell(g),β), xrange,label="Effective Pot @ Beta=$β")
plot!(classical,xrange,label="Classical Pot")
plot!(ylim=(0,1.25)) # same as Fig 3 scale

savefig("MAPI-detailed.png")

xrange=-2:0.01:2 # For the plots; this to agree with FeynmanKleinertPRA
plot(size=(1200,600),fmt=:png) # force aspect ratio;
plot!(x->exp(-U(x)*β)/Zq,xrange,label="Quantum PDF")
plot!(x->exp(-classical(x)*β)/Zc,xrange,label="Classical PDF")

savefig("PDF.png")


text:
module cantilever_examples

using FinEtools
using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
using FinEtoolsDeforLinear: FEMMDeforLinear
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
using LinearAlgebra: norm
using Statistics: mean
using SparseArrays
using DelimitedFiles
using Interpolations
using UnicodePlots

function neohookean_h8()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    tmag = 0.2*phun("MPa");# Magnitude of the traction
    nincr = 8
    utol = 10e-7;
    graphics = ~true;
    maxdu_tol = W/1e7;
    maxiter = 9
    tolerance = W / 1000
    traction_vector = [0.0, 0.0, -tmag]

    fens, fes = H8block(L, W, H, 32, 9, 9)

    l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1)
    e2 = FDataDict("node_list"=>l1, "component"=>2)
    e3 = FDataDict("node_list"=>l1, "component"=>3)

    bfes = meshboundary(fes)
    el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
        v .= time .* traction_vector
        return v
    end
    fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
    t1 = FDataDict("femm"=>FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2))), "traction_vector"=>fi)

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

    # Rmout = fill(0.0, 3, 3)
    # rv = vec([-0.56 -0.1361 0.35])
    # rotmat3!(Rmout, rv)
    # mcsys = CSys(Rmout)
    # femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    region1 = FDataDict("femm"=>femm)
    modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1], "traction_bcs"=>[t1], "essential_bcs"=>[e1, e2, e3])

    modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
    modeldata["maxdu_tol"] = maxdu_tol;
    modeldata["maxiter"] = maxiter;
    modeldata["line_search"]  = true;
    modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
        @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
    end
    Ux = FFlt[]; lambdas = FFlt[]
    modeldata["increment_observer"] = (lambda, modeldata) -> begin
        push!(Ux, mean(modeldata["un1"].values[movel1,3]));
        push!(lambdas, lambda);
    end

    modeldata = nonlinearstatics(modeldata);
    @show Ux / phun("mm"), lambdas

    pl = lineplot(Ux / phun("mm"), lambdas)
    display(pl)

    vtkexportmesh("neohookean_h8.vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
end # function neohookean_h8

function neohookean_h20()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    tmag = 0.2*phun("MPa");# Magnitude of the traction
    nincr = 8
    utol = 10e-7;
    graphics = ~true;
    maxdu_tol = W/1e7;
    maxiter = 9
    tolerance = W / 1000
    traction_vector = [0.0, 0.0, -tmag]

    fens, fes = H20block(L, W, H, 6, 3, 3)

    l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>1)
    e2 = FDataDict("node_list"=>l1, "component"=>2)
    e3 = FDataDict("node_list"=>l1, "component"=>3)

    bfes = meshboundary(fes)
    el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
    setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
        v .= time .* traction_vector
        return v
    end
    fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
    t1 = FDataDict("femm"=>FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2))), "traction_vector"=>fi)

    movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

    # Rmout = fill(0.0, 3, 3)
    # rv = vec([-0.56 -0.1361 0.35])
    # rotmat3!(Rmout, rv)
    # mcsys = CSys(Rmout)
    # femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), mcsys, m)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    region1 = FDataDict("femm"=>femm)
    modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1], "traction_bcs"=>[t1], "essential_bcs"=>[e1, e2, e3])

    modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
    modeldata["maxdu_tol"] = maxdu_tol;
    modeldata["maxiter"] = maxiter;
    modeldata["line_search"]  = true;
    modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
        @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
    end
    Ux = FFlt[]; lambdas = FFlt[]
    modeldata["increment_observer"] = (lambda, modeldata) -> begin
        push!(Ux, mean(modeldata["un1"].values[movel1,3]));
        push!(lambdas, lambda);
    end

    modeldata = nonlinearstatics(modeldata);
    @show Ux / phun("mm"), lambdas

    pl = lineplot(Ux / phun("mm"), lambdas)
    display(pl)

    vtkexportmesh("neohookean_h20.vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
end # function neohookean_h8

function allrun()
    println("#####################################################")
    println("# neohookean_h8 ")
    neohookean_h8()
    println("#####################################################")
    println("# neohookean_h20 ")
    neohookean_h20()
    return true
end # function allrun

end # module cantilever_examples


text:
using Spec, Test

function myreverse end

@pre_spec function myreverse(x)
    # Test that our input has the right methods defined on it to be 'reversed'.
    @assert hasmethod(similar,   Tuple{typeof(x)})
    @assert hasmethod(eachindex, Tuple{typeof(x)})
    @assert hasmethod(getindex,  Tuple{typeof(x), Int})
end

@post_spec function myreverse(x::T) where {T}
    @assert __result__ isa T
    for i in eachindex(x, __result__)
        @assert x[i] == __result__[end-i+1]
    end
    @assert myreverse(__result__) == x
end

function myreverse(x)
    out = similar(x)
    for i ∈ eachindex(x, out)
        out[end-i+1] = x[i]
    end
    out
end

@test_throws MethodError myreverse("hi")
@test_throws AssertionError @validated myreverse("hi")

f(x) = x + 1
@post_spec f(x) = @assert isfinite(__result__)

@test isnan(f(NaN))
@test_throws AssertionError @validated f(NaN)


text:
# # Spatio-Temporal Prediction Examples
# In this page we are simply running files from the
# `examples` folder of the `TimeseriesPrediction` package.

# This is how you can (programmatically) find this folder:
# ```julia
# using TimeseriesPrediction
# exdir = dirname(dirname(pathof(TimeseriesPrediction)))*"/examples"
# ```

# ## Temporal Prediction: Kuramoto-Sivashinsky
# *(this requires `FFTW` to be installed)*

include("1Dfield_temporalprediction.jl")

# ## Cross Prediction: Barkley Model

include("2Dfield_crossprediction.jl")

# ## Temporal Prediction: Periodic Nonlinear Barkley Model

include("2Dfield_temporalprediction.jl")


text:
#!/usr/bin/env julia

# println("\nDANS abbrevs_dict.jl")

# Inspiré du gem abbrev de Ruby
#   http://ruby-doc.org/stdlib-2.5.0/libdoc/abbrev/rdoc/Abbrev.html
#   https://github.com/ruby/ruby/blob/trunk/lib/abbrev.rb
#
# Given a set of strings, calculate the set of unambiguous abbreviations for
# those strings, and return a hash where the keys are all the possible
# abbreviations and the values are the full strings.
#
# Thus, given +words+ is "car" and "cone", the keys pointing to "car" would
# be "ca" and "car", while those pointing to "cone" would be "co", "con", and
# "cone".
#
#   abbrevsdict(["car", "cone"])
#   #=> {"ca"=>"car", "con"=>"cone", "co"=>"cone", "car"=>"car", "cone"=>"cone"}
#
# The optional +pattern+ parameter is a pattern or a string. Only input
# strings that match the pattern or start with the string are included in the
# output hash.
#
#   abbrevsdict(%w{car box cone crab}, /b/)
#   #=> {"box"=>"box", "bo"=>"box", "b"=>"box", "crab" => "crab"}
#
#   abbrevsdict(%w{car box cone}, 'ca')
#   #=> {"car"=>"car", "ca"=>"car"}
#
# TODO: en faire un Pkg Julia officiel
# TODO: accepter une collection de String au sens large
#       (e.g une Base.KeyIterator)
# TODO: créer un module AbbrevsDict et methode AbbrevsDict.abbrevsdict()
#
# HIST
# - xx/06/2018 création par améliorer ArgParse
# - 15/10/2018 table devient table contient des String plutot que des Symbol
#               (pour julia07)
function abbrevsdict(words::Array{String,1}, pat=r".*")
    if isa(pat, String)
        pat = Regex(pat)
    end
    # table = Dict{String, Symbol}()
    table = Dict{String, String}()
    seen = Dict{String, Int}()
    for word in words
        length(word)==0 && next
        for len in length(word):-1:1
            abbrev = word[1:len]
            # if !ismatch(pat, abbrev)
            #     continue
            # end
            if !occursin(pat, abbrev)
                continue
            end
            # Without Base.xxx, could raise error if get is redefined somewhere
            seen[abbrev] = 1 + Base.get(seen, abbrev, 0)
            if seen[abbrev] == 1
                # table[abbrev] = Symbol(word)
                table[abbrev] = word
            elseif seen[abbrev] == 2
                delete!(table, abbrev)
            else
                break
            end
        end
    end
    for word in words
        if !occursin(pat, word)
            continue
        end
        # table[word] = Symbol(word)
        table[word] = word
    end
    table
end
# function abbrevsdict(words::Array{String,1}, pat=r".*")
#     if isa(pat, String)
#         pat = Regex(pat)
#     end
#     table = Dict{String, Symbol}()
#     seen = Dict{String, Int}()
#     for word in words
#         length(word)==0 && next
#         for len in length(word):-1:1
#             abbrev = word[1:len]
#             @show abbrev, typeof(abbrev)
#             if !occursin(pat, abbrev)
#                 continue
#             end
#             # Without Base.xxx, could raise error if get is redefined somewhere
#             seen[abbrev] = 1 + Base.get(seen, abbrev, 0)
#             if seen[abbrev] == 1
#                 table[abbrev] = Symbol(word)
#             elseif seen[abbrev] == 2
#                 delete!(table, abbrev)
#             else
#                 break
#             end
#         end
#     end
#     @show words
#     # for word in words
#     #     # if !ismatch(pat, word)
#     #     #     continue
#     #     # end
#     #     if !occursin(pat, word)
#     #         continue
#     #     end
#     #     table[word] = word
#     # end
#     @show table
#     table
# end

function test_abbrevsdict()
    @show typeof(ARGS)
    @show ARGS
    @show abbrevsdict(ARGS)
    @show abbrevsdict(["car", "cone"])
    @show abbrevsdict(["car", "cone"], "e")
    @show abbrevsdict(["car", "cone"], r"[q-z]")
end

if basename(@__FILE__) == basename(PROGRAM_FILE)
    test_abbrevsdict()
    # main()
end


text:
struct IsometricCamera
  Position::space
  Direction::space
  Sensorx::space
  Sensory::space
  Maxx::Float64
  Maxy::Float64
  function IsometricCamera(Position::space, Direction::space, Sensorx::space, Sensory::space, Maxx::Int64, Maxy::Int64)
    new(Position, normalize(Direction), Sensorx/Maxx, Sensory/Maxy, Maxx/2.0, Maxy/2.0)
  end
end

function los(Cam::IsometricCamera,X::Int64,Y::Int64)::Tuple{space,space}
  ((X-Cam.Maxx)*Cam.Sensorx+(Y-Cam.Maxy)*Cam.Sensory+Cam.Position,Cam.Direction)
end

struct PinholeCamera
  Position::space
  Direction::space
  Sensorx::space
  Sensory::space
  Maxx::Float64
  Maxy::Float64
  FocalLength::Float64
  function PinholeCamera(Position::space, Direction::space, Sensorx::space, Sensory::space, Maxx::Int64, Maxy::Int64, Fl::Float64)
    new(Position, normalize(Direction), Sensorx/Maxx, Sensory/Maxy, Maxx/2.0, Maxy/2.0,Fl)
  end
end # struct

function los(Cam::PinholeCamera,X::Int64,Y::Int64)::Tuple{space,space}
  light_passes_lens = (X-Cam.Maxx)*Cam.Sensorx+(Y-Cam.Maxy)*Cam.Sensory+Cam.Position
  (light_passes_lens,normalize(light_passes_lens-(Cam.Position-Cam.FocalLength*Cam.Direction)))
end

x, y = (1080,1920)
#Höhe Breite
#Heigth Width
isocam = IsometricCamera(space(25.0,0.0,0.0),
                         space(-1.0,0,0),
                         x/600*35*normalize(space(0.0,1.0,0.0)),
                         y/600*35*normalize(space(0.0,0.0,1.0)),x,y);

pincam = PinholeCamera(space(25.0,0.0,0.0),
                         space(-1.0,0,0),
                         x/600*35*normalize(space(0.0,1.0,0.0)),
                         y/600*35*normalize(space(0.0,0.0,1.0)),x,y,90.0);
#frame = Array{RayResults}(undef,x,y);
myPic = Array{RGB}(undef,x,y);


text:
#####
##### Divergence operators
#####

"""
    divᶜᶜᶜ(i, j, k, grid, u, v, w)

Calculates the divergence ∇·𝐔 of a vector field 𝐔 = (u, v, w),

    1/V * [δxᶜᵃᵃ(Ax * u) + δxᵃᶜᵃ(Ay * v) + δzᵃᵃᶜ(Az * w)],

which will end up at the cell centers `ccc`.
"""
@inline function divᶜᶜᶜ(i, j, k, grid, u, v, w)
    return 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, u) +
                                      δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, v) +
                                      δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶜᶠ, w))
end

"""
    div_xyᶜᶜᵃ(i, j, k, grid, u, v)

Returns the discrete `div_xy = ∂x u + ∂y v` of velocity field `u, v` defined as

```
1 / Azᶜᶜᵃ * [δxᶜᵃᵃ(Δyᵃᶜᵃ * u) + δyᵃᶜᵃ(Δxᶜᵃᵃ * v)]
```

at `i, j, k`, where `Azᶜᶜᵃ` is the area of the cell centered on (Center, Center, Any) --- a tracer cell,
`Δy` is the length of the cell centered on (Face, Center, Any) in `y` (a `u` cell),
and `Δx` is the length of the cell centered on (Center, Face, Any) in `x` (a `v` cell).
`div_xyᶜᶜᵃ` ends up at the location `cca`.
"""
@inline function div_xyᶜᶜᶜ(i, j, k, grid, u, v)
    return 1 / Azᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u) +
                                       δyᵃᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v))
end

# Convention
 index_left(i, ::Center) = i
 index_left(i, ::Face)   = i - 1
index_right(i, ::Center) = i + 1
index_right(i, ::Face)   = i

@inline Base.div(i, j, k, grid::AbstractGrid, loc, q_west, q_east, q_south, q_north, q_bottom, q_top) =
    1 / volume(i, j, k, grid, loc...) * (δx_Ax_q(i, j, k, grid, loc, q_west, q_east) +
                                         δy_Ay_q(i, j, k, grid, loc, q_south, q_north) +
                                         δz_Az_q(i, j, k, grid, loc, q_bottom, q_top))

@inline function δx_Ax_q(i, j, k, grid, (LX, LY, LZ), qᵂ, qᴱ)
    iᵂ = index_left(i, LX)
    Axᵂ = Ax(iᵂ, j, k, grid, LX, LY, LZ)

    iᴱ = index_right(i, LX)
    Axᴱ = Ax(iᴱ, j, k, grid, LX, LY, LZ)

    return Axᴱ * qᴱ - Axᵂ * qᵂ
end

@inline function δy_Ay_q(i, j, k, grid, (LX, LY, LZ), qˢ, qᴺ)
    jˢ = index_left(j, LY)
    Ayˢ = Ay(i, jˢ, k, grid, LX, LY, LZ)

    jᴺ = index_right(j, LY)
    Ayᴺ = Ay(i, jᴺ, k, grid, LX, LY, LZ)

    return Ayᴺ * qᴺ - Ayˢ * qˢ
end

@inline function δz_Az_q(i, j, k, grid, (LX, LY, LZ), qᴮ, qᵀ)
    kᴮ = index_left(k, LZ)
    Azᴮ = Az(i, j, kᴮ, grid, LX, LY, LZ)

    kᵀ = index_right(k, LZ)
    Azᵀ = Az(i, j, kᵀ, grid, LX, LY, LZ)

    return Azᵀ * qᵀ - Azᴮ * qᴮ
end

# And flat!

@inline δx_Ax_q(i, j, k, grid::XFlatGrid, args...) = zero(grid)
@inline δy_Ay_q(i, j, k, grid::YFlatGrid, args...) = zero(grid)
@inline δz_Az_q(i, j, k, grid::ZFlatGrid, args...) = zero(grid)


text:
export temp_filename, temp_file, temp_dir, cleanup_path

@fixture function temp_filename(;extension::Union(String, Nothing)=nothing, create::Union(String,Bool)=false)
  const filename = "$(tempname())$(extension!=nothing ? ("."*extension) : "")"
  if create==false
    # Do nothing
  elseif create==true
    touch(filename)
  else
    open(filename, "w") do f
      write(f, create)
    end
  end
  yield_fixture(filename)

  if isdir(filename)
    run(`rm -fr $filename`)
  elseif ispath(filename)
    rm(filename)
  end
end

function temp_file(fn::Function; extension::Union(String, Nothing)=nothing, content::Union(String,Nothing)=nothing, mode::String="w")
  if content!=nothing
    return temp_filename((filename) -> open(fn, filename, mode), extension=extension, create=content)
  else
    return temp_filename((filename) -> open(fn, filename, mode), extension=extension)
  end
end

temp_dir(fn::Function) = temp_filename() do filename
  mkdir(filename)
  fn(filename)
end

@fixture function cleanup_path(path::String; ignore_missing=false)
  yield_fixture()
  if isdir(path)
    run(`rm -rf $path`)
  elseif ignore_missing==false || ispath(path)
    rm(path)
  end
end


text:
using Libdl, Base.Sys

depsfile = joinpath(dirname(@__FILE__), "deps.jl")

if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(libpath)
    open(depsfile,"w") do f
        print(f,"const libhiop = ")
        show(f, libpath)
        println(f)
    end
end

# Name of shared library
libhiop_name = string("libhiop", ".", Libdl.dlext)
paths_to_try = String[]

# First, look into JULIA_HIOP_LIBRARY_PATH
if haskey(ENV, "JULIA_HIOP_LIBRARY_PATH")
    push!(paths_to_try, joinpath(ENV["JULIA_HIOP_LIBRARY_PATH"], libhiop_name))
end
# And look also in LD_LIBRARY_PATH
push!(paths_to_try, libhiop_name)

global found_hiop = false
for path in paths_to_try
    println(path)
    d = Libdl.dlopen_e(path)
    if d != C_NULL
        global found_hiop = true
        # Store pointer to Hiop
        write_depsfile(path)
        break
    end
end

if !found_hiop
    error("Could not load HiOp shared library. ",
            "Make sure it is in your LD_LIBRARY_PATH.")
end


text:
export trystatic, tostatic, ⩢, staticmod

# TODO: tryconvert(T, x) should return convert(T, x) unless that would throw
# an error. tryconvert(StaticInteger{0}, x) would be equivalent to the current
# trystatic(x, static(0))

"""
trystatic(x, y1, y2, ...)
x ⩢ y1 ⩢ y2 ...

Test if a number `x` is equal to any of the numbers `y1`, `y2`, ..., and in that
case return `static(y)`. Otherwise, or if `x` is already a `Static` number, `x is
returned unchanged.

The inferred return type will typically be a small `Union`, which Julia
can handle efficiently.

This function can be used to call specialized methods for certain input values.
For example, `f(x, y ⩢ 0)` will call `f(x, y)` if `y` is nonzero, but
`f(x, static(0))` if y is zero. This is useful if it enables optimizations that
outweigh the cost of branching.

NOTE: When the list of y-values is longer than one, y1, y2, ... must be `Static`
numbers, or inferrence will not work. (In which case `trystatic` is not more
efficient than `static(x)`.)
"""
@inline function trystatic(x::Number, y::Number)
    if isnan(y)
        isnan(x) ? static(NaN) : x
    else
        x==y ? static(y) : x
    end
end

@inline trystatic(x::Number, y::Number, ys::Number...) = trystatic(trystatic(x, y), ys...)

@inline ⩢(x, y) = trystatic(x,y)

@inline trystatic(x::Number, t::T) where {T<:Tuple} = trystatic(x, static.(t)...)

#The numbers `y1`, `y2`, ..., or the range `r` should be such that they can be
#computed at inference. I.e. they should be constructed using literals, `Static`
#numbers, and other constants that are deducible from types types.

"""
trystatic(x, r)
Tests if an integer `x` is in the range `r`, and if so, returns a `Static`
integer from the range. (Otherwise, `x` is returned unchanged.)

NOTE: The range must be completely static, or inferrence will not work.
"""
@inline trystatic(x) = x
@inline trystatic(x::StaticInteger, r::OrdinalRange{<:Integer, <:Integer}) = x
@inline trystatic(x::Integer, r::OrdinalRange{<:Integer, <:Integer}) =
    trystatic(x::Integer, LengthStepRange(static(zeroth(r)), static(step(r)), static(length(r))))
@inline trystatic(x::Integer, r::LengthStepRange{<:Integer, <:StaticInteger, <:StaticInteger}) = x in r ? tostatic(x, r) : x
@inline trystatic(x::StaticInteger, r::LengthStepRange{<:Integer, <:StaticInteger, <:StaticInteger}) = x # disambig
@inline trystatic(x::Integer, r::LengthUnitRange{<:Integer, <:StaticInteger, <:StaticInteger}) = x in r ? tostatic(x, r) : x
@inline trystatic(x::StaticInteger, r::LengthUnitRange{<:Integer, <:StaticInteger, <:StaticInteger}) = x # disambig

#@inline tostatic(x::Integer, r::StepRange) = tostatic(x, LengthStepRange(static(zeroth(r)), static(step(r)), static(lenght(r)))
#@inline tostatic(x::Integer, r::UnitRange) = tostatic(x, LengthUnitRange(static(zeroth(r)), static(lenght(r)))

"""
tostatic(x, r)
Returns a `Static` integer, equal to `x` from the range `r`. If no element in
`r` is equal to `x`, then the behaviour of this function is undefined.
"""
@generated function tostatic(x::Integer, r::LengthStepRange{<:Integer, StaticInteger{Z}, StaticInteger{S}, StaticInteger{L}}) where {Z, S, L}
    if x == StaticInteger
        quote
            x
        end
    else
        quote
            Base.@_inline_meta
            $(tostaticexpr(Z, S, L))
        end
    end
end
@generated function tostatic(x::Integer, r::LengthUnitRange{<:Integer, StaticInteger{Z}, StaticInteger{L}}) where {Z, L}
    if x == StaticInteger
        quote
            x
        end
    else
        quote
            Base.@_inline_meta
            $(tostaticexpr(Z, 1, L))
        end
    end
end

# Simpler variation:
#@inline tostatic(x::Static, r) = x
#@inline tostatic(x::Integer, r) = ntuple(i->static(r[i]), length(r))[x-first(r)

function tostaticexpr(z, s, l, blk=identity, x=:x)
    if l<=1
        blk(:( StaticInteger{$(z+s)}() ))
    else
        mid = l÷2
        :( $(s>0 ? :(<=) : :(>=))($x, $(z + s*mid)) ? $(tostaticexpr(z, s, mid, blk, x)) : $(tostaticexpr(z+mid*s, s, l-mid, blk, x)) )
    end
end

# do-constructs are probably not worth the effort
# @generated function tostatic(f, x::Integer, r::LengthUnitRange{<:Integer, StaticInteger{Z}, StaticInteger{L}}) where {Z, L}
#     quote
#         Base.@_inline_meta
#         $(tostaticexpr(Z, 1, L, y->:(f($y))))
#     end
# end
# @generated function tostatic(f, x::Integer, r::LengthStepRange{<:Integer, StaticInteger{Z}, StaticInteger{S}, StaticInteger{L}}) where {Z, S, L}
#     quote
#         Base.@_inline_meta
#         $(tostaticexpr(Z, S, L, y->:(f($y))))
#     end
# end

# TODO: Make cleaner macro interface!
macro tostatic(x::Symbol, start::Int, stop::Int, blk=x)
    xx = esc(x)
    xblk = esc(blk)
    quote
        if $xx isa Static
            $xblk
        else
            $(StaticNumbers.tostaticexpr(start-1, 1, stop-start+1, y->:( let $xx=$y; $xblk; end), xx))
        end
    end
end

"""
`staticmod(x,y)` returns `static(mod(x,y))` if `y` is `Static`.

If `y` is not `Static`, then `staticmod(x,y)` is the same as `mod(x,y)`.
"""
staticmod(x, y) = mod(x,y)
@inline staticmod(x::Integer, y::StaticInteger) = tostatic(mod(x,y), LengthUnitRange(static(-1), y))

# do-constructs. Remove?
# staticmod(f, x, y) = f(mod(x,y))
# @inline staticmod(f, x::Integer, y::StaticInteger) = tostatic(f, mod(x,y), LengthUnitRange(static(-1), y))


text:
module GameOfLife
export simulate, WorldState

include("World.jl")

using Plots, Images
gr()

function simulate(;max_days::Int = 100, initial_state::Array = GameOfLife.default_state(), speed::Int = 90)
    if !(0 < speed < 101)
        throw("speed must be between 1 to 100")
    end

    world = WorldState(max_days, initial_state)
    for (days, generation) in zip(world, 1:max_days)
        p = plot(Gray.(world.current_state .== 0), framestyle = :none)
        title!(string("Generation ", generation))
        display(p)
        sleep(1 - (speed * 0.01))
    end
end

end

max_stars_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_issues_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_forks_repo_name: VasanthManiVasi/Conways-Game-of-Life


text:
module TestCopulas

using Copulas
using Base.Test

## test missing functions
##-----------------------

u1 = [0.4]
u2 = [0.5]

cop = Copulas.ClaytonPC(3.0)
@test_throws Exception Copulas.pdf(cop, u1, u2)
@test_throws Exception Copulas.cdf(cop, u1, u2)
@test_throws Exception Copulas.hfun(cop, u1, u2)
@test_throws Exception Copulas.vfun(cop, u1, u2)
@test_throws Exception Copulas.hinv(cop, u1, u2)
@test_throws Exception Copulas.vinv(cop, u1, u2)

## test extensions to Float64
##---------------------------

u1 = [0.4, 0.5]
u2 = [0.5, 0.7]

cop = Copulas.ClaytonPC_Cpp([3.0])

expOut = Copulas.pdf(cop, [u1[1], u1[1]], u2)
@test Copulas.pdf(cop, u1[1], u2) == expOut
expOut = Copulas.pdf(cop, u1, [u2[1], u2[1]])
@test Copulas.pdf(cop, u1, u2[1]) == expOut
expOut = Copulas.pdf(cop, [u1[1]], [u2[1]])
@test Copulas.pdf(cop, u1[1], u2[1]) == expOut

expOut = Copulas.cdf(cop, [u1[1], u1[1]], u2)
@test Copulas.cdf(cop, u1[1], u2) == expOut
expOut = Copulas.cdf(cop, u1, [u2[1], u2[1]])
@test Copulas.cdf(cop, u1, u2[1]) == expOut
expOut = Copulas.cdf(cop, [u1[1]], [u2[1]])
@test Copulas.cdf(cop, u1[1], u2[1]) == expOut

expOut = Copulas.hfun(cop, [u1[1], u1[1]], u2)
@test Copulas.hfun(cop, u1[1], u2) == expOut
expOut = Copulas.hfun(cop, u1, [u2[1], u2[1]])
@test Copulas.hfun(cop, u1, u2[1]) == expOut
expOut = Copulas.hfun(cop, [u1[1]], [u2[1]])
@test Copulas.hfun(cop, u1[1], u2[1]) == expOut

expOut = Copulas.vfun(cop, [u1[1], u1[1]], u2)
@test Copulas.vfun(cop, u1[1], u2) == expOut
expOut = Copulas.vfun(cop, u1, [u2[1], u2[1]])
@test Copulas.vfun(cop, u1, u2[1]) == expOut
expOut = Copulas.vfun(cop, [u1[1]], [u2[1]])
@test Copulas.vfun(cop, u1[1], u2[1]) == expOut

expOut = Copulas.hinv(cop, [u1[1], u1[1]], u2)
@test Copulas.hinv(cop, u1[1], u2) == expOut
expOut = Copulas.hinv(cop, u1, [u2[1], u2[1]])
@test Copulas.hinv(cop, u1, u2[1]) == expOut
expOut = Copulas.hinv(cop, [u1[1]], [u2[1]])
@test Copulas.hinv(cop, u1[1], u2[1]) == expOut

expOut = Copulas.vinv(cop, [u1[1], u1[1]], u2)
@test Copulas.vinv(cop, u1[1], u2) == expOut
expOut = Copulas.vinv(cop, u1, [u2[1], u2[1]])
@test Copulas.vinv(cop, u1, u2[1]) == expOut
expOut = Copulas.vinv(cop, [u1[1]], [u2[1]])
@test Copulas.vinv(cop, u1[1], u2[1]) == expOut

end


text:
using Plots
include_ddp()

# Model
include_model("quadrotor")

function fd(model::Quadrotor{Midpoint, FixedTime}, x, u, w, h, t)
	return view(x, 1:model.n) + h * f(model, view(x, 1:model.n) + 0.5 * h * f(model, view(x, 1:model.n), view(u, 1:model.m), w), view(u, 1:model.m), w)
end

n = model.n
m = model.m

# Time
T = 101
h = 0.05

# Initial conditions, controls, disturbances
x1 = zeros(model.n)
x1[3] = 1.0

xT = copy(x1)
# xT[1] = 5.0
# xT[2] = 0.0
xT[3] = 5.0

u_ref = -1.0 * model.mass * model.g[3] / 4.0 * ones(model.m)
ū = [u_ref + 0.001 * randn(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)
# x̄ = linear_interpolation(x1, xT, T)
# plot(hcat(x̄...)')

# Objective
Q = h * [(t < T ? Diagonal([1.0 * ones(6); 1.0 * ones(6)])
        : Diagonal([0.0 * ones(6); 0.0 * ones(6)])) for t = 1:T]
q = h * [-2.0 * Q[t] * xT for t = 1:T]

R = h * [Diagonal(1.0 * ones(model.m)) for t = 1:T-1]
r = h * [-2.0 * R[t] * u_ref * 0.0 for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => zeros(model.m), :uu => 10.0 * ones(model.m), :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		c .= [ul - u; u - uu]
	elseif t == T
		xT = cons.con[T].info[:xT]
		c .= x - xT
	else
		nothing
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	con_tol = 1.0e-3,
	ρ_init = 1.0, ρ_scale = 10.0)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# Trajectories
plot(hcat(u...)', linetype = :steppost)
plot(hcat(x...)[1:3, :]', linetype = :steppost)

# Visualize
# include(joinpath(pwd(), "models/visualize.jl"))
# vis = Visualizer()
# render(vis)
# # open(vis)
# visualize!(vis, model, x̄, Δt = h)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using CairoMakie, AbstractPlotting, GeoMakie

lons = LinRange(-179.5, 179.5, 360)
lats = LinRange(-89.5, 89.5, 180)

field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

source = LonLat()
dest = WinkelTripel()

points, faces = GeoMakie.triangulated_grid(lons, lats)

tpoints = Proj4.transform.(source, dest, points)

scene = mesh(
    tpoints, faces;
    color = GeoMakie.img2colorvec(field),
    shading = false,
    show_axis = false, scale_plot = false, # don't show the normal axis and
                                           # ensure the aspect ratio is correct
    resolution = (1200, 600)               # just a nice resolution for this projection
)

geoaxis!(scene, extrema(lons), extrema(lats); crs = (src = source, dest = dest,))

coastlines!(scene; crs = (src = source, dest = dest))

save("geomakie-readme.png", scene) # works
save("geomakie-readme.svg", scene) # works
save("geomakie-readme.pdf", scene) # works but is really large