text:
m = InfiniteModel();
@infinite_parameter(m, t in [-Inf, Inf])
@infinite_parameter(m, x[1:2] in [-Inf, Inf])

@testset "Default Monte Carlo Sampling" begin
    # test univariate Monte Carlo sampling
    @testset "mc_sampling (univariate)" begin
        (supports, weights) = generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, 0., 2., Val(mc_sampling))
        @test length(supports) == 5
        @test length(weights) == 5
        @test weights == 0.4 * ones(5)
        @test all(supports .>= 0.)
        @test all(supports .<= 2.)
        (supports, weights) = generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, -Inf, 0, Val(mc_sampling))
        @test length(supports) == 5
        @test length(weights) == 5
        (supports, weights) = generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, 0, Inf, Val(mc_sampling))
        @test length(supports) == 5
        @test length(weights) == 5
    end

    # test multivariate Monte Carlo sampling
    @testset "mc_sampling (multivariate)" begin
        lb = convert(JuMPC.SparseAxisArray, [0., 0.])
        ub = convert(JuMPC.SparseAxisArray, [1., 4.])
        (supports, weights) = generate_supports_and_coeffs(
                                  IntervalSet(-Inf, Inf), x, 5, lb, ub, Val(mc_sampling))
        @test length(supports) == 5
        @test length(weights) == 5
        @test weights == 0.8 * ones(5)
        @test all([i >= lb for i in supports])
        @test all([i <= ub for i in supports])
    end
end

@testset "Quadrature Methods" begin
    # test Gauss-Legendre method for bounded interval
    @testset "gauss_legendre" begin
        (supports, weights) = generate_supports_and_coeffs(
                               IntervalSet(-Inf, Inf), t, 5, 1., 5., Val(gauss_legendre))
        (expect_supports, expect_weights) = FGQ.gausslegendre(5)
        expect_supports = expect_supports * 2 .+ 3
        expect_weights = expect_weights * 2
        @test supports == expect_supports
        @test weights == expect_weights
    end

    # test Gauss-Hermite method for infinite interval
    # at the compiling step FGQ.gausshermite creates huge allocation
    @testset "gauss_hermite" begin
        @test_throws ErrorException generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, 0., Inf, Val(gauss_hermite))
        @test_throws ErrorException generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, -Inf, 0., Val(gauss_hermite))
        @test_throws ErrorException generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, Inf, Inf, Val(gauss_hermite))
        @test_throws ErrorException generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, 0., 1., Val(gauss_hermite))
        (supports, weights) = generate_supports_and_coeffs(
                             IntervalSet(-Inf, Inf), t, 5, -Inf, Inf, Val(gauss_hermite))
        @test length(supports) == 5
        (expect_supports, expect_weights) = FGQ.gausshermite(5)
        expect_weights = expect_weights .* exp.(expect_supports.^2)
        @test length(expect_weights) == 5
        @test supports == expect_supports
        @test weights == expect_weights
    end

    # test Gauss-Laguerre method for semi-infinite interval
    @testset "gauss_laguerre" begin
        @test_throws ErrorException generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, -Inf, Inf, Val(gauss_laguerre))
        @test_throws ErrorException generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, 0., 1., Val(gauss_laguerre))
        (supports, weights) = generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, 1., Inf, Val(gauss_laguerre))
        (expect_supports, expect_weights) = gausslaguerre(5)
        expect_weights = expect_weights .* exp.(expect_supports)
        expect_supports = expect_supports .+ 1.
        @test supports == expect_supports
        @test weights == expect_weights
        (supports, weights) = generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 5, -Inf, 1., Val(gauss_laguerre))
        (expect_supports, expect_weights) = gausslaguerre(5)
        expect_weights = expect_weights .* exp.(expect_supports)
        expect_supports = -expect_supports .+ 1.
        @test supports == expect_supports
        @test weights == expect_weights
    end
end

@testset "trapezoid Rule" begin
    (supports, weights) = generate_supports_and_coeffs(IntervalSet(-Inf, Inf), t, 3, 0., 2., Val(trapezoid))
    @test length(supports) == 3
    @test length(weights) == 3
    @test weights == [0.5, 1., 0.5]
    @test supports == [0., 1., 2.]
end

# test the function that transform (semi-)infinite domain to finite domain
# with the default transform function
@testset "Infinite transform for (semi-)infinite domain" begin
    @testset "infinite_transform" begin
        @test_throws ErrorException infinite_transform(IntervalSet(-Inf, Inf), t, 5, 0., 1.)
        # Infinite domain
        (supports, weights) = infinite_transform(IntervalSet(-Inf, Inf), t, 5, -Inf, Inf,
                                                 Val(gauss_legendre))
        (expect_supports, expect_weights) = gausslegendre(5)
        expect_weights = expect_weights .* MEM._default_dx.(expect_supports,
                                                            -Inf, Inf)
        expect_supports = MEM._default_x.(expect_supports, -Inf, Inf)
        @test supports == expect_supports
        @test weights == expect_weights
        # Lower bounded semi-infinite domain
        (supports, weights) = infinite_transform(IntervalSet(-Inf, Inf), t, 5, 0., Inf,
                                                 Val(gauss_legendre))
        (expect_supports, expect_weights) = gausslegendre(5)
        expect_supports = (expect_supports .+ 1) ./ 2
        expect_weights = expect_weights ./ 2
        expect_weights = expect_weights .* MEM._default_dx.(expect_supports,
                                                            0., Inf)
        expect_supports = MEM._default_x.(expect_supports, 0., Inf)
        @test supports == expect_supports
        @test weights == expect_weights
        # Upper bounded semi-infinite domain
        (supports, weights) = infinite_transform(IntervalSet(-Inf, Inf), t, 5, -Inf, 0.,
                                                 Val(gauss_legendre))
        (expect_supports, expect_weights) = gausslegendre(5)
        expect_supports = (expect_supports .+ 1) ./ 2
        expect_weights = expect_weights ./ 2
        expect_weights = expect_weights .* MEM._default_dx.(expect_supports,
                                                            -Inf, 0.)
        expect_supports = MEM._default_x.(expect_supports, -Inf, 0.)
        @test supports == expect_supports
        @test weights == expect_weights
    end
end

@testset "Data Generation" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in [0., 1.])
    @infinite_parameter(m, y[1:2] in [0., 1.])
    dist1 = Normal(0., 1.)
    dist2 = MvNormal([0., 0.], [1. 0.;0. 10.])
    @infinite_parameter(m, β in dist1)
    @infinite_parameter(m, ω[1:2] in dist2)
    @testset "generate_measure_data (univariate)" begin
        @test_throws ErrorException generate_measure_data(x, 10, 1., 2.)
        measure_data = generate_measure_data(x, 10, 0.3, 0.7, eval_method = gauss_legendre)
        (expect_supports, expect_weights) = generate_supports_and_coeffs(IntervalSet(0., 1.), x, 10, 0.3, 0.7, Val(gauss_legendre))
        @test measure_data.supports == expect_supports
        @test measure_data.coefficients == expect_weights
        measure_data = generate_measure_data(β, 10, -1., Inf)
        @test length(measure_data.supports) == 10
        @test all([i >= -1. for i in measure_data.supports])
        measure_data = generate_measure_data(β, 10, -Inf, 1.)
        @test length(measure_data.supports) == 10
        @test all([i <= 1. for i in measure_data.supports])
    end

    @testset "generate_measure_data (multivariate)" begin
        lb = convert(JuMPC.SparseAxisArray, [0.3, 0.2])
        ub = convert(JuMPC.SparseAxisArray, [0.5, 0.9])
        measure_data = generate_measure_data(y, 10, lb, ub)
        @test length(measure_data.supports) == 10
        @test length(measure_data.coefficients) == 10
        @test sum(measure_data.coefficients) == 0.14
        warn = "Truncated distribution for multivariate distribution is " *
               "not supported. Lower bounds and upper bounds are ignored."
        @test_logs (:warn, warn) generate_measure_data(ω, 10,
                                  convert(JuMPC.SparseAxisArray, [-1., -10.]),
                                  convert(JuMPC.SparseAxisArray, [1., 10.]))
    end

    @testset "generate_supports_and_coeffs for unextended types" begin
        @test_throws ErrorException generate_supports_and_coeffs(
                                    BadSet(), x, 10, 0., 1., Val(mc_sampling))
    end
end


text:
using Base: @propagate_inbounds
using ..TopOpt.TopOptProblems: getdh, getE, getν, getdensity, gettypes
using LinearAlgebra: norm
using Ferrite: getdim, getrefshape, value, getlowerdim

"""
Generate element stiffness matrices
"""
function make_Kes_and_fes(problem::TrussProblem, quad_order=1)
    make_Kes_and_fes(problem, quad_order, Val{:Static})
end

function make_Kes_and_fes(problem::TrussProblem, ::Type{Val{mat_type}}) where mat_type
    make_Kes_and_fes(problem, 1, Val{mat_type})
end

function make_Kes_and_fes(problem::TrussProblem{xdim, T}, quad_order, ::Type{Val{mat_type}}) where {xdim, T, mat_type}
    dh = getdh(problem)
    Es = getE(problem)
    # ν = getν(problem)
    # ρ = getdensity(problem)
    As = getA(problem)

    # * Shape functions and quadrature rule
    interpolation_space = Ferrite.default_interpolation(getcelltype(problem.truss_grid.grid))
    # Lagrange{ξdim, refshape, geom_order}()
    ξdim = getdim(interpolation_space)
    refshape = getrefshape(dh.field_interpolations[1])
    quadrature_rule = QuadratureRule{ξdim, refshape}(quad_order)
    cellvalues = GenericCellScalarValues(T, quadrature_rule, interpolation_space; xdim=xdim)

    # * A Line element's faces are not meaningful in truss problems
    # placeholder to make type right
    facevalues = FaceScalarValues(QuadratureRule{ξdim-1, refshape}(quad_order), interpolation_space)

    # * Calculate element stiffness matrices
    n_basefuncs = getnbasefunctions(cellvalues)
    Kesize = xdim*n_basefuncs
    MatrixType, VectorType = gettypes(T, Val{mat_type}, Val{Kesize})
    Kes, weights = _make_Kes_and_weights(dh, Tuple{MatrixType, VectorType}, Val{n_basefuncs}, Val{xdim*n_basefuncs},
        Es, As, quadrature_rule, cellvalues)

    # ! distributed load, not used in a truss problem
    # dloads = _make_dloads(weights, problem, facevalues)

    return Kes, weights, cellvalues, facevalues #dloads,
end

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

"""
    GenericCellScalarValues{ξdim,xdim,T<:Real,refshape<:AbstractRefShape} <: CellScalarValues{xdim,T,refshape}

`GenericCellScalarValues` is a generalization of the `Ferrite.CellScalarValues` to separate the reference domain
dimension `ξdim` and the node coordinate dimension `xdim`. While in most of the solid mechanics cases `ξdim = xdim`,
in a truss element (line element with nodes in 2D or 3D), `xdim = 2` or `3` and `ξdim = 1`.

**Arguments:**
* `ξdim` : reference domain dimension
* `xdim` : node coordinate dimension
"""
struct GenericCellScalarValues{ξdim,xdim,T,refshape} <: CellValues{xdim,T,refshape}
    N::Matrix{T}
    dNdx::Matrix{Vec{xdim,T}}
    dNdξ::Matrix{Vec{ξdim,T}}
    detJdV::Vector{T}
    M::Matrix{T}
    dMdξ::Matrix{Vec{ξdim,T}}
    qr_weights::Vector{T}
end

function GenericCellScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol)
    TrussCellScalarValues(Float64, quad_rule, func_interpol, geom_interpol)
end

