text:
# This API is imported from MathProgBase and is meant to be transitional.
# A more ideal API would pass expressions (with vector-valued nodes) directly
# let the solver call out to AD tools.

# NLP = Nonlinear Programming.
# The non-MOI-like name is intentional because this part of the API is not
# MOI-like.

"""
    AbstractNLPEvaluator

Abstract supertype for the callback object used in `NLPBlock`.
"""
abstract type AbstractNLPEvaluator end

"""
    NLPBlock()

Holds the `NLPBlockData` that represents a set of nonlinear constraints, and
optionally a nonlinear objective.
"""
struct NLPBlock <: AbstractModelAttribute end

"""
    NLPBlockDual(N)
    NLPBlockDual()

The Lagrange multipliers on the constraints from the `NLPBlock` in result `N`.
If `N` is omitted, it is 1 by default.
"""
struct NLPBlockDual <: AbstractModelAttribute
    N::Int
end
NLPBlockDual() = NLPBlockDual(1)
_result_index_field(attr::NLPBlockDual) = attr.N

is_set_by_optimize(::NLPBlockDual) = true

"""
    NLPBlockDualStart()

An initial assignment of the Lagrange multipliers on the constraints from the
`NLPBlock` that the solver may use to warm-start the solve.
"""
struct NLPBlockDualStart <: AbstractModelAttribute end

"""
    NLPBoundsPair(lower,upper)

A struct holding a pair of lower and upper bounds. `-Inf` and `Inf`
can be used to indicate no lower or upper bound, respectively.
"""
struct NLPBoundsPair
    lower::Float64
    upper::Float64
end

"""
    struct NLPBlockData
        constraint_bounds::Vector{NLPBoundsPair}
        evaluator::AbstractNLPEvaluator
        has_objective::Bool
    end

A `struct` encoding a set of nonlinear constraints of the form ``lb \\le g(x)
\\le ub`` and, if `has_objective == true`, a nonlinear objective function
``f(x)``. `constraint_bounds` holds the pairs of ``lb`` and ``ub`` elements.
Nonlinear objectives *override* any objective set by using the
`ObjectiveFunction` attribute. The `evaluator` is a callback object that is used
to query function values, derivatives, and expression graphs. If
`has_objective == false`, then it is an error to query properties of
the objective function, and in Hessian-of-the-Lagrangian queries, `σ` must be
set to zero. Throughout the evaluator, all variables are ordered according to
ListOfVariableIndices(). """
struct NLPBlockData
    constraint_bounds::Vector{NLPBoundsPair}
    evaluator::AbstractNLPEvaluator
    has_objective::Bool
end

"""
    initialize(d::AbstractNLPEvaluator, requested_features::Vector{Symbol})

Must be called before any other methods. The vector `requested_features`
lists features requested by the solver. These may include `:Grad` for gradients
of ``f``, `:Jac` for explicit Jacobians of ``g``, `:JacVec` for
Jacobian-vector products, `:HessVec` for Hessian-vector
and Hessian-of-Lagrangian-vector products, `:Hess` for explicit Hessians and
Hessian-of-Lagrangians, and `:ExprGraph` for expression graphs.
"""
function initialize end

"""
    features_available(d::AbstractNLPEvaluator)

Returns the subset of features available for this problem instance, as a
list of symbols in the same format as in `initialize`.
"""
function features_available end

"""
    eval_objective(d::AbstractNLPEvaluator, x)

Evaluate the objective ``f(x)``, returning a scalar value.
"""
function eval_objective end

"""
    eval_constraint(d::AbstractNLPEvaluator, g, x)

Evaluate the constraint function ``g(x)``, storing the result in the vector `g` which must be of the
appropriate size.
"""
function eval_constraint end

"""
    eval_objective_gradient(d::AbstractNLPEvaluator, g, x)

Evaluate ``\\nabla f(x)`` as a dense vector, storing the result in the vector
`g` which must be of the appropriate size.
"""
function eval_objective_gradient end

"""
    jacobian_structure(d::AbstractNLPEvaluator)::Vector{Tuple{Int64,Int64}}

Returns the sparsity structure of the Jacobian matrix
``J_g(x) = \\left[ \\begin{array}{c} \\nabla g_1(x) \\\\ \\nabla g_2(x) \\\\ \\vdots \\\\ \\nabla g_m(x) \\end{array}\\right]``
where ``g_i`` is the ``i\\text{th}`` component of ``g``. The sparsity structure
is assumed to be independent of the point ``x``. Returns a vector of tuples,
`(row, column)`, where each indicates the position of a structurally nonzero element.
These indices are not required to be sorted and can contain duplicates, in which
case the solver should combine the corresponding elements by adding them together.
"""
function jacobian_structure end

"""
    hessian_lagrangian_structure(d::AbstractNLPEvaluator)::Vector{Tuple{Int64,Int64}}

Returns the sparsity structure of the Hessian-of-the-Lagrangian matrix
``\\nabla^2 f + \\sum_{i=1}^m \\nabla^2 g_i`` as a vector of tuples, where
each indicates the position of a structurally nonzero element. These indices are
not required to be sorted and can contain duplicates, in which case the solver
should combine the corresponding elements by adding them together. Any mix of
lower and upper-triangular indices is valid. Elements `(i,j)` and
`(j,i)`, if both present, should be treated as duplicates.
"""
function hessian_lagrangian_structure end

"""
    eval_constraint_jacobian(d::AbstractNLPEvaluator, J, x)

Evaluates the sparse Jacobian matrix
``J_g(x) = \\left[ \\begin{array}{c} \\nabla g_1(x) \\\\ \\nabla g_2(x) \\\\ \\vdots \\\\ \\nabla g_m(x) \\end{array}\\right]``.
The result is stored in the vector `J` in the same order as the indices returned
by `jacobian_structure`.
"""
function eval_constraint_jacobian end

"""
    eval_constraint_jacobian_product(d::AbstractNLPEvaluator, y, x, w)

Computes the Jacobian-vector product ``J_g(x)w``, storing the result in the vector `y`.
"""
function eval_constraint_jacobian_product end

"""
    eval_constraint_jacobian_transpose_product(d::AbstractNLPEvaluator, y, x, w)

Computes the Jacobian-transpose-vector product ``J_g(x)^Tw``, storing the result
in the vector `y`.
"""
function eval_constraint_jacobian_transpose_product end

"""
    eval_hessian_lagrangian_prod(d::AbstractNLPEvaluator, h, x, v, σ, μ)

Given scalar weight `σ` and vector of constraint weights `μ`,
computes the Hessian-of-the-Lagrangian-vector product
``\\left(\\sigma\\nabla^2 f(x) + \\sum_{i=1}^m \\mu_i \\nabla^2 g_i(x)\\right)v``,
storing the result in the vector `h`.
"""
function eval_hessian_lagrangian_product end

"""
    eval_hessian_lagrangian(d::AbstractNLPEvaluator, H, x, σ, μ)

Given scalar weight `σ` and vector of constraint weights `μ`,
computes the sparse Hessian-of-the-Lagrangian matrix
``\\sigma\\nabla^2 f(x) + \\sum_{i=1}^m \\mu_i \\nabla^2 g_i(x)``,
storing the result in the vector `H` in the same order as the indices
returned by `hessian_lagrangian_structure`.
"""
function eval_hessian_lagrangian end

"""
    objective_expr(d::AbstractNLPEvaluator)

Returns an expression graph for the objective function as a standard Julia `Expr`
object. All sums and products are flattened out as simple `Expr(:+,...)` and
`Expr(:*,...)` objects. The symbol `x` is used as a placeholder for the
vector of decision variables. No other undefined symbols are permitted;
coefficients are embedded as explicit values. For example, the expression
``x_1+\\sin(x_2/\\exp(x_3))`` would be represented as the Julia object
`:(x[1] + sin(x[2]/exp(x[3])))`. Each integer index is wrapped in a [`VariableIndex`](@ref).
See the [Julia manual](https://docs.julialang.org/en/release-0.6/manual/metaprogramming/)
for more information on the structure of `Expr` objects.
There are currently no restrictions on recognized functions; typically these
will be built-in Julia functions like `^`, `exp`, `log`, `cos`, `tan`, `sqrt`,
etc., but modeling interfaces may choose to extend these basic functions.
"""
function objective_expr end

"""
    constraint_expr(d::AbstractNLPEvaluator, i)

Returns an expression graph for the ``i\\text{th}`` constraint in the same
format as described above, with an additional comparison operator indicating
the sense of and bounds on the constraint. The right-hand side of the comparison
must be a constant; that is, `:(x[1]^3 <= 1)` is allowed, while
`:(1 <= x[1]^3)` is not valid. Double-sided constraints are allowed, in which
case both the lower bound and upper bounds should be constants; for example,
`:(-1 <= cos(x[1]) + sin(x[2]) <= 1)` is valid.
"""
function constraint_expr end


