text:
using ReachabilityAnalysis, StructuralDynamicsODESolvers

function load_heat1d(; file)
    vars = matread(file)

    K = vars["K"]
    C = vars["C"]

    return K, C
end

function initial_distribution(n=100, tol=0.1)
    xs = range(0, 1, length=n+2)[2:end-1]
    c0 = @. sin(pi*xs) + 0.5*sin(3*pi*xs)
    r0 = tol * c0
    X0 = Hyperrectangle(c0, r0)
    return X0
end

function solve_heat1d(; kwargs...)
    prob, alg, NSTEPS = _solve_heat1d(; kwargs...)
    return solve(prob, alg=alg, NSTEPS=NSTEPS)
end

function _solve_heat1d(; file,
                         δ=1e-5,
                         NSTEPS=10_000,
                         model=Forward(inv=true),
                         X0=nothing,
                         alg=BOX(δ=δ, approx_model=model))

    ## load model, we convert to a dense matrix to take the inverse in the discretization
    K, C = load_heat1d(file=file)
    A = -inv(Matrix(C)) * K |> Matrix
    n = size(A, 1)

    if isnothing(X0)
        X0 = initial_distribution(n)
    end

    prob = @ivp(x' = A*x, x(0) ∈ X0)

    return prob, alg, NSTEPS
end

maxtemp_heat1d(sol; nodes=1:dim(sol), kwargs...) = maxtemp(sol, nodes=nodes)
maxtemp_heat1d(; nodes=nothing, kwargs...) = maxtemp(solve_heat1d(kwargs...), nodes=nodes)

mintemp_heat1d(sol; nodes=1:dim(sol), kwargs...) = mintemp(sol, nodes=nodes)
mintemp_heat1d(; nodes=nothing, kwargs...) = mintemp(solve_heat1d(kwargs...), nodes=nodes)

# dado el numero de puntos, discretiza el dominio [0, 1] devuelve el perfil de temperaturas analitico
# en los tiempos dados en un tiempo fijo
function perfil_temp_teo(npuntos, Tampl, t::AbstractFloat)
    rho       = 1. ;
    cSpHe     = 1. ;
    kCond     = 1. ;

    xs = range(0, 1, length=npuntos)
    alpha = kCond / ( rho * cSpHe ) ;
    Tout = Tampl * ( exp(-(  pi )^2 * alpha * t ) * sin.(pi * xs) + exp(-(3*pi)^2 * alpha * t ) * 0.5 * sin.( 3 * pi * xs ) );
    return xs, Tout
end

# para un nodo dado, me da la temperatura en funcion del tiempo en los instantes
# de tiempo dados por el vector t
# xs : posicion
# Tampl = 1 : perfil de temperaturas con suma de senos
# Tampl = 1.1 : lo anterior pero con un 10 % mas
function evolucion_temp_teo(xs::AbstractFloat, Tampl, t::AbstractVector)
    rho       = 1. ;
    cSpHe     = 1. ;
    kCond     = 1. ;

    alpha = kCond / ( rho * cSpHe ) ;
    Tout = Tampl * ( exp.(-(  pi )^2 * alpha * t ) * sin.(pi * xs) + exp.(-(3*pi)^2 * alpha * t ) * 0.5 * sin.( 3 * pi * xs ) );
    return Tout
end

# Solucion numerica con Backward Euler
# use tol = 0.1 for 10% increase in the initial uncertainty
# use tol = -0.1 for 10% decrease in the initial uncertainty
function solve_heat1d_implicit_euler(; file, δ=1e-5, tol=0.1, NSTEPS=10_000, U₀=nothing)

    K, C = load_heat1d(file=file)

    n = size(K, 1)
    M = zeros(n, n)
    R = zeros(n)
    sys = SecondOrderAffineContinuousSystem(M, C, K, R)
    X0 = initial_distribution(n, abs(tol))

    ## extremal value
    x₀ = initial_state_singleton(X0, tol, U₀)

    prob = InitialValueProblem(sys, (x₀, x₀))

    alg = BackwardEuler(Δt=δ)
    return solve(prob, alg, NSTEPS=NSTEPS)
end

function initial_state_singleton(X0, tol, U₀)
    if isnothing(U₀)
        x₀ = X0.center + sign(tol) * X0.radius
    elseif U₀ == :sample
        x₀ = hcat(sample(X0, 1)...) # sample(X0, 1)
    elseif U₀ == :zigzag₊
        x₀ = [X0.center[i] + (-1)^i*X0.radius[i] for i in 1:n]
    elseif U₀ == :zigzag₋
        x₀ = [X0.center[i] + (-1)^(i-1)*X0.radius[i] for i in 1:n]
    end
    return x₀
end

## no-op
ReachabilityAnalysis.convexify(R::ReachabilityAnalysis.AbstractReachSet) = R

function _gradient_cpa(sol, t::AbstractFloat)
    X = sol(t) |> convexify
    n = dim(X)  # for 100 nodes, n = 99

    ST = Interval{Float64, IA.Interval{Float64}}
    grad = Vector{ST}(undef, n+1)

    dx = 1/n

    ## condiciones de borde (cero en los bordes)
    R = overapproximate(Projection(X, 1:1), Interval)
    grad[1] = set(R).dat * (1/dx) |> Interval
    L = overapproximate(Projection(X, n:n), Interval)
    grad[n+1] = -set(L).dat * (1/dx) |> Interval

    for i in 1:n-1
        L = overapproximate(Projection(X, i:i), Interval)
        R = overapproximate(Projection(X, i+1:i+1), Interval)
        grad[i+1] = Interval((set(R).dat - set(L).dat) * (1/dx))
    end
    dt = tspan(X)

    return ReachSet(CartesianProductArray(grad), dt)
end

function _gradient(sol, t::AbstractFloat)
    Tg = _gradient_cpa(sol, t)

    arr = array(set(Tg))
    n = dim(sol) + 1
    xs = range(0, 1, length=n+1)
    dx = [Interval(xs[i], xs[i+1]) for i in 1:n]
    u = UnionSetArray([dxi × ai for (dxi, ai) in zip(dx, arr)])
end

## function para calcular el gradiente de una sol numerica
## TODO improve using single loop
function gradiente_numerica(sol, dx)
    U = sol.U
    N = length(U) # numero de pasos temporales
    out = [vcat(x[1], diff(x, dims=1)[:], -x[end]) .* (1/dx) for x in U]
    return out
end

struct SolBE{S, G}
    displacement::S
    gradient::G
end

# example uso _varias_sol_be:
#
#   num_paso = 100
#   n_curvas = 100
#   tol = 0.1
#   δ = 1E-5
#   xs, out = _varias_sol_be(δ=δ, tol=tol, num_paso = num_paso, n_curvas=n_curvas)
#   figg = plot()
#   for i in 1:length(out)-2
#       plot!(figg, xs, out[i], lab="")
#   end
#   plot!(figg, xs, out[end-1], c=:magenta, lab="")
#   plot!(figg, xs, out[end], c=:magenta, lab="")
#   fig
function varias_sol(; file, δ, tol, NSTEPS, NCURVAS)
    ## SM = Vector{Matrix{Float64}}
    ## SV = Vector{Vector{Float64}}
    gradientes = [] # Vector{SolBE{SM, SV}}()
    temperaturas = []

    ## dx = 1/(99 + 1)

    ## calcular sol random
    ## separo el primer paso para calcular dx
    sol = solve_heat1d_implicit_euler(file=file, δ=δ, tol=tol, NSTEPS=NSTEPS, U₀=:sample)
    dx = 1/(length(sol.U[1])+1)
    grad = gradiente_numerica(sol, dx)
    push!(gradientes, grad) # SolBE(sol.U, gr))
    push!(temperaturas, sol)

    ## pasos siguientes
    for k in 2:NCURVAS
        sol = solve_heat1d_implicit_euler(file=file, δ=δ, tol=tol, NSTEPS=NSTEPS, U₀=:sample)
        grad = gradiente_numerica(sol, dx)
        push!(gradientes, grad) # SolBE(sol.U, gr))
        push!(temperaturas, sol)
    end

    ## calcular sol extremas
    for t in [tol, -tol]
        sol = solve_heat1d_implicit_euler(file=file, δ=δ, tol=t, NSTEPS=NSTEPS, U₀=nothing)
        grad = gradiente_numerica(sol, dx)
        push!(gradientes, grad) # SolBE(sol.U, gr))
    end

    for Uopt in [:zigzag₊, :zigzag₋]
        sol = solve_heat1d_implicit_euler(file=file, δ=δ, tol=tol, NSTEPS=NSTEPS, U₀=Uopt)
        grad = gradiente_numerica(sol, dx)
        push!(gradientes, grad) # SolBE(sol.U, gr))
        push!(temperaturas, sol)
    end

    return temperaturas, gradientes
end

## Codigo Para analisis de gradientes
function LGG_dirs_solucion(n)
    ## n = statedim(prob)

    M = Diagonal(ones(n))
    M = [Vector(M[:, i]) for i in 1:size(M, 1)]

    return vcat(M, .-M)

end

function LGG_dirs_gradiente(n)
    B = Bidiagonal(ones(n), -ones(n-1), :U)
    B = [Vector(B[:, i]) for i in 2:size(B, 2)]
    vini = zeros(n); vini[1] = 1.0
    vend = zeros(n); vend[end] = -1.0
    insert!(B, 1, vini);
    push!(B, vend);

    return vcat(B, .-B)
end

## grad_lgg = solucion con LGG usando LGG_grad_dirs
function gradiente_LGG_graddirs(; grad_lgg, L=1.0, paso_temporal)
    n = Int(0.5*length(grad_lgg[1].sf))
    Δx = L/n
    xs = range(0, L, step=Δx)
    α = 1/Δx
    intervs_xs = [IA.Interval(xs[i], xs[i+1]) for i in 1:n]
    return [ReachSet(Interval(-grad_lgg[paso_temporal].sf[i+n] * α, grad_lgg[paso_temporal].sf[i] * α), intervs_xs[i]) for i in 1:n] |> Flowpipe
end

## solu_lgg = solucion con LGG usando LGG_solu_dirs (direcciones cartesianas)
function gradiente_LGG_soludirs(; solu_lgg, LGG_grad_dirs, L=1.0, paso_temporal)
    n = Int(0.5*length(LGG_grad_dirs))
    Δx = L/n
    xs = range(0, L, step=Δx)
    α = 1/Δx
    intervs_xs = [IA.Interval(xs[i], xs[i+1]) for i in 1:n]
    return [ReachSet(Interval(-α*ρ(LGG_grad_dirs[i+n], solu_lgg[paso_temporal]), α*ρ(LGG_grad_dirs[i], solu_lgg[paso_temporal])), intervs_xs[i]) for i in 1:n] |> Flowpipe
end

# ===============================================================
# Methods to plot flowpipe "snapshots"
# ===============================================================

# no-op
ReachabilityAnalysis.convexify(R::ReachabilityAnalysis.AbstractReachSet) = R

function flowpipe_snapshot(sol, t)
    n = dim(sol)
    solt = sol(t) |> convexify
    X = [set(overapproximate(Projection(solt, i:i), Interval)) for i in 1:n]
    R = ReachSet(CartesianProductArray(X), tspan(solt))
    _flowpipe_snapshot(R)
end

# cartesian product array of intervals
function _flowpipe_snapshot(X::ReachSet{N, CartesianProductArray{N, Interval{N, IntervalArithmetic.Interval{N}}}}) where {N}
    n = dim(X)
    xs = range(0, 1, length=n+2)

    fp_espacial = array(set(X))
    n_x = length(fp_espacial)
    fpx1 = [Interval(xs[k+1], xs[k+1]) × fp_espacial[k] for k in 1:n_x]

    T0 = 0.0 # temperature at the endpoints
    aux = [Interval(xs[1], xs[1]) × Interval(T0[1], T0[1])]
    for k = 1:n_x
        push!(aux, fpx1[k])
    end
    push!(aux, Interval(xs[end], xs[end]) × Interval(T0[end], T0[end]))
    fpx2 = UnionSetArray([ConvexHull(aux[k], aux[k+1]) for k in 1:length(aux)-1])

    return fpx2
end

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


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

"""
    create_lifecycle_policy(description, execution_role_arn, policy_details, state)
    create_lifecycle_policy(description, execution_role_arn, policy_details, state, params::Dict{String,<:Any})

Creates a policy to manage the lifecycle of the specified Amazon Web Services resources.
You can create up to 100 lifecycle policies.

# Arguments
- `description`: A description of the lifecycle policy. The characters ^[0-9A-Za-z _-]+ are
  supported.
- `execution_role_arn`: The Amazon Resource Name (ARN) of the IAM role used to run the
  operations specified by the lifecycle policy.
- `policy_details`: The configuration details of the lifecycle policy.
- `state`: The desired activation state of the lifecycle policy after creation.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Tags"`: The tags to apply to the lifecycle policy during creation.
"""
function create_lifecycle_policy(
    Description,
    ExecutionRoleArn,
    PolicyDetails,
    State;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "POST",
        "/policies",
        Dict{String,Any}(
            "Description" => Description,
            "ExecutionRoleArn" => ExecutionRoleArn,
            "PolicyDetails" => PolicyDetails,
            "State" => State,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_lifecycle_policy(
    Description,
    ExecutionRoleArn,
    PolicyDetails,
    State,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "POST",
        "/policies",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "Description" => Description,
                    "ExecutionRoleArn" => ExecutionRoleArn,
                    "PolicyDetails" => PolicyDetails,
                    "State" => State,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_lifecycle_policy(policy_id)
    delete_lifecycle_policy(policy_id, params::Dict{String,<:Any})

Deletes the specified lifecycle policy and halts the automated operations that the policy
specified.

# Arguments
- `policy_id`: The identifier of the lifecycle policy.

"""
function delete_lifecycle_policy(
    policyId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dlm(
        "DELETE",
        "/policies/$(policyId)/";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_lifecycle_policy(
    policyId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "DELETE",
        "/policies/$(policyId)/",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_lifecycle_policies()
    get_lifecycle_policies(params::Dict{String,<:Any})

Gets summary information about all or the specified data lifecycle policies. To get
complete information about a policy, use GetLifecyclePolicy.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"policyIds"`: The identifiers of the data lifecycle policies.
- `"resourceTypes"`: The resource type.
- `"state"`: The activation state.
- `"tagsToAdd"`: The tags to add to objects created by the policy. Tags are strings in the
  format key=value. These user-defined tags are added in addition to the Amazon Web
  Services-added lifecycle tags.
- `"targetTags"`: The target tag for a policy. Tags are strings in the format key=value.
"""
function get_lifecycle_policies(; aws_config::AbstractAWSConfig=global_aws_config())
    return dlm("GET", "/policies"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
end
function get_lifecycle_policies(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dlm(
        "GET", "/policies", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    get_lifecycle_policy(policy_id)
    get_lifecycle_policy(policy_id, params::Dict{String,<:Any})

Gets detailed information about the specified lifecycle policy.

# Arguments
- `policy_id`: The identifier of the lifecycle policy.

"""
function get_lifecycle_policy(policyId; aws_config::AbstractAWSConfig=global_aws_config())
    return dlm(
        "GET",
        "/policies/$(policyId)/";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_lifecycle_policy(
    policyId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "GET",
        "/policies/$(policyId)/",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags for the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dlm(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(tags, resource_arn)
    tag_resource(tags, resource_arn, params::Dict{String,<:Any})

Adds the specified tags to the specified resource.

# Arguments
- `tags`: One or more tags.
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.

"""
function tag_resource(Tags, resourceArn; aws_config::AbstractAWSConfig=global_aws_config())
    return dlm(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    Tags,
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Tags" => Tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes the specified tags from the specified resource.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the resource.
- `tag_keys`: The tag keys.

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dlm(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_lifecycle_policy(policy_id)
    update_lifecycle_policy(policy_id, params::Dict{String,<:Any})

Updates the specified lifecycle policy.

# Arguments
- `policy_id`: The identifier of the lifecycle policy.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: A description of the lifecycle policy.
- `"ExecutionRoleArn"`: The Amazon Resource Name (ARN) of the IAM role used to run the
  operations specified by the lifecycle policy.
- `"PolicyDetails"`: The configuration of the lifecycle policy. You cannot update the
  policy type or the resource type.
- `"State"`: The desired activation state of the lifecycle policy after creation.
"""
function update_lifecycle_policy(
    policyId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dlm(
        "PATCH",
        "/policies/$(policyId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_lifecycle_policy(
    policyId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dlm(
        "PATCH",
        "/policies/$(policyId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

ext: jl


text:
@doc doc"""
    Stiefel{n,k,T} <: Manifold

The Stiefel manifold consists of all $n\times k$, $n\geq k$ orthonormal matrices, i.e.

````math
\mathcal M = \{ x \in \mathbb F^{n\times k} : x^{\mathrm{H}}x = I_k \},
````

where $\mathbb F \in \{\mathbb R, \mathbb C\}$,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and
$I_n \in \mathbb R^{n\times n}$ denotes the $k \times k$ identity matrix.

The tangent space at a point $x\in\mathcal M$ is given by

````math
T_x\mathcal M = \{ v \in \mathbb{F}^{n\times k} : x^{\mathrm{H}}v + v^{\mathrm{H}}x=0_n\},
````

where $0_n$ is the $k\times k$ zero matrix.

The metric is either inherited from $\mathbb R^{n,k}$ for the real-valued case
or the one inherited from interpreting the complex valued entries in the Gaussian
plane $\mathbb R^2$ and then over all entries as before, i.e. the latter
may be called an Hermitian metric in the complex-valued matrices.

The manifold is named after
[Eduard L. Stiefel](https://en.wikipedia.org/wiki/Eduard_Stiefel) (1909–1978).

# Constructor
    Stiefel(n,k,F=ℝ)

Generate the (real-valued) Stiefel manifold of $n\times k$ dimensional orthonormal matrices.
"""
struct Stiefel{n,k,F} <: Manifold end

Stiefel(n::Int, k::Int, F::AbstractNumbers = ℝ) = Stiefel{n,k,F}()

@doc doc"""
    check_manifold_point(M::Stiefel, x; kwargs...)

Check whether `x` is a valid point on the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$,
i.e. that it has the right [`AbstractNumbers`](@ref) type and $x^{\mathrm{H}}x$
is (approximatly) the identity, where $\cdot^{\mathrm{H}}$ is the complex conjugate
transpose. The settings for approximately can be set with `kwargs...`.
"""
function check_manifold_point(M::Stiefel{n,k,T}, x; kwargs...) where {n,k,T}
    if (T === ℝ) && !(eltype(x) <: Real)
        return DomainError(
            eltype(x),
            "The matrix $(x) is not a real-valued matrix, so it does noe lie on the Stiefel manifold of dimension ($(n),$(k)).",
        )
    end
    if (T === ℂ) && !(eltype(x) <: Real) && !(eltype(x) <: Complex)
        return DomainError(
            eltype(x),
            "The matrix $(x) is neiter real- nor complex-valued matrix, so it does noe lie on the complex Stiefel manifold of dimension ($(n),$(k)).",
        )
    end
    if any(size(x) != representation_size(M))
        return DomainError(
            size(x),
            "The matrix $(x) is does not lie on the Stiefel manifold of dimension ($(n),$(k)), since its dimensions are wrong.",
        )
    end
    c = x' * x
    if !isapprox(c, one(c); kwargs...)
        return DomainError(
            norm(c - one(c)),
            "The point $(x) does not lie on the Stiefel manifold of dimension ($(n),$(k)), because x'x is not the unit matrix.",
        )
    end
end

@doc doc"""
    check_tangent_vector(M::Stiefel, x, v; kwargs...)

Check whether `v` is a valid tangent vector at `x` on the [`Stiefel`](@ref)
`M`=$\operatorname{St}(n,k)$, i.e. the [`AbstractNumbers`](@ref) fits and
it (approximtly) holds that $x^{\mathrm{H}}v + v^{\mathrm{H}}x = 0$, where
`kwargs...` is passed to the `isapprox`.
"""
function check_tangent_vector(M::Stiefel{n,k,T}, x, v; kwargs...) where {n,k,T}
    mpe = check_manifold_point(M, x)
    mpe === nothing || return mpe
    if (T === ℝ) && !(eltype(v) <: Real)
        return DomainError(
            eltype(v),
            "The matrix $(v) is not a real-valued matrix, so it can not be a tangent vector to the Stiefel manifold of dimension ($(n),$(k)).",
        )
    end
    if (T === ℂ) && !(eltype(v) <: Real) && !(eltype(v) <: Complex)
        return DomainError(
            eltype(v),
            "The matrix $(v) is neiter real- nor complex-valued matrix, so it can not bea tangent vectorto the complex Stiefel manifold of dimension ($(n),$(k)).",
        )
    end
    if any(size(v) != representation_size(M))
        return DomainError(
            size(v),
            "The matrix $(v) is does not lie in the tangent space of $(x) on the Stiefel manifold of dimension ($(n),$(k)), since its dimensions are wrong.",
        )
    end
    if !isapprox(x' * v + v' * x, zeros(k, k); kwargs...)
        return DomainError(
            norm(x' * v + v' * x),
            "The matrix $(v) is does not lie in the tangent space of $(x) on the Stiefel manifold of dimension ($(n),$(k)), since x'v + v'x is not the zero matrix.",
        )
    end
end

@doc doc"""
    exp(M, x, v)

Compute the exponential map on the [`Stiefel`](@ref)`{n,k,T}`() manifold `M`
eminating from `x` into tangent direction `v`.

$\operatorname{exp}_{x} v = \begin{pmatrix}
   x\\v
 \end{pmatrix}
 \exp
 \left(
 \begin{pmatrix} x^{\mathrm{H}}v & - v^{\mathrm{H}}v\\
 I_n & x^{\mathrm{H}}v\end{pmatrix}
 \right)
\begin{pmatrix}  \exp( -x^{\mathrm{H}}v) \\ 0_n\end{pmatrix}$

where $\exp$ denotes matrix exponential,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k$ and
$0_k$ are the identity matrix and the zero matrix of dimension $k \times k$, respectively.
"""
exp(::Stiefel, ::Any...)

function exp!(M::Stiefel{n,k}, y, x, v) where {n,k}
    return copyto!(
        y,
        [x v] *
        exp([x'v -v' * v; one(zeros(eltype(x), k, k)) x' * v]) *
        [exp(-x'v); zeros(eltype(x), k, k)],
    )
end

@doc doc"""
    inner(M::Stiefel, x, v, w)

Compute the inner product for two tangent vectors `v`, `w` from the
tangent space of `x` on the [`Stiefel`](@ref) manifold `M`. The formula reads

````math
(v,w)_x = \operatorname{trace}(v^{\mathrm{H}}w),
````
i.e. the [`EuclideanMetric`](@ref) from the embedding restricted to the tangent
space. For the complex-valued case this is the Hermitian metric, to be precise.
"""
inner(::Stiefel, x, v, w) = dot(v, w)

@doc doc"""
    inverse_retract(M::Stiefel, x, y, ::PolarInverseRetraction)

Compute the inverse retraction based on a singular value decomposition
for two points `x`, `y` on the [`Stiefel`](@ref) manifold `M` and return
the resulting tangent vector in `v`. This follows the folloing approach:
From the Polar retraction we know that

````math
\operatorname{retr}_x^{-1}(y) = ys - q
````

if such a symmetric positive definite $k\times k$ matrix exists. Since $ys-q$ is
also a tangent vector at $x$ we obtain

````math
x^{\mathrm{H}}ys + s(x^{\mathrm{H}}y)^{\mathrm{H}} + 2I_k = 0,
````
which can either be solved by a Lyapunov approach or a continuous-time
algebraic Riccati equation as described in [^KanekoFioriTanaka2013]

This implementation follows the Lyapunov approach.

    inverse_retract(M, x, y, ::QRInverseRetraction)

Compute the inverse retraction based on a qr decomposition
for two points `x`, `y` on the [`Stiefel`](@ref) manifold `M` and return
the resulting tangent vector in `v`. The computation follows Algorithm 1
in [^KanekoFioriTanaka2013].

[^KanekoFioriTanaka2013]:
    > T. Kaneko, S. Fiori, T. Tanaka: "Empirical Arithmetic Averaging over the
    > Compact Stiefel Manifold", IEEE Transactions on Signal Processing, 2013,
    > doi: [10.1109/TSP.2012.2226167](https://doi.org/10.1109/TSP.2012.2226167).
"""
inverse_retract(::Stiefel, ::Any...)

function inverse_retract!(::Stiefel, v, x, y, ::PolarInverseRetraction)
    A = x' * y
    H = -2 * one(x' * x)
    B = lyap(A, H)
    return copyto!(v, y * B - x)
end
function inverse_retract!(::Stiefel{n,k}, v, x, y, ::QRInverseRetraction) where {n,k}
    A = x' * y
    R = zeros(typeof(one(eltype(x)) * one(eltype(y))), k, k)
    for i = 1:k
        b = zeros(i)
        b[i] = 1
        b[1:(end-1)] = -transpose(R[1:(i-1), 1:(i-1)]) * A[i, 1:(i-1)]
        R[1:i, i] = A[1:i, 1:i] \ b
    end
    return copyto!(v, y * R - x)
end

function isapprox(M::Stiefel, x, v, w; kwargs...)
    return isapprox(sqrt(inner(M, x, zero_tangent_vector(M, x), v - w)), 0; kwargs...)
end
isapprox(M::Stiefel, x, y; kwargs...) = isapprox(norm(x - y), 0; kwargs...)

@doc doc"""
    manifold_dimension(M::Stiefel)

Return the dimension of the [`Stiefel`](@ref) manifold `M`=$\operatorname{St}(n,k,𝔽)$.
The dimension is given by

````math
\dim \mathrm{Stiefel}(n, k, ℝ) &= nk - \frac{1}{2}k(k+1)\\
\dim \mathrm{Stiefel}(n, k, ℂ) &= 2nk - k^2\\
\dim \mathrm{Stiefel}(n, k, ℍ) &= 4nk - k(2k-1)
````
"""
manifold_dimension(::Stiefel{n,k,ℝ}) where {n,k} = n * k - div(k * (k + 1), 2)
manifold_dimension(::Stiefel{n,k,ℂ}) where {n,k} = 2 * n * k - k * k
manifold_dimension(::Stiefel{n,k,ℍ}) where {n,k} = 4 * n * k - k * (2k - 1)

@doc doc"""
    project_tangent(M, x, v)

Project `v` onto the tangent space of `x` to the [`Stiefel`](@ref) manifold `M`.
The formula reads

````math
\operatorname{proj}_{\mathcal M}(x,v) = v - x \operatorname{Sym}(x^{\mathrm{H}}v),
````

where $\operatorname{Sym}(y)$ is the symmetrization of $y$, e.g. by
$\operatorname{Sym}(y) = \frac{y^{\mathrm{H}}+y}{2}$.
"""
project_tangent(::Stiefel, ::Any...)

project_tangent!(::Stiefel, w, x, v) = copyto!(w, v - x * Symmetric(x' * v))

@doc doc"""
    retract(M, x, v, ::PolarRetraction)

Compute the SVD-based retraction [`PolarRetraction`](@ref) on the
[`Stiefel`](@ref) manifold `M`. With $USV = x + v$ the retraction reads
````math
\operatorname{retr}_x(v) = U\bar{V}^\mathrm{H}.
````

    retract(M, x, v, ::QRRetraction )

Compute the QR-based retraction [`QRRetraction`](@ref) on the
[`Stiefel`](@ref) manifold `M`. With $QR = x + v$ the retraction reads
````math
\operatorname{retr}_x(v) = QD,
````
where D is a $n\times k$ matrix with
````math
D = \operatorname{diag}\bigl(\operatorname{sgn}(R_{ii}+0,5)_{i=1}^k \bigr),
````
where $\operatorname{sgn}(x) = \begin{cases}
1 & \text{ for } x > 0,\\
0 & \text{ for } x = 0,\\
-1& \text{ for } x < 0.
\end{cases}$
"""
retract(::Stiefel, ::Any...)

function retract!(::Stiefel, y, x, v, ::PolarRetraction)
    s = svd(x + v)
    return mul!(y, s.U, s.Vt)
end
function retract!(::Stiefel, y, x, v, ::QRRetraction)
    qrfac = qr(x + v)
    d = diag(qrfac.R)
    D = Diagonal(sign.(sign.(d .+ 0.5)))
    return copyto!(y, Matrix(qrfac.Q) * D)
end

@doc doc"""
    representation_size(M::Stiefel)

Returns the representation size of the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$,
i.e. `(n,k)`, which is the matrix dimensions.
"""
@generated representation_size(::Stiefel{n,k}) where {n,k} = (n, k)

@doc doc"""
    zero_tangent_vector(M::Stiefel, x)

Returns the zero tangent vector from the tangent space at `x` on the [`Stiefel`](@ref)
`M`=$\operatorname{St}(n,k)$, i.e. an `(n,k)` zero matrix.
"""
zero_tangent_vector(::Stiefel, ::Any...)

zero_tangent_vector!(::Stiefel, v, x) = fill!(v, 0)


text:
"""
    sfmodel_MoMTest(<keyword arguments>)

Provide method-of-moments (MoM) tests on the distribution assumptions of the
composed error (v-u or v+u) of stochastic frontier models. Also provide MoM
estimates of the model parameters as well as the inefficiency (JLMS) and
efficiency (BC) index. Return various results in a distionary.

# Arguments

- `sfdist(::Vararg)`: the distribution assumption of the one-sided stochastic
  variable (aka inefficiency term) of the model. Currently the
  choices include `half` (or `h`) and `exponential` (or `expo`, `e`)
- `sftype(::Vararg)`: whether the model is a `production` (or `prod`) frontier
  or a `cost` frontier.
- `@depvar(::Vararg)`: the dependent variable from the `data`.
- `@frontier(::Vararg)`: a list of variables, separated by commas, in the frontier function.
- `data::DataFrames`: the data in the DataFrames format.
- `omega (or ω)::Union{Real, Vector, Tuple}=<a vector of numbers>`:
  The `ω` parameter used in the test. It could take the form of a single value
  (e.g., `ω=1`) or a vector or tuple (e.g., `ω=(0.5, 1, 2)`). Default is
  `ω=1.0`.
- level::Real=<number>: The significance level (default=0.05) of the bias-corrected
  confidence intervals. If `level`>0.5, it is automatically transformed to
  `1-level`, such that `level=0.05` and `level=0.95` both return 95%
  confidene intervals at the 5% significance level.
- `verbose::Bool`: whether print results on screen. Default is `=true`.
- `testonly::Bool`: whether to print only the test results. Default is
  `false`, which print both the test and the estimation results on
  the screen.

# Remarks

- Currently, the two-sided error (aka `v`) of the model is assumed to
  follow a normal distribuion in the test.
- Reference: Chen, Y.T., & Wang, H.J. (2012). Centered-Residuals-Based Moment
  Estimator and Test for Stochastic Frontier Models. Econometric Reviews,
  31(6), 625-653.

# Examples
```julia-repl
julia> df = DataFrame(CSV.File("demodata.csv"))
julia> df[!, :_cons] .=1.0;
julia> res = sfmodel_MoMTest(sftype(prod), sfdist(expo),
                             @depvar(y), @frontier(_cons, x1, x2, x3),
                             data=df, ω=(0.5,1,2))
julia> res.list  # a list of the saved results
julia> res.jlms  # the JLMS inefficiency index
```
"""
function sfmodel_MoMTest(arg::Vararg; data::DataFrame, ω::Union{Vector, Real, Tuple}=1.0,
        omega::Union{Vector, Real, Tuple}=1.0, level::Real=0.05,
        verbose::Bool=true, testonly::Bool=false)

    _dicMoM = Dict{Symbol,Any}()

    for k in (:type, :dist, :depvar, :frontier)
        _dicMoM[k] = nothing
    end

    for d in :($(arg))
        _dicMoM[d[1]] = d[2]
    end

    # --- check syntax; in house so that no need using global dic ----

    for k in (:type, :dist, :depvar, :frontier)
        if  (_dicMoM[k] === nothing)
            throw("For the method of moment test/estimation, the `$k` equation is missing in sfmodel_MoMTest().")
        end
    end

    if length(keys(_dicMoM)) > 4  # should have only 4 at this moment
       throw("For the MoM model specification, only `sftype`, `sfdist`, `@depvar`, and `@frontier` are needed. You apparently provided more than these.")
    end

    if typeof(data) != DataFrame
       throw("`data=` needs to be a DataFrame.")
    end

    #-- retrieve parameters -------

    if !(1 in ω) || (length(ω) > 1)
        ω = ω
    elseif !(1 in omega) || (length(omega) > 1)
        ω = omega
    end

    s = uppercase(String(_dicMoM[:dist][1])[1:1])  # s = "H"(half), s = "E"(exponential)
    if _dicMoM[:type][1] == :cost
        _porc = -1.0
    else
        _porc = 1.0
    end

    ((level > 0.0) && (level < 1.0)) || throw("The significance level (`level`) should be between 0 and 1.")
    if level > 0.5
       level = 1-level  # 0.95 -> 0.05
    end
    α = level

    ##### print the title ############

     if verbose

        println("\n****************************************")
        println("** Moment Based Tests and Estimations **")
        println("****************************************")

     end

    ############################### Start Calculation #####################################################

    if (s == "H") || (s == "E")

        yvar = data[:, _dicMoM[:depvar]]
        xvar = data[:, _dicMoM[:frontier]]

        sT = size(yvar,1)

        #---- pick the constant ----------

        cnsname, cnspos, cnsnum = pickConsNameFromDF(xvar)
        if cnsnum == 0
           printstyled("\nThere is no intercept (no constant variable) in the model. Are you sure?\n"; color = :green)
        elseif cnsnum == 1
           cnspos = Int(cnspos[1])
        elseif cnsnum > 1
           printstyled("\nThere are multiple constant varibles in the model resulting in multicollinearity. The estimation is abort.\n"; color = :red)
           throw("Model specification errors.")
        else
           throw("Something is not right about the specification of constant variables in the model")
        end

        yvar  =  convert(Array{Float64}, Matrix(yvar))
        xvar  =  convert(Array{Float64}, Matrix(xvar))

        β_ols  = xvar \ yvar  # OLS estiamte, uses a pivoted QR factorization;
        resid  = yvar - xvar*β_ols
        sse    = sum((resid).^2)
        ssd    = sqrt(sse/(size(resid,1)-size(xvar,2))) # sample standard deviation; σ² = (1/(N-K))* Σ ϵ^2
        ll_ols = sum(normlogpdf.(0, ssd, resid)) # ols log-likelihood
        sk_ols = sum((resid).^3) / ((ssd^3)*(size(resid,1))) # skewnewss of ols residuals
        std_ols = sqrt.(diag((ssd^2)*inv(xvar'*xvar)))
        Resid   = _porc*resid  # prod: v-(u-E(u)),  cost: v+(u-E(u))

        CResid = Resid.-mean(Resid)

        dm2 = CResid.^2
        m2 = mean(dm2)
        dm3 = CResid.^3
        m3 = mean(dm3)

        if s == "H"
            minfo1 = "normal and half-normal"
            λ₂ = 1-2/π
            λ₃ = (1-4/π) * sqrt(2/π)
        elseif s == "E"
            minfo1 = "normal and exponential"
            λ₂ = 1
            λ₃ = -2
        end

        σᵤ² = ((m3/λ₃)^2)^(1/3)
        σᵥ² = m2 - λ₂ * σᵤ²

        if cnsnum == 1
            if s == "H"
                β_ols[cnspos] = β_ols[cnspos] .+  _porc*sqrt(σᵤ²)*sqrt(2/π)
            elseif s == "E"
                β_ols[cnspos] = β_ols[cnspos] .+  _porc*sqrt(σᵤ²)
            end

            v_adj  = m2 * (inv(xvar' * xvar))
            std_ols[cnspos] = Real(sqrt(v_adj[cnspos, cnspos]))
        end

        if s == "H"
            Resid2 = Resid .- _porc*sqrt(σᵤ²)*sqrt(2/π)
        elseif s == "E"
            Resid2 = Resid .- _porc*sqrt(σᵤ²)
        end

        Grad2 = CResid.^2 .- (σᵥ²+λ₂*σᵤ²)
        Grad3 = CResid.^3 .- λ₃*(σᵤ²^1.5).-3*(σᵥ²+λ₂*σᵤ²)*(CResid)
        HessE = [-1 -λ₂ ; 0  -1.5*λ₃*sqrt(σᵤ²)]

        Grad22=Grad2.^2
        Grad33=Grad3.^2
        Grad23=Grad2.*Grad3
        mGrad22 = mean(Grad22)
        mGrad33 = mean(Grad33)
        mGrad23 = mean(Grad23)

        mMat = [mGrad22 mGrad23 ; mGrad23 mGrad33]
        vcovM  = (1/sT)*inv(HessE)*mMat*inv(HessE)'
        myvstd = sqrt(vcovM[1,1])
        myustd = sqrt(vcovM[2,2])

        nofom = length(ω)
        mycos = zeros(nofom)
        mysin = zeros(nofom)
        testcs = zeros(nofom)

        for i in 1:nofom
            Ecoswv = exp(-0.5*(ω[i]^2)*σᵥ²)
            EDcoswv= -0.5*(ω[i]^2)*Ecoswv

            if s == "H"  #The normal-half normal model
                μ = sqrt(σᵤ²)*sqrt(2/π)
                Dmu= sqrt(2/π)/(2*sqrt(σᵤ²))
                Ecoswu = exp(-0.5*(ω[i]^2)*σᵤ²)
                Edcoswu = -0.5*(ω[i]^2)*Ecoswu
                t = ω[i]*sqrt(σᵤ²/2)
                erfiw  = erfi(t)
                Esinwu = Ecoswu*erfiw
                EDsinwu =Edcoswu*erfiw+ω[i]/sqrt(2*π*σᵤ²)
            elseif s == "E" #The normal-exponential model
                μ = sqrt(σᵤ²)
                Dmu = 1/(2*sqrt(σᵤ²))
                Ecoswu = 1/(1+σᵤ²*ω[i]^2)
                Edcoswu = -(ω[i]^2)/((1+σᵤ²*ω[i]^2)^2)
                Esinwu = ω[i]*sqrt(σᵤ²)/(1+σᵤ²*(ω[i]^2))
                EDsinwu = ω[i]/(2*sqrt(σᵤ²)*(1+σᵤ²*ω[i]^2))-ω[i]^3*sqrt(σᵤ²)/(1+σᵤ²*(ω[i]^2))^2
            end

            hc = Ecoswu*cos(ω[i]*μ)+Esinwu*sin(ω[i]*μ)
            Mc = Ecoswv*hc
            Dhc= Edcoswu*cos(ω[i]*μ)+EDsinwu*sin(ω[i]*μ)-Ecoswu*sin(ω[i]*μ)*ω[i]*Dmu+Esinwu*cos(ω[i]*μ)*ω[i]*Dmu
            Dc = [-EDcoswv*hc , -Ecoswv*Dhc]

            hs = Ecoswu*sin(ω[i]*μ)-Esinwu*cos(ω[i]*μ)
            Ms = Ecoswv*hs
            Dhs = Edcoswu*sin(ω[i]*μ)-EDsinwu*cos(ω[i]*μ)+Ecoswu*cos(ω[i]*μ)*ω[i]*Dmu+Esinwu*sin(ω[i]*μ)*ω[i]*Dmu
            Ds = [-EDcoswv*hs , -Ecoswv*Dhs ]

            TFc=cos.(ω[i]*CResid).-Mc
            TFs=sin.(ω[i]*CResid).-Ms
            Gradc=TFc.+Ms*ω[i]*CResid
            Grads=TFs.-Mc*ω[i]*CResid

            invHess = inv(HessE')
            Xi = [Gradc,Grads] - ([Dc Ds]' * invHess' * [Grad2,Grad3])
            xic = Xi[1]
            xis = Xi[2]

            mTFc = mean(TFc)
            mTFs = mean(TFs)
            mTFcs = [mTFc, mTFs]

            # mat CTestcNH = sT*mTFcNH*inv(xicNH'*xicNH/sT)*mTFcNH
            xic2 = xic.^2
            xic2sum = 1/mean(xic2)
            mycosine = sT * mTFc * xic2sum * mTFc
            mycos[i] = mycosine

            # mat CTestsNH = sT*mTFsNH*inv(xisNH'*xisNH/sT)*mTFsNH
            xis2 = xis.^2
            xis2sum = 1/mean(xis2)
            mysine = sT * mTFs * xis2sum * mTFs
            mysin[i] = mysine

            # mat CTestcsNH = sT*mTFcsNH'*inv(xicsNH'*xicsNH/sT)*mTFcsNH
            xics2 = xic.*xis
            inv11 = sum(xic2)
            inv22 =sum(xis2)
            inv12 = sum(xics2)
            invmat = [inv11 inv12 ; inv12 inv22]/sT

            CTestcs = sT*mTFcs'*inv(invmat)*mTFcs
            Testcs = CTestcs[1,1]
            testcs[i] = Testcs

        end
    end

    # ********** calculate log-likelihood value ***********

    if s == "H"
        σ²  = σᵤ² + σᵥ²
        μₛ  = ( - σᵤ² .* Resid2) ./ σ²
        σₛ² = (σᵥ² * σᵤ²) / σ²

        llike = sum( .- 0.5 * log(σ²)
                     .+ normlogpdf.( Resid2 ./ sqrt(σ²))
                     .+ normlogcdf.(μₛ ./ sqrt(σₛ²))
                     .- normlogcdf(0) )
    elseif s == "E"
        σᵤ = sqrt(σᵤ²)
        σᵥ = sqrt(σᵥ²)

        llike = sum( -log(σᵤ) .+ normlogcdf.(-(Resid2 ./ σᵥ)  .- (σᵥ/σᵤ))
                    .+ Resid2 ./ σᵤ .+ σᵥ²/(2*σᵤ²) )
    end

    #***** JLMS & BC ********************************

    if s == "H"
        myeps = Resid .- (_porc)*sqrt(σᵤ²)*sqrt(2/π)  # prod: v-(u-E(u)) - E(u) = v-u;  cost: v+(u-E(u)) + E(u) = v+u
        mustar  = -σᵤ²*myeps/(σᵤ² + σᵥ²)
        sig2star = σᵤ²*σᵥ²/(σᵤ² + σᵥ²)
        rat1 = mustar/sqrt(sig2star)
        _jlms = sqrt(sig2star)*normpdf.(rat1)./normcdf.(rat1) .+ (mustar)
        _bc  = exp.(-mustar .+ 0.5*sig2star).*normcdf.(rat1 .- sqrt(sig2star))./normcdf.(rat1)
    elseif s == "E"
        myeps = Resid .- (_porc)*sqrt(σᵤ²) # prod: v-(u-E(u)) - E(u) = v-u;  cost: v+(u-E(u)) + E(u) = v+u
        mustar = -myeps .- σᵥ²/sqrt(σᵤ²)
        rat1 = mustar/sqrt(σᵥ²)
        _jlms = sqrt(σᵥ²)*normpdf.(rat1)./normcdf.(rat1) .+ mustar
        _bc = exp.(-mustar .+ 0.5*σᵥ²).*normcdf.(rat1 .- sqrt(σᵥ²))./normcdf.(rat1)
    end

    _jlmsM = mean(_jlms)
    _bcM   = mean(_bc)

    #***** critical values **************************

    crit1 = chisqinvcdf(1,1-0.01)
    crit2 = chisqinvcdf(1,1-0.05)
    crit3 = chisqinvcdf(1,1-0.10)

    #***** prepare for printing test results***********

    omega_list = [i for i in ω]
    test_res = hcat(omega_list, mysin, mycos)
    test_header = (["ω", "sine", "cosine"])

    test_res2 = vcat(reshape(test_header, 1, :), test_res)

    #***********************************

   if verbose

    if s == "H"
        print( "\n* Null Hypothesis: v is "); printstyled("normal"; bold=true, color=:yellow); print(" AND u is"); printstyled(" half-normal"; bold=true, color=:yellow); println(".")
    elseif s == "E"
        print( "\n* Null Hypothesis: v is "); printstyled("normal"; bold=true, color=:yellow); print(" AND u is"); printstyled(" exponential"; bold=true, color=:yellow); println(".")
    end
    println(" ")

    println("  Test Statistics (χ² distribution)")

    #* make PrettyTable for sine and cosine test

    pretty_table(test_res;
                 header = test_header,
                 header_crayon = crayon"yellow bold",
                 formatters = ft_printf("%4.5f", 2:3))
    println("  Note: Chen and Wang (2012 EReviews) indicates that cosine test with ω=1 has good overall performance.\n")

    println("\n  Critical Values (χ²(1))")
    pretty_table([crit1 crit2 crit3];
                 header = ["1%","5%","10%"],
                 header_crayon = crayon"yellow bold",
                 formatters = ft_printf("%4.5f", 1:3))

  end # if_verbose

    #***** prepare for printing estimation results *******

    nofpar = length(_dicMoM[:frontier]) + 2
        var_list  = cat([String(i) for i in _dicMoM[:frontier]], ["σᵥ²", "σᵤ²"], dims=1)
        coef_list = cat(β_ols, [σᵥ², σᵤ²], dims=1)
        std_list  = cat(std_ols, [myvstd, myustd], dims=1)
        t_stats = cat(zeros(nofpar-2), ["n.a.","n.a."], dims=1)
        pv_list = cat(zeros(nofpar-2), ["n.a.","n.a."], dims=1)
        ci_down = zeros(nofpar)
        ci_up = zeros(nofpar)
        tt = cquantile(Normal(0,1), α/2)
        for i in 1:nofpar
            if i <= nofpar-2
                t_stats[i] = coef_list[i] / std_list[i]
                pv_list[i] = pvalue(TDist(sT - nofpar), t_stats[i]; tail=:both)
                ci_down[i] = coef_list[i] - tt * std_list[i]
                ci_up[i] = coef_list[i] + tt * std_list[i]
            else
                ci_down[i] = (sT-1)/cquantile(Chisq(sT-1),α/2) * coef_list[i]
                ci_up[i] = (sT-1)/cquantile(Chisq(sT-1),1-(α/2)) * coef_list[i]
            end
        end

    est_res = hcat(var_list, coef_list, std_list, t_stats, pv_list, ci_down, ci_up)
    est_header = ([" ", "Coef.", "Std. Err.", "z", "P>|z|", "$(Int(100-100α))%CI_l", "$(Int(100-100α))%CI_u"])

    est_res2 = vcat(reshape(est_header, 1, :), est_res)  # for save in dic

   #******* make PrettyTable for estimation results *******

   if verbose && !testonly

    println("\n\n* Method of Moments Estimates of the Model (Chen and Wang 2012 EReviews)")

    print("  ** Model type: "); printstyled(minfo1; color=:yellow); println()
    print("  ** The constant variable (for intercept) in the model: ");
          if  cnsnum == 0
            printstyled("n.a.\n"; color = :yellow)
          else
            printstyled("$(cnsname[1]).\n"; color = :yellow)
          end
    print("  ** Number of observations: "); printstyled(sT; color=:yellow); println()
    print("  ** Log-likelihood value: "); printstyled(round(llike; digits=5); color=:yellow); println()
    println()

    pretty_table(est_res;
                 header = est_header,
                 header_crayon = crayon"yellow bold",
                 formatters = ft_printf("%4.5f", 2:7))
    println("  Note: CI of σᵥ² and σᵤ² is calculated based on the χ² distribution.\n\n")

    # --- print additional information -----

    printstyled("***** Additional Information *********\n"; color=:cyan)

    print("* OLS (frontier-only) log-likelihood: "); printstyled(round(ll_ols; digits=5); color=:yellow); println("")
    print("* Skewness of OLS residuals: "); printstyled(round(sk_ols; digits=5); color=:yellow); println("")
    print("* The sample mean of the JLMS inefficiency index: "); printstyled(round(_jlmsM; digits=5); color=:yellow); println("")
    print("* The sample mean of the BC efficiency index: "); printstyled(round(_bcM; digits=5); color=:yellow); println("\n")

    println("* Use `name.list` to see saved results (keys and values) where `name` is the return specified in `name = sfmodel_MoMTest(..)`. Values may be retrieved using the keys. For instance:")
    println("   ** `name.MoM_loglikelihood`: the log-likelihood value of the model;")
    println("   ** `name.jlms`: Jondrow et al. (1982) inefficiency index;")
    println("   ** `name.bc`: Battese and Coelli (1988) efficiency index;")
    println("* Use `keys(name)` to see available keys.")

    printstyled("**************************************\n\n\n"; color=:cyan)

   end  # if_verbose

 #* ########### create a dictionary and make a tuple for return ########### *#

      _dicRES = OrderedDict{Symbol, Any}()
      _dicRES[:n_observations]  = sT
      _dicRES[:MoMtest]         = [test_res2][1]
      _dicRES[:table]           = [est_res2][1]
      _dicRES[:coeff]           = coef_list
      _dicRES[:std_err]         = std_list
      _dicRES[:jlms]            = _jlms
      _dicRES[:bc]              = _bc
      _dicRES[:MoM_loglikelihood] = llike #! hjw: see SFloglikefun
      _dicRES[:OLS_loglikelihood] = ll_ols
      _dicRES[:OLS_resid_skew]    = sk_ols
      _dicRES[:depvar]        = _dicMoM[:depvar]
      _dicRES[:frontier]      = _dicMoM[:frontier]
      _dicRES[:type]          = _dicMoM[:type]
      _dicRES[:dist]          = _dicMoM[:dist]
      _dicRES[:PorC]          = _porc

     #* ----- Delete optional keys that have value nothing,

     #* ----- Create a NamedTuple from the dic as the final output;
     #* -----     put the dic in the tuple.

         _ntRES = NamedTuple{Tuple(keys(_dicRES))}(values(_dicRES))
         _ntRES = (; _ntRES..., list    = _dicRES)

     #* ---- Create a gloal dictionary for sf_predict ----

  #* ############  make returns  ############ *#

      return _ntRES

end


text:
type AdversarialArm <: BanditArm
  t::Int64
  active_start::Int64
  active_end::Int64
end

function draw(arm::AdversarialArm)
  if arm.active_start <= arm.t <= arm.active_end
    1
  else
    0
  end
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
using ExactWrightFisher, Random, Test, Distributions
using FeynmanKacParticleFilters

@testset "Testing the WF particle filtering functions" begin

    Random.seed!(0);
    α_vec = [1.2, 1.4, 1.3]
    K = length(α_vec)
    Pop_size_WF3 = 10
    Nparts = 100
    time_grid_WF3 = range(0, stop = 1, length = 10)
    wfchain = Wright_Fisher_K_dim_exact_trajectory([0.2, 0.4, 0.4], time_grid_WF3, α_vec)
    wfobs_WF3 = [rand(Multinomial(Pop_size_WF3, wfchain[:,k])) for k in 1:size(wfchain,2)] |> l -> hcat(l...)
    data_WF3 = Dict(zip(time_grid_WF3 , [wfobs_WF3[:,t] for t in 1:size(wfobs_WF3,2)]))

    logGt = FeynmanKacParticleFilters.create_potential_functions(data_WF3, DualOptimalFiltering.multinomial_logpotential)
    Mt = DualOptimalFiltering.create_transition_kernels_WF(data_WF3, α_vec)
    RS(W) = rand(Categorical(W), length(W))

    # @test data_WF3[time_grid_WF3[2]] == [2, 4, 4]
    @test_nowarn Mt[time_grid_WF3[2]](α_vec / sum(α_vec))
    @test_nowarn pf_adaptive = FeynmanKacParticleFilters.generic_particle_filtering_adaptive_resampling_logweights(Mt, logGt, Nparts, RS)

    data_WF3 = Dict(zip(time_grid_WF3 , [wfobs_WF3[:,t:t]' for t in 1:size(wfobs_WF3,2)]))

    @test_nowarn FeynmanKacParticleFilters.create_potential_functions(data_WF3, DualOptimalFiltering.multinomial_logpotential)

end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using Fake
using Test

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


text:
using OpenSAFT, Test

@testset "ogSAFT" begin
    test_system = ogSAFT(["methanol"])
    Temp        = 298.15
    Vol         = 1e-4
    z           = [1.]
    #@test typeof(test_system) == OpenSAFT.ogSAFT
    @testset "params" begin
        @test test_system.params.segment.values[1] ≈ 1.6 rtol = 1e-8
        @test test_system.params.sigma.values[1] ≈ 3.203e-10 rtol = 1e-8
        @test test_system.params.epsilon.values[1] ≈ 163.15 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 2964 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 0.053 rtol = 1e-8
        @test test_system.params.n_sites[Set(["methanol"])]["e"] ≈ 1 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_seg(test_system, z, Vol, Temp) ≈ -0.540464960274975 rtol = 1e-8
        @test OpenSAFT.a_chain(test_system, z, Vol, Temp) ≈ -0.22445414690633522 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -4.695916382523908 rtol = 1e-8
    end
end

@testset "CKSAFT" begin
    test_system = CKSAFT(["carbon dioxide","2-propanol"])
    Temp        = 298.15
    Vol         = 1e-4
    z           = [0.5,0.5]
    @test typeof(test_system) == OpenSAFT.CKSAFT
    @testset "params" begin
        @test test_system.params.segment[Set(["2-propanol"])] ≈ 3.249 rtol = 1e-8
        @test test_system.params.sigma[Set(["2-propanol"])] ≈ 3.0430915331057595e-10 rtol = 1e-8
        @test test_system.params.sigma[Set(["carbon dioxide", "2-propanol"])] ≈ 3.1070590057871705e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["2-propanol"])] ≈ 202.94 rtol = 1e-8
        @test test_system.params.epsilon[Set(["carbon dioxide", "2-propanol"])] ≈ 198.51779876711507 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["2-propanol"]), "H"), (Set(["2-propanol"]), "e")])] ≈ 2670.0 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["2-propanol"]), "H"), (Set(["2-propanol"]), "e")])] ≈ 0.021 rtol = 1e-8
        @test test_system.params.n_sites[Set(["2-propanol"])]["e"] ≈ 1 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_seg(test_system, z, Vol, Temp) ≈ -1.2395529662948277 rtol = 1e-8
        @test OpenSAFT.a_chain(test_system, z, Vol, Temp) ≈ -0.7747586154084931 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -1.7937079004096872 rtol = 1e-8
    end
end

@testset "SAFTVRSW" begin
    test_system = system(["water","ethane"], "SAFTVRSW")
    Temp        = 298.15
    Vol         = 1e-4
    z           = create_z(test_system, [0.5, 0.5])
    @test typeof(test_system) == OpenSAFT.SAFTVRSW
    @testset "params" begin
        @test test_system.params.segment[Set(["ethane"])] ≈ 1.33 rtol = 1e-8
        @test test_system.params.sigma[Set(["water"])] ≈ 3.036e-10 rtol = 1e-8
        @test test_system.params.sigma[Set(["water","ethane"])] ≈ 3.6345e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["water"])] ≈ 253.3 rtol = 1e-8
        @test test_system.params.epsilon[Set(["water","ethane"])] ≈ 238.6248939234966 rtol = 1e-8
        @test test_system.params.lambda[Set(["water"])] ≈ 1.8 rtol = 1e-8
        @test test_system.params.lambda[Set(["water","ethane"])] ≈ 1.595600082542303 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["water"]), "H"), (Set(["water"]), "e")])] ≈ 1366.0 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["water"]), "H"), (Set(["water"]), "e")])] ≈ 1.028e-30 rtol = 1e-8
        @test test_system.params.n_sites[Set(["water"])]["e"] ≈ 2 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_mono(test_system, z, Vol, Temp) ≈ -1.4659622048306407 rtol = 1e-8
        @test OpenSAFT.a_chain(test_system, z, Vol, Temp) ≈ 0.022703334973543182 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -0.5091186885233859 rtol = 1e-8
    end
end

@testset "softSAFT" begin
    test_system = system(["methanol"], "softSAFT")
    Temp        = 298.15
    Vol         = 1e-4
    z           = create_z(test_system, [1.])
    @test typeof(test_system) == OpenSAFT.softSAFT
    @testset "params" begin
        @test test_system.params.segment[Set(["methanol"])] ≈ 1.491 rtol = 1e-8
        @test test_system.params.sigma[Set(["methanol"])] ≈ 3.375e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["methanol"])] ≈ 220.4 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 3213.0 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 4.847e-27 rtol = 1e-8
        @test test_system.params.n_sites[Set(["methanol"])]["e"] ≈ 1 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_LJ(test_system, z, Vol, Temp) ≈ -1.299697047509115 rtol = 1e-8
        @test OpenSAFT.a_chain(test_system, z, Vol, Temp) ≈ 0.33553325545339646 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -3.7490459421490447 rtol = 1e-8
    end
end

@testset "PCSAFT" begin
    test_system = system(["butane", "ethanol"], "PCSAFT")
    Temp        = 298.15
    Vol         = 1e-4
    z           = create_z(test_system, [0.5, 0.5])
    @test typeof(test_system) == OpenSAFT.PCSAFT
    @testset "params" begin
        @test test_system.params.segment[Set(["ethanol"])] ≈ 2.3827 rtol = 1e-8
        @test test_system.params.sigma[Set(["ethanol"])] ≈ 3.1771e-10 rtol = 1e-8
        @test test_system.params.sigma[Set(["butane", "ethanol"])] ≈ 3.44285e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["ethanol"])] ≈ 198.24 rtol = 1e-8
        @test test_system.params.epsilon[Set(["butane", "ethanol"])] ≈ 204.31368602729677 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["ethanol"]), "H"), (Set(["ethanol"]), "e")])] ≈ 2653.4 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["ethanol"]), "H"), (Set(["ethanol"]), "e")])] ≈ 0.032384 rtol = 1e-8
        @test test_system.params.n_sites[Set(["ethanol"])]["e"] ≈ 1 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_hc(test_system, z, Vol, Temp) ≈ 3.114823074155765 rtol = 1e-8
        @test OpenSAFT.a_disp(test_system, z, Vol, Temp) ≈ -6.090736624622517 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -2.121606453473655 rtol = 1e-8
    end
end

@testset "sPCSAFT" begin
    test_system = system(["pentane", "methanol"], "sPCSAFT")
    Temp        = 298.15
    Vol         = 1e-4
    z           = create_z(test_system, [0.5, 0.5])
    @test typeof(test_system) == OpenSAFT.sPCSAFT
    @testset "params" begin
        @test test_system.params.segment[Set(["methanol"])] ≈ 2.7921 rtol = 1e-8
        @test test_system.params.sigma[Set(["methanol"])] ≈ 2.651e-10 rtol = 1e-8
        @test test_system.params.sigma[Set(["pentane", "methanol"])] ≈ 3.21195e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["methanol"])] ≈ 186.6 rtol = 1e-8
        @test test_system.params.epsilon[Set(["pentane", "methanol"])] ≈ 207.7063311504972 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 2090.2 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 0.146 rtol = 1e-8
        @test test_system.params.n_sites[Set(["methanol"])]["e"] ≈ 1 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_hc(test_system, z, Vol, Temp) ≈ 3.568650770403549 rtol = 1e-8
        @test OpenSAFT.a_disp(test_system, z, Vol, Temp) ≈ -6.994181358803752 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -1.7525112985184315 rtol = 1e-8
    end
end

@testset "SAFTVRMie" begin
    test_system = system(["methanol"], "SAFTVRMie")
    Temp        = 298.15
    Vol         = 1e-4
    z           = create_z(test_system, [1.])
    @test typeof(test_system) == OpenSAFT.SAFTVRMie
    @testset "params" begin
        @test test_system.params.segment[Set(["methanol"])] ≈ 1.67034 rtol = 1e-8
        @test test_system.params.sigma[Set(["methanol"])] ≈ 3.2462e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["methanol"])] ≈ 307.69 rtol = 1e-8
        @test test_system.params.lambdaA[Set(["methanol"])] ≈ 6.0 rtol = 1e-8
        @test test_system.params.lambdaR[Set(["methanol"])] ≈ 19.235 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 2062.1 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["methanol"]), "H"), (Set(["methanol"]), "e")])] ≈ 1.0657e-28 rtol = 1e-8
        @test test_system.params.n_sites[Set(["methanol"])]["e"] ≈ 2 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_mono(test_system, z, Vol, Temp) ≈ -1.7176380421592838 rtol = 1e-8
        @test OpenSAFT.a_chain(test_system, z, Vol, Temp) ≈ -0.030259270092795967 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -3.1565551121889293 rtol = 1e-8
    end
end

@testset "SAFTVRQMie" begin
    test_system = system(["helium"], "SAFTVRQMie")
    Temp        = 298.15
    Vol         = 1e-4
    z           = create_z(test_system, [1.])
    @test typeof(test_system) == OpenSAFT.SAFTVRQMie
    @testset "params" begin
        @test test_system.params.segment[Set(["helium"])] ≈ 1. rtol = 1e-8
        @test test_system.params.sigma[Set(["helium"])] ≈ 2.549e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["helium"])] ≈ 10.952 rtol = 1e-8
        @test test_system.params.lambdaA[Set(["helium"])] ≈ 6.0 rtol = 1e-8
        @test test_system.params.lambdaR[Set(["helium"])] ≈ 13.0 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_mono(test_system, z, Vol, Temp) ≈ 0.12253715358076675 rtol = 1e-8
    end
end

@testset "SAFTgammaMie" begin
    test_system = system(gc("ethanol"), "SAFTgammaMie")
    Temp        = 298.15
    Vol         = exp10(-3.5)
    z           = create_z(test_system, [1.])
    @test typeof(test_system) == OpenSAFT.SAFTgammaMie
    @testset "params" begin
        @test test_system.params.segment[Set(["CH2OH"])] ≈ 2.0 rtol = 1e-8
        @test test_system.params.shapefactor[Set(["CH2OH"])] ≈ 0.58538 rtol = 1e-8
        @test test_system.params.sigma[Set(["CH2OH"])] ≈ 3.4054000000000005e-10 rtol = 1e-8
        @test test_system.params.sigma[Set(["CH2OH","CH3"])] ≈ 3.7413000000000003e-10 rtol = 1e-8
        @test test_system.params.epsilon[Set(["CH2OH"])] ≈ 407.22 rtol = 1e-8
        @test test_system.params.epsilon[Set(["CH2OH","CH3"])] ≈ 333.2 rtol = 1e-8
        @test test_system.params.lambda_a[Set(["CH2OH"])] ≈ 6.0 rtol = 1e-8
        @test test_system.params.lambda_a[Set(["CH2OH","CH3"])] ≈ 6.0 rtol = 1e-8
        @test test_system.params.lambda_r[Set(["CH2OH"])] ≈ 22.699 rtol = 1e-8
        @test test_system.params.lambda_r[Set(["CH2OH","CH3"])] ≈ 18.406912409694556 rtol = 1e-8
        @test test_system.params.epsilon_assoc[Set([(Set(["CH2OH"]), "H"), (Set(["CH2OH"]), "e1")])] ≈ 2097.9 rtol = 1e-8
        @test test_system.params.bond_vol[Set([(Set(["CH2OH"]), "H"), (Set(["CH2OH"]), "e1")])] ≈ 6.2309e-29 rtol = 1e-8
        @test test_system.params.n_sites[Set(["CH2OH"])]["e1"] ≈ 2 rtol = 1e-8
    end
    # We can try to extract specific terms, although I don't know how useful it will be
    @testset "equations" begin
        @test OpenSAFT.a_mono(test_system, z, Vol, Temp) ≈ -1.151043781769667 rtol = 1e-8
        @test OpenSAFT.a_chain(test_system, z, Vol, Temp) ≈ -0.1255227354789658 rtol = 1e-8
        @test OpenSAFT.a_assoc(test_system, z, Vol, Temp) ≈ -1.9386416653191778 rtol = 1e-8
    end
end


text:
"""
    SurfaceFluxes

## Interface
  - [`surface_conditions`](@ref) computes
    - Monin-Obukhov length
    - Potential temperature flux (if not given) using Monin-Obukhov theory
    - transport fluxes using Monin-Obukhov theory
    - friction velocity/temperature scale/tracer scales
    - exchange coefficients

## References

@article{nishizawa2018surface,
  title={A Surface Flux Scheme Based on the Monin-Obukhov Similarity for Finite Volume Models},
  author={Nishizawa, S and Kitamura, Y},
  journal={Journal of Advances in Modeling Earth Systems},
  volume={10},
  number={12},
  pages={3159--3175},
  year={2018},
  publisher={Wiley Online Library}
}

"""
module SurfaceFluxes

using NLsolve

using ..Thermodynamics
using DocStringExtensions
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet: molmass_ratio, grav
using CLIMAParameters.SubgridScale: von_karman_const
using StaticArrays

include("UniversalFunctions.jl")
using .UniversalFunctions
const UF = UniversalFunctions

abstract type SurfaceFluxesModel end

struct FVScheme end
struct DGScheme end

export surface_conditions, exchange_coefficients

"""
    SurfaceFluxConditions{FT}

Surface flux conditions, returned from `surface_conditions`.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SurfaceFluxConditions{FT}
    L_MO::FT
    VDSE_flux_star::FT
    flux::AbstractVector
    x_star::AbstractVector
    K_exchange::AbstractVector
end

function Base.show(io::IO, sfc::SurfaceFluxConditions)
    println(io, "----------------------- SurfaceFluxConditions")
    println(io, "L_MO           = ", sfc.L_MO)
    println(io, "VDSE_flux_star = ", sfc.VDSE_flux_star)
    println(io, "flux           = ", sfc.flux)
    println(io, "x_star         = ", sfc.x_star)
    println(io, "K_exchange     = ", sfc.K_exchange)
    println(io, "-----------------------")
end

"""
    surface_conditions
Surface conditions given
 - `x_initial` initial guess for solution (`L_MO, u_star, θ_star, ϕ_star, ...`)
 - `x_ave` volume-averaged value for variable `x`
 - `x_s` surface value for variable `x`
 - `z_0` roughness length for variable `x`
 - `F_exchange` flux at the top for variable `x`
 - `VDSE_scale` virtual dry static energy (i.e., basic potential temperature)
 - `Δz` layer thickness (not spatial discretization)
 - `z` coordinate axis
 - `a` free model parameter with prescribed value of 4.7
 - `VDSE_flux_star_given` potential temperature flux (optional)

If `VDSE_flux_star` is not given, then it is computed by iteration
of equations 3, 17, and 18 in Nishizawa2018.
"""
function surface_conditions(
    param_set::AbstractEarthParameterSet,
    x_initial::AbstractVector,
    x_ave::AbstractVector,
    x_s::AbstractVector,
    z_0::AbstractVector,
    F_exchange::AbstractVector,
    VDSE_scale::FT,
    qt_bar::FT,
    Δz::FT,
    z::FT,
    a::FT,
    scheme,
    VDSE_flux_star_given::Union{Nothing, FT} = nothing,
) where {FT <: AbstractFloat, AbstractEarthParameterSet}

    n_vars = length(x_initial) - 1
    @assert length(x_initial) == n_vars + 1
    @assert length(x_ave) == n_vars
    @assert length(x_s) == n_vars
    @assert length(z_0) == n_vars
    @assert length(F_exchange) == n_vars
    local sol
    u, θ, qt = x_initial[1], x_initial[2], x_initial[3]
    let param_set = param_set
        function f!(F, x_all)
            L_MO, x_vec = x_all[1], x_all[2:end]
            u, θ = x_vec[1], x_vec[2]
            VDSE_flux_star =
                VDSE_flux_star_given == nothing ? -u * θ : VDSE_flux_star_given
            flux = VDSE_flux_star
            F[1] = L_MO - monin_obukhov_length(param_set, u, VDSE_scale, flux)
            for i in 1:n_vars
                ϕ = x_vec[i]
                transport = i == 1 ? MomentumTransport() : HeatTransport()
                L_MO = monin_obukhov_length(param_set, u, VDSE_scale, flux)
                uf = Businger(param_set, L_MO)
                F[i + 1] =
                    ϕ - compute_physical_scale(
                        uf,
                        Δz,
                        z_0[i],
                        x_ave[i],
                        x_s[i],
                        transport,
                        scheme,
                    )
            end
        end
        sol = nlsolve(f!, x_initial, autodiff = :forward)
    end
    if converged(sol)
        L_MO, x_star = sol.zero[1], collect(sol.zero[2:end])
        u_star, θ_star = x_star[1], x_star[2]
    else
        #println("Warning: Unconverged Surface Fluxes")
        L_MO, x_star = sol.zero[1], sol.zero[2:end]
        u_star, θ_star = x_star[1], x_star[2]
    end

    _grav::FT = grav(param_set)
    _von_karman_const::FT = von_karman_const(param_set)
    VDSE_flux_star = -u_star^3 * VDSE_scale / (_von_karman_const * _grav * L_MO)
    flux = -u_star * x_star

    K_exchange = exchange_coefficients(
        param_set,
        z,
        F_exchange,
        a,
        x_star,
        VDSE_scale,
        L_MO,
    )

    C_exchange =
        get_flux_coefficients(param_set, z, a, x_star, VDSE_scale, L_MO, z_0)

    return SurfaceFluxConditions(L_MO, VDSE_flux_star, flux, x_star, C_exchange)
end

"""
    compute_physical_scale
    (3), (17), (18), (19)
"""
function compute_physical_scale(
    uf::AbstractUniversalFunction{FT},
    Δz,
    z_0,
    x_ave,
    x_s,
    transport,
    ::FVScheme,
) where {FT}
    _von_karman_const::FT = von_karman_const(uf.param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    R_z0 = 1 - z_0 / Δz
    temp1 = log(Δz / z_0)
    temp2 = -Psi(uf, Δz / uf.L, transport)
    temp3 = z_0 / Δz * Psi(uf, z_0 / uf.L, transport)
    temp4 = R_z0 * (psi(uf, z_0 / uf.L, transport) - 1)
    Σterms = temp1 + temp2 + temp3 + temp4
    return _π_group⁻¹ * _von_karman_const / Σterms * (x_ave - x_s)
end

function compute_physical_scale(
    uf::AbstractUniversalFunction{FT},
    Δz,
    z_0,
    x_ave,
    x_s,
    transport,
    ::DGScheme,
) where {FT}
    _von_karman_const::FT = von_karman_const(uf.param_set)
    _π_group = FT(UF.π_group(uf, transport))
    _π_group⁻¹ = (1 / _π_group)
    temp1 = log(Δz / z_0)
    temp2 = -psi(uf, Δz / uf.L, transport)
    temp3 = psi(uf, z_0 / uf.L, transport)
    Σterms = temp1 + temp2 + temp3
    return _π_group⁻¹ * _von_karman_const / Σterms * (x_ave - x_s)
end

### Generic terms

"""
    monin_obukhov_length(param_set, u, θ, flux)

Monin-Obukhov length. Eq. 3
"""
function monin_obukhov_length(
    param_set::AbstractEarthParameterSet,
    u_star,
    VDSE_scale,
    flux,
)
    FT = typeof(u_star)
    _grav::FT = grav(param_set)
    _von_karman_const::FT = von_karman_const(param_set)
    return -u_star^3 * VDSE_scale / (_von_karman_const * _grav * flux)
end

"""
    exchange_coefficients(z, F_exchange, a, x_star, VDSE_scale, L_MO)

Computes exchange transfer coefficients
  - `K_exchange` exchange coefficients
"""
function exchange_coefficients(
    param_set,
    z,
    F_exchange,
    a,
    x_star,
    VDSE_scale,
    L_MO,
)
    N = length(F_exchange)
    FT = typeof(z)
    _von_karman_const::FT = von_karman_const(param_set)
    K_exchange = zeros(length(x_star))
    uf = Businger(param_set, L_MO)
    for i in 1:N
        transport = i == 1 ? MomentumTransport() : HeatTransport()
        phi_t = phi(uf, z / L_MO, transport)
        _π_group = FT(UF.π_group(uf, transport))
        K_exchange[i] =
            -F_exchange[i] * _von_karman_const * z / _π_group /
            (x_star[i] * phi_t) # Eq. 19 in
    end
    return K_exchange
end

function get_flux_coefficients(param_set, z, a, x_star, VDSE_scale, L_MO, z0)
    N = length(x_star)
    FT = typeof(z)
    _von_karman_const::FT = von_karman_const(param_set)
    C = zeros(length(x_star))
    uf = Businger(param_set, L_MO)
    psi_m = psi(uf, z / L_MO, MomentumTransport())
    psi_h = psi(uf, z / L_MO, HeatTransport())
    for i in 1:N
        logζ_ψ_m = log(z / z0[i] - psi_m)
        logζ_ψ_h = log(z / z0[i] - psi_h)
        if i == 1
            C[i] = _von_karman_const^2 / logζ_ψ_m^2
        else
            C[i] = _von_karman_const^2 / logζ_ψ_m / logζ_ψ_h
        end
    end
    return (C)

end

end # SurfaceFluxes module

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using LogicCircuits
using LogicCircuits: JuiceTransformer, dimacs_comments, zoo_version

using LazyArtifacts
using Lerche: Lerche, Lark, Transformer, @rule, @inline_rule

include("jpc_io.jl")
include("psdd_io.jl")
include("spn_io.jl")
include("clt_io.jl")
include("ensemble_io.jl")
include("plot.jl")

# if no logic circuit file format is given on read, infer file format from extension

function file2pcformat(file)
    if endswith(file,".gz")
        # duplicate code from `LogicCircuits.file2logicformat` -- not sure how to make this nicer
        file_inner, _ = splitext(file)
        format_inner = file2pcformat(file_inner)
        GzipFormat(format_inner)
    elseif endswith(file,".jpc")
        JpcFormat()
    elseif endswith(file,".psdd")
        PsddFormat()
    elseif endswith(file,".spn")
        SpnFormat()
    else
        # try a logic circuit format
        LogicCircuits.file2logicformat(file)
    end
end

"""
    Base.read(file::AbstractString, ::Type{C}) where C <: ProbCircuit

Reads circuit from file; uses extension to detect format type, for example ".psdd" for PSDDs.
"""
Base.read(file::AbstractString, ::Type{C}) where C <: ProbCircuit =
    read(file, C, file2pcformat(file))

Base.read(files::Tuple{AbstractString,AbstractString}, ::Type{C}) where C <: StructProbCircuit =
    read(files, C, (file2pcformat(files[1]), VtreeFormat()))

"""
    Base.write(file::AbstractString, circuit::ProbCircuit)

Writes circuit to file; uses file name extention to detect file format.
"""
Base.write(file::AbstractString, circuit::ProbCircuit) =
    write(file, circuit, file2pcformat(file))

"""
    Base.write(files::Tuple{AbstractString,AbstractString}, circuit::StructProbCircuit)

Saves circuit and vtree to file.
"""
Base.write(files::Tuple{AbstractString,AbstractString},
           circuit::StructProbCircuit) =
    write(files, circuit, (file2pcformat(files[1]), VtreeFormat()))

#  when asked to parse/read as `ProbCircuit`, default to `PlainProbCircuit`

Base.parse(::Type{ProbCircuit}, args...) =
    parse(PlainProbCircuit, args...)

Base.read(io::IO, ::Type{ProbCircuit},  args...) =
    read(io, PlainProbCircuit,  args...)

Base.read(io::IO, ::Type{ProbCircuit}, f::GzipFormat) =
    # avoid method ambiguity
    read(io, PlainProbCircuit,  f)

# copy read/write API for tuples of files

function Base.read(files::Tuple{AbstractString, AbstractString}, ::Type{C}, args...) where C <: StructProbCircuit
    open(files[1]) do io1
        open(files[2]) do io2
            read((io1, io2), C, args...)
        end
    end
end

function Base.write(files::Tuple{AbstractString,AbstractString},
                    circuit::StructProbCircuit, args...)
    open(files[1], "w") do io1
        open(files[2], "w") do io2
            write((io1, io2), circuit, args...)
        end
    end
end


text:
include("parserMomkpZL.jl")
include("parserMomkpPG.jl")
include("displayGraphic.jl")
include("vOptMomkp.jl")
include("VNS.jl")
include("MOKP.jl")
include("tools.jl")
include("nsga_2.jl")
using CPLEX
function main(name::String, flag::Int, exact::Bool)
    verbose = false
    nGenerations = 10000000
    nIndividus = 100 #size(prob.C)[2]
    probCrossover = 0.9
    probMutation = 0.1
    timeAllowed = 10
    if flag == 1
        fname = "instancesZL/" * name
        momkpZL = readInstanceMOMKPformatZL(verbose, fname)
    else
        fname = "instancesPG/" * name
        momkpZL = readInstanceMOMKPformatPG(verbose, fname)
    end
    prob = _bi01IP(momkpZL.P[1:2, :], momkpZL.W, momkpZL.ω)
    #println("--------------------------------------")
    #println(name)
    if exact
        if !isfile("YN/" * name * ".DAT")
            # we must solve with cplex first
            solverSelected = CPLEX.Optimizer
            YN, XE = vSolveBi01IP(solverSelected, prob.C, prob.A, prob.b)
            Yn_int = Vector{Vector{Int}}(undef, length(YN))
            for j = 1:length(YN)
                Yn_int[j] = [0, 0]
                Yn_int[j][1] = round(Int, YN[j][1])
                Yn_int[j][2] = round(Int, YN[j][2])
            end
            TOOL_write_YN(Yn_int, name)
            YN = Yn_int
        else
            YN = TOOL_read_YN(name)
        end
        # displayGraphics(fname, YN, "black", L"y \in Y_N")
    end
    #VNS
    #println("VNS")
    #print("temps gen pop : ")
    pop = initPopEpsilon(10, prob)
    VNS_sol, real_time = GVNS(pop, 2, timeAllowed, 1, prob, 2, 2)
    YN_VNS = get_YN_VNS(VNS_sol)
    # displayGraphics(fname, YN_VNS, "green", L"y \in Y_{VNS}")
    # NSGA-II
    #println("NSGA")
    NSGA_sol = nsga2(nIndividus, nGenerations, probCrossover, probMutation, prob, real_time)
    YN_NSGA = get_YN_nsga(NSGA_sol)

    # displayGraphics(fname, YN_NSGA, "red", L"y \in Y_{NSGA}")
    if exact
        #println("hypervolumes (VNS,exact,ratio) ", compareHV(YN_VNS, YN))
        #println("hypervolumes (NSGA-II,exact,ratio) ", compareHV(YN_NSGA, YN))
    else
        #println("hypervolumes (VNS,NSGA-II,ratio) ", compareHV(YN_VNS, YN_NSGA))
    end
    v1 = round(compareHV(YN_VNS, YN)[3], digits = 4)
    v2 = round(compareHV(YN_NSGA, YN)[3], digits = 4)
    real_time = round(digits = 1, real_time)
    #println("(|YN_VNS|,|YN_NSGA|,|YN|)", (length(YN_VNS), length(YN_NSGA), exact ? length(YN) : 0))
    println("\\hline")
    println("$name & $(length(YN)) & $(length(YN_VNS)) & ", v1, " & $(length(YN_NSGA)) & ", v2, " & ", real_time, "\\\\")
    #println("--------------------------------------")
end

#main("knapsack.100.2", 1, true)
#main("knapsack.100.3", 1, true)
#main("knapsack.100.4", 1, true)
#main("knapsack.250.2", 1, true)
#main("knapsack.250.3", 1, true)
#main("knapsack.250.4",1,true)
#main("knapsack.500.2", 1, true)
#main("knapsack.500.3",1,true)
#main("knapsack.500.4",1,true)
#main("knapsack.750.2",1,true)
#main("knapsack.750.3",1,true)
#main("knapsack.750.4",1,true)
#main("set1/ZL28.DAT", 2, true)
#main("set1/ZL100.DAT", 2, true)
#main("set1/ZL105.DAT", 2, true)
#main("set1/ZL250.DAT", 2, true)
#main("set1/ZL500.DAT",2,true)
#main("set1/ZL750.DAT",2,true)
#main("set2/A1.DAT", 2, true)
#main("set2/A2.DAT", 2, true)
#main("set2/A3.DAT", 2, true)
#main("set2/A4.DAT", 2, true)
#main("set2/D1.DAT", 2, true)
#main("set2/D2.DAT", 2, true)
#main("set2/D3.DAT", 2, true)
#main("set2/D4.DAT", 2, true)
#main("set2/kp28.DAT", 2, true)
#main("set2/kp28-2.DAT", 2, true)
main("set2/kp28c1W-c2ZL.DAT", 2, true)
main("set2/kp28cW-WZL.DAT", 2, true)
main("set2/kp28W-Perm.DAT", 2, true)
main("set2/kp28W-ZL.DAT", 2, true)
main("set2/kp28W.DAT", 2, true)
main("set3/W7B1-rnd1-1800.DAT", 2, true)
main("set3/W7B1-rnd1-3000.DAT", 2, true)
main("set3/W7B1-tube1-1800.DAT", 2, true)
main("set3/W7B1-tube1-3000.DAT", 2, true)
main("set3/W7B1-tube1-asyn.DAT", 2, true)
main("set3/W7B1-tube2-1800.DAT", 2, true)
main("set3/Wcollage-tube.DAT", 2, true)


text:
#  this function takes in a fits_name and returns
#  a dictionary with the exact same member names.

# TO FINISH:   Implement gcount loop to create multi_dimensional
#              return dict objec.
#              multi_D_dict = fill(return_dictionary, gcount_value)
#
#              Loop through return types as well to define the return type.
#              Currently, everything is Float32.
#
#              Figure out how to dereference a PyObject chararry in order
#              to pass string arrays.
#              It might be a PyCall issue.
#
#
# TEST CASE is 1D
# fits_name = "/Users/alex/Desktop/make_pros/old/completeness_test_results.fits"
# fits_name = "/Users/alex/Desktop/final_files/s0927/s0927_analysis.fits"
import FITSIO

function open_fits(fits_name)

include("/Users/alex/CODES/Julia_Custom_jl/alex/src/get_fits_members.jl")
include("/Users/alex/CODES/Julia_Custom_jl/alex/src/remove_char.jl")

# open the PyObject
data=FITSIO.fits_open_data(fits_name)

# get the fits members
members=get_fits_members(fits_name)

# pass the ttype_names to itterate over
# in order to create the custom return object
iters = members[:ttype_names]

#######################################################################################################
#######################################################################################################
#                                                                                ############
#                                                                                ############
#                                   make return object                                     ############
#                                                                                ############
#                                                                                ############
#######################################################################################################
#######################################################################################################
                   new_dict = Dict{Symbol,Any}()                                           ############
                        println("making new dictionary with members:")                     ############
                                                                                ############
                                                                                ############
                        # create new dictionary with each field empty                      ############
                        for sym in iters                                                   ############
                        println(sym)                                                       ############

                        # temp=Array(Float64, 6898)
                        # fits_read_col(data, Float32, 11, 1, 1, temp)
                        fits_read_col(data, String, 11, 1, 1, temp)

# function fits_read_col{T}(f::FITSFile, ::Type{T}, colnum::Int64, firstrow::Int64, firstelem::Int64, data::Array{T})

                         new_dict[symbol(sym)] = temp                                      ############

                        end                                                                ############
#######################################################################################################
#######################################################################################################

return new_dict

end


text:
@assert supertype(Int32) == Signed
@assert supertype(Signed) == Integer
@assert supertype(Integer) == Real
@assert supertype(Symbol) == Any
@assert supertype(DataType) == Type
@assert supertype(Union) == Type
@assert supertype(UnionAll) == Type

@assert Float64 <: Real
@assert isa(Array, Type)

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
struct Estep
    N
    method
    quadpts
    atol
    dentype
    Estep(;N = 500, method = :BAEM, quadpts = 21, atol = 1e-5, dentype = :DN) = new(N, method, quadpts, atol, dentype)
end

struct Mstep
    N
    method # :Newton or BFGS()
    atol
    Mstep(;N = 10, method = BFGS(), atol = 1e-5) = new(N, method, atol)
end

"""
Run EM cycle to estimate the model item parameters from data and Q matrix.

# Graded Response Model (Samejima, 1969) -> `graded()`

# 3 parameter logistic model -> `guessing()`

# Logistic exponential model -> `LPE()`

"""
function estimate(D::AbstractDataFrame, ndims::Int64, model::IRTmodel; group = fill(1, size(D, 1)), E::Estep = Estep(), M::Mstep = Mstep())
    if typeof(model) <: IRTmodel
        model = fill(model, size(D, 2))
    end
    Q = explanatoryMIRT(ndims, size(D, 2))
    U = combine(groupby([D DataFrame(__group = group)], All()), nrow => :__count)
    _estimate(U, Q, model, E::Estep, M::Mstep)
end

# Internal function to run EM cycle for the general models.
# This function assume to be passed the multidimensional Q matrix.
function _estimate(D, Q, model, E::Estep, M::Mstep)

    # Initialization step
    U = D[:, Not([:__count, :__group])]
    F = D[:, :__count]
    G = D[:, :__group]
    U = mapcols(c -> rescore(c), U) # rescoreは，オプションで選択可能にする。エラー検知もできるとよい。
    pind, G = takepersonindex(G)
    nind = Dict(zip(axes(U, 1), F))
    iind = [takeitemindex(@view U[i,:]) for i in axes(U, 1)]
    P, I = size(U)
    lnL = initializematrix(P, E.quadpts^size(Q, 2))
    post = initializematrix(P, E.quadpts^size(Q, 2))
    pars = [init(model[i], Q[i,:], D[:, Not([:__count, :__group])][:, i]) for i in 1:I]
    [setgroupindex!(pars[j], U[:,j], G) for j in 1:I]
    ghsize = length(unique(G))
    gh = [initializenode(size(Q, 2), E.quadpts) for i in 1:ghsize]
    bgh = initializenode(size(Q, 2), E.quadpts) # for reference
    pmoments = [Moments(size(Q, 2)) for _ in unique(G)]
    # Start EM Cycle
    if E.method == :BAEM
        for iter in 1:E.N
            print("EMcycle = ", @sprintf "%4.0f : " iter)
            old = copy.(pars)
            if any(checitemparameters.(pars))
                error("\nStop the estimation.\n")
            end
            N, r = BAEMcycle(lnL, post, pind, nind, iind, gh, bgh, U, pars, Q; Mstep = M)
            updateDN!(pind, nind, post, gh, pmoments)
            @fastmath @inbounds [rescale!(pars[i], pmoments[1].μ, pmoments[1].Σ, Q[i,:]) for i in axes(pars, 1)]
            flag = checkconvergence(old, pars; atol = E.atol)
            print(@sprintf " MOMENTS μ = %1.3f, σ² = %1.3f" pmoments[1].μ[1] pmoments[1].Σ[1])
            if all(flag)
                print("\nConverged.\n\n")
                break
            end
            print("\r")
        end
    end
    # Report estimated parameters, SE, and the other information.
    return pars, pmoments
end


text:
using LaserTag
using Plots
using POMDPToolbox
using ContinuousPOMDPTreeSearchExperiments
using Plots
using POMCPOW
using POMDPs

ro = MoveTowards()
solver = POMCPOWSolver(tree_queries=10_000,
                       criterion=MaxUCB(20.0),
                       final_criterion=MaxTries(),
                       max_depth=100,
                       enable_action_pw=false,
                       k_observation=4.0,
                       alpha_observation=1/20,
                       estimate_value=FORollout(ro),
                       check_repeat_act=false,
                       check_repeat_obs=false,
                       init_N=InevitableInit(),
                       init_V=InevitableInit(),
                       rng=MersenneTwister(13)
                      )

p = gen_lasertag(rng=MersenneTwister(4))

policy = solve(solver, p)

action(policy, initial_state_distribution(p))

# blink(policy)


text:
function BigFraction(arg0::BigInteger)
    return BigFraction((BigInteger,), arg0)
end

function BigFraction(arg0::BigInteger, arg1::BigInteger)
    return BigFraction((BigInteger, BigInteger), arg0, arg1)
end

function BigFraction(arg0::jdouble)
    return BigFraction((jdouble,), arg0)
end

function BigFraction(arg0::jdouble, arg1::jdouble, arg2::jint)
    return BigFraction((jdouble, jdouble, jint), arg0, arg1, arg2)
end

function BigFraction(arg0::jdouble, arg1::jint)
    return BigFraction((jdouble, jint), arg0, arg1)
end

function BigFraction(arg0::jint)
    return BigFraction((jint,), arg0)
end

function BigFraction(arg0::jint, arg1::jint)
    return BigFraction((jint, jint), arg0, arg1)
end

function BigFraction(arg0::jlong)
    return BigFraction((jlong,), arg0)
end

function BigFraction(arg0::jlong, arg1::jlong)
    return BigFraction((jlong, jlong), arg0, arg1)
end

function abs(obj::BigFraction)
    return jcall(obj, "abs", BigFraction, ())
end

function add(obj::BigFraction, arg0::BigFraction)
    return jcall(obj, "add", BigFraction, (BigFraction,), arg0)
end

function add(obj::BigFraction, arg0::BigInteger)
    return jcall(obj, "add", BigFraction, (BigInteger,), arg0)
end

function add(obj::BigFraction, arg0::jint)
    return jcall(obj, "add", BigFraction, (jint,), arg0)
end

function add(obj::BigFraction, arg0::jlong)
    return jcall(obj, "add", BigFraction, (jlong,), arg0)
end

function big_decimal_value(obj::BigFraction)
    return jcall(obj, "bigDecimalValue", BigDecimal, ())
end

function big_decimal_value(obj::BigFraction, arg0::jint)
    return jcall(obj, "bigDecimalValue", BigDecimal, (jint,), arg0)
end

function big_decimal_value(obj::BigFraction, arg0::jint, arg1::jint)
    return jcall(obj, "bigDecimalValue", BigDecimal, (jint, jint), arg0, arg1)
end

function byte_value(obj::Number)
    return jcall(obj, "byteValue", jbyte, ())
end

function compare_to(obj::BigFraction, arg0::BigFraction)
    return jcall(obj, "compareTo", jint, (BigFraction,), arg0)
end

function divide(obj::BigFraction, arg0::BigFraction)
    return jcall(obj, "divide", BigFraction, (BigFraction,), arg0)
end

function divide(obj::BigFraction, arg0::BigInteger)
    return jcall(obj, "divide", BigFraction, (BigInteger,), arg0)
end

function divide(obj::BigFraction, arg0::jint)
    return jcall(obj, "divide", BigFraction, (jint,), arg0)
end

function divide(obj::BigFraction, arg0::jlong)
    return jcall(obj, "divide", BigFraction, (jlong,), arg0)
end

function double_value(obj::BigFraction)
    return jcall(obj, "doubleValue", jdouble, ())
end

function equals(obj::BigFraction, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function float_value(obj::BigFraction)
    return jcall(obj, "floatValue", jfloat, ())
end

function get_denominator(obj::BigFraction)
    return jcall(obj, "getDenominator", BigInteger, ())
end

function get_denominator_as_int(obj::BigFraction)
    return jcall(obj, "getDenominatorAsInt", jint, ())
end

function get_denominator_as_long(obj::BigFraction)
    return jcall(obj, "getDenominatorAsLong", jlong, ())
end

function get_field(obj::BigFraction)
    return jcall(obj, "getField", BigFractionField, ())
end

function get_numerator(obj::BigFraction)
    return jcall(obj, "getNumerator", BigInteger, ())
end

function get_numerator_as_int(obj::BigFraction)
    return jcall(obj, "getNumeratorAsInt", jint, ())
end

function get_numerator_as_long(obj::BigFraction)
    return jcall(obj, "getNumeratorAsLong", jlong, ())
end

function get_reduced_fraction(::Type{BigFraction}, arg0::jint, arg1::jint)
    return jcall(BigFraction, "getReducedFraction", BigFraction, (jint, jint), arg0, arg1)
end

function hash_code(obj::BigFraction)
    return jcall(obj, "hashCode", jint, ())
end

function int_value(obj::BigFraction)
    return jcall(obj, "intValue", jint, ())
end

function long_value(obj::BigFraction)
    return jcall(obj, "longValue", jlong, ())
end

function multiply(obj::BigFraction, arg0::BigFraction)
    return jcall(obj, "multiply", BigFraction, (BigFraction,), arg0)
end

function multiply(obj::BigFraction, arg0::BigInteger)
    return jcall(obj, "multiply", BigFraction, (BigInteger,), arg0)
end

function multiply(obj::BigFraction, arg0::jint)
    return jcall(obj, "multiply", BigFraction, (jint,), arg0)
end

function multiply(obj::BigFraction, arg0::jlong)
    return jcall(obj, "multiply", BigFraction, (jlong,), arg0)
end

function negate(obj::BigFraction)
    return jcall(obj, "negate", BigFraction, ())
end

function percentage_value(obj::BigFraction)
    return jcall(obj, "percentageValue", jdouble, ())
end

function pow(obj::BigFraction, arg0::BigInteger)
    return jcall(obj, "pow", BigFraction, (BigInteger,), arg0)
end

function pow(obj::BigFraction, arg0::jdouble)
    return jcall(obj, "pow", jdouble, (jdouble,), arg0)
end

function pow(obj::BigFraction, arg0::jint)
    return jcall(obj, "pow", BigFraction, (jint,), arg0)
end

function pow(obj::BigFraction, arg0::jlong)
    return jcall(obj, "pow", BigFraction, (jlong,), arg0)
end

function reciprocal(obj::BigFraction)
    return jcall(obj, "reciprocal", BigFraction, ())
end

function reduce(obj::BigFraction)
    return jcall(obj, "reduce", BigFraction, ())
end

function short_value(obj::Number)
    return jcall(obj, "shortValue", jshort, ())
end

function signum(obj::BigFraction)
    return jcall(obj, "signum", jint, ())
end

function subtract(obj::BigFraction, arg0::BigFraction)
    return jcall(obj, "subtract", BigFraction, (BigFraction,), arg0)
end

function subtract(obj::BigFraction, arg0::BigInteger)
    return jcall(obj, "subtract", BigFraction, (BigInteger,), arg0)
end

function subtract(obj::BigFraction, arg0::jint)
    return jcall(obj, "subtract", BigFraction, (jint,), arg0)
end

function subtract(obj::BigFraction, arg0::jlong)
    return jcall(obj, "subtract", BigFraction, (jlong,), arg0)
end

function to_string(obj::BigFraction)
    return jcall(obj, "toString", JString, ())
end


text:
function input_ui()

    o1 = Observable("")
    graph = vbox(node(:div, "Raster resistance map or network/graph: ",
                      attributes = Dict(:style => "margin-top: 12px")),
                filepicker(value = o1))

    o2 = Observable(false)
    graph_is_res = checkbox("Data represents resistances instead of conductances",
                            value = o2)

    input = vbox(graph,
                 vskip(1em),
                 graph_is_res)

	# Get file path
	s = Scope()
	s.dom = input
    @private s["filepath"] = o1
    @private s["check"] = o2
    s
end

function pairwise_input_ui()

    o1 = Observable("")
    additional_input = node(:div, tachyons_css, "Pairwise Mode Options") |> class"f4 lh-title"
    focal = vbox(node(:div, "Select focal node locations from file: ",
                      attributes = Dict(:style => "margin-top: 12px")),
                 filepicker(value = o1))

    pair = vbox(additional_input,
                 focal)

    s = Scope()
    s.dom = pair
    @private s["focal"] = o1

    s
end

node


text:
module  TestUtils

using PopGenCore
using PopGenSims
using Test

cats = @nancycats ;
cats2 = copy(cats)

@testset "append PopData" begin
    tmp = append(cats, cats2)
    @test tmp.metadata.samples == 2 * cats.metadata.samples
    @test length(tmp.genodata.name) == 2 * length(cats.genodata.name)
    @test tmp isa PopData
    append!(cats2, cats)
    @test cats2.metadata.samples == 2 * cats.metadata.samples
    @test length(cats2.genodata.name) == 2 * length(cats.genodata.name)
    @test cats2 isa PopData
end

@testset "allele pool" begin
    @test PopGenSims.allele_pool(cats.genodata.genotype[1:30]) |> length == 56
    @test PopGenSims.allele_pool(cats.genodata.genotype[1:30]) |> typeof <: NTuple
    a,b = PopGenSims.allele_pool(cats)
    @test eltype(a) == String
    @test length(a) == 9
    @test typeof(b) <: Dict
    @test length(b) == length(a)
    c = PopGenSims.simulate_sample(b, a, ploidy = 2)
    @test length(c) == length(a) == length(b)
    @test all(length.(c) .== 2)
    c3 = PopGenSims.simulate_sample(b, a, ploidy = 3)
    @test all(length.(c3) .== 3)
end

end # module


text:
#####
##### `LAPACK.trsyl!`
#####

function ChainRules.frule(
    (_, _, _, ΔA, ΔB, ΔC),
    ::typeof(LAPACK.trsyl!),
    transa::AbstractChar,
    transb::AbstractChar,
    A::AbstractMatrix{T},
    B::AbstractMatrix{T},
    C::AbstractMatrix{T},
    isgn::Int,
) where {T<:BlasFloat}
    C, scale = LAPACK.trsyl!(transa, transb, A, B, C, isgn)
    Y = (C, scale)
    ΔAtrans = transa === 'T' ? transpose(ΔA) : (transa === 'C' ? ΔA' : ΔA)
    ΔBtrans = transb === 'T' ? transpose(ΔB) : (transb === 'C' ? ΔB' : ΔB)
    mul!(ΔC, ΔAtrans, C, -1, scale)
    mul!(ΔC, C, ΔBtrans, -isgn, true)
    ΔC, scale2 = LAPACK.trsyl!(transa, transb, A, B, ΔC, isgn)
    rmul!(ΔC, inv(scale2))
    ∂Y = Composite{typeof(Y)}(ΔC, Zero())
    return Y, ∂Y
end


text:
"""
Creates a switched system from a multitude of single-mode systems.

# Arguments
- `system_dirs::Array{String,1}`: directories containing the results produced by single mode systems.
- `exp_dir::String`: location to create the switched-system directory.
"""
function create_switched_system(system_dirs::Array{String,1}, exp_dir::String)

    # Create a new system dir with the tags
    new_dirname = "$exp_dir/switched-system"
    !isdir(new_dirname) && mkpath(new_dirname)

    # Copy the region info from the first system
    first_dir = system_dirs[1]
    cp("$first_dir/regions.bson", "$new_dirname/regions.bson", force=true)

    # Load all of the transition matrices into memory
    minPr_mats = []
    maxPr_mats = []

    for system_dir in system_dirs
        res_mats = load_transition_matrices("$system_dir/transition_mats.bson")
        push!(minPr_mats, res_mats["minPr"])
        push!(maxPr_mats, res_mats["maxPr"])
    end

    mm_minPr, mm_maxPr = create_switched_system_matrix(minPr_mats, maxPr_mats)

    # Save and return directory
    bson("$new_dirname/transition_mats.bson", Dict("minPr" => mm_minPr, "maxPr" => mm_maxPr))

    return new_dirname, mm_minPr, mm_maxPr
end

"""
Creates a switched system matrix from a list of min and max probabilities of transition.
"""
function create_switched_system_matrix(minPr_mats, maxPr_mats)
    num_modes = length(minPr_mats)
    mat_size = size(minPr_mats[1], 1)
    mm_minPr = hcat([[0.0 for i=1:mat_size*num_modes] for j=1:mat_size]...)
    mm_maxPr = hcat([[0.0 for i=1:mat_size*num_modes] for j=1:mat_size]...)

    for i=1:num_modes
        mm_minPr[i:num_modes:end, :] = minPr_mats[i]
        mm_maxPr[i:num_modes:end, :] = maxPr_mats[i]
    end

    return mm_minPr, mm_maxPr
end


text:
"docstring"
function function_name(arg)
    for value in 0:arg
        if value > 3
            return true
        elseif value < 4
            return nothing
        else
            return false
        end
    end
end


text:
module TestExpressionsFindColumnNames

import Test: @testset, @test
import Volcanito: find_column_names, ColumnName

test_cases = (
    (
        input = :(a + b + sin(c)),
        output = (
            ColumnName(:a, false),
            ColumnName(:b, false),
            ColumnName(:c, false),
        ),
    ),
    (
        input = :(x + sin(y) + cos(z^2) + $(Expr(:$, :a))),
        output = (
            ColumnName(:x, false),
            ColumnName(:y, false),
            ColumnName(:z, false),
        ),
    ),
    (
        input = :(1 + exp(x) + y),
        output = (
            ColumnName(:x, false),
            ColumnName(:y, false),
        ),
    ),
    (
        input = :(1 + exp(x) + y + $(Expr(:$, :z))),
        output = (
            ColumnName(:x, false),
            ColumnName(:y, false),
        ),
    ),
    (
        input = :(1 + exp(x, e = z) + y + $(Expr(:$, :z))),
        output = (
            ColumnName(:x, false),
            ColumnName(:y, false),
            ColumnName(:z, false),
        ),
    ),
)

@testset "find_column_names" begin
    for case in test_cases
        @test sort([find_column_names(case.input)...]) == sort([case.output...])
    end
end

end


text:
module ConvNets

end