function GenericCellScalarValues(::Type{T}, quad_rule::QuadratureRule{ξdim,shape}, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol; xdim=ξdim) where {ξdim,T,shape<:Ferrite.AbstractRefShape}
    @assert getdim(func_interpol) == getdim(geom_interpol)
    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
    n_qpoints = length(getweights(quad_rule))
    # * Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol)
    N    = fill(zero(T)           * T(NaN), n_func_basefuncs, n_qpoints)
    dNdx = fill(zero(Vec{xdim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Vec{ξdim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
    # * Geometry interpolation
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M    = fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{ξdim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)
    for (qp, ξ) in enumerate(quad_rule.points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
        end
        for i in 1:n_geom_basefuncs
            dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
        end
    end
    detJdV = fill(T(NaN), n_qpoints)
    GenericCellScalarValues{ξdim,xdim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
end

# common values
using Ferrite: getnbasefunctions, getngeobasefunctions, getnquadpoints
getn_scalarbasefunctions(cv::GenericCellScalarValues) = size(cv.N, 1)
@propagate_inbounds shape_gradient(cv::GenericCellScalarValues, q_point::Int, base_func::Int) = cv.dNdx[base_func, q_point]
@propagate_inbounds getdetJdV(cv::GenericCellScalarValues, q_point::Int) = cv.detJdV[q_point]

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

"""
    pinv(::Tensor{order, dim})

Compute the pseudo-inverse of a Vector tensor.
"""
function LinearAlgebra.pinv(t::Vec{dim, T}) where {dim, T}
    LinearAlgebra.Transpose{T, Vec{dim, T}}(t / sum(t.^2))
end

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

"""
    Kes, weights = _make_Kes_and_weights(dof_handler, Tuple{MatrixType, VectorType}, Val{n_basefuncs}, Val{dim*n_basefuncs}, Es, As, quadrature_rule, cellvalues)

`weights` : a vector of `xdim*n_basefuncs` vectors, element_id => self-weight load vector, in truss elements, they are all zeros.
"""
function _make_Kes_and_weights(
    dh::DofHandler{xdim, N, T},
    ::Type{Tuple{MatrixType, VectorType}},
    ::Type{Val{n_basefuncs}},
    ::Type{Val{Kesize}},
    Es::Vector{T}, As::Vector{T},
    quadrature_rule, cellvalues::GenericCellScalarValues) where {xdim, N, T, MatrixType <: StaticArray, VectorType, n_basefuncs, Kesize}
    nel = getncells(dh.grid)
    Kes = Symmetric{T, MatrixType}[]
    sizehint!(Kes, nel)
    # body_force = ρ .* g # Force per unit volume
    weights = [zeros(VectorType) for i in 1:nel]
    Ke_e = zeros(T, xdim, xdim)
    fe = zeros(T, Kesize)
    Ke_0 = Matrix{T}(undef, Kesize, Kesize)

    celliterator = CellIterator(dh)
    for (k, cell) in enumerate(celliterator)
        Ke_0 .= 0
        truss_reinit!(cellvalues, cell, As[k])
        # fe = weights[k]
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for b in 1:n_basefuncs
                ∇ϕb = shape_gradient(cellvalues, q_point, b)
                # ϕb = shape_value(cellvalues, q_point, b)
                for d2 in 1:xdim
                    # self weight force calculation
                    # fe = @set fe[(b-1)*dim + d2] += ϕb * body_force[d2] * dΩ
                    for a in 1:n_basefuncs
                        ∇ϕa = shape_gradient(cellvalues, q_point, a)
                        # TODO specialized KroneckerDelta struct to make dotdot more efficient
                        Ke_e .= Es[k] * ∇ϕa ⊗ ∇ϕb * dΩ
                        for d1 in 1:xdim
                            #if dim*(b-1) + d2 >= dim*(a-1) + d1
                            Ke_0[xdim*(a-1) + d1, xdim*(b-1) + d2] += Ke_e[d1,d2]
                            #end
                        end
                    end
                end
            end
        end
        # weights[k] = fe
        if MatrixType <: SizedMatrix # Work around because full constructor errors
            push!(Kes, Symmetric(SizedMatrix{Kesize,Kesize,T}(Ke_0)))
        else
            push!(Kes, Symmetric(MatrixType(Ke_0)))
        end
    end
    return Kes, weights
end

@inline function truss_reinit!(cv::GenericCellScalarValues{ξdim,xdim,T}, ci::CellIterator{xdim,N,T}, crossec::T) where {ξdim,xdim,N,T}
    Ferrite.check_compatible_geointerpolation(cv, ci)
    truss_reinit!(cv, ci.coords, crossec)
end

"""
Reinit a cell for a truss element, using the nodal coordinates `x`, cross section `crossec`
"""
function truss_reinit!(cv::GenericCellScalarValues{ξdim,xdim,T}, x::AbstractVector{Vec{xdim,T}}, crossec::T) where {ξdim,xdim,T}
    n_geom_basefuncs = getngeobasefunctions(cv)
    n_func_basefuncs = getn_scalarbasefunctions(cv)
    @assert length(x) == n_geom_basefuncs
    isa(cv, CellVectorValues) && (n_func_basefuncs *= xdim)

    @inbounds for i in 1:length(cv.qr_weights)
        w = cv.qr_weights[i]
        dxdξ = zero(Tensor{1,xdim})
        for j in 1:n_geom_basefuncs
            # in a truss element, x_j ∈ R, dMdξ_j ∈ R, ξ ∈ R
            # cv.dMdξ[j, i] is a 1-1 tensor here
            dxdξ += x[j] * cv.dMdξ[j, i][1]
        end
        # detJ = √(J' J), J = dxdξ ∈ R(n x 1)
        detJ = norm(dxdξ)
        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
        cv.detJdV[i] = detJ * w * crossec
        Jinv = pinv(dxdξ)
        for j in 1:n_func_basefuncs
            # cv.dNdξ[j, i] is a 1-1 tensor here
            cv.dNdx[j, i] = cv.dNdξ[j, i][1] * Jinv'
        end
    end
end

ext: jl


text:
using Documenter, GeoParams, Plots
push!(LOAD_PATH,"../src/")

@info "Making documentation..."
makedocs(
  sitename="GeoParams.jl",
  authors = "Boris Kaus and contributors",
  modules = [GeoParams],
  format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"), # easier local build
  pages = [
        "Home" => "index.md",
        "User Guide" => Any[
            "Nondimensionalization" =>  "man/nondimensionalize.md",
            "Material Parameters" =>  "man/materialparameters.md",
            "Density" =>  "man/density.md",
            "Creep laws" =>  "man/creeplaws.md",
            "Heat Capacity" =>  "man/heatcapacity.md",
            "Conductivity" =>  "man/conductivity.md",
            "Latent heat" =>  "man/latentheat.md",
            "Radioactive heat" =>  "man/radioactiveheating.md",
            "Shear heating" =>  "man/shearheating.md",
            "Gravity" =>  "man/gravity.md",
            "Partial Melting"  => "man/melting.md",
            "Zircon Ages"  => "man/zirconages.md",
            "Phase Diagrams" => "man/phasediagrams.md",
            "Seismic Velocity" => "man/seismicvelocity.md",
        ],
        "Plotting"  => "man/plotting.md",
        "List of functions"  => "man/listfunctions.md",
        "Contributing" => "man/contributing.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaGeodynamics/GeoParams.jl.git",
    devbranch = "main"
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ObjectMetricStatus indicates the current value of a metric describing a kubernetes object (for example, hits-per-second on an Ingress object).

    IoK8sApiAutoscalingV2beta2ObjectMetricStatus(;
        current=nothing,
        describedObject=nothing,
        metric=nothing,
    )

    - current::IoK8sApiAutoscalingV2beta2MetricValueStatus : current contains the current value for the given metric
    - describedObject::IoK8sApiAutoscalingV2beta2CrossVersionObjectReference
    - metric::IoK8sApiAutoscalingV2beta2MetricIdentifier : metric identifies the target metric by name and selector
"""
mutable struct IoK8sApiAutoscalingV2beta2ObjectMetricStatus <: SwaggerModel
    current::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricValueStatus } # spec name: current
    describedObject::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2CrossVersionObjectReference } # spec name: describedObject
    metric::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricIdentifier } # spec name: metric

    function IoK8sApiAutoscalingV2beta2ObjectMetricStatus(;current=nothing, describedObject=nothing, metric=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2ObjectMetricStatus, Symbol("current"), current)
        setfield!(o, Symbol("current"), current)
        validate_property(IoK8sApiAutoscalingV2beta2ObjectMetricStatus, Symbol("describedObject"), describedObject)
        setfield!(o, Symbol("describedObject"), describedObject)
        validate_property(IoK8sApiAutoscalingV2beta2ObjectMetricStatus, Symbol("metric"), metric)
        setfield!(o, Symbol("metric"), metric)
        o
    end
end # type IoK8sApiAutoscalingV2beta2ObjectMetricStatus

const _property_map_IoK8sApiAutoscalingV2beta2ObjectMetricStatus = Dict{Symbol,Symbol}(Symbol("current")=>Symbol("current"), Symbol("describedObject")=>Symbol("describedObject"), Symbol("metric")=>Symbol("metric"))
const _property_types_IoK8sApiAutoscalingV2beta2ObjectMetricStatus = Dict{Symbol,String}(Symbol("current")=>"IoK8sApiAutoscalingV2beta2MetricValueStatus", Symbol("describedObject")=>"IoK8sApiAutoscalingV2beta2CrossVersionObjectReference", Symbol("metric")=>"IoK8sApiAutoscalingV2beta2MetricIdentifier")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2ObjectMetricStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2ObjectMetricStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2ObjectMetricStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2ObjectMetricStatus)
    (getproperty(o, Symbol("current")) === nothing) && (return false)
    (getproperty(o, Symbol("describedObject")) === nothing) && (return false)
    (getproperty(o, Symbol("metric")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricStatus }, name::Symbol, val)
end


text:
using Test
using AbstractLattices

function (∨)(a::Number,b::Number)
    return max(a,b)
end

function (∧)(a::Number,b::Number)
    return min(a,b)
end

@test 5 ∧ 10 == 5
@test 5 ∨ 10 == 10


text:
#!/usr/bin/env julia

function main()
    if length(ARGS) > 0
        n = parse(Int64, ARGS[1])
        accum = 0
        for i in 1:(n-1)
            if i % 3 == 0 || i % 5 == 0
                accum += i
            end
        end
        @printf("The sum of the numbers up to %d divisible by 3 or 5 is %d\n",
                n, accum)
    else
        error("this program needs a single numeric argument")
    end
end

main()


text:
using Nextion
using Test

@testset "button" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    nexButtonPlus = NexButton(nexSerial, Name("b0"), cid=ComponentID(1))
    nexNumber0 = NexNumber(nexSerial, Name("n0"), cid=ComponentID(2))
    nexButtonMinus = NexButton(nexSerial, Name("b1"), cid=ComponentID(3))
    nexButtonEnter = NexButton(nexSerial, Name("b2"), cid=ComponentID(4))

    send(nexSerial, "page pg_but")

    nexButtonPlus.backcolor = Colour.GREEN
    nexButtonMinus.backcolor = Colour.RED

    #nexNumber0.value = value
    value = nexNumber0.value
    while true
        sleep(0.2)
        new_value = nexNumber0.value
        if new_value != value
            value = new_value
            println(value)
        end
    end

    close(nexSerial)

    @test true

end


text:
# multiply.jl

#export LinearMapAA_test_vmul # testing
#export LinearMapAA_test_mul # testing
using LinearMapsAA
using Test: @test, @testset

"""
    LinearMapAA_test_vmul(A::LinearMapAM)
    LinearMapAA_test_vmul(A::LinearMapAO)
tests for multiply by vector, scalar, and Matrix
"""
function LinearMapAA_test_vmul(A::LinearMapAM)
    B = Matrix(A)
    (M,N) = size(A)

    u = rand(M)
    v = rand(N)

    @testset "A*v" begin
        Bv = B * v
        Bpu = B' * u

        y = A * v
        x = A' * u
        @test isapprox(Bv, y)
        @test isapprox(Bpu, x)

        mul!(y, A, v)
        mul!(x, A', u)
        @test isapprox(Bv, y)
        @test isapprox(Bpu, x)
    end

    #= nah
    @testset "u'*A" begin
        uB = u' * B
        x = u' * A
        @test isapprox(uB, x)
        mul!(x, u', A)
        @test isapprox(uB, x)
    end
    =#

    @testset "A*X" begin
        X = rand(N,2)
        BX = B * X
        Y = A * X
        @test Y isa LinearMapAM
        Y = Matrix(Y)
        @test isapprox(BX, Y)
        Y .= 0
        # mul!(Y, A, X) # doesn't work because A*X is AM
        # @test isapprox(BX, Y)
    end

    @testset "X*A" begin
        X = rand(2,M)
        XB = X * B
        Y = X * A
        @test Y isa LinearMapAM
        Y = Matrix(Y)
        @test isapprox(XB, Y)
        Y .= 0
        # mul!(Y, X, A) # doesn't work because X*A is AM
        # @test isapprox(XB, Y)
    end

    @testset "5-arg" begin
        y1 = randn(M)
        y2 = copy(y1)
        mul!(y1, A, v, 2, 3)
        mul!(y2, B, v, 2, 3)
        @test isapprox(y1, y2)

        x1 = randn(N)
        x2 = copy(x1)
        mul!(x1, A', u, 4, 3)
        mul!(x2, B', u, 4, 3)
        @test isapprox(x1, x2)
    end

    @testset "*s" begin
        s = 4.2
        C = s * A
        @test isapprox(Matrix(C), s * B)
        C = A * s
        @test isapprox(Matrix(C), B * s)
    end

#=
    s = 4.2
    C = copy(A)
    lmul!(s, C)
    @test isapprox(s * B * v, C * v)

    C = copy(A)
    rmul!(C, s)
    @test isapprox(s * B * v, C * v)
=#

    true
end

function LinearMapAA_test_vmul(A::LinearMapAO)
    B = Matrix(A)
    (M,N) = size(A)

    u = rand(M)
    v = rand(N)
    u = reshape(u, A._odim)
    v = reshape(v, A._idim)

    @testset "A*v" begin
        Bv = reshape(B * v[:], A._odim)
        Bpu = reshape(B' * u[:], A._idim)

        y = A * v
        x = A' * u
        @test isapprox(Bv, y)
        @test isapprox(Bpu, x)

        mul!(y, A, v)
        mul!(x, A', u)
        @test isapprox(Bv, y)
        @test isapprox(Bpu, x)

        A1 = redim(A ; idim = (N,))
        y1 = A1 * vec(v)
        @test y1 ≈ y
    end

    #= nah
    @testset "u'*A" begin
        uB = reshape(u[:]' * B, A._idim)
        x = u' * A
        @test isapprox(uB, x)
        mul!(x, u', A)
        @test isapprox(uB, x)
    end
    =#

    @testset "A*X" begin
        K = 2
        X = rand(A._idim..., K)
        BX = reshape(B * reshape(X, :, K), A._odim..., K)
        Y = A * X
        @test Y isa AbstractArray
        @test isapprox(BX, Y)

        Y .= 0
        mul!(Y, A, X)
        @test isapprox(BX, Y)

        # 5-arg
        Y1 = copy(Y)
        Y2 = copy(Y)
        mul!(Y1, A, X, 2, 3)
        mul!(reshape(Y2,:,K), B, reshape(X,:,K), 2, 3)
        Y2 = reshape(Y2, A._odim..., K)
        @test isapprox(Y1, Y2)
    end

    # todo: should eventually test complex X*A as well

    @testset "X*A" begin
        K = 2
        X = rand(A._odim..., K)
        XB = reshape(X, :, K)' * B
        XB = reshape(XB', A._idim..., K)
        Y = X * A
        @test Y isa AbstractArray
        @test isapprox(XB, Y)

        Y .= 0
        mul!(Y, X, A)
        @test isapprox(XB, Y)

        # 5-arg
        Y1 = copy(Y)
        Y2 = copy(Y)
        mul!(Y1, X, A, 2, 3)
           mul!(reshape(Y2,:,K), B', reshape(X,:,K), 2, 3)
        Y2 = reshape(Y2, A._idim..., K)
           @test isapprox(Y1, Y2)
    end

    @testset "*s" begin
        s = 4.2
        C = s * A
        @test isapprox(Matrix(C), s * B)
        C = A * s
        @test isapprox(Matrix(C), B * s)
    end

    true
end

# object multiplication
function LinearMapAA_test_mul( ;
    A::LinearMapAO = LinearMapAA(rand(6,4) ; odim=(2,3), idim=(1,4)),
)
    (M,N) = size(A)

    @testset "O*O" begin
        @test A'*A isa LinearMapAO
    end

    @testset "O*M" begin
        B = LinearMapAA(rand(N,3)) # AM
        O = redim(A ; idim=(N,))
        @test O*B isa LinearMapAO
    end

    @testset "M*O" begin
        B = LinearMapAA(rand(3,M)) # AM
        O = redim(A ; odim=(M,))
        @test B*O isa LinearMapAO
    end

    true
end


text:
export TabularApproximator

"""
    TabularApproximator(table<:AbstractArray)

For `table` of 1-d, it will create a [`V_APPROXIMATOR`](@ref). For `table` of 2-d, it will create a [`QApproximator`].

!!! warning
    For `table` of 2-d, the first dimension is action and the second dimension is state.
"""
struct TabularApproximator{N,T<:AbstractArray} <: AbstractApproximator
    table::T
    function TabularApproximator(table::T) where {T<:AbstractArray}
        n = ndims(table)
        n <= 2 || throw(ArgumentError("the dimention of table must be <= 2"))
        new{n,T}(table)
    end
end

"""
    TabularApproximator(; n_state, n_action = nothing, init = 0.0)
"""
function TabularApproximator(; n_state, n_action = nothing, init = 0.0)
    table = isnothing(n_action) ? fill(init, n_state) : fill(init, n_action, n_state)
    TabularApproximator(table)
end

(app::TabularApproximator{1})(s::Int) = @views app.table[s]

(app::TabularApproximator{2})(s::Int) = @views app.table[:, s]
(app::TabularApproximator{2})(s::Int, a::Int) = app(s)[a]

function RLBase.update!(app::TabularApproximator{1}, correction::Pair)
    s, e = correction
    app.table[s] += e
end

function RLBase.update!(app::TabularApproximator{2}, correction::Pair)
    (s, a), e = correction
    app.table[a, s] += e
end

function RLBase.update!(Q::TabularApproximator{2}, correction::Pair{Int,Vector{Float64}})
    s, errors = correction
    for (a, e) in enumerate(errors)
        Q.table[a, s] += e
    end
end

RLBase.ApproximatorStyle(::TabularApproximator{1}) = VApproximator()
RLBase.ApproximatorStyle(::TabularApproximator{2}) = QApproximator()


text:
################################################################################
# Dynamics Violation
################################################################################

mutable struct DynamicsViolation{SVd,T}
	N::Int
	vio::SVd
	max::T
end

function DynamicsViolation(N::Int)
	vio = zeros(N-1)
	max = 0.0
	TYPE = typeof.((vio,max))
	return DynamicsViolation{TYPE...}(N,vio,max)
end

function dynamics_violation(model::AbstractGameModel, pdtraj::PrimalDualTraj)
    N = pdtraj.probsize.N
    dyn_vio = DynamicsViolation(N)
    for k = 1:N-1
        dyn_vio.vio[k] = maximum(abs.(dynamics_residual(model, pdtraj, k)))
    end
	dyn_vio.max = maximum(dyn_vio.vio)
    return dyn_vio
end

function dynamics_violation(model::AbstractGameModel, pdtraj::PrimalDualTraj, i::Int)
    N = pdtraj.probsize.N
	pi = pdtraj.probsize.pz[i]

    dyn_vio = DynamicsViolation(N)
    for k = 1:N-1
        dyn_vio.vio[k] = maximum(abs.(dynamics_residual(model, pdtraj, k)[pi]))
    end
	dyn_vio.max = maximum(dyn_vio.vio)
    return dyn_vio
end

################################################################################
# Control Violation
################################################################################

mutable struct ControlViolation{SVd,T}
	N::Int
	vio::SVd
	max::T
end

function ControlViolation(N::Int)
	vio = zeros(N-1)
	max = 0.0
	TYPE = typeof.((vio,max))
	return ControlViolation{TYPE...}(N,vio,max)
end

function control_violation(game_con::GameConstraintValues, pdtraj::PrimalDualTraj)
	N = pdtraj.probsize.N
	u_vio = ControlViolation(N)
	for conval in game_con.control_conval
		TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
		TrajectoryOptimization.max_violation!(conval)
		u_vio.vio[conval.inds] = max.(u_vio.vio[conval.inds], conval.c_max)
	end
	u_vio.max = maximum(u_vio.vio)
    return u_vio
end

function control_violation(game_con::GameConstraintValues, pdtraj::PrimalDualTraj, i::Int)
	N = pdtraj.probsize.N
	pi = pdtraj.probsize.pu[i]

	u_vio = ControlViolation(N)
	for conval in game_con.control_conval
		TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
		TrajectoryOptimization.max_violation!(conval)
		c_max = [max(0.0, maximum(v[pi])) for v in conval.vals]
		u_vio.vio[conval.inds] = max.(u_vio.vio[conval.inds], c_max)
	end
	u_vio.max = maximum(u_vio.vio)
    return u_vio
end

################################################################################
# State Violation
################################################################################

mutable struct StateViolation{SVd,T}
	N::Int
	vio::SVd
	max::T
end

function StateViolation(N::Int)
	vio = zeros(N)
	max = 0.0
	TYPE = typeof.((vio,max))
	return StateViolation{TYPE...}(N,vio,max)
end

function state_violation(game_con::GameConstraintValues, pdtraj::PrimalDualTraj)
	N = pdtraj.probsize.N
	p = game_con.probsize.p
	x_vio = StateViolation(N)
	for i = 1:p
		for conval in game_con.state_conval[i]
			TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
			TrajectoryOptimization.max_violation!(conval)
			x_vio.vio[conval.inds] = max.(x_vio.vio[conval.inds], conval.c_max)
		end
	end
	x_vio.max = maximum(x_vio.vio)
    return x_vio
end

function state_violation(game_con::GameConstraintValues, pdtraj::PrimalDualTraj, i::Int)
	N = pdtraj.probsize.N
	p = game_con.probsize.p

	x_vio = StateViolation(N)
	for conval in game_con.state_conval[i]
		TrajectoryOptimization.evaluate!(conval, pdtraj.pr)
		TrajectoryOptimization.max_violation!(conval)
		# @show size(conval.vals)
		# @show size(conval.vals[1])
		# @show conval.vals[1]
		# ned to identify the constraint that apply to player i
		# c_max = [max(0.0, maximum(v)) for v in conval.vals]
		# x_vio.vio[conval.inds] = max.(x_vio.vio[conval.inds], c_max)
		x_vio.vio[conval.inds] = max.(x_vio.vio[conval.inds], conval.c_max)
	end
	x_vio.max = maximum(x_vio.vio)
    return x_vio
end

################################################################################
# Optimality Violation
################################################################################

mutable struct OptimalityViolation{SVd,T}
	N::Int
	vio::SVd
	max::T
end

function OptimalityViolation(N::Int)
	vio = zeros(N)
	max = 0.0
	TYPE = typeof.((vio,max))
	return OptimalityViolation{TYPE...}(N,vio,max)
end

function optimality_violation(core::NewtonCore)
	N = core.probsize.N
	p = core.probsize.p
	o_vio = OptimalityViolation(N)
	stamp = VStamp()
	for i = 1:p
		for k = 1:N
			stampify!(stamp, :opt, i, :x, 1, k)
			valid(stamp, N, p) ? o_vio.vio[k] = max(o_vio.vio[k], maximum(abs.(core.res_sub[stamp]))) : nothing
			stampify!(stamp, :opt, i, :u, i, k)
			valid(stamp, N, p) ? o_vio.vio[k] = max(o_vio.vio[k], maximum(abs.(core.res_sub[stamp]))) : nothing
		end
	end
	o_vio.max = maximum(o_vio.vio)
    return o_vio
end

function optimality_violation(core::NewtonCore, i::Int)
	N = core.probsize.N
	p = core.probsize.p
	o_vio = OptimalityViolation(N)
	stamp = VStamp()
	for k = 1:N
		stampify!(stamp, :opt, i, :x, 1, k)
		valid(stamp, N, p) ? o_vio.vio[k] = max(o_vio.vio[k], maximum(abs.(core.res_sub[stamp]))) : nothing
		stampify!(stamp, :opt, i, :u, i, k)
		valid(stamp, N, p) ? o_vio.vio[k] = max(o_vio.vio[k], maximum(abs.(core.res_sub[stamp]))) : nothing
	end
	o_vio.max = maximum(o_vio.vio)
    return o_vio
end


text:
module PATHSolver

using ForwardDiff
using SparseArrays
using Random
using LinearAlgebra

const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("PATHSolver not properly installed. Please re-build PATHSolver.")
end

export solveMCP, solveLCP, options

const cached_J = [convert(SparseMatrixCSC{Cdouble, Cint}, zeros(0, 0))]
const cached_J_filled = Ref(false)

const status =
   [  :Solved,                          # 1 - solved
      :StationaryPointFound,            # 2 - stationary point found
      :MajorIterationLimit,             # 3 - major iteration limit
      :CumulativeMinorIterationLimit,   # 4 - cumulative minor iteration limit
      :TimeLimit,                       # 5 - time limit
      :UserInterrupt,                   # 6 - user interrupt
      :BoundError,                      # 7 - bound error (lb is not less than ub)
      :DomainError,                     # 8 - domain error (could not find a starting point)
      :InternalError                    # 9 - internal error
  ]

count_nonzeros(M::AbstractSparseMatrix) = nnz(M)
count_nonzeros(M::AbstractMatrix) = count(x -> x != 0, M) # fallback for dense matrices

###############################################################################
# wrappers for callback functions
###############################################################################
# static int (*f_eval)(int n, double *z, double *f);
# static int (*j_eval)(int n, int nnz, double *z, int *col_start, int *col_len,
      # int *row, double *data);
function f_user_wrap(n::Cint, z_ptr::Ptr{Cdouble}, f_ptr::Ptr{Cdouble})
  z = unsafe_wrap(Array{Cdouble}, z_ptr, Int(n), own=false)
  f = unsafe_wrap(Array{Cdouble}, f_ptr, Int(n), own=false)
  f .= user_f(z)
  return Cint(0)
end

function j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
                     col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint},
                     row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})

  z = unsafe_wrap(Array{Cdouble}, z_ptr, Int(n), own=false)
  J::SparseMatrixCSC{Cdouble, Cint} = user_j(z)
  if nnz(J) > expected_nnz
    println("nnz(J) = ", nnz(J))
    println("expected_nnz = ", expected_nnz)
    error("Evaluated jacobian has more nonzero entries than were initially provided in solveMCP(). Try solveMCP(..., nnz=n^2).")
  end
  load_sparse_matrix(J, n, expected_nnz, col_start_ptr, col_len_ptr, row_ptr, data_ptr)
  return Cint(0)
end

function cached_j_user_wrap(n::Cint, expected_nnz::Cint, z_ptr::Ptr{Cdouble},
                     col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint},
                     row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})
  if !(cached_J_filled[])
    load_sparse_matrix(cached_J[], n, expected_nnz, col_start_ptr, col_len_ptr, row_ptr, data_ptr)
    cached_J_filled[] = true
  end
  return Cint(0)
end

function load_sparse_matrix(J::SparseMatrixCSC, n::Cint, expected_nnz::Cint,
                            col_start_ptr::Ptr{Cint}, col_len_ptr::Ptr{Cint},
                            row_ptr::Ptr{Cint}, data_ptr::Ptr{Cdouble})
  # Transfer data from the computed jacobian into the sparse format that PATH
  # expects. Fortunately, PATH uses a compressed-sparse-column storage which
  # is compatible with Julia's default SparseMatrixCSC format.

  # col_start in PATH corresponds to J.colptr[1:end-1]
  col_start = unsafe_wrap(Array{Cint}, col_start_ptr, Int(n), own=false)
  # col_len in PATH corresponds to diff(J.colptr)
  col_len = unsafe_wrap(Array{Cint}, col_len_ptr, Int(n), own=false)
  # row in PATH corresponds to rowvals(J)
  row = unsafe_wrap(Array{Cint}, row_ptr, Int(expected_nnz), own=false)
  # data in PATH corresponds to nonzeros(J)
  data = unsafe_wrap(Array{Cdouble}, data_ptr, Int(expected_nnz), own=false)

  @inbounds for i in 1:n
    col_start[i] = J.colptr[i]
    col_len[i] = J.colptr[i + 1] - J.colptr[i]
  end

  rv = rowvals(J)
  nz = nonzeros(J)
  num_nonzeros = nnz(J)
  @inbounds for i in 1:num_nonzeros
    row[i] = rv[i]
    data[i] = nz[i]
  end
end

###############################################################################
# solveMCP
###############################################################################

# solveMCP without z0, without j_eval
function solveMCP(f_eval::Function,
                  lb::AbstractVector{T}, ub::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  nnz=-1) where {T <: Number, S <: String}

  z0 = (lb + ub) ./ 2
  return solveMCP(f_eval, lb, ub, z0, var_name, con_name; nnz=nnz)
end

# solveMCP without z0, with j_eval
function solveMCP(f_eval::Function, j_eval::Function,
                  lb::AbstractVector{T}, ub::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  nnz=-1) where {T <: Number, S <: String}

  z0 = (lb + ub) ./ 2
  return solveMCP(f_eval, j_eval, lb, ub, z0, var_name, con_name; nnz=nnz)
end

# solveMCP with z0, without j_eval
function solveMCP(f_eval::Function,
                  lb::AbstractVector{T}, ub::AbstractVector{T}, z0::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  nnz=-1) where {T <: Number, S <: String}

  j_eval = x -> ForwardDiff.jacobian(f_eval, x)
  return solveMCP(f_eval, j_eval, lb, ub, z0, var_name, con_name; nnz=nnz)
end

# Full implementation of solveMCP  / solveMCP with z0, with j_eval
function solveMCP(f_eval::Function, j_eval::Function,
                  lb::AbstractVector{T}, ub::AbstractVector{T}, z0::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  nnz=-1) where {T <: Number, S <: String}

  if length(var_name)==0
    var_name = C_NULL
  end

  if length(con_name)==0
    con_name = C_NULL
  end

  global user_f = f_eval
  global user_j = j_eval
  f_user_cb = @cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
  j_user_cb = @cfunction(j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))

  n = length(lb)
  z = max.(lb, min.(ub, z0))
  f = zeros(n)

  if nnz == -1
    # overestimating number of nonzeros in Jacobian
    nnz = max( count_nonzeros(j_eval(lb)), count_nonzeros(j_eval(ub)) )
    nnz = max( nnz, count_nonzeros(j_eval(z)) )
    for i in 1:2
      z_rand = max.(lb, min.(ub, rand(Float64, n)))
      nnz_rand = count_nonzeros(j_eval(z_rand))
      nnz = max( nnz, nnz_rand )
    end
    nnz = min( 2 * nnz, n * n )
  end

  t = ccall( (:path_main, libpath47julia), Cint,
          (Cint, Cint,
           Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Ptr{Cchar}}, Ptr{Ptr{Cchar}},
           Ptr{Cvoid}, Ptr{Cvoid}),
           n, nnz, z, f, lb, ub, var_name, con_name, f_user_cb, j_user_cb)

  remove_option_file()
  return status[t], z, f
end

###############################################################################
# solveLCP, LCP functions
###############################################################################

# solveLCP without z, without M
function solveLCP(f_eval::Function,
                  lb::AbstractVector{T}, ub::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  lcp_check=false) where {T <: Number, S <: String}

    return solveLCP(f_eval, lb, ub, copy(lb), var_name, con_name, lcp_check=lcp_check)
end

# solveLCP with z0, without M
function solveLCP(f_eval::Function,
                  lb::AbstractVector{T}, ub::AbstractVector{T}, z0::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  lcp_check=false) where {T <: Number, S <: String}

  J = ForwardDiff.jacobian(f_eval, lb)
  if lcp_check
      Jr = ForwardDiff.jacobian(f_eval, rand(Float64, size(lb)))
      if opnorm(J-Jr, 1) > 1e-8
          error("The problem does not seem linear. Use `solveMCP()` instead.")
      end
  end

  return solveLCP(f_eval, J, lb, ub, z0, var_name, con_name)
end

# solveLCP without z, with M
function solveLCP(f_eval::Function, M::AbstractMatrix,
                  lb::AbstractVector{T}, ub::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  lcp_check=false) where {T <: Number, S <: String}

    return solveLCP(f_eval, M, lb, ub, copy(lb), var_name, con_name, lcp_check=lcp_check)
end

# Full implmentation of solveLCP / solveLCP with z0, with M
function solveLCP(f_eval::Function, M::AbstractMatrix,
                  lb::AbstractVector{T}, ub::AbstractVector{T}, z0::AbstractVector{T},
                  var_name::AbstractVector{S}=String[],
                  con_name::AbstractVector{S}=String[];
                  lcp_check=false) where {T <: Number, S <: String}

  if length(var_name)==0
      var_name = C_NULL
  end

  if length(con_name)==0
      con_name = C_NULL
  end

  if lcp_check
      J = ForwardDiff.jacobian(f_eval, lb)
      if opnorm(J-M, 1) > 1e-8
          # warn("The user supplied Jacobian does not match with the result by FowardDiff.jacobian(). It proceeds with the user supplied Jacobian.")
          error("The user supplied Jacobian does not match with the result by FowardDiff.jacobian().")
      end
  end

  global user_f = f_eval
  cached_J[] = M
  cached_J_filled[] = false
  f_user_cb = @cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
  j_user_cb = @cfunction(cached_j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))

  n = length(lb)
  z = max.(lb, min.(ub, z0))
  f = zeros(n)

  nnzs = count_nonzeros(M)

  t = ccall( (:path_main, libpath47julia), Cint,
          (Cint, Cint,
           Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Ptr{Cchar}}, Ptr{Ptr{Cchar}},
           Ptr{Cvoid}, Ptr{Cvoid}),
           n, nnzs, z, f, lb, ub, var_name, con_name, f_user_cb, j_user_cb)

  remove_option_file()
  return status[t], z, f
end

###############################################################################
# handling PATH options
###############################################################################

function remove_option_file()
  if isfile("path.opt")
    rm("path.opt")
  end
end

function options(;kwargs...)
  opt_file = open("path.opt", "w")
  println(opt_file, "* Generated by PATHSolver.jl. Do not edit.")
  for (key, value) in kwargs
    println(opt_file, key, " ", value)
  end
  close(opt_file)
end

end # Module


text:
clf()
plot(X0a.tt[1:10:div(end,10)]*10, X0a.yy[1:10:div(end,10)]*10, label="c")
plot(X0b.tt[1:100:end], X0b.yy[1:100:end], label="b")
plot(X0a.tt[1:1:end], X0a.yy[1:1:end], label="a")
legend()


text:
export CUSOLVERError

struct CUSOLVERError <: Exception
    code::cusolverStatus_t
    msg::AbstractString
end
Base.show(io::IO, err::CUSOLVERError) = print(io, "CUSOLVERError(code $(err.code), $(err.msg))")

function CUSOLVERError(code::cusolverStatus_t)
    msg = status_message(code)
    return CUSOLVERError(code, msg)
end

function status_message(status)
    if status == CUSOLVER_STATUS_SUCCESS
        return "the operation completed successfully"
    elseif status == CUSOLVER_STATUS_NOT_INITIALIZED
        return "the library was not initialized"
    elseif status == CUSOLVER_STATUS_ALLOC_FAILED
        return "the resource allocation failed"
    elseif status == CUSOLVER_STATUS_INVALID_VALUE
        return "an invalid value was used as an argument"
    elseif status == CUSOLVER_STATUS_ARCH_MISMATCH
        return "an absent device architectural feature is required"
    elseif status == CUSOLVER_STATUS_EXECUTION_FAILED
        return "the GPU program failed to execute"
    elseif status == CUSOLVER_STATUS_INTERNAL_ERROR
        return "an internal operation failed"
    elseif status == CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED
        return "the matrix type is not supported."
    else
        return "unknown status"
    end
end

macro check(solver_func)
    quote
        local err::cusolverStatus_t
        err = $(esc(solver_func::Expr))
        if err != CUSOLVER_STATUS_SUCCESS
            throw(CUSOLVERError(err))
        end
        err
    end
end


text:
@testset "Variational Student-T Processes" begin
    N = 50

    x = rand(N)
    X = collect(reshape(x, :, 1))
    y = rand(N)
    ν = 3.0
    k = SqExponentialKernel()
    l = GaussianLikelihood()
    vi = AnalyticVI()
    m = VStP(x, y, k, l, vi, ν)
    @test_throws ErrorException VStP(x, y, k, l, vi, 0.5)
    @test_throws ErrorException VStP(x, y, k, l, QuadratureVI(), 0.5)
    @test_nowarn println(m)
    AGP.compute_kernel_matrices!(m)
    @test_nowarn AGP.local_prior_updates!(m, collect(eachrow(X)))
    @test AGP.objective(m) == ELBO(m)

    train!(m, 20)
end


text:
export meanSquaredError, meanAverageError
@doc doc"""
    meanSquaredError(M,x,y)
Computes the (mean) squared error between the two
[`MPoint`](@ref)s `x` and `y` on the ([`Power`](@ref)) manifold `M`.
"""
function meanSquaredError(M::mT,x::P,y::P) where {mT <: Manifold, P <: MPoint}
    return distance(M,x,y)^2
end
function meanSquaredError(M::Power,x::P,y::P) where {mT <: Manifold, P <: MPoint}
    return 1/prod(M.powerSize) * sum( distance.( Ref(M.manifold), getValue(x), getValue(y) ).^2 )
end
@doc doc"""
    meanSquaredError(M,x,y)
Computes the (mean) squared error between the two
[`MPoint`](@ref)s `x` and `y` on the ([`Power`](@ref)) manifold `M`.
"""
function meanAverageError(M::mT,x::P,y::P) where {mT <: Manifold, P <: MPoint}
    return distance(M,x,y)
end
function meanAverageError(M::Power,x::P,y::P) where {mT <: Manifold, P <: MPoint}
    return 1/prod(M.powerSize) * sum( distance.( Ref(M.manifold), getValue(x), getValue(y) ) )
end


text:
function greet(name)
    "Hello, $(name)!"
end


text:
using Pkg;Pkg.activate(@__DIR__)
using Comrade
using PyCall
using Plots

using ComradeAHMC

load_ehtim()
@pyimport eht_dmc as ed

# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
obs.add_scans()

obsavg = obs.avg_coherent(0.0, scan_avg=true)

img = ehtim.image.make_empty(512, μas2rad(250.0), obs.ra, obs.dec, obs.rf, obs.source)
img = img.add_gauss(1.0, [μas2rad(40.0), μas2rad(30.0), π/3, 0.0, 0.0])

obslist = obsavg.split_obs()

gainoff = 0.0
gainp = Dict("AA"=> 0.01,
             "AP" => 0.1,
             "AZ" => 0.1,
             "JC" => 0.1,
             "LM" => 0.2,
             "PV" => 0.1,
             "SM" => 0.1,
             "SR" => 0.1)
obsim = img.observe_same(ehtim.obsdata.merge_obs(obslist), ttype="fast",  gainp=gainp, gain_offset=gainoff, ampcal=false, phasecal=false, stabilize_scan_phase=true, stabilize_scan_amp=true)
# now get cal table
ctable = ehtim.calibrating.self_cal.self_cal(obsim, img, caltable=true)
ehtim.caltable.save_caltable(ctable, obsim, datadir=joinpath(@__DIR__ ,"CaltableTest"))

dvis = extract_vis(obsim)
dlcamp = extract_lcamp(obsim)
dcphase = extract_cphase(obsim)
damp = extract_amp(obsim)
st = scantable(dvis)

gcache = Comrade.GainCache(st)

struct Test{G}
    gcache::G
end

function Test(st::Comrade.ScanTable)
    gcache = Comrade.GainCache(st)
    return Test{typeof(gcache)}(gcache)
end

function (mod::Test)(θ)
    (;f, σ, τ, ξ,  gamp, gphase) = θ
    g = exp.(gamp).*cis.(gphase)
    m = f*rotated(stretched(Gaussian(), σ*τ, σ), ξ)
    return GainModel(mod.gcache, g, m)
end

# define the priors
distamp = (AA = Normal(0.0, 0.1),
           AP = Normal(0.0, 0.1),
           LM = Normal(0.0, 0.2),
           AZ = Normal(0.0, 0.1),
           JC = Normal(0.0, 0.1),
           PV = Normal(0.0, 0.1),
           #SM = LogNormal(0.0, 0.1)
         )
distphase = (AA = Normal(0.0, 1e-4),
             AP = Normal(0.0, π),
             LM = Normal(0.0, π),
             AZ = Normal(0.0, 1.0*π),
             JC = Normal(0.0, 1.0*π),
             PV = Normal(0.0, π),
             #SM = Normal(0.0, 1.0*π),
             )
prior = (
          f = Uniform(0.9, 1.1),
          σ = Uniform(μas2rad(10.0), μas2rad(30.0)),
          τ = Uniform(0.1, 1.0),
          ξ = Uniform(-π/2, π/2),
          gamp = Comrade.GainPrior(distamp, st),
          gphase = Comrade.GainPrior(distphase, st)
        )
# Now form the posterior

mms  = Test(st)

cllklhd = RadioLikelihood(dlcamp, dcphase)
clpost = Posterior(cllklhd, prior, mms)

lklhd = RadioLikelihood(dvis)#damp, dcphase)
post = Posterior(lklhd, prior, mms)
tpost = asflat(post)
# We will use HMC to sample the posterior.
# First to reduce burn in we use pathfinder

ndim = dimension(tpost)
f = OptimizationFunction(tpost, GalacticOptim.AutoForwardDiff())
x0 = Comrade.HypercubeTransform.inverse(tpost, rand(post.prior))

prob = OptimizationProblem(f, x0, nothing, lb=fill(-5.0, ndim), ub=fill(5.0, ndim))
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxiters=1_000_000)

prob = GalacticOptim.OptimizationProblem(f, 0.1*randn(ndim), nothing)
sol = solve(prob, LBFGS(); g_tol=1e-2, maxiters=2_000)

@info sol.minimum

xopt = transform(tpost, sol)

residual(mms(xopt), dvis)

metric = DenseEuclideanMetric(ndim)
chain, stats = sample(post, AHMC(;metric), 7_000; nadapts=5_000, init_params=xopt)


text:
module CLIMA

include(joinpath("Utilities","TicToc","TicToc.jl"))
include(joinpath("Utilities","ParametersType","ParametersType.jl"))
include(joinpath("Common","PlanetParameters","PlanetParameters.jl"))
include(joinpath("Utilities","RootSolvers","RootSolvers.jl"))
include(joinpath("Utilities","VariableTemplates","VariableTemplates.jl"))
include(joinpath("Common","MoistThermodynamics","MoistThermodynamics.jl"))
include(joinpath("Atmos","Parameterizations","CloudPhysics","MicrophysicsParameters.jl"))
include(joinpath("Atmos","Parameterizations","CloudPhysics","Microphysics.jl"))
include(joinpath("Atmos","Parameterizations","SurfaceFluxes","SurfaceFluxes.jl"))
include(joinpath("Atmos","Parameterizations","TurbulenceConvection","TurbulenceConvection.jl"))
include(joinpath("Atmos","Parameterizations","SubgridScaleTurbulence","SubgridScaleParameters.jl"))
include(joinpath("Arrays","MPIStateArrays.jl"))
include(joinpath("Mesh","Mesh.jl"))
include(joinpath("DGmethods","SpaceMethods.jl"))
include(joinpath("DGmethods","DGmethods.jl"))
include(joinpath("Ocean","ShallowWater","ShallowWaterModel.jl"))
include(joinpath("Ocean","HydrostaticBoussinesq","HydrostaticBoussinesqModel.jl"))
include(joinpath("Ocean","HydrostaticBoussinesq","SimpleBoxProblem.jl"))
include(joinpath("DGmethods_old","DGBalanceLawDiscretizations.jl"))
include(joinpath("LinearSolvers","LinearSolvers.jl"))
include(joinpath("LinearSolvers","GeneralizedConjugateResidualSolver.jl"))
include(joinpath("LinearSolvers","GeneralizedMinimalResidualSolver.jl"))
include(joinpath("LinearSolvers","ColumnwiseLUSolver.jl"))
include(joinpath("ODESolvers","ODESolvers.jl"))
include(joinpath("ODESolvers","GenericCallbacks.jl"))
include(joinpath("Atmos","Model","AtmosModel.jl"))
include(joinpath("Diagnostics","Diagnostics.jl"))
include(joinpath("InputOutput","VTK","VTK.jl"))
include(joinpath("Driver","Driver.jl"))
include(joinpath("Atmos","Model","courant.jl"))

end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# Reference frame transformations

export get_hci2eci_matrix

"""
Compute the rotation matrix for the HCI -> ECI transformation

jd is the Julian date
"""
function get_hci2eci_matrix(jd::Float64)

    T_JC = (jd - 2451545.0)/36525.

    eps = 0.409092802283074 - 0.000226966106587847 * T_JC -
    2.8623399732707e-9 * T_JC^2 +
    8.79645943005142e-9 * T_JC^3

    # Rotation matrix from HCI to ECI
    return [1. 0 0;
            0 cos(eps) -sin(eps);
            0 sin(eps) cos(eps)]
end


text:
@testset "broadcast" begin
    @testset "sin.(x)" begin
        @testset "rrule" begin
            x = rand(3, 3)
            y, pullback = rrule(broadcast, sin, x)
            @test y == sin.(x)
            (dself, dsin, dx) = pullback(One())
            @test dself == NO_FIELDS
            @test dsin == DoesNotExist()
            @test extern(dx) == cos.(x)

            x̄, ȳ = rand(), rand()
            ∂x = pullback(ȳ)[3]
            @test isequal(
                extern(ChainRules.accumulate(x̄, ∂x)),
                x̄ .+ ȳ .* cos.(x)
            )

            x̄, ȳ = Zero(), rand(3, 3)
            ∂x = pullback(ȳ)[3]
            @test extern(extern(accumulate(x̄, ∂x))) == ȳ .* cos.(x)
        end
        @testset "frule" begin
            x = rand(3, 3)
            y, pushforward = frule(broadcast, sin, x)
            @test y == sin.(x)

            ẏ = pushforward(NamedTuple(), NamedTuple(), One())
            @test extern(ẏ) == cos.(x)
        end
    end
end


text:
get_filtered_rules(genotype::String) = filter_rules(get_rules(genotype))

function get_rules(genotype::String)
    rules = []

    for (i, _) in enumerate(genotype)
        rule = _get_rule_at(genotype, i)

        if rule !== nothing
            push!(rules, rule)
        end
    end

    return rules
end

function filter_rules(rules)
    if length(rules) <= 1 return rules end

    rules = unique(rules)

    all_chosen_rules = []

    currently_chosen_rules = []

    # add best-matching rules for axiom network

    for node in AXIOM_NODES
        append!(currently_chosen_rules, find_best_matches(node, rules))
    end

    currently_chosen_rules = unique(currently_chosen_rules)
    currently_chosen_rules = _append_non_overlapping!(all_chosen_rules, currently_chosen_rules)

    # repeatedly add the best matching predecessor for each successor

    while 0 < length(currently_chosen_rules)
        last_chosen_rules = currently_chosen_rules
        currently_chosen_rules = []

        for chosen_rule in last_chosen_rules
            append!(currently_chosen_rules, find_best_matches(chosen_rule.successor_1, rules))
            append!(currently_chosen_rules, find_best_matches(chosen_rule.successor_2, rules))
        end

        currently_chosen_rules = unique(currently_chosen_rules)
        setdiff!(currently_chosen_rules, all_chosen_rules)
        currently_chosen_rules = _append_non_overlapping!(all_chosen_rules, currently_chosen_rules)
    end

    all_chosen_rules = unique(all_chosen_rules)

    return all_chosen_rules
end

function _append_non_overlapping!(existing_rules, rules_to_add)
    rules_added = []

    for to_add in rules_to_add
        exists_overlapping_rule = false

        for existing in existing_rules
            if _overlap(to_add, existing)
                exists_overlapping_rule = true
                break
            end
        end

        if !exists_overlapping_rule
            push!(existing_rules, to_add)
            push!(rules_added, to_add)
        end
    end

    return rules_added
end

function _overlap(rule_1, rule_2)
    rule_1.start_index <= rule_2.start_index <= rule_1.end_index ||
        rule_2.start_index <= rule_1.start_index <= rule_2.end_index
end

find_best_matches(node::Node, rules) = find_best_matches(node.string, rules)
function find_best_matches(node::String, rules)
    best_matches = []
    best_match_length = 0

    for rule in rules
        if _is_match(rule, node)
            if best_match_length == length(rule)
                push!(best_matches, rule)
            elseif best_match_length < length(rule)
                best_match_length = length(rule)
                best_matches = [rule]
            end
        end
    end

    return best_matches
end

_is_match(rule::Rule, node) = _is_match(rule.predecessor, node)
_is_match(rule, node) = startswith(node, rule)

function _get_rule_at(genotype, position)
    start = position

    len_genotype = length(genotype)

    if len_genotype <= position + 2 return nothing end

    #  valid rule starts with 11
    if genotype[position:position + 1] != "11" return nothing end

    position += 2
    if len_genotype < position return nothing end

    # read predecessor
    predecessor = _segment_at(genotype, position)
    position += 2*length(predecessor) + 1
    if len_genotype < position return nothing end

    # read successor 1
    successor_1 = _segment_at(genotype, position)
    position += 2*length(successor_1) + 1
    if len_genotype < position return nothing end

    # read successor 2
    successor_2 = _segment_at(genotype, position)
    position += 2*length(successor_2) + 1
    if len_genotype < position return nothing end

    # read bits
    if len_genotype < position + 5 return nothing end
    bits = genotype[position:position + 5]

    return Rule(predecessor, successor_1, successor_2,
                _to_bool(bits[1]), _to_bool(bits[2]),
                _to_bool(bits[3]), _to_bool(bits[4]),
                _to_bool(bits[5]), _to_bool(bits[6]),
                start, position + 5)
end

function _segment_at(genotype, start)
    segment = ""

    for i in start:2:(length(genotype) - 1)
        if genotype[i] == '1'
            return segment
        end

        segment = segment * genotype[i + 1]
    end

    return segment
end

_to_bool(x::Char) = x == '1'

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
const ConcreteFloat = Union{Float64,Float32,Float16,BigFloat}
const ConcreteInt = Union{Bool,Int8,Int16,Int32,Int64,Int128,BigInt}

abstract type AbstractParticles{T,N} <: Real end
"""
    struct Particles{T, N} <: AbstractParticles{T, N}

This type represents uncertainty using a cloud of particles.
# Constructors:
- `Particles()`
- `Particles(N::Integer)`
- `Particles(d::Distribution)`
- `Particles(N::Integer, d::Distribution)`
- `Particles(N::Integer, d::Distribution; permute, systematic)`
- `Particles(v::Array{T,1} where T)`
- `Particles(m::Array{T,2} where T)`
"""
struct Particles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
end

"""
    struct StaticParticles{T, N} <: AbstractParticles{T, N}

See `?Particles` for help. The difference between `StaticParticles` and `Particles` is that the `StaticParticles` store particles in a static vecetor. This makes runtimes much shorter, but compile times longer. See the documentation for some benchmarks.
"""
struct StaticParticles{T,N} <: AbstractParticles{T,N}
    particles::SArray{Tuple{N}, T, 1, N}
end

"""
Particles with weights.
To weight the particles `p`, modify the field `p.logweights`. You can resample the particles using `resample!(p)`, where each particles is resampled with a probability proportional to its weight.
"""
struct WeightedParticles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
    # weights::Vector{T}
    logweights::Vector{T}
end
function WeightedParticles{T,N}(v::AbstractVector) where {T,N}
    # weights = fill(1/N, N)
    logweights = fill(-log(N), N)
    WeightedParticles{T,N}(v,logweights)
end

const MvParticles = Vector{<:AbstractParticles} # This can not be AbstractVector since it causes some methods below to be less specific than desired
const MvWParticles = Vector{<:WeightedParticles}


text:
## The first set of functions in this file are baseline functions that
##  modify a matrix representing the Green's function, operating at
##  complex frequency

"""

    GFPECGG!(GscaGG::AbstractArray{Complex{FT}, 2},
             startidx1::Integer, startidx2::Integer,
             freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
             posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
             latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
             alphafield::FT, alphasource::FT, numdims::Integer) where FT<:AbstractFloat

Fill `GscaGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 scattering
Green's function interaction tensor above a perfect electrically conducting (PEC) plane
(coplanar with the xy-plane) from the vacuum Green's function via image theory, at frequency
`freq` and wavevector `blochk` (a 3-element vector), corresponding to source position
`possource` and field position `posfield` (both 3-element vectors) described by Gaussian
basis functions parameterized by real atomic polarizabilities `alphasource` and `alphafield`
respectively, using the appropriate function for `numdims` periodic dimensions (with lattice
parameters `latticevecs` for the real lattice vectors and `reciprocalvecs` for the
reciprocal lattice vectors as appropriate). `numdims` must be 0, 1, or 2.

"""
function GFPECGG!(GscaGG::AbstractArray{Complex{FT}, 2},
                  startidx1::Integer, startidx2::Integer,
                  freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
                  alphafield::FT, alphasource::FT, numdims::Integer) where FT<:AbstractFloat

    if (numdims == 0)
        GFPECGG0!(GscaGG, startidx1, startidx2, freq, posfield, possource, alphafield, alphasource);
    elseif (numdims == 1)
        GFPECGG1Ewald!(GscaGG, startidx1, startidx2, freq, blochk, posfield, possource, latticevecs, reciprocalvecs, alphafield, alphasource);
    elseif (numdims == 2)
        GFPECGG2Ewald!(GscaGG, startidx1, startidx2, freq, blochk, posfield, possource, latticevecs, reciprocalvecs, alphafield, alphasource);
    else
        error("GFPECGG!(): Periodic dimensionality numdims (currently ", numdims,
              ") must be 0 (no periodicity), 1, or 2");
    end

end

"""

    GFPECGG0!(GscaGG::AbstractArray{Complex{FT}, 2},
              startidx1::Integer, startidx2::Integer,
              freq::Union{FT, Complex{FT}},
              posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
              alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Fill `GscaGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 scattering
Green's function interaction tensor above a perfect electrically conducting (PEC) plane
(coplanar with the xy-plane) from the vacuum Green's function via image theory, at frequency
`freq` without spatial periodicity, corresponding to source position `possource` and field
position `posfield` (both 3-element vectors) described by Gaussian basis functions
parameterized by real atomic polarizabilities `alphasource` and `alphafield` respectively.

"""
function GFPECGG0!(GscaGG::AbstractArray{Complex{FT}, 2},
                   startidx1::Integer, startidx2::Integer,
                   freq::Union{FT, Complex{FT}},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (posfield[3] < zero(eltype(posfield)))
        error("GFPECGG0!(): can only handle field positions z (currently ",
              posfield[3], ") larger than 0");
    end
    if (possource[3] < zero(eltype(possource)))
        error("GFPECGG0!(): can only handle source positions z' (currently ",
              possource[3], ") larger than 0");
    end

    possource[3] *= -1;
    GFVACGG0!(GscaGG, startidx1, startidx2, freq,
              posfield, possource, alphafield, alphasource);

    for jj=1:2
        for ii=1:3
            GscaGG[ii, jj] *= -1;
        end
    end

    possource[3] *= -1;

end

"""

    GFPECGG1Ewald!(GscaGG::AbstractArray{Complex{FT}, 2},
                   startidx1::Integer, startidx2::Integer,
                   freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   latticevecs::AbstractArray{FT, 2},
                   reciprocalvecs::AbstractArray{FT, 2},
                   alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Fill `GscaGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 scattering
Green's function interaction tensor above a perfect electrically conducting (PEC) plane
(coplanar with the xy-plane) from the vacuum Green's function via image theory, at frequency
`freq` and wavevector `blochk` (a 3-element vector), corresponding to source position
`possource` and field position `posfield` (both 3-element vectors) described by Gaussian
basis functions parameterized by real atomic polarizabilities `alphasource` and `alphafield`
respectively, in 1 periodic dimension (with lattice parameters `latticevecs` for the real
lattice vectors and `reciprocalvecs` for the reciprocal lattice vectors as appropriate).
This is done via Ewald summation over Gaussian basis functions in 1 dimension.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFPECGG1Ewald!(GscaGG::AbstractArray{Complex{FT}, 2},
                        startidx1::Integer, startidx2::Integer,
                        freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                        posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                        latticevecs::AbstractArray{FT, 2},
                        reciprocalvecs::AbstractArray{FT, 2},
                        alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (posfield[3] < zero(eltype(posfield)))
        error("GFPECGG1Ewald!(): can only handle field positions z (currently ",
              posfield[3], ") larger than 0");
    end
    if (possource[3] < zero(eltype(possource)))
        error("GFPECGG1Ewald!(): can only handle source positions z' (currently ",
              possource[3], ") larger than 0");
    end
    if (abs(latticevecs[3, 1])/norm(latticevecs[:, 1]) > eps())
        error("GFPECGG1Ewald!(): 1D-periodic lattice vector must be",
              " perpendicular to e_z = [0, 0, 1]");
    end

    possource[3] *= -1;
    GFVACGG1Ewald!(GscaGG, startidx1, startidx2, freq, blochk, posfield, possource,
                   latticevecs, reciprocalvecs, alphafield, alphasource);

    for jj=1:2
        for ii=1:3
            GscaGG[ii, jj] *= -1;
        end
    end

    possource[3] *= -1;

end

"""

    GFPECGG2Ewald!(GscaGG::AbstractArray{Complex{FT}, 2},
                   startidx1::Integer, startidx2::Integer,
                   freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                   posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                   latticevecs::AbstractArray{FT, 2},
                   reciprocalvecs::AbstractArray{FT, 2},
                   alphafield::FT, alphasource::FT) where FT<:AbstractFloat

Fill `GscaGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 scattering
Green's function interaction tensor above a perfect electrically conducting (PEC) plane
(coplanar with the xy-plane) from the vacuum Green's function via image theory, at frequency
`freq` and wavevector `blochk` (a 3-element vector), corresponding to source position
`possource` and field position `posfield` (both 3-element vectors) described by Gaussian
basis functions parameterized by real atomic polarizabilities `alphasource` and `alphafield`
respectively, in 2 periodic dimensions (with lattice parameters `latticevecs` for the real
lattice vectors and `reciprocalvecs` for the reciprocal lattice vectors as appropriate).
This is done via Ewald summation over Gaussian basis functions in 2 dimensions.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFPECGG2Ewald!(GscaGG::AbstractArray{Complex{FT}, 2},
                        startidx1::Integer, startidx2::Integer,
                        freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                        posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                        latticevecs::AbstractArray{FT, 2},
                        reciprocalvecs::AbstractArray{FT, 2},
                        alphafield::FT, alphasource::FT) where FT<:AbstractFloat

    if (posfield[3] < zero(eltype(posfield)))
        error("GFPECGG2Ewald!(): can only handle field positions z (currently ",
              posfield[3], ") larger than 0");
    end
    if (possource[3] < zero(eltype(possource)))
        error("GFPECGG2Ewald!(): can only handle source positions z' (currently ",
              possource[3], ") larger than 0");
    end
    if (abs(latticevecs[3, 1])/norm(latticevecs[:, 1]) > eps() ||
        abs(latticevecs[3, 2])/norm(latticevecs[:, 2]) > eps())
        error("GFPECGG2Ewald!(): 2D-periodic lattice vectors must be perpendicular",
              " to e_z = [0, 0, 1]");
    end

    possource[3] *= -1;
    GFVACGG2Ewald!(GscaGG, startidx1, startidx2, freq, blochk, posfield, possource,
                   latticevecs, reciprocalvecs, alphafield, alphasource);

    for jj=1:2
        for ii=1:3
            GscaGG[ii, jj] *= -1;
        end
    end

    possource[3] *= -1;

end

## The second set of functions in this file consists of aliases for
##  functions with fewer arguments, useful for compact geometries

"""

    GFPECGG!(GscaGG::AbstractArray{Complex{FT}, 2},
             startidx1::Integer, startidx2::Integer,
             freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
             possource::AbstractArray{FT, 1}, alphafield::FT,
             alphasource::FT) where {FT<:AbstractFloat}

Fill `GscaGG[startidx1:startidx1+2, startidx2:startidx2+2]` in-place with the 3-by-3 scattering
Green's function interaction tensor above a perfect electrically conducting (PEC) plane
(coplanar with the xy-plane) from the vacuum Green's function via image theory, at frequency
`freq`, corresponding to source position `possource` and field position `posfield` (both
3-element vectors) described by Gaussian basis functions parameterized by real atomic
polarizabilities `alphasource` and `alphafield` respectively, for nonperiodic geometries.

"""
GFPECGG!(GscaGG::AbstractArray{Complex{FT}, 2},
         startidx1::Integer, startidx2::Integer,
         freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
         possource::AbstractArray{FT, 1}, alphafield::FT,
         alphasource::FT) where {FT<:AbstractFloat} =
GFPECGG0!(GscaGG, startidx1, startidx2, freq, posfield, possource, alphafield, alphasource);

## The third set of functions in this file consists of aliases for not
##  in-place function versions of the first and second sets
## (Note that all return types are 3-by-3 arrays with Complex{Float64}
##  elements, independent of arguments)

"""

    GFPECGG(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                 posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                 latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
                 alphafield::FT, alphasource::FT, numdims::Integer) where FT<:AbstractFloat

Return the 3-by-3 scattering Green's function interaction tensor above a perfect
electrically conducting (PEC) plane (coplanar with the xy-plane) from the vacuum Green's
function via image theory, at frequency `freq` and wavevector `blochk` (a 3-element vector),
corresponding to source position `possource` and field position `posfield` (both 3-element
vectors) described by Gaussian basis functions parameterized by real atomic polarizabilities
`alphasource` and `alphafield` respectively, using the appropriate function for `numdims`
periodic dimensions (with lattice parameters `latticevecs` for the real lattice vectors and
`reciprocalvecs` for the reciprocal lattice vectors as appropriate). `numdims` must be 0, 1,
or 2.

"""
function GFPECGG(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                 posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                 latticevecs::AbstractArray{FT, 2}, reciprocalvecs::AbstractArray{FT, 2},
                 alphafield::FT, alphasource::FT, numdims::Integer) where FT<:AbstractFloat

    GscaGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFPECGG!(GscaGG, 1, 1, freq, blochk, posfield, possource,
             latticevecs, reciprocalvecs, alphafield, alphasource, numdims);
    return GscaGG;

end

"""

    GFPECGG0(freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
             possource::AbstractArray{FT, 1}, alphafield::FT,
             alphasource::FT) where FT<:AbstractFloat

Return the 3-by-3 scattering Green's function interaction tensor above a perfect
electrically conducting (PEC) plane (coplanar with the xy-plane) from the vacuum Green's
function via image theory, at frequency `freq` without spatial periodicity, corresponding to
source position `possource` and field position `posfield` (both 3-element vectors) described
by Gaussian basis functions parameterized by real atomic polarizabilities `alphasource` and
`alphafield` respectively.

"""
function GFPECGG0(freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
                  possource::AbstractArray{FT, 1}, alphafield::FT,
                  alphasource::FT) where FT<:AbstractFloat

    GscaGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFPECGG0!(GscaGG, 1, 1, freq, posfield, possource, alphafield, alphasource);
    return GscaGG;

end

"""

    GFPECGG1Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  latticevecs::AbstractArray{FT, 2},
                  reciprocalvecs::AbstractArray{FT, 2}, alphafield::FT,
                  alphasource::FT) where FT<:AbstractFloat

Return the 3-by-3 scattering Green's function interaction tensor above a perfect
electrically conducting (PEC) plane (coplanar with the xy-plane) from the vacuum Green's
function via image theory, at frequency `freq` and wavevector `blochk` (a 3-element vector),
corresponding to source position `possource` and field position `posfield` (both 3-element
vectors) described by Gaussian basis functions parameterized by real atomic polarizabilities
`alphasource` and `alphafield` respectively, in 1 periodic dimension (with lattice parameters
`latticevecs` for the real lattice vectors and `reciprocalvecs` for the reciprocal lattice
vectors as appropriate). This is done via Ewald summation over Gaussian basis functions in 1
dimension.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFPECGG1Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                       posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                       latticevecs::AbstractArray{FT, 2},
                       reciprocalvecs::AbstractArray{FT, 2}, alphafield::FT,
                       alphasource::FT) where FT<:AbstractFloat

    GscaGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFPECGG1Ewald!(GscaGG, 1, 1, freq, blochk, posfield, possource,
                   latticevecs, reciprocalvecs, alphafield, alphasource);
    return GscaGG;

end

"""

    GFPECGG2Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                  posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                  latticevecs::AbstractArray{FT, 2},
                  reciprocalvecs::AbstractArray{FT, 2}, alphafield::FT,
                  alphasource::FT) where FT<:AbstractFloat

Return the 3-by-3 scattering Green's function interaction tensor above a perfect
electrically conducting (PEC) plane (coplanar with the xy-plane) from the vacuum Green's
function via image theory, at frequency `freq` and wavevector `blochk` (a 3-element vector),
corresponding to source position `possource` and field position `posfield` (both 3-element
vectors) described by Gaussian basis functions parameterized by real atomic polarizabilities
`alphasource` and `alphafield` respectively, in 2 periodic dimensions (with lattice
parameters `latticevecs` for the real lattice vectors and `reciprocalvecs` for the reciprocal
lattice vectors as appropriate). This is done via Ewald summation over Gaussian basis
functions in 2 dimensions.

!!! warning

    The Ewald summation convergence parameters are hard-coded into this function, though
    they should work for most common cases.

"""
function GFPECGG2Ewald(freq::Union{FT, Complex{FT}}, blochk::AbstractArray{FT, 1},
                       posfield::AbstractArray{FT, 1}, possource::AbstractArray{FT, 1},
                       latticevecs::AbstractArray{FT, 2},
                       reciprocalvecs::AbstractArray{FT, 2}, alphafield::FT,
                       alphasource::FT) where FT<:AbstractFloat

    GscaGG = Array{Complex{FT}, 2}(undef, 3, 3);
    GFPECGG2Ewald!(GscaGG, 1, 1, freq, blochk, posfield, possource,
                   latticevecs, reciprocalvecs, alphafield, alphasource);
    return GscaGG;

end

"""

    GFPECGG(freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
            possource::AbstractArray{FT, 1}, alphafield::FT,
            alphasource::FT) where {FT<:AbstractFloat}

Return the 3-by-3 scattering Green's function interaction tensor above a perfect
electrically conducting (PEC) plane (coplanar with the xy-plane) from the vacuum Green's
function via image theory, at frequency `freq`, corresponding to source position `possource`
and field position `posfield` (both 3-element vectors) described by Gaussian basis functions
parameterized by real atomic polarizabilities `alphasource` and `alphafield` respectively,
for nonperiodic geometries.

"""
GFPECGG(freq::Union{FT, Complex{FT}}, posfield::AbstractArray{FT, 1},
         possource::AbstractArray{FT, 1}, alphafield::FT,
         alphasource::FT) where {FT<:AbstractFloat} =
GFPECGG0(freq, posfield, possource, alphafield, alphasource);


text:
"""
	FrustumMap(viewport::Matrix4,projection::Matrix4,modelview::Matrix4)

`mutable struct` to contain the current values of matrices of viewing pipeline.

The constructor `FrustumMap()` returns the `viewport`, `projection`, `modelview` matrices and their inverses.

# Example
```
function projectPoint(map::FrustumMap,p3::Point3d)
	p4=(map.viewport * (map.projection * (map.modelview * Point4d(p3[1],p3[2],p3[3],1.0))))
	return Point3d(p4[1]/p4[4],p4[2]/p4[4],p4[3]/p4[4])
end
```
"""
mutable struct FrustumMap

	viewport::Matrix4
	projection::Matrix4
	modelview::Matrix4

	inv_viewport::Matrix4
	inv_projection::Matrix4
	inv_modelview::Matrix4

	# constructor
	function FrustumMap(viewport,projection::Matrix4,modelview::Matrix4)
		x,y,w,h = viewport
		viewport_T=Matrix4(
			w/2.0,   0.0,   0.0, x+w/2.0,
			  0.0, h/2.0,   0.0, y+h/2.0,
			  0.0,   0.0, 1/2.0,   1/2.0,
			  0.0,   0.0,   0.0,     1.0)
		new(viewport_T,projection,modelview,inv(viewport_T),inv(projection),inv(modelview))
	end

end

"""
	projectPoint(map::FrustumMap,p3::Point3d)

Return in 3D screen coordinates the *point of projection*, given in 3D world coordinates.
"""
function projectPoint(map::FrustumMap,p3::Point3d)
	p4=(map.viewport * (map.projection * (map.modelview * Point4d(p3[1],p3[2],p3[3],1.0))))
	return Point3d(p4[1]/p4[4],p4[2]/p4[4],p4[3]/p4[4])
end

"""
	nprojectPoint(map::FrustumMap,x::Float32,y::Float32, z::Float32)

Return in 3D world coordinates a `GL.Point3d`, given in 3D screen coordinates.
Use the inverse matrices of viewing transformation, contained in `map::FrustumMap`.
"""
function unprojectPoint(map::FrustumMap,x::Float64,y::Float64, z::Float64)
	p4 = (map.inv_modelview * (map.inv_projection * (map.inv_viewport * Point4d(x,y,z, 1.0))))
	return Point3d(p4[1]/p4[4],p4[2]/p4[4],p4[3]/p4[4])
end


text:
# A Julia-native gamma function. Pretty darn accurate and fast, and using this
# instead of SpecialFunctions.gamma would mean that the dependency can be
# removed. That doesn't seem worth it to me at this point in time, but putting
# this here in case.

const _G  = 607/128
const _C1 = 0.99999999999999709182
const _C  = (57.156235665862923517,
            -59.597960355475491248,
             14.136097974741747174,
            -0.49191381609762019978,
            .33994649984811888699e-4,
            .46523628927048575665e-4,
           -.98374475304879564677e-4,
            .15808870322491248884e-3,
           -.21026444172410488319e-3,
            .21743961811521264320e-3,
           -.16431810653676389022e-3,
            .84418223983852743293e-4,
           -.26190838401581408670e-4,
            .36899182659531622704e-5)
const _Cl = length(_C)
const SQRT_2_PI = sqrt(2*pi)

# Lanczos approximation.
@inline function _gamma_pos(_z::T) where{T}
  T_one    = one(T)
  T_half   = T_one/(T_one+T_one)
  z_half   = _z-T_half
  z_g_half = z_half+_G
  out = 0.99999999999999709182
  zpj = _z
  @inbounds for c in _C
    out += c/zpj
    zpj += T_one
  end
  out*SQRT_2_PI*exp(z_half*log(z_g_half) - z_g_half)
end

function _gamma(x)
  iszero(x) && return Inf
  (isinteger(x) && x < zero(x)) && return NaN
  x > zero(x) && return _gamma_pos(x)
  flx = floor(x)
  xf  = x - flx
  g   = _gamma_pos(xf)
  _x  = xf-one(x)
  for _ in 1:Int(-flx)
    g  /= _x
    _x -= one(x)
  end
  g
end


text:
include("calculus.jl")
include("expralgebra.jl")

function parsefoc(foc1,vlist,dlist,plist)
    foc=subs(addindex!(subs(foc1,plist)),dlist)
    @assert foc.head==:vcat || foc.head==:vect
    list = :([])
    for i = 1:length(foc.args)
        foc.args[i],list= getexpectation(foc.args[i],list,length(list.args)+1)
    end

    foc2 = deepcopy(foc)
    subs!(foc2,Dict(zip([Expr(:ref,:Expect,i) for i = 1:length(list.args)],list.args)))
    filter!((k,v)->k.args[2]<1,dlist)
    j2  = jacobian(foc2,[Expr(:ref,v,0) for v in names(Policy(vlist))])
    foc,foc2,j2,list
end

function getexpectation(x,list,ieq)
    if isa(x,Expr)
        if x.head==:call && x.args[1]==:Expect
            push!(list.args,:ProbWeights*x.args[2])
            x.head = :ref
            x.args[2] = ieq
        else
            for i = 1:length(x.args)
                x.args[i],list = getexpectation(x.args[i],list,ieq)
            end
        end
    end
    return x,list
end

function addpweights!(x,nP)
  if typeof(x) == Expr
    if (x.head == :call) && (x.args[1] == :*) && (x.args[2]==:ProbWeights)
      return Expr(:call,:+,[subs(x.args[3],:j=>j)*:(M.ProbWeights[i,$j]) for j = 1:nP]...)
    #   return sum([subs(x.args[3],:j=>j)*:(M.future.P[i,$j]) for j = 1:nP])
    else
      for i = 1:length(x.args)
        x.args[i]=addpweights!(x.args[i],nP)
      end
    end
  end
  return x
end

addpweights(x,nP) = addpweights(deepcopy(x),nP)

function inserthardloc(ex::Expr,vlist::Vector{ModVar})
    list2 = Dict()
    cnt = 0
    for s ∈ Endogenous(vlist)
        cnt+=1
        push!(list2,Expr(:ref,s.name,-1)=>:(M.S[i,$cnt]))
    end
    for s ∈ Stochastic(vlist)
        cnt+=1
        push!(list2,Expr(:ref,s.name,0)=>:(M.S[i,$cnt]))
        push!(list2,Expr(:ref,s.name,1)=>:(M.SP[i+(j-1)*length(M.G),$cnt]))
    end
    cnt=0
    for s ∈ Future(State(vlist,false),false)
        cnt+=1
        if vtype(s)==Aggregate
            push!(list2,Expr(:ref,s.name,s.val.target[2])=>:(M.X[i,$cnt]))
        else
            push!(list2,Expr(:ref,s.name,0)=>:(M.X[i,$cnt]))
        end
    end
    cnt=0
    for s ∈ Future(vlist)
        cnt+=1
        push!(list2,Expr(:ref,s.name,1)=>:(M.XP[i+(j-1)*length(M.G),$cnt]))
    end
    return subs(ex,list2)
end

function buildfunc(ex::Expr,targ,t=0)
    F = Expr(:for,:(i=1:length(M.G)),Expr(:block))
    for i = 1:length(ex.args)
        push!(F.args[2].args,:($targ[i,$(i+t)] = $(ex.args[i])))
    end
    return F
end

function buildJ(vJ)
    ex = Expr(:block)
    for i = 1:length(vJ.args)
        push!(ex.args,:(M.temp.J[$i] = $(vJ.args[i])))
    end
    push!(ex.args,:(return))
    return ex
end


text: include("projectednewton.jl")


text:
module QMPS
using Yao
using Yao.ConstGate: SWAP
using BitBasis: packbits

using StatsBase
using StatsBase: mean

using LinearAlgebra
using KrylovKit

using QuAlgorithmZoo
PauliGate{T} = Union{XGate{T}, YGate{T}, ZGate{T}}

include("Adam.jl")
include("Core.jl")
include("AbstractModel.jl")
include("gradient.jl")
include("correlation.jl")
include("ansatz/ansatz.jl")
end

max_stars_repo_name: GiggleLiu/QuantumMPS

max_issues_repo_name: GiggleLiu/QuantumMPS

max_forks_repo_name: GiggleLiu/QuantumMPS


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

JLLWrappers.@generate_main_file_header("MozillaCACerts")
JLLWrappers.@generate_main_file("MozillaCACerts", UUID("14a3606d-f60d-562e-9121-12d972cd8159"))
end  # module MozillaCACerts_jll


text:
# evaluate right bound of scan_func
function get_right_endpoint(
    theta_init::Vector{Float64}, # initial point of parameters
    scan_func::Function, # h(theta) function for predictions or parameters
    loss_func::Function, # lambda(theta) - labmbda_min - delta_lambda
    method::Val{:CICO_ONE_PASS}; # function works only for method ONE_PASS;

    theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
        (-Inf, Inf), length(theta_init)
        ),
    scan_bound::Float64 = 9.0,
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 1e-3, # i do not know how to use it
    # good results in :LN_NELDERMEAD, :LN_COBYLA, :LN_PRAXIS,
    # errors in :LN_BOBYQA, :LN_SBPLX, :LN_NEWUOA
    local_alg::Symbol = :LN_NELDERMEAD,
    # options for local fitter :max_iter
    max_iter::Int = 10^5,
    ftol_abs::Float64 = 1e-3,
    kwargs...
    )
    # dim of the theta vector
    n_theta = length(theta_init)

    # checking arguments
    # methods which are not supported
    if local_alg in [:LN_BOBYQA, :LN_SBPLX, :LN_NEWUOA]
        @warn "Using local_alg = :"*String(local_alg)*" may result in wrong output."
    end
    # when using :LN_NELDERMEAD initial parameters should not be zero
    if local_alg == :LN_NELDERMEAD
        zeroParameter = [ isapprox(theta_init[i], 0., atol=1e-2) for i in 1:n_theta]
        if any(zeroParameter)
            @warn "Close-to-zero parameters found when using :LN_NELDERMEAD."
            show(findall(zeroParameter))
        end
    end

    # optimizer
    local_opt = Opt(local_alg, n_theta)
    ftol_abs!(local_opt, scan_tol) #ftol_abs

    # Constraints function
    out_of_bound::Bool = false
    function constraints_func(x, g)
        # this part is necessary to understand the difference between
        # "stop out of bounds" and "stop because of function call error"
        try
            loss = loss_func(x)
        catch e
            @warn "Error when call loss_func($x)"
            throw(e)
        end

        if (loss < 0.) && (scan_func(x) > scan_bound)
            out_of_bound = true
            throw(ForcedStop("Out of the scan bound but in ll constraint."))
        #elseif isapprox(loss, 0., atol=loss_tol)
            #@warn "loss_tol reached... but..."
            #return loss
        end

        return loss
    end

    # constrain optimizer
    opt = Opt(:LN_AUGLAG, n_theta)
    ftol_abs!(opt, scan_tol)
    max_objective!(
        opt,
        (x, g) -> scan_func(x)
        )
    local_optimizer!(opt, local_opt)
    maxeval!(opt, max_iter)

    # inequality constraints
    inequality_constraint!(
        opt,
        constraints_func,
        loss_tol
    )
    [ inequality_constraint!(
        opt,
        (x, g) -> x[i] - theta_bounds[i][2],
        0.
    ) for i in 1:n_theta ]
    [ inequality_constraint!(
        opt,
        (x, g) -> theta_bounds[i][1] - x[i],
        0.
    ) for i in 1:n_theta ]

    # start optimization
    (optf, optx, ret) = optimize(opt, theta_init)

    if (ret == :FORCED_STOP && !out_of_bound)
        pp = ProfilePoint[]
        res = (nothing, pp, :LOSS_ERROR_STOP)
    elseif ret == :MAXEVAL_REACHED
        pp = ProfilePoint[]
        res = (nothing, pp, :MAX_ITER_STOP)
    elseif (ret == :FORCED_STOP && out_of_bound) # successfull result
        pp = ProfilePoint[]
        res = (nothing, pp, :SCAN_BOUND_REACHED)
    elseif ret == :FTOL_REACHED # successfull result
        loss = loss_func(optx)
        pp = [ ProfilePoint(optf, loss, optx, ret, nothing) ]
        res = (optf, pp, :BORDER_FOUND_BY_SCAN_TOL)
    else
        # this part is not normally reached, just for case
        throw(ErrorException("No interpretation of the optimization results."))
        # do not throw
        #pp = ProfilePoint[]
        #res = (nothing, pp, :UNKNOWN_STOP)
    end

    return res
end # of bound_right

function get_right_endpoint(
    theta_init::Vector{Float64}, # initial point of parameters
    theta_num::Int, # number of parameter to scan
    loss_func::Function, # lambda(theta) - labmbda_min - delta_lambda
    method::Val{:CICO_ONE_PASS}; # function works only for method ONE_PASS;

    theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
        (-Inf, Inf), length(theta_init)
        ),
    scan_bound::Float64 = 9.0,
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 1e-3,
    local_alg::Symbol = :LN_NELDERMEAD,
    kwargs... # options for local fitter
    )
    # checking arguments
    if theta_num > length(theta_init)
        throw(DomainError(theta_num, "theta_num exceed theta dimention"))
    end

    function scan_func(theta::Vector{Float64})
        theta[theta_num]
    end

    get_right_endpoint(
        theta_init,
        scan_func,
        loss_func,
        method;

        theta_bounds = theta_bounds,
        scan_bound = scan_bound,
        scan_tol = scan_tol,
        loss_tol = loss_tol,
        local_alg = local_alg,
        kwargs... # options for local fitter
    )
end


text:
module TestDataTypes

using Mutabilities
using Test

@testset begin
    x = meltproperties(1 + 2im)
    @test (x.re, x.im) == (1, 2)
    x.re *= 100
    @test freezeproperties(x) === 100 + 2im

    @test sprint(show, meltproperties(1 + 2im)) == "meltproperties(1 + 2im)"
    @test sprint(show, "text/plain", meltproperties(1 + 2im)) == "mutable handle to 1 + 2im"
end

end  # module