text:
using LinearAlgebra
using LDLFactorizations
using SparseArrays
using BenchmarkTools

n = 3
A = zeros(3,3)
A[1,1] = 9.0

B = 9.0*Matrix(I,n,n)

@benchmark cholesky($B)

function get_svd_sqrt(A::AbstractMatrix)
    x = svd(A)
    return Diagonal(sqrt.(x.S))*x.V'
end

function get_eig_sqrt(A::AbstractMatrix)
    x = eigen(A)
    return Diagonal(sqrt.(x.values))*x.vectors'
end

@benchmark get_svd_sqrt($B)
@benchmark get_eig_sqrt($B)
@benchmark ldlt(sparse($B); shift = 0.0, check = true, perm=Int[])

cholesky(B).U
C = SparseMatrixCSC(B)
D = sparse(B)
E = sparse(A)

F = ldlt(sparse(B); shift = 0.0, check = true, perm=Int[])

Array(sparse(cholesky(C)))

:

ext: jl


text:
@rule MatrixDirichlet(:out, Marginalisation) (m_a::PointMass, ) = begin
    return MatrixDirichlet(mean(m_a))
end

@rule MatrixDirichlet(:out, Marginalisation) (q_a::PointMass, ) = begin
    return MatrixDirichlet(mean(q_a))
end


text:
begin
    #= none:2 =# Core.@doc "    eval_interp(m::Module, ex)\n\nevaluate the interpolation operator in `ex` inside given module `m`.\n" function eval_interp(m::Module, ex)
            ex isa Expr || return ex
            if ex.head === :$
                x = ex.args[1]
                if x isa Symbol && isdefined(m, x)
                    return Base.eval(m, x)
                else
                    return ex
                end
            end
            return Expr(ex.head, map((x->begin
                                eval_interp(m, x)
                            end), ex.args)...)
        end
    #= none:20 =# Core.@doc "    eval_literal(m::Module, ex)\n\nEvaluate the literal values and insert them back to the expression.\nThe literal value can be checked via [`is_literal`](@ref).\n" function eval_literal(m::Module, ex)
            ex isa Expr || return ex
            if ex.head === :call && all(is_literal, ex.args[2:end])
                return Base.eval(m, ex)
            end
            return Expr(ex.head, map((x->begin
                                eval_literal(m, x)
                            end), ex.args)...)
        end
    replace_symbol(x::Symbol, name::Symbol, value) = begin
            if x === name
                value
            else
                x
            end
        end
    replace_symbol(x, ::Symbol, value) = begin
            x
        end
    function replace_symbol(ex::Expr, name::Symbol, value)
        Expr(ex.head, map((x->begin
                        replace_symbol(x, name, value)
                    end), ex.args)...)
    end
    #= none:41 =# Core.@doc "    subtitute(ex::Expr, old=>new)\n\nSubtitute the old symbol `old` with `new`.\n" function subtitute(ex::Expr, replace::Pair)
            (name, value) = replace
            return replace_symbol(ex, name, value)
        end
    #= none:51 =# Core.@doc "    name_only(ex)\n\nRemove everything else leaving just names, currently supports\nfunction calls, type with type variables, subtype operator `<:`\nand type annotation `::`.\n\n# Example\n\n```julia\njulia> using Expronicon\n\njulia> name_only(:(sin(2)))\n:sin\n\njulia> name_only(:(Foo{Int}))\n:Foo\n\njulia> name_only(:(Foo{Int} <: Real))\n:Foo\n\njulia> name_only(:(x::Int))\n:x\n```\n" function name_only(#= none:76 =# @nospecialize(ex))
            ex isa Symbol && return ex
            ex isa QuoteNode && return ex.value
            ex isa Expr || error("unsupported expression $(ex)")
            ex.head in [:call, :curly, :<:, :(::), :where, :function, :kw, :(=), :->] && return name_only(ex.args[1])
            ex.head === :. && return name_only(ex.args[2])
            ex.head === :module && return name_only(ex.args[2])
            error("unsupported expression $(ex)")
        end
    #= none:86 =# Core.@doc "    rm_lineinfo(ex)\n\nRemove `LineNumberNode` in a given expression.\n\n!!! tips\n\n    the `LineNumberNode` inside macro calls won't be removed since\n    the `macrocall` expression requires a `LineNumberNode`. See also\n    [issues/#9](https://github.com/Roger-luo/Expronicon.jl/issues/9).\n" function rm_lineinfo(ex)
            let
                cache_1 = nothing
                return_1 = nothing
                x_1 = ex
                if x_1 isa Expr
                    if begin
                                if cache_1 === nothing
                                    cache_1 = Some((x_1.head, x_1.args))
                                end
                                x_2 = cache_1.value
                                x_2 isa Tuple{Symbol, var2} where var2<:AbstractArray
                            end && (x_2[1] == :macrocall && (begin
                                        x_3 = x_2[2]
                                        x_3 isa AbstractArray
                                    end && ((ndims(x_3) === 1 && length(x_3) >= 2) && begin
                                            x_4 = x_3[1]
                                            x_5 = x_3[2]
                                            x_6 = (SubArray)(x_3, (3:length(x_3),))
                                            true
                                        end)))
                        return_1 = let line = x_5, name = x_4, args = x_6
                                Expr(:macrocall, name, line, map(rm_lineinfo, args)...)
                            end
                        $(Expr(:symbolicgoto, Symbol("##final#746_1")))
                    end
                    if begin
                                x_7 = cache_1.value
                                x_7 isa Tuple{var1, var2} where {var2<:AbstractArray, var1}
                            end && (begin
                                    x_8 = x_7[1]
                                    x_9 = x_7[2]
                                    x_9 isa AbstractArray
                                end && ((ndims(x_9) === 1 && length(x_9) >= 0) && begin
                                        x_10 = (SubArray)(x_9, (1:length(x_9),))
                                        true
                                    end))
                        return_1 = let args = x_10, head = x_8
                                Expr(head, map(rm_lineinfo, filter((x->begin
                                                    !(x isa LineNumberNode)
                                                end), args))...)
                            end
                        $(Expr(:symbolicgoto, Symbol("##final#746_1")))
                    end
                end
                return_1 = let
                        ex
                    end
                $(Expr(:symbolicgoto, Symbol("##final#746_1")))
                (error)("matching non-exhaustive, at #= none:98 =#")
                $(Expr(:symboliclabel, Symbol("##final#746_1")))
                return_1
            end
        end
    #= none:105 =# Base.@kwdef struct PrettifyOptions
            rm_lineinfo::Bool = true
            flatten_blocks::Bool = true
            rm_nothing::Bool = true
            rm_single_block::Bool = true
            alias_gensym::Bool = true
        end
    #= none:113 =# Core.@doc "    prettify(ex; kw...)\n\nPrettify given expression, remove all `LineNumberNode` and\nextra code blocks.\n\n# Options (Kwargs)\n\nAll the options are `true` by default.\n\n- `rm_lineinfo`: remove `LineNumberNode`.\n- `flatten_blocks`: flatten `begin ... end` code blocks.\n- `rm_nothing`: remove `nothing` in the `begin ... end`.\n- `rm_single_block`: remove single `begin ... end`.\n- `alias_gensym`: replace `##<name>#<num>` with `<name>_<id>`.\n\n!!! tips\n\n    the `LineNumberNode` inside macro calls won't be removed since\n    the `macrocall` expression requires a `LineNumberNode`. See also\n    [issues/#9](https://github.com/Roger-luo/Expronicon.jl/issues/9).\n" function prettify(ex; kw...)
            prettify(ex, PrettifyOptions(; kw...))
        end
    function prettify(ex, options::PrettifyOptions)
        ex isa Expr || return ex
        ex = if options.alias_gensym
                alias_gensym(ex)
            else
                ex
            end
        for _ = 1:10
            curr = prettify_pass(ex, options)
            ex == curr && break
            ex = curr
        end
        return ex
    end
    function prettify_pass(ex, options::PrettifyOptions)
        ex = if options.rm_lineinfo
                rm_lineinfo(ex)
            else
                ex
            end
        ex = if options.flatten_blocks
                flatten_blocks(ex)
            else
                ex
            end
        ex = if options.rm_nothing
                rm_nothing(ex)
            else
                ex
            end
        ex = if options.rm_single_block
                rm_single_block(ex)
            else
                ex
            end
        return ex
    end
    #= none:158 =# Core.@doc "    flatten_blocks(ex)\n\nRemove hierachical expression blocks.\n" function flatten_blocks(ex)
            ex isa Expr || return ex
            ex.head === :block || return Expr(ex.head, map(_flatten_blocks, ex.args)...)
            has_block = any(ex.args) do x
                    x isa Expr && x.head === :block
                end
            if has_block
                return flatten_blocks(_flatten_blocks(ex))
            end
            return Expr(ex.head, map(flatten_blocks, ex.args)...)
        end
    function _flatten_blocks(ex)
        ex isa Expr || return ex
        ex.head === :block || return Expr(ex.head, map(flatten_blocks, ex.args)...)
        args = []
        for stmt = ex.args
            if stmt isa Expr && stmt.head === :block
                for each = stmt.args
                    push!(args, flatten_blocks(each))
                end
            else
                push!(args, flatten_blocks(stmt))
            end
        end
        return Expr(:block, args...)
    end
    #= none:193 =# Core.@doc "    rm_nothing(ex)\n\nRemove the constant value `nothing` in given expression `ex`.\n" function rm_nothing(ex)
            let
                cache_2 = nothing
                return_2 = nothing
                x_11 = ex
                if x_11 isa Expr
                    if begin
                                if cache_2 === nothing
                                    cache_2 = Some((x_11.head, x_11.args))
                                end
                                x_12 = cache_2.value
                                x_12 isa Tuple{Symbol, var2} where var2<:AbstractArray
                            end && (x_12[1] == :block && (begin
                                        x_13 = x_12[2]
                                        x_13 isa AbstractArray
                                    end && ((ndims(x_13) === 1 && length(x_13) >= 0) && begin
                                            x_14 = (SubArray)(x_13, (1:length(x_13),))
                                            true
                                        end)))
                        return_2 = let args = x_14
                                Expr(:block, filter((x->begin
                                                x !== nothing
                                            end), args)...)
                            end
                        $(Expr(:symbolicgoto, Symbol("##final#760_1")))
                    end
                    if begin
                                x_15 = cache_2.value
                                x_15 isa Tuple{var1, var2} where {var2<:AbstractArray, var1}
                            end && (begin
                                    x_16 = x_15[1]
                                    x_17 = x_15[2]
                                    x_17 isa AbstractArray
                                end && ((ndims(x_17) === 1 && length(x_17) >= 0) && begin
                                        x_18 = (SubArray)(x_17, (1:length(x_17),))
                                        true
                                    end))
                        return_2 = let args = x_18, head = x_16
                                Expr(head, map(rm_nothing, args)...)
                            end
                        $(Expr(:symbolicgoto, Symbol("##final#760_1")))
                    end
                end
                return_2 = let
                        ex
                    end
                $(Expr(:symbolicgoto, Symbol("##final#760_1")))
                (error)("matching non-exhaustive, at #= none:199 =#")
                $(Expr(:symboliclabel, Symbol("##final#760_1")))
                return_2
            end
        end
    function rm_single_block(ex)
        let
            cache_3 = nothing
            return_3 = nothing
            x_19 = ex
            if x_19 isa Expr
                if begin
                            if cache_3 === nothing
                                cache_3 = Some((x_19.head, x_19.args))
                            end
                            x_20 = cache_3.value
                            x_20 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_20[1] == :(=) && (begin
                                    x_21 = x_20[2]
                                    x_21 isa AbstractArray
                                end && (ndims(x_21) === 1 && length(x_21) >= 0)))
                    return_3 = let
                            ex
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_22 = cache_3.value
                            x_22 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_22[1] == :-> && (begin
                                    x_23 = x_22[2]
                                    x_23 isa AbstractArray
                                end && (ndims(x_23) === 1 && length(x_23) >= 0)))
                    return_3 = let
                            ex
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_24 = cache_3.value
                            x_24 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_24[1] == :quote && (begin
                                    x_25 = x_24[2]
                                    x_25 isa AbstractArray
                                end && ((ndims(x_25) === 1 && length(x_25) >= 0) && begin
                                        x_26 = (SubArray)(x_25, (1:length(x_25),))
                                        true
                                    end)))
                    return_3 = let xs = x_26
                            ex
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_27 = cache_3.value
                            x_27 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_27[1] == :block && (begin
                                    x_28 = x_27[2]
                                    x_28 isa AbstractArray
                                end && (length(x_28) === 1 && (begin
                                            cache_4 = nothing
                                            x_29 = x_28[1]
                                            x_29 isa Expr
                                        end && (begin
                                                if cache_4 === nothing
                                                    cache_4 = Some((x_29.head, x_29.args))
                                                end
                                                x_30 = cache_4.value
                                                x_30 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                            end && (x_30[1] == :quote && (begin
                                                        x_31 = x_30[2]
                                                        x_31 isa AbstractArray
                                                    end && ((ndims(x_31) === 1 && length(x_31) >= 0) && begin
                                                            x_32 = (SubArray)(x_31, (1:length(x_31),))
                                                            true
                                                        end))))))))
                    return_3 = let xs = x_32
                            ex
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_33 = cache_3.value
                            x_33 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_33[1] == :try && (begin
                                    x_34 = x_33[2]
                                    x_34 isa AbstractArray
                                end && (length(x_34) === 4 && (begin
                                            cache_5 = nothing
                                            x_35 = x_34[1]
                                            x_35 isa Expr
                                        end && (begin
                                                if cache_5 === nothing
                                                    cache_5 = Some((x_35.head, x_35.args))
                                                end
                                                x_36 = cache_5.value
                                                x_36 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                            end && (x_36[1] == :block && (begin
                                                        x_37 = x_36[2]
                                                        x_37 isa AbstractArray
                                                    end && ((ndims(x_37) === 1 && length(x_37) >= 0) && (begin
                                                                x_38 = (SubArray)(x_37, (1:length(x_37),))
                                                                x_34[2] === false
                                                            end && (x_34[3] === false && (begin
                                                                        cache_6 = nothing
                                                                        x_39 = x_34[4]
                                                                        x_39 isa Expr
                                                                    end && (begin
                                                                            if cache_6 === nothing
                                                                                cache_6 = Some((x_39.head, x_39.args))
                                                                            end
                                                                            x_40 = cache_6.value
                                                                            x_40 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                                                        end && (x_40[1] == :block && (begin
                                                                                x_41 = x_40[2]
                                                                                x_41 isa AbstractArray
                                                                                end && ((ndims(x_41) === 1 && length(x_41) >= 0) && begin
                                                                                x_42 = (SubArray)(x_41, (1:length(x_41),))
                                                                                true
                                                                                end)))))))))))))))
                    return_3 = let try_stmts = x_38, finally_stmts = x_42
                            Expr(:try, Expr(:block, rm_single_block.(try_stmts)...), false, false, Expr(:block, rm_single_block.(finally_stmts)...))
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_43 = cache_3.value
                            x_43 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_43[1] == :try && (begin
                                    x_44 = x_43[2]
                                    x_44 isa AbstractArray
                                end && (length(x_44) === 3 && (begin
                                            cache_7 = nothing
                                            x_45 = x_44[1]
                                            x_45 isa Expr
                                        end && (begin
                                                if cache_7 === nothing
                                                    cache_7 = Some((x_45.head, x_45.args))
                                                end
                                                x_46 = cache_7.value
                                                x_46 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                            end && (x_46[1] == :block && (begin
                                                        x_47 = x_46[2]
                                                        x_47 isa AbstractArray
                                                    end && ((ndims(x_47) === 1 && length(x_47) >= 0) && (begin
                                                                x_48 = (SubArray)(x_47, (1:length(x_47),))
                                                                x_49 = x_44[2]
                                                                cache_8 = nothing
                                                                x_50 = x_44[3]
                                                                x_50 isa Expr
                                                            end && (begin
                                                                    if cache_8 === nothing
                                                                        cache_8 = Some((x_50.head, x_50.args))
                                                                    end
                                                                    x_51 = cache_8.value
                                                                    x_51 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                                                end && (x_51[1] == :block && (begin
                                                                            x_52 = x_51[2]
                                                                            x_52 isa AbstractArray
                                                                        end && ((ndims(x_52) === 1 && length(x_52) >= 0) && begin
                                                                                x_53 = (SubArray)(x_52, (1:length(x_52),))
                                                                                true
                                                                            end)))))))))))))
                    return_3 = let try_stmts = x_48, catch_stmts = x_53, catch_var = x_49
                            Expr(:try, Expr(:block, rm_single_block.(try_stmts)...), catch_var, Expr(:block, rm_single_block.(catch_stmts)...))
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_54 = cache_3.value
                            x_54 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_54[1] == :try && (begin
                                    x_55 = x_54[2]
                                    x_55 isa AbstractArray
                                end && (length(x_55) === 4 && (begin
                                            cache_9 = nothing
                                            x_56 = x_55[1]
                                            x_56 isa Expr
                                        end && (begin
                                                if cache_9 === nothing
                                                    cache_9 = Some((x_56.head, x_56.args))
                                                end
                                                x_57 = cache_9.value
                                                x_57 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                            end && (x_57[1] == :block && (begin
                                                        x_58 = x_57[2]
                                                        x_58 isa AbstractArray
                                                    end && ((ndims(x_58) === 1 && length(x_58) >= 0) && (begin
                                                                x_59 = (SubArray)(x_58, (1:length(x_58),))
                                                                x_60 = x_55[2]
                                                                cache_10 = nothing
                                                                x_61 = x_55[3]
                                                                x_61 isa Expr
                                                            end && (begin
                                                                    if cache_10 === nothing
                                                                        cache_10 = Some((x_61.head, x_61.args))
                                                                    end
                                                                    x_62 = cache_10.value
                                                                    x_62 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                                                end && (x_62[1] == :block && (begin
                                                                            x_63 = x_62[2]
                                                                            x_63 isa AbstractArray
                                                                        end && ((ndims(x_63) === 1 && length(x_63) >= 0) && (begin
                                                                                x_64 = (SubArray)(x_63, (1:length(x_63),))
                                                                                cache_11 = nothing
                                                                                x_65 = x_55[4]
                                                                                x_65 isa Expr
                                                                                end && (begin
                                                                                if cache_11 === nothing
                                                                                cache_11 = Some((x_65.head, x_65.args))
                                                                                end
                                                                                x_66 = cache_11.value
                                                                                x_66 isa Tuple{Symbol, var2} where var2<:AbstractArray
                                                                                end && (x_66[1] == :block && (begin
                                                                                x_67 = x_66[2]
                                                                                x_67 isa AbstractArray
                                                                                end && ((ndims(x_67) === 1 && length(x_67) >= 0) && begin
                                                                                x_68 = (SubArray)(x_67, (1:length(x_67),))
                                                                                true
                                                                                end))))))))))))))))))
                    return_3 = let try_stmts = x_59, catch_stmts = x_64, catch_var = x_60, finally_stmts = x_68
                            Expr(:try, Expr(:block, rm_single_block.(try_stmts)...), catch_var, Expr(:block, rm_single_block.(catch_stmts)...), Expr(:block, rm_single_block.(finally_stmts)...))
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_69 = cache_3.value
                            x_69 isa Tuple{Symbol, var2} where var2<:AbstractArray
                        end && (x_69[1] == :block && (begin
                                    x_70 = x_69[2]
                                    x_70 isa AbstractArray
                                end && (length(x_70) === 1 && begin
                                        x_71 = x_70[1]
                                        true
                                    end)))
                    return_3 = let stmt = x_71
                            stmt
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
                if begin
                            x_72 = cache_3.value
                            x_72 isa Tuple{var1, var2} where {var2<:AbstractArray, var1}
                        end && (begin
                                x_73 = x_72[1]
                                x_74 = x_72[2]
                                x_74 isa AbstractArray
                            end && ((ndims(x_74) === 1 && length(x_74) >= 0) && begin
                                    x_75 = (SubArray)(x_74, (1:length(x_74),))
                                    true
                                end))
                    return_3 = let args = x_75, head = x_73
                            Expr(head, map(rm_single_block, args)...)
                        end
                    $(Expr(:symbolicgoto, Symbol("##final#772_1")))
                end
            end
            return_3 = let
                    ex
                end
            $(Expr(:symbolicgoto, Symbol("##final#772_1")))
            (error)("matching non-exhaustive, at #= none:207 =#")
            $(Expr(:symboliclabel, Symbol("##final#772_1")))
            return_3
        end
    end
    #= none:235 =# Core.@doc "    rm_annotations(x)\n\nRemove type annotation of given expression.\n" function rm_annotations(x)
            x isa Expr || return x
            if x.head == :(::)
                if length(x.args) == 1
                    return gensym("::$(x.args[1])")
                else
                    return x.args[1]
                end
            elseif x.head in [:(=), :kw]
                return rm_annotations(x.args[1])
            else
                return Expr(x.head, map(rm_annotations, x.args)...)
            end
        end
    #= none:255 =# Core.@doc "    alias_gensym(ex)\n\nReplace gensym with `<name>_<id>`.\n\n!!! note\n    Borrowed from [MacroTools](https://github.com/FluxML/MacroTools.jl).\n" alias_gensym(ex) = begin
                alias_gensym!(Dict{Symbol, Symbol}(), Dict{Symbol, Int}(), ex)
            end
    function alias_gensym!(d::Dict{Symbol, Symbol}, count::Dict{Symbol, Int}, ex)
        if is_gensym(ex)
            haskey(d, ex) && return d[ex]
            name = Symbol(gensym_name(ex))
            id = get(count, name, 0) + 1
            d[ex] = Symbol(name, :_, id)
            count[name] = id
            return d[ex]
        end
        ex isa Expr || return ex
        args = map(ex.args) do x
                alias_gensym!(d, count, x)
            end
        return Expr(ex.head, args...)
    end
end


text:
using CodeTracer
using Documenter

makedocs(;
    modules=[CodeTracer],
    authors="Robert Feldt",
    repo="https://github.com/robertfeldt/CodeTracer.jl/blob/{commit}{path}#L{line}",
    sitename="CodeTracer.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://robertfeldt.github.io/CodeTracer.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/robertfeldt/CodeTracer.jl",
)


text:
function cheb2ultra(c_cheb::AbstractVector{T},λ::T,plan::ChebyshevUltrasphericalPlan{BACKWARD,T}) where T<:AbstractFloat
    M,N,nM₀,αN,K = getconstants(plan)
    i₁,i₂,j₁,j₂ = getindices(plan)
    p₁,p₂,rp,c₁,c₂,um,vm,θ,tempsin,tempsin2,tempsinλ,tempsinλm,tempmindices,cnλ,cnmλ,w,anλ,c_cheb2,pr = getplan(plan)

    # Initialize c_cheb2
    @inbounds for i=1:N+1 c_cheb2[i] = c_cheb[i] end
    @inbounds for i=N+2:2N+1 c_cheb2[i] = 0 end

    # Initialize values at Clenshaw-Curtis points
    v_ultra = zero(c_cheb2)

    # Perform transposed DCT-I
    applyTN!(c_cheb2,p₁)

    # Scale values by Clenshaw-Curtis weights
    @inbounds for i=1:2N+1 c_cheb2[i] *= w[i] end

    k=0
    for κ=1:K
        k = κ
        copyto!(cnmλ,cnλ)
        copyto!(tempsinλm,tempsinλ)

        j₁[k] < nM₀ && (k-=1; break)

        @inbounds for i=i₁[k]:i₁[k+1]-1
            forward_recurrence!(pr,j₂[k],θ[i],tempsin2[2N+2-i],tempsin2[i],rp)
            cc2i = c_cheb2[i]
            for j=1:j₂[k] v_ultra[j] += cc2i*pr[j] end
        end

        @inbounds for i=i₂[k+1]+1:i₂[k]
            forward_recurrence!(pr,j₂[k],θ[i],tempsin2[2N+2-i],tempsin2[i],rp)
            cc2i = c_cheb2[i]
            for j=1:j₂[k] v_ultra[j] += cc2i*pr[j] end
        end

        for m=0:M-1
            # Compute u_m(θ) and v_m(θ)
            compute_umvm!(um,vm,λ,tempsinλm,m,θ,i₁[k+1]:i₂[k+1])

            # Multiply point-wise by u_m(θ) and v_m(θ) for valid indices
            init_c₁c₂!(c₁,c₂.parent,um,vm,c_cheb2,i₁[k+1],i₂[k+1])

            # Apply planned DCT-I and DST-I in-place
            applyTN!(c₁,p₁);applyUN!(c₂,p₂)

            # Compute diagonal 2N-scaling multiplied by local coefficients and zero out excess
            @inbounds for j=j₁[k]:j₂[k] v_ultra[j] += cnmλ[j]*(c₁[j]+c₂.parent[j]) end

            # Update C_{n,m}^λ by recurrence in m
            cst = (λ+m)/2/(m+1)*(m+1-λ)
            @inbounds for j=1:2N+1 cnmλ[j] *= cst/(j+λ+m) end

            # Update sin^{λ+m}(θ)
            @inbounds for j=1:2N+1 tempsinλm[j] *= tempsin[j] end
        end
    end

    # Finish off recurrence
    @inbounds for i=i₁[k+1]:i₂[k+1]
        forward_recurrence!(pr,j₂[k+1],θ[i],tempsin2[2N+2-i],tempsin2[i],rp)
        cc2i = c_cheb2[i]
        for j=1:j₂[k+1] v_ultra[j] += cc2i*pr[j] end
    end

    # Scale by orthonormality constants
    ret = zero(c_cheb)
    @inbounds for i=1:N+1 ret[i] = v_ultra[i]/anλ[i] end

    ret
end


text:
using PSPlot

x = 1:10
y = 2x
plot(x, y, label="ll")
xlabel("xx")
ylabel("yy")
legend()
printfig("output.eps")


text:
using ADFPEPS
using ADFPEPS:bulk,HORIZONTAL_RULES
using Random
using Test
using OMEinsum

@testset "ADFPEPS:bulk" begin
    Random.seed!(100)
    T = rand(2,2,4,2,2)
    # @show bulk(T)
    @show
end


text:
using PPT
using Base.Test

const testfile = joinpath(dirname(@__FILE__), "literate_org_tangled_tests.jl")
if isfile(testfile)
    include(testfile)
else
    error("PPT not properly installed. Please run Pkg.build(\"PPT\") then restart Julia.")
end


text:
function JuMP.set_objective_sense(m::InnerBilevelModel, sense::MOI.OptimizationSense)
    JuMP.set_objective_sense(mylevel_model(m), sense)
end
function JuMP.set_objective(m::InnerBilevelModel, sense::MOI.OptimizationSense,
    f::JuMP.AbstractJuMPScalar)
    level_f = replace_variables(f, bilevel_model(m), mylevel_var_list(m), level(m))
    JuMP.set_objective(mylevel_model(m), sense, level_f)
end
JuMP.objective_sense(m::InnerBilevelModel) = JuMP.objective_sense(mylevel_model(m))
function JuMP.objective_function_type(m::InnerBilevelModel)
    tp = JuMP.objective_function_type(mylevel_model(m))
    return bilevel_type(m, tp)
end
bilevel_type(::InnerBilevelModel, ::Type{JuMP.VariableRef}) = BilevelVariableRef
function bilevel_type(::InnerBilevelModel, ::Type{JuMP.GenericAffExpr{C, V}}
) where {C, V}
    JuMP.GenericAffExpr{C, BilevelVariableRef}
end
function bilevel_type(::InnerBilevelModel, ::Type{JuMP.GenericQuadExpr{C, V}}
) where {C, V}
    JuMP.GenericAffExpr{C, BilevelVariableRef}
end
# JuMP.objective_function(m::InnerBilevelModel) = mylevel_obj_function(m)
function JuMP.objective_function(m::InnerBilevelModel)
    f = JuMP.objective_function(mylevel_model(m))
    reverse_replace_variable(f, m)
end
function JuMP.objective_function(m::InnerBilevelModel, FT::Type)
    f = JuMP.objective_function(mylevel_model(m), FT)
    f2 = reverse_replace_variable(f, m)
    f2 isa FT || error("The objective function is not of type $FT, show $(typeof(f2))")
    return f2
end

function JuMP.relative_gap(bm::BilevelModel)::Float64
    _check_solver(bm)
    return MOI.get(bm.solver, MOI.RelativeGap())
end
function JuMP.relative_gap(bm::UpperModel)::Float64
    return JuMP.relative_gap(bm.m)
end
function JuMP.dual_objective_value(bm::BilevelModel; result::Int=1)::Float64
    _check_solver(bm)
    return MOI.get(bm.solver, MOI.DualObjectiveValue(result))
end
function JuMP.objective_bound(bm::BilevelModel)::Float64
    _check_solver(bm)
    return MOI.get(bm.solver, MOI.ObjectiveBound())
end
function JuMP.objective_bound(bm::UpperModel)::Float64
    return JuMP.objective_bound(bm.m)
end

function JuMP.set_objective(m::BilevelModel, sense::MOI.OptimizationSense,
    f::JuMP.AbstractJuMPScalar)
    bilevel_obj_error()
end
JuMP.objective_sense(m::BilevelModel) = JuMP.objective_sense(m.upper)# bilevel_obj_error()
JuMP.objective_function_type(model::BilevelModel) = bilevel_obj_error()
JuMP.objective_function(model::BilevelModel) = bilevel_obj_error()
function JuMP.objective_function(model::BilevelModel, FT::Type)
    bilevel_obj_error()
end

bilevel_obj_error() = error("There is no objective for BilevelModel use Upper(.) and Lower(.)")

function JuMP.objective_value(model::BilevelModel)
    _check_solver(model)
    return MOI.get(model.solver, MOI.ObjectiveValue())
end
function JuMP.objective_value(model::UpperModel)
    return JuMP.objective_value(model.m)
end
function JuMP.objective_value(model::LowerModel)
    return lower_objective_value(model.m)
end

function lower_objective_value(model::BilevelModel; result::Int = 1)
    f = JuMP.objective_function(Lower(model))
    # Evaluate the lower objective expression
    return JuMP.value(f, v -> JuMP.value(v, result = result))
    # return JuMP.value(v -> inner_ref_to_value(Lower(model), v), f)
end

function JuMP.set_objective_coefficient(
    ::UpperModel, variable::BilevelVariableRef, coeff::Real)
    level_var = variable.model.var_upper[variable.idx]
    model = JuMP.owner_model(level_var)
    JuMP.set_objective_coefficient(
        model, level_var, coeff)
end
function JuMP.set_objective_coefficient(
    ::LowerModel, variable::BilevelVariableRef, coeff::Real)
    level_var = variable.model.var_lower[variable.idx]
    model = JuMP.owner_model(level_var)
    JuMP.set_objective_coefficient(
        model, level_var, coeff)
end


text:
T=3
tsname="joDCT"
@testset "$tsname" begin
for t=1:T # start test loop
    m=4^t
    v1=rand(m)
    mv1=rand(m,m)
    v2=rand(m,m)
    vv2=vec(v2)
    mv2=rand(m,m,m)
    vmv2=reshape(mv2,m*m,m)

    A1=joDCT(m)
    A2=joDCT(m,m)
    verbose && println("$tsname ($m[,$m]) - planned")
    @testset "$m [x $m]" begin
        @test isadjoint(joDCT(m))[1]
        @test islinear(joDCT(m))[1]
        @test isadjoint(joDCT(m,m))[1]
        @test islinear(joDCT(m,m))[1]

        @test norm(A1*v1-dct(v1))<joTol
        @test norm(A1\v1-idct(v1))<joTol
        @test norm(adjoint(A1)*v1-idct(v1))<joTol
        @test norm((adjoint(A1)*A1)*v1-v1)<joTol
        @test norm(A1*mv1-dct(mv1,1))<joTol
        @test norm(A1\mv1-idct(mv1,1))<joTol
        @test norm(adjoint(A1)*mv1-idct(mv1,1))<joTol
        @test norm((adjoint(A1)*A1)*mv1-mv1)<joTol

        @test norm(A2*vv2-vec(dct(v2)))<joTol
        @test norm(A2\vv2-vec(idct(v2)))<joTol
        @test norm(adjoint(A2)*vv2-vec(idct(v2)))<joTol
        @test norm((adjoint(A2)*A2)*vv2-vv2)<joTol
        @test norm(vec(A2*vmv2)-vec(dct(mv2,1:2)))<joTol
        @test norm(vec(A2\vmv2)-vec(idct(mv2,1:2)))<joTol
        @test norm(vec(adjoint(A2)*vmv2)-vec(idct(mv2,1:2)))<joTol
        @test norm((adjoint(A2)*A2)*vmv2-vmv2)<joTol
    end

    A1=joDCT(m;planned=false)
    A2=joDCT(m,m;planned=false)
    verbose && println("$tsname ($m[,$m]) - not planned")
    @testset "$m [x $m]" begin
        @test isadjoint(joDCT(m;planned=false))[1]
        @test islinear(joDCT(m;planned=false))[1]
        @test isadjoint(joDCT(m,m;planned=false))[1]
        @test islinear(joDCT(m,m;planned=false))[1]

        @test norm(A1*v1-dct(v1))<joTol
        @test norm(A1\v1-idct(v1))<joTol
        @test norm(adjoint(A1)*v1-idct(v1))<joTol
        @test norm((adjoint(A1)*A1)*v1-v1)<joTol
        @test norm(A1*mv1-dct(mv1,1))<joTol
        @test norm(A1\mv1-idct(mv1,1))<joTol
        @test norm(adjoint(A1)*mv1-idct(mv1,1))<joTol
        @test norm((adjoint(A1)*A1)*mv1-mv1)<joTol

        @test norm(A2*vv2-vec(dct(v2)))<joTol
        @test norm(A2\vv2-vec(idct(v2)))<joTol
        @test norm(adjoint(A2)*vv2-vec(idct(v2)))<joTol
        @test norm((adjoint(A2)*A2)*vv2-vv2)<joTol
        @test norm(vec(A2*vmv2)-vec(dct(mv2,1:2)))<joTol
        @test norm(vec(A2\vmv2)-vec(idct(mv2,1:2)))<joTol
        @test norm(vec(adjoint(A2)*vmv2)-vec(idct(mv2,1:2)))<joTol
        @test norm((adjoint(A2)*A2)*vmv2-vmv2)<joTol
    end

end # end test loop
end


text:
# ## Getting started

using MLJ, RDatasets, PrettyPrinting, Random
MLJ.color_off() # hide

# We start by generating a 2D cloud of points

Random.seed!(3203)
X = randn(20, 2)
y = vcat(-ones(10), ones(10))

# which we can visualise

using PyPlot
figure(figsize=(8,6))

ym1 = y .== -1
ym2 = .!ym1
plot(X[ym1, 1], X[ym1, 2], ls="none", marker="o")
plot(X[ym2, 1], X[ym2, 2], ls="none", marker="x")

savefig(joinpath(@OUTPUT, "ISL-lab-9-g1.svg")) # hide

# \figalt{Toy points}{ISL-lab-9-g1.svg}

# let's wrap the data as a table:

X = MLJ.table(X)
y = categorical(y);

# and fit a SVM classifier

@load SVC pkg=LIBSVM

svc_mdl = SVC()
svc = machine(svc_mdl, X, y)

fit!(svc);

# As usual we can check how it performs

ypred = predict(svc, X)
misclassification_rate(ypred, y)

# Not bad.

# ### Basic tuning

# As usual we could tune the model, for instance the penalty encoding the tradeoff between margin width and misclassification:

rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = predict(mtm, X)
misclassification_rate(ypred, y)

# You could also change the kernel etc.
PyPlot.close_figs() # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
push!(LOAD_PATH, string(pwd(),"/src"))

using ModifiedStokesSolver
using FastGaussQuadrature
using Base.Test

panelorder1 = 32
panelorder2 = 160

glpoints1, glweights1 = gausslegendre(panelorder1)
glpoints2, glweights2 = gausslegendre(panelorder2)

@testset "Upsampling" begin
    U = ModifiedStokesSolver.upsampling_matrix(panelorder1, panelorder2)
    f(x) = @. 1/(x^2+2)
    fu = U * f(glpoints1)
    err = fu - f(glpoints2)
    relerr_inf = norm(err, Inf) / norm(fu, Inf)
    @test relerr_inf < 2e-14
end


text:
using Documenter, Jot, Pkg

makedocs(
  sitename = "Jot Documentation",
  format = Documenter.HTML(
      prettyurls = get(ENV, "CI", nothing) == "true"
  ),
  pages = [
    "Home" => "index.md",
    "Manual" => [
      "Guide" => "Guide.md",
      "Examples" => "Examples.md",
      "Managing Lambdas" => "Managing_lambdas.md",
      "Labels" => "Labels.md",
    ],
    "API" => [
      "Functions" => "Functions.md",
      "Types" => "Types.md",
    ]
  ],
)

deploydocs(
  repo = "github.com/harris-chris/Jot.jl.git",
  devbranch = "main",
)


text:
using Base.Test
using DataFrames

let
	# @assert isequal(xtab([1, 2, 2, 2, 3]), xtab([1, 2, 3],[1, 3, 1]))
	@assert isequal(xtabs([1, 2, 2, 2, 3]), [2 => 3, 3 => 1, 1 => 1])
end


text:
using BenchmarkTools, AssociatedLegendrePolynomials, Markdown

dotune  = "--tune" in ARGS
saveres = "--save" in ARGS

const PARAMS_PATH = joinpath(dirname(@__FILE__), "params.json")
const SUITE   = BenchmarkGroup()
const ASSERTS = BenchmarkGroup()

const LMAX = 700
const MMAX = 350
const NORMS = (LegendreUnitNorm(),
               LegendreSphereNorm(),
               LegendreUnitCoeff{Float64}(LMAX),
               LegendreSphereCoeff{Float64}(LMAX))

const xs = 0.5
const x0 = fill(0.5)
const x1 = collect(range(-1, 1, length=100))
const x2 = collect(reshape(x1, 10, 10))

##############################################################################
# Benchmarks which should be compared revision-to-revision
##############################################################################

@info "Loading benchmarks..."
# Run suite of benchmarks on each normalization style
for norm in NORMS
    nname = string(typeof(norm))
    gnorm = addgroup!(SUITE, nname)

    # Organize by output shape --- single value at (lmax,mmax), vector of values
    # ({(l,mmax) | l <= lmax}), and matrix of values ({(l,m) | l <= lmax and m <= l}).
    for (outdim, osz) in zip(("dim0", "dim1", "dim2"), ((), (LMAX,), (LMAX,MMAX)))
        gname = "out" * outdim
        # Then use true scalar and 0 through 2-dimensional inputs
        for (indim, x) in zip(("scalar", "dim0", "dim1", "dim2"), (xs, x0, x1, x2))
            bname = "in" * indim
            Λ = zeros(size(x)..., (osz .+ 1)...)
            gnorm[(gname, bname)] = @benchmarkable legendre!($norm, $Λ, $LMAX, $MMAX, $x)
        end
    end
end

##############################################################################
# (Optionally tune and) Run the benchmarks
##############################################################################

if dotune
    @info "Tuning benchmarks..."
    tune!(SUITE, verbose=true)
    BenchmarkTools.save(PARAMS_PATH, params(SUITE))
end
if isfile(PARAMS_PATH)
    @info "Loaded benchmarking parameters." PARAMS_PATH
    loadparams!(SUITE, BenchmarkTools.load(PARAMS_PATH)[1], :evals)
end

@info "Running benchmarks..."
res = run(SUITE, verbose=true)
# Truncate down to the minimums immediately. Make the saved file smaller.
res = minimum(res)

##############################################################################
# Assertions on performance characteristics which sould be true for any
# given revision. Uses results of the above benchmarks without comparison to
# another revision.
##############################################################################

@info "Making intra-benchmark comparisons..."
for norm in NORMS
    nname = string(typeof(norm))
    gn = addgroup!(ASSERTS, nname)
    resn = res[nname]

    # Make sure scalar arguments perform better than an equivalent 0-dimensional input
    for outdim in ("dim0", "dim1", "dim2")
        gname = "out" * outdim
        gn[(gname, "inscalar vs indim0")] =
            judge(resn[(gname, "inscalar")], resn[(gname, "indim0")])
    end
end

if saveres
    branch = strip(read(`git symbolic-ref --short HEAD`, String))
    if branch == "master"
        branch = strip(read(`git describe --tags`, String))
    end
    # TODO: BenchmarkTools.load() doesn't work with the judged results, so save the
    #       benchmarks and assertions separately.
    benchfile   = "bench-" * branch * ".json"
    assertsfile = "bench-" * branch * "-asserts.json"
    @info "Saving benchmarks to `bench-$branch.json`" pwd=pwd()
    BenchmarkTools.save(benchfile, res)
    BenchmarkTools.save(assertsfile, ASSERTS)
end

##############################################################################
# Summary printing
##############################################################################

function printtable(group)
    for k in sort!(collect(keys(group)))
        g = group[k]
        md = Markdown.MD()
        push!(md.content, Markdown.parse("## " * k))
        tab = Markdown.Table(Any[], Any[:l, :l])
        push!(tab.rows, Any["Name", "Trial"])
        for l in sort!(leaves(g))
            n, t = l
            push!(tab.rows, Any[n[1], repr(t)])
        end
        push!(md.content, tab)
        show(md)
        print("\n\n")
    end
    print("\n")
end

println("---")
println("# BENCHMARKS")
printtable(res)
println()

println("---")
println("# ASSERTIONS")
printtable(ASSERTS)
println()


text:
using Krylov
# convex
# ACG = readdlm("NewResults/convpbsCGlin.txt")
# ACR = readdlm("NewResults/convpbsCRlin.txt")

# nonconvex
ACG = readdlm("NewResults/nonconvpbsCGlin.txt")
ACR = readdlm("NewResults/nonconvpbsCRlin.txt")

Problems = ACG[2:end,1]

function perfprofileOpt()
    stats = Dict{String, Any}("CR" => ["#f" "#g" "#Hv"], "CG" => ["#f" "#g" "#Hv"])
    k = 1
    for p in Problems
        k += 1
        statCR = ACR[k,12] ? ACR[k,8:10] : -ACR[k,8:10]
        stats["CR"] = vcat(values(stats["CR"]), statCR')
        statCG = ACG[k,12] ? ACG[k,8:10] : -ACG[k,8:10]
        stats["CG"] = vcat(values(stats["CG"]), statCG')
    end
    return stats
end

stats = perfprofileOpt()
tf = font(17) # titlefont
f = font(17)

pb_type = "nonconvex problems"
algo_used = "CR and CG"

p = performance_profile(hcat([p[2:end, 1] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("#f (", pb_type, ")"),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues, w=2, style=:auto) #color=:black) # Profile for #f
# Plots.xlabel!("Within this factor of the best (log₂ scale)")
savefig(p, string("profil_f_", pb_type, ".pdf"))

p = performance_profile(hcat([p[2:end, 2] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("#g (", pb_type, ")"),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues, w=2, style=:auto) # Profile for #g
# Plots.xlabel!("Within this factor of the best (log₂ scale)")
savefig(p, string("profil_g_", pb_type, ".pdf"))

p = performance_profile(hcat([p[2:end, 3] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("#Hv (", pb_type, ")"),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues, w=2, style=:auto) # Profile for #Hv
# Plots.xlabel!("Within this factor of the best (log₂ scale)")
Plots.savefig(p, string("profil_Hv_", pb_type, ".pdf"))

p = performance_profile(hcat([p[2:end, 1]+p[2:end, 2]+p[2:end, 3] for p in values(stats)]...),
                            collect(String, [string(s) for s in keys(stats)]),
                            title=string("#f + #g + #Hv (", pb_type, ")"),
                            titlefont = tf, legendfont = f, guidefont = f,
                            legend=:bottomright, palette=:blues, w=2, style=:auto) # Profile for #f + #g + #Hv
# Plots.xlabel!("Within this factor of the best (log₂ scale)")
savefig(p, string("profil_f+g+Hv_", pb_type, ".pdf"))


text:
#define module WvDerivativeTi containing relevant functions
#to calculate the derivative of the RBM wavefunction
#for a given state and for translation invariance symmetry
module WvDerivativeTi

#Return a matrix that is used to sort the N independent weights in the matrix W_ij
function cod(y,L_x,L_y)

  #define vectors for useful storing of indexes
  A = Array{Array{Int64}}(undef,L_x)
  B = Array{Array{Int64}}(undef,L_y)
  C = Array{Array{Int64}}(undef,L_x*L_y)

  for i=1:L_y
    B[i] = collect(1+y*L_x*L_y+(i-1)*L_x:i*L_x+y*L_x*L_y)
  end

  for i=1:L_x
      A[i] = Int64[]
  end

  q = 1
  while q <= L_x
      for i=1:L_y
          append!(A[q],circshift(B[i],-(q-1)))
      end

      q += 1

  end

  for i=1:L_x*L_y

      C[i] = Int64[]
  end

  k = 1
  while k < L_x*L_y
      for h=0:L_y-1
          for s=1:L_x
              append!(C[k],circshift(A[s],-L_x*h))
              k = k + 1
          end
      end
  end

  Z =   Array{Int64}(undef,L_x*L_y,L_x*L_y)

  for i=1:L_x*L_y
      for j=1:L_x*L_y
          Z[i,j] = C[i][j]
      end
  end

  return Z

end

#Define useful objects for the implementation of symmetry
#and efficient calculation of observables for the translation symmetric RBM.
function List(Cod::Array{Int64},ii::Int64,nspins::Int64)

    K = Array{Int64}(undef,nspins)

    for t=1:nspins
        K[t]= findfirst(isequal(ii),Cod[:,t])
    end

    return K

end

#Return a matrix used to efficiently calculate the derivatives
#of the wavefunction for the translation symmetric RBM
function transl_inv_tool(nspins::Int64, alpha::Int64,L_x::Int64,L_y::Int64)

  #useful matrices for filling in the filter matrix
  Cod = Array{Array{Int64}}(undef,alpha)
  ss = Array{Array{Int64}}(undef,alpha)

  #Matrix storing the filter index matrix for transaltion invariance
  filter_ind = Array{Int64}(undef,nspins,nspins*alpha)

  for h=1:alpha
  Cod[h] = cod(h-1,L_x,L_y)
  ss[h] = Cod[h][:,1]
  end

  for h=1:alpha
      for i=1:nspins
          filter_ind[:,i+(h-1)*nspins] = List(Cod[h],ss[h][i],nspins)
      end
  end

  return filter_ind

end

#Return a vector used to efficiently calculate the derivatives
#of the wavefunction for the translation symmetric RBM
function ti_tool(nspins::Int64, L_x::Int64,L_y::Int64,ind_h::Int64,ind_i::Int64)

  Cod = cod(ind_h - 1,L_x,L_y)
  ss = Cod[:,1]

  return List(Cod,ss[ind_i],nspins)

end

#Calculate the derivative of the wavefunction respect to the variational parameters for a given state
function Var_derivatives(state::Array{Int64,1},ENWfDer::Array{Complex{Float64},2},Lt::Array{Complex{Float64},1},FilterInd::Array{Int64,2},nspins::Int64,alpha::Int64, ind::Int64)

    #Fill ENWfDer with the derivative of the wavefunction w.r.t. b and W
    k = 1
    for h=1:alpha
        r_ = (h-1)*nspins; @views t_ = tanh.(Lt[(r_+1):(h*nspins)])
        @views    ENWfDer[ind,1 + h] = sum(t_)
        for i=1:nspins
                @views    lt = state[FilterInd[:,i+r_]]'*t_
                    ENWfDer[ind,1 + alpha + k] = lt
                    k += 1
        end
    end

end

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
export bezsurf

"""
	bezsurf(npts, mptsg, b[], udpts, wdpts)

Calculate a Bezier surface over the `npts` by `mpts` control net points in `b`.

The method evaluates `udpts` by `wdpts` points of the curve and gives back the tuple `(P, EV, FV)`.
Here:
 - `P` is the vector of points memorized by columns;
 - `EV` is the `1`-dimensional cellular complex of the edges;
 - `FV` is the `2`-dimensional cellular complex of the faces.

This structure as it is formed could be directly imported in Plasm module by typing:
```julia
julia> P, EV, FV = bezier(npts, b[], dpts);
julia> Plasm.view(Plasm.mkpol(P, FV))
```

---

# Arguments

- `npts::Int64`: the number of the control net rows in `u` direction.
- `mpts::Int64`: the number of the control net rows in `w` direction;
- `b::Array{Float64, 2}`: 2-dimensional Array containing the control polygon vertices in the x, y, z coordinates for columns.
- `udpts::Int64`: number of data net rows in `u` direction.
- `wdpts::Int64`: number of data net rows in `w` direction.

---

# Examples

```jldoctest
julia> b = [
		  0.0 0.0 0.0 1.0 1.0 1.0 2.0 2.0 2.0;
		  0.0 2.0 4.0 0.0 2.0 4.0 0.0 2.0 4.0;
		  2.0 0.0 2.0 0.0 1.0 0.0 2.0 0.0 2.0
		]
3×9 Array{Float64,2}:
 0.0  0.0  0.0  1.0  1.0  1.0  2.0  2.0  2.0
 0.0  2.0  4.0  0.0  2.0  4.0  0.0  2.0  4.0
 2.0  0.0  2.0  0.0  1.0  0.0  2.0  0.0  2.0

julia> bezsurf(3, 3, b, 5, 5)[1]
3×25 Array{Float64,2}:
 0.0  0.5   1.0  1.5   2.0  0.0   0.5       1.0     1.5       2.0   0.0  0.5     1.0   1.5     2.0  0.0   0.5       1.0     1.5       2.0   0.0  0.5   1.0  1.5   2.0
 0.0  0.0   0.0  0.0   0.0  1.0   1.0       1.0     1.0       1.0   2.0  2.0     2.0   2.0     2.0  3.0   3.0       3.0     3.0       3.0   4.0  4.0   4.0  4.0   4.0
 2.0  1.25  1.0  1.25  2.0  1.25  0.921875  0.8125  0.921875  1.25  1.0  0.8125  0.75  0.8125  1.0  1.25  0.921875  0.8125  0.921875  1.25  2.0  1.25  1.0  1.25  2.0

julia> bezsurf(3, 3, b, 5, 5)[2]
44-element Array{Array{Int64,1},1}:
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 5]
 [5, 6]
 [6, 7]
 [7, 8]
 [8, 9]
 [9, 10]
 [10, 11]
 ⋮
 [12, 17]
 [13, 18]
 [14, 19]
 [15, 20]
 [16, 21]
 [17, 22]
 [18, 23]
 [19, 24]
 [20, 25]

julia> bezsurf(3, 3, b, 5, 5)[3]
16-element Array{Array{Int64,1},1}:
 [1, 2, 6, 7]
 [2, 3, 7, 8]
 [3, 4, 8, 9]
 [4, 5, 9, 10]
 [6, 7, 11, 12]
 [7, 8, 12, 13]
 [8, 9, 13, 14]
 [9, 10, 14, 15]
 [11, 12, 16, 17]
 [12, 13, 17, 18]
 [13, 14, 18, 19]
 [14, 15, 19, 20]
 [16, 17, 21, 22]
 [17, 18, 22, 23]
 [18, 19, 23, 24]
 [19, 20, 24, 25]
```

---

_By Elia Onofri_
"""

function bezsurf(npts::Int64, mpts::Int64, b::Array{Float64,2}, udpts::Int64, wdpts::Int64)::Tuple{Array{Float64,2}, Array{Array{Int64,1},1}, Array{Array{Int64,1},1}}

	@assert npts*mpts == length(b[1,:]) ("ERROR: there are not npts * mpts = $(npts) * $(mpts) = $(npts * mpts) points of the net in b[]")

	ustep::Float64 = 1 / (udpts-1)
	wstep::Float64 = 1 / (wdpts - 1)
	n::Int64 = npts - 1
	m::Int64 = mpts - 1

	#=== x, y, z components of B ===#
	x::Array{Float64,2} = zeros(npts, mpts)
	y::Array{Float64,2} = zeros(npts, mpts)
	z::Array{Float64,2} = zeros(npts, mpts)

	for i = 1 : npts
		for j = 1 : mpts
			ij = (i - 1) * mpts + j
			x[i, j] = b[1, ij]
			y[i, j] = b[2, ij]
			z[i, j] = b[3, ij]
		end
	end

	#=== SYNC BLOCK BEGIN ===#
	#=== all the following @async are totally indipendent blocks ===#
	@sync begin

		#=== N = C * D ===#
		C::Array{Float64,2} = zeros(2, 2)
		D::Array{Float64,2} = zeros(2, 2)

		@async begin
			C = zeros(npts, npts)
			D = zeros(npts, npts)
			for i = 0 : n
				for j = 0 : (n - i)
					if (n-j-i)%2==0
						C[i+1, j+1] = binomial(n-j, n-j-i)
					else
						C[i+1, j+1] = 0.0-binomial(n-j, n-j-i)
					end
				end
			end

			for i = 1 : npts
				D[i, i] = abs(C[i, 1])
			end
		end

		#=== M = E * F ===#
		E::Array{Float64,2} = zeros(2, 2)
		F::Array{Float64,2} = zeros(2, 2)

		@async begin
			E = zeros(mpts, mpts)
			F = zeros(mpts, mpts)
			for i = 0 : m
				for j = 0 : (m - i)
					if (m - j - i) % 2 == 0
						E[i+1, j+1] = binomial(m-j, m-j-i)
					else
						E[i+1, j+1] = 0.0 - binomial(m-j, m-j-i)
					end
				end
			end

			for i = 1 : mpts
				F[i, i] = abs(E[i, 1])
			end
		end

		#=== U = [u^n ... u 1] ===#
		U::Array{Float64,2} = zeros(2, 2)

		@async begin
			U = ones(udpts, npts)
			u::Float64 = 0.0
			for i = 1 : udpts
				for j = 0 : n-1  #Per j=0 il valore è 1 ed è già impostato
					U[i, n-j] = U[i, npts-j] * u
				end
				u = u + ustep
			end
		end

		#=== W = [w^m ... w 1] ===#
		W::Array{Float64,2} = zeros(2, 2)

		@async begin
			W = ones(wdpts, mpts)
			w::Float64 = 0.0
			for i = 1 : wdpts
				for j = 0 : m-1  #Per j=0 il valore è 1 ed è già impostato
					W[i, m-j] = W[i, mpts-j] * w
				end
				w = w + wstep
			end
		end

		#=== Edge Vector ===#
		EV::Array{Array{Int64,1},1} = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

		@async begin
			for i = 1 : ((udpts * wdpts) - 1)
				app = Array{Int64}(2)
				app[1] = i
				app[2] = i+1
				push!(EV, app)
			end

			for i = 1 : (wdpts * (udpts - 1))
				app = Array{Int64}(2)
				app[1] = i
				app[2] = i + wdpts
				push!(EV, app)
			end
		end

		#=== Faces Vector ===#
		FV::Array{Array{Int64,1},1} = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

		@async begin
			for j = 0 : (udpts - 2)
				for i = 1 : wdpts - 1
					app = Array{Int64}(4)
					app[1] = j * wdpts + i
					app[2] = j * wdpts + i + 1
					app[3] = (j + 1) * wdpts + i
					app[4] = (j + 1) * wdpts + i + 1
					push!(FV, app)
				end
			end
		end

	end
	#=== SYNC BLOCK END ===#

	#=== Points Vector ===#
	UCD::Array{Float64,2} = U * C * D
	FEW::Array{Float64,2} = transpose(W * E * F)
	xP::Array{Float64,2} = transpose(UCD * x * FEW)
	yP::Array{Float64,2} = transpose(UCD * y * FEW)
	zP::Array{Float64,2} = transpose(UCD * z * FEW)

	P::Array{Float64,2} = zeros(3, udpts * wdpts)

	for i = 1 : udpts
		for j = 1 : wdpts
			ij::Int64 = (i-1) * udpts + j
			P[1, ij] = xP[i, j]
			P[2, ij] = yP[i, j]
			P[3, ij] = zP[i, j]
		end
	end

	return (P, EV, FV)
end


text:
# This file was generated, do not modify it. # hide
X3 = select(X2, [:Lag1, :Lag2])
clf = machine(LogisticClassifier(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
mean(ŷ .== y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
module CUSOLVER

import CUDAdrv: CUDAdrv, CuContext, CuStream_t, CuPtr, PtrOrCuPtr, CU_NULL
import CUDAapi

using ..CuArrays
using ..CuArrays: libcusolver, active_context, _getindex, unsafe_free!

using LinearAlgebra
using SparseArrays

import Base.one
import Base.zero
import CuArrays.CUSPARSE.CuSparseMatrixCSR
import CuArrays.CUSPARSE.CuSparseMatrixCSC
import CuArrays.CUSPARSE.cusparseMatDescr_t

include("libcusolver_types.jl")
include("error.jl")

const _dense_handles = Dict{CuContext,cusolverDnHandle_t}()
const _dense_handle = Ref{cusolverDnHandle_t}(C_NULL)
const _sparse_handles = Dict{CuContext,cusolverSpHandle_t}()
const _sparse_handle = Ref{cusolverSpHandle_t}(C_NULL)

function dense_handle()
    if _dense_handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _dense_handle[] = get!(_dense_handles, active_context[]) do
            context = active_context[]
            handle = cusolverDnCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusolverDnDestroy(handle))
            handle
        end
    end
    return _dense_handle[]
end

function sparse_handle()
    if _sparse_handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _sparse_handle[] = get!(_sparse_handles, active_context[]) do
            context = active_context[]
            handle = cusolverSpCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusolverSpDestroy(handle))
            handle
        end
    end
    return _sparse_handle[]
end

include("libcusolver.jl")
include("sparse.jl")
include("dense.jl")
include("highlevel.jl")

version() = VersionNumber(cusolverGetProperty(CUDAapi.MAJOR_VERSION),
                          cusolverGetProperty(CUDAapi.MINOR_VERSION),
                          cusolverGetProperty(CUDAapi.PATCH_LEVEL))

end


text:
function BenchPerfConfigSweeps.load(resultdir)
    sweep = BenchmarkConfigSweeps.load(resultdir)
    perfs = Vector{Union{GroupPerf,Nothing}}(undef, length(sweep.results))
    fill!(perfs, nothing)
    for (i, r) in enumerate(sweep.results)
        if r !== nothing
            perfpath = BenchPerf.perfpath(resultpath(resultdir, i))
            perfs[i] = BenchPerf.loadperf(perfpath)
        end
    end
    replacetrials!(sweep.results, perfs)
    return sweep
    # return BenchPerfSweepResult(sweep, perfs)
end

# struct BenchPerfSweepResult
#     sweep::SweepResult
#     perfs::Vector{Union{GroupPerf,Nothing}}
# end

function replacetrials!(results::AbstractVector, perfs::AbstractVector)
    @assert axes(results) == axes(perfs)
    for (r, p) in zip(results, perfs)
        r === nothing && continue
        @assert p !== nothing
        replacetrials!(r, p)
    end
    return results
end

function replacetrials!(bench::BenchmarkGroup, perf::GroupPerf)
    @assert issetequal(keys(bench), keys(perf))
    for (k, bv) in bench
        pv = perf[k]
        bench[k] = replacetrials!(bv, pv)
    end
    return bench
end

replacetrials!(bench::BenchmarkTools.Trial, perf::TrialPerf) = TrialResult(bench, perf)


text:
using JuMP, FlexibilityAnalysis, Gurobi

# Set the covariance matrix for the uncertain parameters
means = [4.0; 5.0]
bea = 1
covar = [2. bea;
         bea 3.]

# Setup the model
m = FlexibilityModel(solver = GurobiSolver(OutputFlag = 0))

# Define variables
@randomvariable(m, x, mean = means[1])
@randomvariable(m, y, mean = means[2])

# Define the constraints
@constraint(m, x + y - 14 <= 0.0)
@constraint(m, x - 2y - 2 <= 0.0)
@constraint(m, -x <= 0.0)
@constraint(m, -y <= 0.0)

# Define the uncertainty set
setuncertaintyset(m, :Ellipsoid, covar)

# Solve
status = solve(m, active_constr = true)
if status == :Optimal
    # Retrieve optimized data
    flexibility_index = getflexibilityindex(m)
    x_opt = getvalue(x)
    y_opt = getvalue(y)
    actives = getactiveconstraints(m)

    # Get the confidence level
    conf_lvl = getconfidencelevel(m)

    # Rank the inequality constraints
    rank_data = rankinequalities(m, active_constr = true)
end

# Estimate the value of SF
SF = findstochasticflexibility(m, use_flexibility_index = true)