text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    _replace_moi_variables(
        nodes::Vector{Nonlinear.Node},
        moi_index_to_consecutive_index::Dict{MOI.VariableIndex,Int},
    )

Return a new `Vector{Nonlinear.Node}` where all occurances of
`NODE_MOI_VARIABLE` are replaced by `NODE_VARIABLE` that is 1-indexed and
ordered.
"""
function _replace_moi_variables(
    nodes::Vector{Nonlinear.Node},
    moi_index_to_consecutive_index::Dict{MOI.VariableIndex,Int},
)
    new_nodes = Vector{Nonlinear.Node}(undef, length(nodes))
    for (i, node) in enumerate(nodes)
        if node.type == Nonlinear.NODE_MOI_VARIABLE
            new_nodes[i] = Nonlinear.Node(
                Nonlinear.NODE_VARIABLE,
                moi_index_to_consecutive_index[MOI.VariableIndex(node.index)],
                node.parent,
            )
        else
            new_nodes[i] = node
        end
    end
    return new_nodes
end

@enum(Linearity, CONSTANT, LINEAR, PIECEWISE_LINEAR, NONLINEAR)

"""
    _classify_linearity(
        nodes::Vector{Nonlinear.Node},
        adj::SparseArrays.SparseMatrixCSC,
        subexpression_linearity::Vector{Linearity},
    )

Classify the nodes in a tree as constant, linear, or nonlinear with respect to
the input.
"""
function _classify_linearity(
    nodes::Vector{Nonlinear.Node},
    adj::SparseArrays.SparseMatrixCSC,
    subexpression_linearity::Vector{Linearity},
)
    linearity = Array{Linearity}(undef, length(nodes))
    children_arr = SparseArrays.rowvals(adj)
    for k in length(nodes):-1:1
        node = nodes[k]
        if node.type == Nonlinear.NODE_VARIABLE
            linearity[k] = LINEAR
            continue
        elseif node.type == Nonlinear.NODE_VALUE
            linearity[k] = CONSTANT
            continue
        elseif node.type == Nonlinear.NODE_PARAMETER
            linearity[k] = CONSTANT
            continue
        elseif node.type == Nonlinear.NODE_SUBEXPRESSION
            linearity[k] = subexpression_linearity[node.index]
            continue
        end
        children_idx = SparseArrays.nzrange(adj, k)
        num_constant_children, any_nonlinear = 0, false
        for r in children_idx
            if linearity[children_arr[r]] == NONLINEAR
                any_nonlinear = true
                break
            elseif linearity[children_arr[r]] == CONSTANT
                num_constant_children += 1
            end
        end
        if any_nonlinear
            # If any children are nonlinear, then we're nonlinear...
            linearity[k] = NONLINEAR
            # ...except in the case of ifelse. If the operands are linear then
            # we're piecewise linear.
            op = get(
                Nonlinear.DEFAULT_MULTIVARIATE_OPERATORS,
                node.index,
                nothing,
            )
            if (
                node.type == Nonlinear.NODE_CALL_MULTIVARIATE &&
                op == :ifelse &&
                linearity[children_arr[children_idx[2]]] == LINEAR &&
                linearity[children_arr[children_idx[3]]] == LINEAR
            )
                linearity[k] = PIECEWISE_LINEAR
            end
            continue
        elseif num_constant_children == length(children_idx)
            # If all children are constant, then we're constant.
            linearity[k] = CONSTANT
            continue
        end
        # By this point, some children are constant and some are linear, so if
        # the operator is nonlinear, then we're nonlinear.
        if node.type == Nonlinear.NODE_CALL_UNIVARIATE
            op =
                get(Nonlinear.DEFAULT_UNIVARIATE_OPERATORS, node.index, nothing)
            if op == :+ || op == :-
                linearity[k] = LINEAR
            else
                linearity[k] = NONLINEAR
            end
        elseif node.type == Nonlinear.NODE_CALL_MULTIVARIATE
            op = get(
                Nonlinear.DEFAULT_MULTIVARIATE_OPERATORS,
                node.index,
                nothing,
            )
            if op == :+
                linearity[k] = LINEAR
            elseif op == :-
                linearity[k] = LINEAR
            elseif op == :*
                # Multiplication is linear if there is one non-constant term.
                one_op = num_constant_children == length(children_idx) - 1
                linearity[k] = one_op ? LINEAR : NONLINEAR
            elseif op == :^
                linearity[k] = NONLINEAR
            elseif op == :/
                if linearity[children_arr[children_idx[2]]] == CONSTANT
                    # If the denominator is constant, we're linear.
                    linearity[k] = LINEAR
                else
                    linearity[k] = NONLINEAR
                end
            elseif op == :ifelse
                linearity[k] = NONLINEAR
            else  # User-defined functions
                linearity[k] = NONLINEAR
            end
        elseif node.type == Nonlinear.NODE_LOGIC
            linearity[k] = NONLINEAR
        else
            @assert node.type == Nonlinear.NODE_COMPARISON
            linearity[k] = NONLINEAR
        end
    end
    return linearity
end

"""
    _compute_gradient_sparsity!(
        indices::Coloring.IndexedSet,
        nodes::Vector{Nonlinear.Node},
    )

Compute the sparsity pattern of the gradient of an expression (i.e., a list of
which variable indices are present).
"""
function _compute_gradient_sparsity!(
    indices::Coloring.IndexedSet,
    nodes::Vector{Nonlinear.Node},
)
    for node in nodes
        if node.type == Nonlinear.NODE_VARIABLE
            push!(indices, node.index)
        elseif node.type == Nonlinear.NODE_MOI_VARIABLE
            error(
                "Internal error: Invalid to compute sparsity if Nonlinear.NODE_MOI_VARIABLE " *
                "nodes are present.",
            )
        end
    end
    return
end

"""
    _compute_hessian_sparsity(
        nodes::Vector{Nonlinear.Node},
        adj,
        input_linearity::Vector{Linearity},
        indexedset::Coloring.IndexedSet,
        subexpression_edgelist::Vector{Set{Tuple{Int,Int}}},
        subexpression_variables::Vector{Vector{Int}},
    )

Compute the sparsity pattern the Hessian of an expression.

 * `input_linearity` is the linearity with respect to the input, computed by
   `_classify_linearity`
 * `subexpression_edgelist` is the edge_list of each subexpression
 * `subexpression_variables` is the list of all variables which appear in a
   subexpression (including recursively).

Idea: consider the (non)linearity of a node *with respect to the output*. The
children of any node which is nonlinear with respect to the output should have
nonlinear interactions, hence nonzeros in the hessian. This is not true in
general, but holds for everything we consider.

A counter example is `f(x, y, z) = x + y * z`, but we don't have any functions
like that. By "nonlinear with respect to the output", we mean that the output
depends nonlinearly on the value of the node, regardless of how the node itself
depends on the input.
"""
function _compute_hessian_sparsity(
    nodes::Vector{Nonlinear.Node},
    adj,
    input_linearity::Vector{Linearity},
    indexedset::Coloring.IndexedSet,
    subexpression_edgelist::Vector{Set{Tuple{Int,Int}}},
    subexpression_variables::Vector{Vector{Int}},
)
    # So start at the root of the tree and classify the linearity wrt the output.
    # For each nonlinear node, do a mini DFS and collect the list of children.
    # Add a nonlinear interaction between all children of a nonlinear node.
    edge_list = Set{Tuple{Int,Int}}()
    nonlinear_wrt_output = fill(false, length(nodes))
    children_arr = SparseArrays.rowvals(adj)
    stack = Int[]
    stack_ignore = Bool[]
    nonlinear_group = indexedset
    if length(nodes) == 1 && nodes[1].type == Nonlinear.NODE_SUBEXPRESSION
        # Subexpression comes in linearly, so append edge_list
        for ij in subexpression_edgelist[nodes[1].index]
            push!(edge_list, ij)
        end
    end
    for k in 2:length(nodes)
        nod = nodes[k]
        if nod.type == Nonlinear.NODE_MOI_VARIABLE
            error(
                "Internal error: Invalid to compute sparsity if Nonlinear.NODE_MOI_VARIABLE " *
                "nodes are present.",
            )
        end
        if nonlinear_wrt_output[k]
            continue # already seen this node one way or another
        elseif input_linearity[k] == CONSTANT
            continue # definitely not nonlinear
        end
        @assert !nonlinear_wrt_output[nod.parent]
        # check if the parent depends nonlinearly on the value of this node
        par = nodes[nod.parent]
        if par.type == Nonlinear.NODE_CALL_UNIVARIATE
            op = get(Nonlinear.DEFAULT_UNIVARIATE_OPERATORS, par.index, nothing)
            if op === nothing || (op != :+ && op != :-)
                nonlinear_wrt_output[k] = true
            end
        elseif par.type == Nonlinear.NODE_CALL_MULTIVARIATE
            op = get(
                Nonlinear.DEFAULT_MULTIVARIATE_OPERATORS,
                par.index,
                nothing,
            )
            if op === nothing
                nonlinear_wrt_output[k] = true
            elseif op in (:+, :-, :ifelse)
                # pass
            elseif op == :*
                # check if all siblings are constant
                sibling_idx = SparseArrays.nzrange(adj, nod.parent)
                if !all(
                    i ->
                        input_linearity[children_arr[i]] == CONSTANT ||
                            children_arr[i] == k,
                    sibling_idx,
                )
                    # at least one sibling isn't constant
                    nonlinear_wrt_output[k] = true
                end
            elseif op == :/
                # check if denominator is nonconstant
                sibling_idx = SparseArrays.nzrange(adj, nod.parent)
                if input_linearity[children_arr[last(sibling_idx)]] != CONSTANT
                    nonlinear_wrt_output[k] = true
                end
            else
                nonlinear_wrt_output[k] = true
            end
        end
        if nod.type == Nonlinear.NODE_SUBEXPRESSION && !nonlinear_wrt_output[k]
            # subexpression comes in linearly, so append edge_list
            for ij in subexpression_edgelist[nod.index]
                push!(edge_list, ij)
            end
        end
        if !nonlinear_wrt_output[k]
            continue
        end
        # do a DFS from here, including all children
        @assert isempty(stack)
        @assert isempty(stack_ignore)
        sibling_idx = SparseArrays.nzrange(adj, nod.parent)
        for sidx in sibling_idx
            push!(stack, children_arr[sidx])
            push!(stack_ignore, false)
        end
        empty!(nonlinear_group)
        while length(stack) > 0
            r = pop!(stack)
            should_ignore = pop!(stack_ignore)
            nonlinear_wrt_output[r] = true
            if nodes[r].type == Nonlinear.NODE_LOGIC ||
               nodes[r].type == Nonlinear.NODE_COMPARISON
                # don't count the nonlinear interactions inside
                # logical conditions or comparisons
                should_ignore = true
            end
            children_idx = SparseArrays.nzrange(adj, r)
            for cidx in children_idx
                push!(stack, children_arr[cidx])
                push!(stack_ignore, should_ignore)
            end
            if should_ignore
                continue
            end
            if nodes[r].type == Nonlinear.NODE_VARIABLE
                push!(nonlinear_group, nodes[r].index)
            elseif nodes[r].type == Nonlinear.NODE_SUBEXPRESSION
                # append all variables in subexpression
                union!(nonlinear_group, subexpression_variables[nodes[r].index])
            end
        end
        for i_ in 1:nonlinear_group.nnz
            i = nonlinear_group.nzidx[i_]
            for j_ in 1:nonlinear_group.nnz
                j = nonlinear_group.nzidx[j_]
                if j > i
                    continue  # Only lower triangle.
                end
                push!(edge_list, (i, j))
            end
        end
    end
    return edge_list
end

"""
    _list_subexpressions(nodes::Vector{Nonlinear.Node})

Returns the list of subexpressions which a given tape depends on directly
"""
function _list_subexpressions(nodes::Vector{Nonlinear.Node})
    indices = Set{Int}(
        n.index for n in nodes if n.type == Nonlinear.NODE_SUBEXPRESSION
    )
    return sort(collect(indices))
end

"""
    _topological_sort(
        starts::Vector{Int},
        subexpressions::Vector{Vector{Nonlinear.Node}},
        subexpression_dependency_graph::Vector{Vector{Int}} =
            Vector{Vector{Int}}(undef, length(subexpressions)),
    )

Return a topologically sorted list of the integer subexpresssion indices that
need to be computed to evaluate `subexpressions[s]` for all `s in starts`.

`starts` should be ordered, and not contain duplicates.

`subexpression_dependency_graph[i]` is a lazily-computed list of "out" edges
from node `i`, in terms of the integer-valued subexpression index (i.e.,
`node.index`). This list should be unique and ordered.

If calling `_topological_sort` a single time, you may omit the
`subexpression_dependency_graph` argument.

However, if calling `_topological_sort` multiple times on the _same_ vector of
subexpresssions, you should create `subexpression_dependency_graph` once (either
as the uninitialized vector, or by explicitly computing the full
`subexpression_dependency_graph`), and pass it in.

## Notes

* It is important to not use recursion here, because expressions may have
  arbitrary levels of nesting!
* This function assumes `subexpressions` is acyclic.
"""
function _topological_sort(
    starts,
    subexpressions::Vector{Vector{Nonlinear.Node}},
    subexpression_dependency_graph::Vector{Vector{Int}} = Vector{Vector{Int}}(
        undef,
        length(subexpressions),
    ),
)
    ordered = Int[]
    in_order = fill(false, length(subexpressions))
    stack = Tuple{Int,Bool}[]
    for s in starts
        if in_order[s]
            continue  # s is already in `ordered`.
        end
        push!(stack, (s, true))
        while !isempty(stack)
            node, needs_checking = pop!(stack)
            if !needs_checking
                # We must be returning to this node for a second time, and we
                # have already checked all of the children. Therefore, we can
                # add it to the set of ordered nodes.
                push!(ordered, node)
                in_order[node] = true
                continue
            elseif in_order[node]
                continue  # This node has already been added to `ordered`.
            end
            # Re-add the node to the stack, but set the `false` flag this time
            # so next time we visit it, it will go on the `ordered` list
            # instead.
            push!(stack, (node, false))
            if !isassigned(subexpression_dependency_graph, node)
                subexpression_dependency_graph[node] =
                    _list_subexpressions(subexpressions[node])
            end
            for child in subexpression_dependency_graph[node]
                if !in_order[child]
                    push!(stack, (child, true))
                end
            end
        end
    end
    return ordered
end

"""
    _order_subexpressions(
        main_expressions::Vector{Vector{Nonlinear.Node}},
        subexpressions::Vector{Vector{Nonlinear.Node}};
    )

Topologically sort the subexpression needed to evaluate `main_expressions`.

Returns two things:

 * A `Vector{Int}` containing the ordered list of subexpression-indices that
   need to be evaluated to compute all `main_expressions`
 * A `Vector{Vector{Int}}`, containing a list of ordered lists of
   subexpression-indices that need to be evaluated to compute
   `main_expressions[i]`.

**Warning:** This doesn't handle cyclic expressions! But this should be fine
because we can't compute them in JuMP anyway.
"""
function _order_subexpressions(
    main_expressions::Vector{Vector{Nonlinear.Node}},
    subexpressions::Vector{Vector{Nonlinear.Node}},
)
    # The graph of node dependencies. Constructed lazily.
    subexpression_dependency_graph =
        Vector{Vector{Int}}(undef, length(subexpressions))
    # Node dependencies of the main expressions.
    starts = Set{Int}()
    individual_sorts = Vector{Int}[]
    for expression in main_expressions
        s = _list_subexpressions(expression)
        union!(starts, s)
        push!(
            individual_sorts,
            _topological_sort(
                s,
                subexpressions,
                subexpression_dependency_graph,
            ),
        )
    end
    full_sort = _topological_sort(
        starts,
        subexpressions,
        subexpression_dependency_graph,
    )
    return full_sort, individual_sorts
end

"""
    _order_subexpressions(
        main_expression::Vector{Nonlinear.Node},
        subexpressions::Vector{Vector{Nonlinear.Node}},
    )

Return a `Vector{Int}` containing the topologically ordered list of
subexpression-indices that need to be evaluated to compute `main_expression`.
"""
function _order_subexpressions(
    main_expression::Vector{Nonlinear.Node},
    subexpressions::Vector{Vector{Nonlinear.Node}},
)
    s = _list_subexpressions(main_expression)
    return _topological_sort(s, subexpressions)
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
@doc raw"""
    y = prox_distance(M,λ,f,x [,p=2])

compute the proximal map $\operatorname{prox}_{\lambda\varphi}$ with
parameter λ of $\varphi(x) = \frac{1}{p}d_{\mathcal M}^p(f,x)$.

# Input
* `M` – a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) $\mathcal M$
* `λ` – the prox parameter
* `f` – a point $f ∈ \mathcal M$ (the data)
* `x` – the argument of the proximal map

# Optional argument
* `p` – (`2`) exponent of the distance.

# Ouput
* `y` – the result of the proximal map of $\varphi$
"""
function prox_distance(M::Manifold, λ, f, x, p::Int=2)
    d = distance(M, f, x)
    if p == 2
        t = λ / (1 + λ)
    elseif p == 1
        if λ < d
            t = λ / d
        else
            t = 1.0
        end
    else
        throw(
            ErrorException(
                "Proximal Map of distance(M,f,x) not implemented for p=$(p) (requires p=1 or 2)",
            ),
        )
    end
    return exp(M, x, log(M, x, f), t)
end
@doc raw"""
    (y1,y2) = prox_TV(M,λ,(x1,x2) [,p=1])

Compute the proximal map $\operatorname{prox}_{\lambda\varphi}$ of
$\varphi(x,y) = d_{\mathcal M}^p(x,y)$ with
parameter `λ`.

# Input
* `M` – a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold)
* `λ` – a real value, parameter of the proximal map
* `(x1,x2)` – a tuple of two points,

# Optional
(default is given in brackets)
* `p` – (1) exponent of the distance of the TV term

# Ouput
* `(y1,y2)` – resulting tuple of points of the
  $\operatorname{prox}_{\lambda\varphi}($ `(x1,x2)` $)$
"""
function prox_TV(M::mT, λ::Number, pointTuple::Tuple{T,T}, p::Int=1) where {mT<:Manifold,T}
    x1 = pointTuple[1]
    x2 = pointTuple[2]
    d = distance(M, x1, x2)
    if p == 1
        t = min(0.5, λ / d)
    elseif p == 2
        t = λ / (1 + 2 * λ)
    else
        throw(
            ErrorException(
                "Proximal Map of TV(M,x1,x2,p) not implemented for p=$(p) (requires p=1 or 2)",
            ),
        )
    end
    return (exp(M, x1, log(M, x1, x2), t), exp(M, x2, log(M, x2, x1), t))
end
@doc raw"""
    ξ = prox_TV(M,λ,x [,p=1])

compute the proximal maps $\operatorname{prox}_{\lambda\varphi}$ of
all forward differences orrucirng in the power manifold array, i.e.
$\varphi(xi,xj) = d_{\mathcal M}^p(xi,xj)$ with `xi` and `xj` are array
elemets of `x` and `j = i+e_k`, where `e_k` is the $k$th unitvector.
The parameter `λ` is the prox parameter.

# Input
* `M` – a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold)
* `λ` – a real value, parameter of the proximal map
* `x` – a point.

# Optional
(default is given in brackets)
* `p` – (1) exponent of the distance of the TV term

# Ouput
* `y` – resulting  point containinf with all mentioned proximal
  points evaluated (in a cylic order).
"""
function prox_TV(M::PowerManifold{𝔽,N,T}, λ, x, p::Int=1) where {𝔽,N<:Manifold,T}
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R).I
    y = copy(x)
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i -> (i == k) ? 1 : 0, d)) #k th unit vector
        for l in 0:1
            for i in R # iterate over all pixel
                if (i[k] % 2) == l
                    J = i.I .+ ek.I #i + e_k is j
                    if all(J .<= maxInd) # is this neighbor in range?
                        j = CartesianIndex(J...) # neigbbor index as Cartesian Index
                        (y[i], y[j]) = prox_TV(M.manifold, λ, (y[i], y[j]), p) # Compute TV on these
                    end
                end
            end # i in R
        end # even odd
    end # directions
    return y
end
@doc raw"""
    ξ = prox_parallel_TV(M,λ,x [,p=1])

compute the proximal maps $\operatorname{prox}_{\lambda\varphi}$ of
all forward differences orrucirng in the power manifold array, i.e.
$\varphi(xi,xj) = d_{\mathcal M}^p(xi,xj)$ with `xi` and `xj` are array
elemets of `x` and `j = i+e_k`, where `e_k` is the $k$th unitvector.
The parameter `λ` is the prox parameter.

# Input
* `M`     – a `PowerManifold` manifold
* `λ`     – a real value, parameter of the proximal map
* `x`     – a point

# Optional
(default is given in brackets)
* `p` – (`1`) exponent of the distance of the TV term

# Ouput
* `y`  – resulting Array of points with all mentioned proximal
  points evaluated (in a parallel within the arrays elements).

*See also* [`prox_TV`](@ref)
"""
function prox_parallel_TV(M::PowerManifold, λ, x::Array{T,1}, p::Int=1) where {T}
    R = CartesianIndices(x[1])
    d = ndims(x[1])
    if length(x) != 2 * d
        throw(
            ErrorException(
                "The number of inputs from the array ($(length(x))) has to be twice the data dimensions ($(d)).",
            ),
        )
    end
    maxInd = Tuple(last(R))
    # create an array for even/odd splitted proxes along every dimension
    y = reshape(deepcopy(x), d, 2)
    x = reshape(x, d, 2)
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i -> (i == k) ? 1 : 0, d)) #k th unit vector
        for l in 0:1 # even odd
            for i in R # iterate over all pixel
                if (i[k] % 2) == l
                    J = i.I .+ ek.I #i + e_k is j
                    if all(J .<= maxInd) # is this neighbor in range?
                        j = CartesianIndex(J...) # neigbbor index as Cartesian Index
                        # parallel means we apply each (direction even/odd) to a seperate copy of the data.
                        (y[k, l + 1][i], y[k, l + 1][j]) = prox_TV(
                            M.manifold, λ, (x[k, l + 1][i], x[k, l + 1][j]), p
                        ) # Compute TV on these
                    end
                end
            end # i in R
        end # even odd
    end # directions
    return y[:] # return as onedimensional array
end
@doc raw"""
    (y1,y2,y3) = prox_TV2(M,λ,(x1,x2,x3),[p=1], kwargs...)

Compute the proximal map $\operatorname{prox}_{\lambda\varphi}$ of
$\varphi(x_1,x_2,x_3) = d_{\mathcal M}^p(c(x_1,x_3),x_2)$ with
parameter `λ`>0, where $c(x,z)$ denotes the mid point of a shortest
geodesic from `x1` to `x3` that is closest to `x2`.

# Input

* `M`          – a manifold
* `λ`          – a real value, parameter of the proximal map
* `(x1,x2,x3)` – a tuple of three points

* `p` – (`1`) exponent of the distance of the TV term

# Optional
`kwargs...` – parameters for the internal [`subgradient_method`](@ref)
    (if `M` is neither `Euclidean` nor `Circle`, since for these a closed form
    is given)

# Output
* `(y1,y2,y3)` – resulting tuple of points of the proximal map
"""
function prox_TV2(
    M::Manifold,
    λ,
    pointTuple::Tuple{T,T,T},
    p::Int=1;
    stopping_criterion::StoppingCriterion=StopAfterIteration(5),
    kwargs...,
) where {T}
    if p != 1
        throw(
            ErrorException(
                "Proximal Map of TV2(M,λ,pT,p) not implemented for p=$(p) (requires p=1) on general manifolds.",
            ),
        )
    end
    PowX = SVector(pointTuple)
    PowM = PowerManifold(M, NestedPowerRepresentation(), 3)
    xR = PowX
    F(x) = 1 / 2 * distance(PowM, PowX, x)^2 + λ * costTV2(PowM, x)
    ∂F(x) = log(PowM, x, PowX) + λ * ∇TV2(PowM, x)
    subgradient_method!(PowM, F, ∂F, xR; stopping_criterion=stopping_criterion, kwargs...)
    return (xR...,)
end
function prox_TV2(::Circle, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T}
    w = @SVector [1.0, -2.0, 1.0]
    x = SVector(pointTuple)
    if p == 1 # Theorem 3.5 in Bergmann, Laus, Steidl, Weinmann, 2014.
        sr_dot_xw = sym_rem(sum(x .* w))
        m = min(λ, abs(sr_dot_xw) / (dot(w, w)))
        s = sign(sr_dot_xw)
        return sym_rem.(x .- m .* s .* w)
    elseif p == 2 # Theorem 3.6 ibd.
        t = λ * sym_rem(sum(x .* w)) / (1 + λ * dot(w, w))
        return sym_rem.(x - t .* w)
    else
        throw(
            ErrorException(
                "Proximal Map of TV2(Circle,λ,pT,p) not implemented for p=$(p) (requires p=1 or 2)",
            ),
        )
    end
end
function prox_TV2(::Euclidean, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T}
    w = @SVector [1.0, -2.0, 1.0]
    x = SVector(pointTuple)
    if p == 1 # Example 3.2 in Bergmann, Laus, Steidl, Weinmann, 2014.
        m = min.(Ref(λ), abs.(x .* w) / (dot(w, w)))
        s = sign.(sum(x .* w))
        return x .- m .* s .* w
    elseif p == 2 # Theorem 3.6 ibd.
        t = λ * sum(x .* w) / (1 + λ * dot(w, w))
        return x .- t .* w
    else
        throw(
            ErrorException(
                "Proximal Map of TV2(Euclidean,λ,pT,p) not implemented for p=$(p) (requires p=1 or 2)",
            ),
        )
    end
end
@doc raw"""
    ξ = prox_TV2(M,λ,x,[p])

compute the proximal maps $\operatorname{prox}_{\lambda\varphi}$ of
all centered second order differences orrucirng in the power manifold array, i.e.
$\varphi(x_k,x_i,x_j) = d_2(x_k,x_i.x_j)$, where $k,j$ are backward and forward
neighbors (along any dimension in the array of `x`).
The parameter `λ` is the prox parameter.

# Input
* `M` – a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold)
* `λ` – a real value, parameter of the proximal map
* `x` – a points.

# Optional
(default is given in brackets)
* `p` – (`1`) exponent of the distance of the TV term

# Ouput
* `y` – resulting point with all mentioned proximal points
  evaluated (in a cylic order).
"""
function prox_TV2(M::PowerManifold{N,T}, λ, x, p::Int=1) where {N,T}
    power_size = power_dimensions(M)
    R = CartesianIndices(power_size)
    d = length(size(x))
    minInd = first(R).I
    maxInd = last(R).I
    y = copy(x)
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i -> (i == k) ? 1 : 0, d)) #k th unit vector
        for l in 0:1
            for i in R # iterate over all pixel
                if (i[k] % 3) == l
                    JForward = i.I .+ ek.I #i + e_k
                    JBackward = i.I .- ek.I # i - e_k
                    if all(JForward .<= maxInd) && all(JBackward .>= minInd)
                        (y[jBackward], y[i], y[jForward]) =
                            prox_TV2(
                                M.manifold,
                                λ,
                                (y[M, JBackward...], y[M, i.I...], y[M, JForward...]),
                                p,
                            ).data # Compute TV on these
                    end
                end # if mod 3
            end # i in R
        end # for mod 3
    end # directions
    return y
end
@doc raw"""
    project_collaborative_TV(M,λ,x [,p=2,q=1])

compute the projection onto collaborative Norm unit (or α-) ball, i.e. of the function

```math
F^q(x) = \sum_{i\in\mathcal G}
  \Bigl( \sum_{j\in\mathcal I_i}
    \sum_{k=1^d} \lVert X_{i,j}\rVert_x^p\Bigr)^\frac{q/p},
```

where $\mathcal G$ is the set of indices for $x\in\mathcal M$ and $\mathcal I_i$
is the set of its forward neighbors.
This is adopted from the paper by Duran, Möller, Sbert, Cremers:
_Collaborative Total Variation: A General Framework for Vectorial TV Models_
(arxiv: [1508.01308](https://arxiv.org/abs/1508.01308)), where the most inner
norm is not on a manifold but on a vector space, see their Example 3 for
details.
"""
function project_collaborative_TV(N::PowerManifold, λ, x, Ξ, p=2.0, q=1.0, α=1.0)
    pdims = power_dimensions(N)
    if length(pdims) == 1
        d = 1
        s = 1
        iRep = (1,)
    else
        d = pdims[end]
        s = length(pdims) - 1
        if s != d
            throw(
                ErrorException(
                    "the last dimension ($(d)) has to be equal to the number of the previous ones ($(s)) but its not.",
                ),
            )
        end
        iRep = (Integer.(ones(d))..., d)
    end
    if q == 1 # Example 3 case 2
        if p == 1
            normΞ = norm.(Ref(N.manifold), x, Ξ)
            return max.(normΞ .- λ, 0.0) ./ ((normΞ .== 0) .+ normΞ) .* Ξ
        elseif p == 2 # Example 3 case 3
            norms = sqrt.(sum(norm.(Ref(N.manifold), x, Ξ) .^ 2; dims=d + 1))
            if length(iRep) > 1
                norms = repeat(norms; inner=iRep)
            end
            # if the norm is zero add 1 to avoid division by zero, also then the
            # nominator is already (max(-λ,0) = 0) so it stays zero then
            return max.(norms .- λ, 0.0) ./ ((norms .== 0) .+ norms) .* Ξ
        else
            throw(ErrorException("The case p=$p, q=$q is not yet implemented"))
        end
    elseif q == Inf
        if p == 2
            norms = sqrt.(sum(norm.(Ref(N.manifold), x, Ξ) .^ 2; dims=d + 1))
            if length(iRep) > 1
                norms = repeat(norms; inner=iRep)
            end
        elseif p == 1
            norms = sum(norm.(Ref(N.manifold), x, Ξ); dims=d + 1)
            if length(iRep) > 1
                norms = repeat(norms; inner=iRep)
            end
        elseif p == Inf
            norms = norm.(Ref(N.manifold), x, Ξ)
        else
            throw(ErrorException("The case p=$p, q=$q is not yet implemented"))
        end
        return (α .* Ξ) ./ max.(Ref(α), norms)
    end # end q
    return throw(ErrorException("The case p=$p, q=$q is not yet implemented"))
end
function project_collaborative_TV(N::PowerManifold, λ, x, Ξ, p::Int, q::Float64=1.0, α=1.0)
    return project_collaborative_TV(N, λ, x, Ξ, Float64(p), q, α)
end
function project_collaborative_TV(N::PowerManifold, λ, x, Ξ, p::Float64, q::Int, α=1.0)
    return project_collaborative_TV(N, λ, x, Ξ, p, Float64(q), α)
end
function project_collaborative_TV(N::PowerManifold, λ, x, Ξ, p::Int, q::Int, α=1.0)
    return project_collaborative_TV(N, λ, x, Ξ, Float64(p), Float64(q), α)
end

ext: jl


text:
using Z3

mutable struct Z3Solver <: NLSolver
    vars::Dict{Symbol,Type}
    hard_clauses::ClauseSet
    soft_clauses::ClauseSet
    function Z3Solver()
        new(Dict(), ClauseSet(), ClauseSet())
    end
end

function add!(s::Z3Solver, cs::ClauseSet)
    s.hard_clauses &= expand(cs)
end
add_soft!(s::Z3Solver, cs::ClauseSet) = s.soft_clauses &= cs
add_vars!(s::Z3Solver, d::Dict{Symbol,Type}) = push!(s.vars, d...)

# ------------------------------------------------------------------------------

function solve(s::Z3Solver; timeout::Int=-1)
    ctx = main_ctx()
    if timeout > 0
        set(ctx, "timeout", timeout*1000)
    end
    z3vars = Dict{Symbol,Z3Expr}()
    for (k,v) in s.vars
        push!(z3vars, k=>real_const(ctx, string(k)))
    end

    z3 = Solver(ctx)
    for cl in s.hard_clauses
        add(z3, Z3Expr(ctx, z3vars, cl))
    end

    elapsed_sec = @elapsed res = check(z3)
    elapsed = Millisecond(round(elapsed_sec*1000))
    if res == Z3.sat
        return NLSat.sat, elapsed, parse_model(get_model(z3))
    elseif res == Z3.unsat
        return NLSat.unsat, elapsed, nothing
    elseif res == Z3.unknown
        reason = reason_unknown(z3)
        if reason == "timeout"
            return NLSat.timeout, elapsed, nothing
        end
        return NLSat.unknown, elapsed, nothing
    end
end


text:
#=
DelayABM method for system of fractional delay differential equations.

```tex
@inproceedings{Bhalekar2011APS,
  title={A PREDICTOR-CORRECTOR SCHEME FOR SOLVING NONLINEAR DELAY DIFFERENTIAL EQUATIONS OF FRACTIONAL ORDER},
  author={Sachin Bhalekar and Varsha Daftardar-Gejji},
  year={2011}
}
```
=#

function solve(FDDESys::FDDESystem, h, ::DelayABM)
    @unpack f, ϕ, α, τ, T = FDDESys
    len = length(ϕ)
    N::Int = round(Int, T/h)
    Ndelay = round(Int, τ/h)
    x1 = zeros(Ndelay+N+1, len)
    x = zeros(Ndelay+N+1, len)
    du = zeros(len)

    # Put the delay term in the array
    for i=1:len
        x[1:Ndelay, i] = ϕ[i]*ones(Ndelay)
    end

    x0 = copy(x[Ndelay, :])

    for i=1:len
        αi = α[i]
        f(du, x0, x[1, :], 0)
        x1[Ndelay+1, i] = x0[i] + h^αi*du[i]/(gamma(αi)*αi)
    end

    for i=1:len
        αi = α[i]
        f(du, x[Ndelay+1, :], x[1, :], 0)
        du1 = copy(du)
        f(du, x0, x[1, :], 0)
        x[Ndelay+1, i] = x0[i] + h^αi*(du1[i] + αi*du[i])/gamma(αi+2)
    end

    M1=zeros(len); N1=zeros(len)
    @fastmath @inbounds @simd for n=1:N
        for i=1:len
            αi = α[i]
            f(du, x0, x[1, :], 0)
            M1[i]=(n^(αi+1)-(n-αi)*(n+1)^αi)*du[i]
            N1[i]=((n+1)^αi-n^αi)*du[i]
        end

        @fastmath @inbounds @simd for j=1:n
            for i=1:len
                αi = α[i]
                f(du, x[Ndelay+j, :], x[j, :], 0)
                M1[i] = M1[i]+((n-j+2)^(αi+1)+(n-j)^(αi+1)-2*(n-j+1)^(αi+1))*du[i]
                N1[i] = N1[i]+((n-j+1)^αi-(n-j)^αi)*du[i]
            end
        end

        for i=1:len
            αi = α[i]
            x1[Ndelay+n+1, i] = x0[i]+h^αi*N1[i]/(gamma(αi)*αi)
            f(du, x[Ndelay+n+1, :], x[n+1, :], 0)
            x[Ndelay+n+1, i] = x0[i]+h^αi*(du[i]+M1[i])/gamma(αi+2)
        end
    end

    xresult = zeros(N-Ndelay+1, len)
    yresult = zeros(N-Ndelay+1, len)

    xresult[N-Ndelay+1, :] .= 0
    yresult[N-Ndelay+1, :] .= 0

    @fastmath @inbounds @simd for n=2*Ndelay+1:N+Ndelay+1
       xresult[n-2*Ndelay, :] = x[n, :]
       yresult[n-2*Ndelay, :] = x[n-Ndelay, :]
    end

    return xresult, yresult
end


text:
mutable struct PerfectBipartiteMatchingHungarianSolver <: PerfectBipartiteMatchingSolver
  weights_matrix::Matrix{Float64}

  function PerfectBipartiteMatchingHungarianSolver()
    return new(zeros(0, 0))
  end
end

function build!(solver::PerfectBipartiteMatchingHungarianSolver, n_arms::Int)
  solver.weights_matrix = zeros(n_arms, n_arms)
  nothing
end

function solve_linear(solver::PerfectBipartiteMatchingHungarianSolver, rewards::Dict{Tuple{Int, Int}, Float64})
  fill!(solver.weights_matrix, 0.0)
  for (k, v) in rewards
    solver.weights_matrix[k[1], k[2]] = v
  end
  solver.weights_matrix = maximum(solver.weights_matrix) .- solver.weights_matrix

  matching = Hungarian.munkres!(solver.weights_matrix)
  indices = CartesianIndices(matching)
  return [(indices[v][1], indices[v][2]) for v in findall(matching .== Hungarian.STAR)]
end

has_lp_formulation(::PerfectBipartiteMatchingHungarianSolver) = false


text:
```
ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется ровно одна внутренняя перегородка в форме прямоугольника.
Робот - в произвольной клетке поля между внешней и внутренней перегородками.

РЕЗУЛЬТАТ: Робот - в исходном положении и по всему периметру внутренней перегородки поставлены маркеры.
```

using HorizonSideRobots
originRobot = Robot(animate = true, "6example.sit")
include("../AbstractType.jl")
include("../functions.jl")

function test1(robot::Robot, borderSide::HorizonSide)
    height = 0
    width = 0
    stepsToBack = 0
    reqBorder = borderSide
    # можно было через while сделать: задать итератор = 4, и при каждом проходе по циклу менять значение side на next_side_pr(side)
    for side in [reqBorder, next_side_pr(reqBorder), next_side_pr(next_side_pr(reqBorder)), anti_side(next_side_pr(reqBorder))]
        while isborder(robot, side)
            if (side == next_side_pr(next_side_pr(reqBorder)))
                height += 1
            elseif (side == reqBorder)
                stepsToBack += 1
            end
            move!(robot, next_side(side))
            putmarker!(robot)
        end
        move!(robot, side)
        putmarker!(robot)
    end
    for _ in 1:(height - stepsToBack)
        move!(robot, next_side_pr(moveSide))
        putmarker!(robot)
    end
end

function findReq(robot::Robot, moveSide::HorizonSide)
    if isborder(robot, moveSide) && !isangle!(robot)
        isReq = false
        steps = 0
        while isborder(robot, moveSide) && !isborder(robot, next_side_pr(moveSide))
            move!(robot, next_side_pr(moveSide))
            steps += 1
        end
        if !isborder(robot, next_side_pr(moveSide))
            isReq = true
        end
        moving_defsteps!(robot, anti_side(next_side_pr(moveSide)), steps)
        return isReq
    else
        return false
    end
end

function moveLikeSnake!(robot, moveSides::NTuple{2,HorizonSide}, borderSide::HorizonSide)::Nothing
    flag = false
    while !(isborder(robot, borderSide)) && !(flag)
        for side in moveSides
            while !flag
                if (!isborder(robot, side))
                    move!(robot,side)
                elseif !(findReq(robot, side))
                    #println("stop there")
                    break
                else
                    flag = true
                    break
                end
            end

            if flag
                test1(robot, side)
                break
            end

            if !(isborder(robot, borderSide))
                move!(robot, borderSide)
            else
                break
            end
        end
    end
end

# master! главная функция программы
function master6!(robot::Robot)::Nothing
    back_path = moving_in_angle!(robot)

    moveLikeSnake!(robot, (Nord, Sud), Ost)

    moving_in_angle!(robot)
    moving_back_to_start!(robot, back_path)
end

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
using FilePathsBase
using JLSO
using LinearAlgebra
using Test

using FilePathsBase.TestPaths

include("testpkg.jl")

@testset "FilePathsBase" begin
    include("mode.jl")
    include("buffer.jl")
    include("system.jl")
    # include("path.jl")

    @static if Sys.isunix()
        # Test that our weird registered path works
        ps = PathSet(TestPkg.posix2test(tmpdir()) / "pathset_root"; symlink=true)

        @testset "$(typeof(ps.root))" begin
            testsets = [
                test_constructor,
                test_registration,
                test_show,
                test_parse,
                test_convert,
                test_components,
                test_parents,
                test_join,
                test_basename,
                test_filename,
                test_extensions,
                test_isempty,
                test_norm,
                test_real,
                test_relative,
                test_abs,
                test_isdir,
                test_isfile,
                test_stat,
                test_size,
                test_modified,
                test_created,
                test_cd,
                test_readpath,
                test_walkpath,
                test_read,
                test_write,
                test_mkdir,
                test_cp,
                test_mv,
                test_symlink,
                test_touch,
                test_tmpname,
                test_tmpdir,
                test_mktmp,
                test_mktmpdir,
                test_download,
                test_issocket,
                # These will also all work for our custom path type,
                # but many implementations won't support them.
                test_isfifo,
                test_ischardev,
                test_isblockdev,
                test_ismount,
                test_isexecutable,
                test_isreadable,
                test_iswritable,
                test_chown,
                test_chmod,
            ]

            # Run all of the automated tests
            test(ps, testsets)

            # TODO: Copy over specific tests that can't be tested reliably from the general case.
        end
    end
end


text:
using OptimalEstimationProject
using Test

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

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
using Pkg
pkg"precompile"

using MPIClusterManagers: mpiexec
using Test

# Code coverage command line options; must correspond to src/julia.h
# and src/ui/repl.c
const JL_LOG_NONE = 0
const JL_LOG_USER = 1
const JL_LOG_ALL = 2
const coverage_opts =
    Dict{Int, String}(JL_LOG_NONE => "none",
                      JL_LOG_USER => "user",
                      JL_LOG_ALL => "all")

function runtests()
    nprocs = clamp(Sys.CPU_THREADS, 2, 4)
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    extra_args = []
    @static if !Sys.iswindows()
        if occursin( "OpenRTE", read(`$mpiexec --version`, String))
            push!(extra_args,"--oversubscribe")
        end
    end

    coverage_opt = coverage_opts[Base.JLOptions().code_coverage]

    jlexec = `$exename --code-coverage=$coverage_opt`
    mpijlexec = `$mpiexec $extra_args -n $nprocs $jlexec`

    @info "Testing: test_cman_julia.jl"
    run(`$jlexec $(joinpath(@__DIR__, "test_cman_julia.jl"))`)

    @info "Testing: test_cman_mpi.jl"
    run(`$mpijlexec $(joinpath(@__DIR__, "test_cman_mpi.jl"))`)

    @info "Testing: test_cman_tcp.jl"
    run(`$mpijlexec $(joinpath(@__DIR__, "test_cman_tcp.jl"))`)
end

runtests()


text:
module TestLalePreprocessing

using Random
using Test
using Lale
using Statistics
using DataFrames: DataFrame, nrow

Random.seed!(1)

const IRIS = getiris()
extra = rand(150,3) |> x->DataFrame(x,:auto)
const X = hcat(IRIS[:,1:4],extra)
const Y = IRIS[:,5] |> Vector

# "KernelCenterer","MissingIndicator","KBinsDiscretizer","OneHotEncoder",
const preprocessors = [
     "FactorAnalysis", "FastICA", "IncrementalPCA",
     "KernelPCA",
     "MiniBatchDictionaryLearning",
     "MiniBatchSparsePCA", "NMF", "PCA",
     "TruncatedSVD",
     "SimpleImputer",
     "Binarizer",# "FunctionTransformer",
     "MaxAbsScaler", "MinMaxScaler", "Normalizer",
     "OrdinalEncoder", "PolynomialFeatures", "PowerTransformer",
     "QuantileTransformer", "RobustScaler"#, #"StandardScaler"
 ]

function fit_test(preproc::String,in::DataFrame,out::Vector)
   _preproc=LaleOp(preproc,"autogen")
   fit!(_preproc,in,out)
   @test _preproc.model != Dict()
   return _preproc
end

function transform_test(preproc::String,in::DataFrame,out::Vector)
   _preproc=LaleOp(preproc,"autogen")
   fit!(_preproc,in,out)
   res = transform(_preproc,in)
   @test size(res)[1] == size(out)[1]
end

@testset "lale preprocessors fit test" begin
   Random.seed!(123)
   for cl in preprocessors
      fit_test(cl,X,Y)
   end
end

@testset "lale preprocessors transform test" begin
   Random.seed!(123)
   for cl in preprocessors
      transform_test(cl,X,Y)
   end
end

function skptest()
    features = X
    labels = Y

    pca = LaleOp("PCA","autogen",n_components=3)
    fit!(pca,features)
    @test transform(pca,features) |> x->size(x,2) == 3

    pca = LaleOp("PCA","autogen")
    fit!(pca,features)
    @test transform(pca,features) |> x->size(x,2) == 7

    svd = LaleOp("TruncatedSVD","autogen",n_components=2)
    fit!(svd,features)
    @test transform(svd,features) |> x->size(x,2) == 2

    ica = LaleOp("FastICA","autogen",n_components=2)
    fit!(ica,features)
    @test transform(ica,features) |> x->size(x,2) == 2

    stdsc = LaleOp("StandardScaler","autogen")
    fit!(stdsc,features)
    @test abs(mean(transform(stdsc,features) |> Matrix)) < 0.00001

    minmax = LaleOp("MinMaxScaler","autogen")
    fit!(minmax,features)
    @test mean(transform(minmax,features) |> Matrix) > 0.30

    vote = VoteEnsemble()
    stack = StackEnsemble()
    best = BestLearner()
    cat = CatFeatureSelector()
    num = NumFeatureSelector()
    disc = CatNumDiscriminator()
    ohe = OneHotEncoder()

    mpipeline = @pipeline (stdsc |> pca) |> best
    fit!(mpipeline,features,labels)
    pred = transform(mpipeline,features)
    @test score(:accuracy,pred,labels) > 50.0

    fpipe = @pipeline ((cat + num) + (num + pca))  |> stack
    fit(fpipe,features,labels)
    @test ((transform(fpipe,features) .== labels) |> sum ) / nrow(features) > 0.50

end
@testset "scikit preprocessor fit/transform test with real data" begin
    Random.seed!(123)
    skptest()
end

end


text:
using Dates: parse_components, default_format
using TimeZones: ParseNextError, _parsesub_tzabbr, _parsesub_offset, _parsesub_time, _parsesub_tzdate, _parsesub_tz

@testset "parse" begin
    @test isequal(
        parse(ZonedDateTime, "2017-11-14 11:03:53 +0100", dateformat"yyyy-mm-dd HH:MM:SS zzzzz"),
        ZonedDateTime(2017, 11, 14, 11, 3, 53, tz"UTC+01"),
    )
    @test isequal(
        parse(ZonedDateTime, "2016-04-11 08:00 UTC", dateformat"yyyy-mm-dd HH:MM ZZZ"),
        ZonedDateTime(2016, 4, 11, 8, tz"UTC"),
    )
    # two-digit time zone
    @test isequal(
        parse(ZonedDateTime, "2000+00", dateformat"yyyyz"),
        ZonedDateTime(2000, tz"UTC"),
    )
    @test_throws ArgumentError parse(ZonedDateTime, "2016-04-11 08:00 EST", dateformat"yyyy-mm-dd HH:MM zzz")
    # test AbstractString
    @test isequal(
        parse(ZonedDateTime, Test.GenericString("2018-01-01 00:00 UTC"), dateformat"yyyy-mm-dd HH:MM ZZZ"),
        ZonedDateTime(2018, 1, 1, 0, tz"UTC"),
    )

end

@testset "tryparse" begin
    @test isequal(
        tryparse(ZonedDateTime, "2013-03-20 11:00:00+04:00", dateformat"y-m-d H:M:SSz"),
        ZonedDateTime(2013, 3, 20, 11, tz"UTC+04"),
    )
    @test isequal(
        tryparse(ZonedDateTime, "2016-04-11 08:00 EST", dateformat"yyyy-mm-dd HH:MM zzz"),
        nothing,
    )
end

@testset "parse components" begin
    local test = ("2017-11-14 11:03:53 +0100", dateformat"yyyy-mm-dd HH:MM:SS zzzzz")
    local expected = [
        Dates.Year(2017),
        Dates.Month(11),
        Dates.Day(14),
        Dates.Hour(11),
        Dates.Minute(3),
        Dates.Second(53),
        tz"UTC+01",
    ]
    @test parse_components(test...) == expected
end

@testset "default format" begin
    @test default_format(ZonedDateTime) === TimeZones.ISOZonedDateTimeFormat
end

@testset "parse constructor" begin
    @test isequal(
        ZonedDateTime("2000-01-02T03:04:05.006+0700"),
        ZonedDateTime(2000, 1, 2, 3, 4, 5, 6, tz"UTC+07")
    )
    @test isequal(
        ZonedDateTime("2000-01-02T03:04:05.006Z"),
        ZonedDateTime(2000, 1, 2, 3, 4, 5, 6, tz"UTC+00")
    )
    @test isequal(
        ZonedDateTime("2018-11-01-0600", dateformat"yyyy-mm-ddzzzz"),
        ZonedDateTime(2018, 11, 1, tz"UTC-06"),
    )

    # Validate that error message contains the original string and the format used
    str = "2018-11-01"
    try
        ZonedDateTime(str)
    catch e
        @test e isa ArgumentError
        @test occursin(str, e.msg)
        @test occursin(string(TimeZones.ISOZonedDateTimeFormat), e.msg)
    end
end

@testset "_parsesub_tzabbr" begin
    empty_msg = "Time zone abbreviation must start with a letter or the less-than (<) character"
    not_closed_msg = "Expected expanded time zone abbreviation end with the greater-than sign (>)"
    three_char_sim_msg = "Time zone abbreviation must be at least three alphabetic characters"
    three_char_exp_msg = "Time zone abbreviation must be at least three characters which are either alphanumeric, the plus sign (+), or the minus sign (-)"

    @test _parsesub_tzabbr("") == ParseNextError(empty_msg, "", 1, 0)
    @test _parsesub_tzabbr("&") == ParseNextError(empty_msg, "&", 1, 1)
    @test _parsesub_tzabbr("FOO") == ("FOO", 4)
    @test _parsesub_tzabbr("FOO&") == ("FOO", 4)
    @test _parsesub_tzabbr("FOO+1") == ("FOO", 4)
    @test _parsesub_tzabbr("<FOO+1>") == ("FOO+1", 8)
    @test _parsesub_tzabbr("<FOO+1") == ParseNextError(not_closed_msg, "<FOO+1", 1, 7)
    @test _parsesub_tzabbr("<") == ParseNextError(not_closed_msg, "<", 1, 2)
    @test _parsesub_tzabbr(">") == ParseNextError(empty_msg, ">", 1, 1)
    @test _parsesub_tzabbr("AB") == ParseNextError(three_char_sim_msg, "AB", 1, 2)
    @test _parsesub_tzabbr("<>") == ParseNextError(three_char_exp_msg, "<>", 2, 1)
    @test _parsesub_tzabbr("αβc") == ParseNextError(empty_msg, "αβc", 1, 1)
end

@testset "_parsesub_offset" begin
    end_of_string_msg = "Expected offset and instead found end of string"
    missing_hours_msg = "Expected offset hour digits"
    hours_range_msg = "Hours outside of expected range [0, 24]"
    minutes_range_msg = "Minutes outside of expected range [0, 59]"
    minutes_digits_msg = "Expected offset minute digits after colon delimiter"
    seconds_range_msg = "Seconds outside of expected range [0, 59]"
    seconds_digits_msg = "Expected offset second digits after colon delimiter"

    @testset "sign" begin
        @test _parsesub_offset("") == ParseNextError(end_of_string_msg, "", 1, 0)
        @test _parsesub_offset("&") == ParseNextError(missing_hours_msg, "&", 1, 1)
        @test _parsesub_offset("-") == ParseNextError("Offset sign (-) is not followed by a value", "-", 1, 1)
        @test _parsesub_offset("+") == ParseNextError("Offset sign (+) is not followed by a value", "+", 1, 1)
    end

    @testset "hours" begin
        @test _parsesub_offset("0") == (0, 2)
        @test _parsesub_offset("0&") == (0, 2)

        @test _parsesub_offset("1")   == ( 1 * 3600, 2)
        @test _parsesub_offset("-1")  == (-1 * 3600, 3)
        @test _parsesub_offset("+1")  == (+1 * 3600, 3)
        @test _parsesub_offset("-24") == (-24 * 3600, 4)
        @test _parsesub_offset("+24") == (+24 * 3600, 4)

        @test _parsesub_offset("-25") == ParseNextError(hours_range_msg, "-25", 2, 3)
        @test _parsesub_offset("+25") == ParseNextError(hours_range_msg, "+25", 2, 3)
    end

    @testset "minutes" begin
        @test _parsesub_offset("0:59")  == ( 59 * 60, 5)
        @test _parsesub_offset("-0:59") == (-59 * 60, 6)
        @test _parsesub_offset("+0:59") == (+59 * 60, 6)

        @test _parsesub_offset("-0:60") == ParseNextError(minutes_range_msg, "-0:60", 4, 5)
        @test _parsesub_offset("+0:60") == ParseNextError(minutes_range_msg, "+0:60", 4, 5)
        @test _parsesub_offset("0:") == ParseNextError(minutes_digits_msg, "0:", 3, 2)
        @test _parsesub_offset("0:-1") == ParseNextError(minutes_digits_msg, "0:-1", 3, 3)
    end

    @testset "seconds" begin
        @test _parsesub_offset("0:0:59")  == ( 59, 7)
        @test _parsesub_offset("-0:0:59") == (-59, 8)
        @test _parsesub_offset("+0:0:59") == (+59, 8)
        @test _parsesub_offset("0:0:59:") == ( 59, 7)

        @test _parsesub_offset("-0:0:60") == ParseNextError(seconds_range_msg, "-0:0:60", 6, 7)
        @test _parsesub_offset("+0:0:60") == ParseNextError(seconds_range_msg, "+0:0:60", 6, 7)
        @test _parsesub_offset("0:0:") == ParseNextError(seconds_digits_msg, "0:0:", 5, 4)
        @test _parsesub_offset("0:0:-1") == ParseNextError(seconds_digits_msg, "0:0:-1", 5, 5)
    end
end

@testset "_parsesub_time" begin
    # Note: `_parsesub_time` primarily makes use of `_parsesub_offset`. Additional tests
    # should be added here if that is no longer the case.
    @test _parsesub_time("1") == (3600, 2)
    @test _parsesub_time("-1") == ParseNextError("Time should not have a sign", "-1", 1, 1)
    @test _parsesub_time("+1") == ParseNextError("Time should not have a sign", "+1", 1, 1)
end

@testset "_parsesub_tzdate" begin
    @testset "Julian day" begin
        # Note: Hard-coding years in tests below make failed testsets easier to read

        # Non-leap year
        for d in [1; 59:61; 365]
            f, i = TimeZones._parsesub_tzdate("J$d")
            @test f(2019) == Date(2019) + Day(d - 1)
            @test dayofyear(f(2019)) == d
        end

        # Leap year
        for (j, d) in zip([1; 59:61; 365], [1; 59; 61:62; 366])
            f, i = TimeZones._parsesub_tzdate("J$j")
            @test f(2020) == Date(2020) + Day(d - 1)
            @test dayofyear(f(2020)) == d
        end
    end

    @testset "zero-based Julian day" begin
        # Note: Hard-coding years in tests below make failed testsets easier to read

        # Non-leap year
        # Note: For non-leap years day 365 is actually the first day of the next year which
        # is why the test for `dayofyear` requires mod 365.
        for d in [0; 58:60; 365]
            f, i = _parsesub_tzdate("$d")
            @test f(2019) == Date(2019) + Day(d)
            @test dayofyear(f(2019)) == (d % 365) + 1
        end

        # Leap year
        for d in [0; 58:60; 365]
            f, i = _parsesub_tzdate("$d")
            @test f(2020) == Date(2020) + Day(d)
            @test dayofyear(f(2020)) == d + 1
        end
    end

    @testset "month, week-of-month, day-of-week" begin
        @testset "First Sunday in October" begin
            f, i = _parsesub_tzdate("M10.1.0")
            @test f(2020) == Date(2020, 10, 4)
            @test f(2019) == Date(2019, 10, 6)
            @test f(2018) == Date(2018, 10, 7)
            @test f(2017) == Date(2017, 10, 1)
            @test i == 8
        end

        @testset "Third Sunday in March" begin
            f, i = _parsesub_tzdate("M3.3.0")
            @test f(2020) == Date(2020, 3, 15)
            @test f(2019) == Date(2019, 3, 17)
            @test f(2018) == Date(2018, 3, 18)
            @test f(2017) == Date(2017, 3, 19)
            @test i == 7
        end

        @testset "Last Saturday in January" begin
            # Note: 2017 - 2020 only have 4 Saturdays in January
            f1, i = _parsesub_tzdate("M1.5.6")
            @test i == 7

            f2, i = _parsesub_tzdate("M1.4.6")
            @test i == 7

            @test f1(2020) == f1(2020) == Date(2020, 1, 25)
            @test f1(2019) == f2(2019) == Date(2019, 1, 26)
            @test f1(2018) == f2(2018) == Date(2018, 1, 27)
            @test f1(2017) == f2(2017) == Date(2017, 1, 28)
        end
    end
end

@testset "_parsesub_tz" begin
    dst_start_date_msg = "Unable to parse daylight saving start date. Expected date and instead found end of string"
    dst_start_time_msg = "Expected daylight saving start time and instead found end of string"
    dst_end_missing_msg = "Expected to find daylight saving end and instead found end of string"
    dst_end_date_msg = "Unable to parse daylight saving end date. Expected date and instead found end of string"

    @testset "empty" begin
        @test _parsesub_tz("") == (FixedTimeZone("UTC"), 1)
    end

    @testset "standard time only" begin
        @test _parsesub_tz("FOO") == ParseNextError("Expected standard offset and instead found end of string", "FOO", 4, 3)
        @test _parsesub_tz("FOO+1") == (FixedTimeZone("FOO", -3600), 6)
        @test _parsesub_tz("<FOO+1>-1") == (FixedTimeZone("FOO+1", 3600), 10)
    end

    @testset "standard/daylight saving time" begin
        tz, i = _parsesub_tz("FOO+0BAR")
        bar, foo = filter(t -> year(t.utc_datetime) == 2020, tz.transitions)
        @test foo == Transition(DateTime(2020, 11, 1, 1), FixedTimeZone("FOO", 0))
        @test bar == Transition(DateTime(2020, 3, 8, 2), FixedTimeZone("BAR", 0, 3600))
        @test i == 9

        tz, i = _parsesub_tz("FOO+0BAR+1")
        bar, foo = filter(t -> year(t.utc_datetime) == 2020, tz.transitions)
        @test foo == Transition(DateTime(2020, 11, 1, 3), FixedTimeZone("FOO", 0))
        @test bar == Transition(DateTime(2020, 3, 8, 2), FixedTimeZone("BAR", 0, -3600))
        @test i == 11
    end

    @testset "transition dates" begin
        @test _parsesub_tz("FOO+0BAR+0,") == ParseNextError(dst_start_date_msg, "FOO+0BAR+0,", 12, 11)
        @test _parsesub_tz("FOO+0BAR+0,0") == ParseNextError(dst_end_missing_msg, "FOO+0BAR+0,0", 13, 12)
        @test _parsesub_tz("FOO+0BAR+0,0,") == ParseNextError(dst_end_date_msg, "FOO+0BAR+0,0,", 14, 13)

        tz, i = _parsesub_tz("FOO+0BAR+0,0,1")
        bar, foo = filter(t -> year(t.utc_datetime) == 2020, tz.transitions)
        @test foo == Transition(DateTime(2020, 1, 2, 2), FixedTimeZone("FOO", 0))
        @test bar == Transition(DateTime(2020, 1, 1, 2), FixedTimeZone("BAR", 0))
        @test i == 15
    end

    @testset "transition times" begin
        @test _parsesub_tz("FOO+0BAR+0,0/") == ParseNextError(dst_start_time_msg, "FOO+0BAR+0,0/", 14, 13)
        @test _parsesub_tz("FOO+0BAR+0,0/0") == ParseNextError(dst_end_missing_msg, "FOO+0BAR+0,0/0", 15, 14)
        @test _parsesub_tz("FOO+0BAR+0,0/0,") == ParseNextError(dst_end_date_msg, "FOO+0BAR+0,0/0,", 16, 15)

        tz, i = _parsesub_tz("FOO+0BAR+0,0/3,1/4")
        bar, foo = filter(t -> year(t.utc_datetime) == 2020, tz.transitions)
        @test foo == Transition(DateTime(2020, 1, 2, 4), FixedTimeZone("FOO", 0))
        @test bar == Transition(DateTime(2020, 1, 1, 3), FixedTimeZone("BAR", 0))
        @test i == 19

        # Times beyond 24-hours
        # - Linux clamps hours when beyond 24.
        # - macOS falls back to UTC hours beyond 167.
        @test _parsesub_tz("FOO+0BAR+0,0/25,1/720") == ParseNextError("Hours outside of expected range [0, 24]", "FOO+0BAR+0,0/25,1/720", 14, 15)

        # - Linux parses +/- times succesfully but clamps negative times to zero.
        # - macOS fails to parse +/- times.
        @test _parsesub_tz("FOO+0BAR+0,0/+1,1/1") == ParseNextError("Daylight saving start time should not have a sign", "FOO+0BAR+0,0/+1,1/1", 14)
        @test _parsesub_tz("FOO+0BAR+0,0/1,1/+1") == ParseNextError("Daylight saving end time should not have a sign", "FOO+0BAR+0,0/1,1/+1", 18)

        @test _parsesub_tz("FOO+0BAR+0,0/-1,1/1") == ParseNextError("Daylight saving start time should not have a sign", "FOO+0BAR+0,0/-1,1/1", 14)
        @test _parsesub_tz("FOO+0BAR+0,0/1,1/-1") == ParseNextError("Daylight saving end time should not have a sign", "FOO+0BAR+0,0/1,1/-1", 18)
    end

    # Example found in the `tzset 3` man page
    @testset "New Zealand example" begin
        tz, i = _parsesub_tz("NZST-12:00:00NZDT-13:00:00,M10.1.0,M3.3.0")
        nzst, nzdt = filter(t -> year(t.utc_datetime) == 2020, tz.transitions)
        @test nzst == Transition(DateTime(2020, 3, 14, 13), FixedTimeZone("NZST", 43200))
        @test nzdt == Transition(DateTime(2020, 10, 3, 14), FixedTimeZone("NZDT", 43200, 3600))
        @test i == 42
    end

    # Daylight saving offset is optional
    @testset "EST/EDT" begin
        tz, i = _parsesub_tz("EST5EDT,M3.2.0,M11.1.0")
        edt, est = filter(t -> year(t.utc_datetime) == 2020, tz.transitions)
        @test est == Transition(DateTime(2020, 11, 1, 6), FixedTimeZone("EST", -18000))
        @test edt == Transition(DateTime(2020, 3, 8, 7), FixedTimeZone("EDT", -18000, 3600))
        @test i == 23
    end

    # Validate the direct specification by comparing the result with a time zone computed
    # from tzdata.
    @testset "equivalent" begin
        # The tzdata time zone only uses consistent rules on/after 2007
        wpg = first(compile("America/Winnipeg", tzdata["northamerica"]))
        consistent_years = t -> year(t.utc_datetime) >= 2007

        tz, i = _parsesub_tz("CST+6CDT+5,M3.2.0/2,M11.1.0/2")
        @test tz.name == "CST/CDT"
        @test tz.name != wpg.name
        @test filter(consistent_years, tz.transitions) == filter(consistent_years, wpg.transitions)
        @test tz.cutoff == wpg.cutoff
        @test i == 30

        # Shorthand versions of the direct time zone specification we created
        expected = tz
        @test first(_parsesub_tz("CST+6CDT+5,M3.2.0,M11.1.0")) == expected
        @test first(_parsesub_tz("CST+6CDT+5")) == expected
    end
end


text:
__precompile__(true)

module Magic

export @μ, @ᛗ

include("mu.jl")
include("futhark.jl")

using .μodule
using .ᛗodule

end # module


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

vjp_algorithm(d::AbstractDensity) = ForwardDiffAD()
vjp_algorithm(d::AbstractPosteriorDensity) = vjp_algorithm(getlikelihood(d))
vjp_algorithm(d::DensityWithShape) = vjp_algorithm(parent(d))
vjp_algorithm(d::TransformedDensity) = vjp_algorithm(parent(d))
vjp_algorithm(d::RenormalizedDensity) = vjp_algorithm(parent(d))

vjp_algorithm(f::Union{LogDensityOf,NegLogDensityOf}) = vjp_algorithm(f.density)

jvp_algorithm(d::AbstractDensity) = ForwardDiffAD()
jvp_algorithm(d::AbstractPosteriorDensity) = jvp_algorithm(getprior(d))
jvp_algorithm(d::DensityWithShape) = jvp_algorithm(parent(d))
jvp_algorithm(d::TransformedDensity) = jvp_algorithm(parent(d))
jvp_algorithm(d::RenormalizedDensity) = jvp_algorithm(parent(d))

jvp_algorithm(f::Union{LogDensityOf,NegLogDensityOf}) = jvp_algorithm(f.density)

struct DensityWithDiff{JVP<:DifferentiationAlgorithm,D<:AbstractDensity,VJP<:DifferentiationAlgorithm} <: AbstractDensity
    vjpalg::VJP
    density::D
    jvpalg::JVP
end

@inline Base.parent(density::DensityWithDiff) = density.density

vjp_algorithm(density::DensityWithDiff) = density.vjpalg
jvp_algorithm(density::DensityWithDiff) = density.jvpalg

Base.:(|)(vjpalg::DifferentiationAlgorithm, density::AbstractDensity) = DensityWithDiff(vjpalg, density, jvp_algorithm(density))
Base.:(|)(vjpalg::DifferentiationAlgorithm, density::DensityWithDiff) = DensityWithDiff(vjpalg, parent(density), jvp_algorithm(density))

Base.:(|)(density::AbstractDensity, jvpalg::DifferentiationAlgorithm) = DensityWithDiff(vjp_algorithm(density), density, jvpalg)
Base.:(|)(density::DensityWithDiff, jvpalg::DifferentiationAlgorithm) = DensityWithDiff(vjp_algorithm(density), parent(density), jvpalg)

function Base.show(io::IO, density::DensityWithDiff)
    vjpalg = vjp_algorithm(density)
    jvpalg = jvp_algorithm(density)
    parent_vjpalg = vjp_algorithm(parent(density))
    parent_jvpalg = jvp_algorithm(parent(density))
    if vjpalg != parent_vjpalg || jvpalg == parent_jvpalg
        Base.show(io, vjpalg)
        print(io, "|")
    end
    Base.show(io, parent(density))
    if jvpalg != parent_jvpalg
        print(io, "|")
        Base.show(io, jvpalg)
    end
end

var_bounds(density::DensityWithDiff) = var_bounds(density.density)

ValueShapes.varshape(density::DensityWithDiff) = varshape(density.density)

ValueShapes.unshaped(density::DensityWithDiff) = DensityWithDiff(density.vjpalg, unshaped(density.density), density.jvpalg)

@inline function eval_logval_unchecked(density::DensityWithDiff, v::Any)
    eval_logval_unchecked(density.density, v)
end

@inline function eval_logval(density::DensityWithDiff, v::Any, T::Type{<:Real})
    eval_logval(density.density, v, T)
end

@inline get_initsrc_from_target(target::DensityWithDiff) = get_initsrc_from_target(target.density)

function bat_transform_impl(target::NoDensityTransform, density::DensityWithDiff, algorithm::DensityIdentityTransform)
    (result = density, trafo = IdentityVT(varshape(density)))
end

function bat_transform_impl(target::AbstractDensityTransformTarget, density::DensityWithDiff, algorithm::TransformAlgorithm)
    transformed_density, trafo = bat_transform_impl(target, parent(density), algorithm)
    (result = DensityWithDiff(vjp_algorithm(density), transformed_density, jvp_algorithm(density)), trafo = trafo)
end


text:
struct UncertaintySetConstraint <: RobustConstraintType
    name::String
end

struct UncertaintySetConst{S} <: RobustConstraint{S}
    name::String
    func::GenericAffExpr
    set::S
end

export UncertaintySetConstraint, UncertaintySetConst


text:
module HighLevelTypes

export @hl, @concretify, tuplejoin

const _hl_types = Dict{Symbol, Tuple}()

@inline tuplejoin(x) = (x,)
@inline tuplejoin(x::Tuple, y::Tuple) = (x..., y...)
@inline tuplejoin(x::Tuple, y) = (x..., y)
@inline tuplejoin(x, y::Tuple) = (x, y...)
@inline tuplejoin(x, y) = (x,y)
@inline tuplejoin(x::Tuple, y, z...) = (x..., tuplejoin(y, z...)...)
@inline tuplejoin(x, y, z...) = (x, tuplejoin(y, z...)...)

function maketypesconcrete(expression)
    for (i,arg) in enumerate(expression.args)
        if isa(arg, Symbol)
            if haskey(_hl_types, arg)
                expression.args[i] = Symbol("_", arg)
            end
        elseif :args in fieldnames(typeof(arg))
            maketypesconcrete(arg)
        end
    end
end

#creates a highlevel type: uses abstract types of attribute fields
macro hl(typeexpr::Expr)
    mutable, nameblock, args = typeexpr.args
    createhltype(mutable, nameblock, args)
end

#uses concrete types for all highlevel types
macro concretify(expr::Expr)
    maketypesconcrete(expr)
    esc(expr)
end

function createhltype(mutable, nameblock, args)
    if !isa(nameblock, Symbol) && nameblock.head == :<:
        parentname = nameblock.args[2]
        nameblock = nameblock.args[1]
    else
        parentname = :NOPARENT
    end

    if isa(nameblock, Symbol)
        name, params = nameblock, []
    elseif nameblock.head == :curly
        name, params = nameblock.args[1], nameblock.args[2:end]
    end

    if parentname == :NOPARENT
        _hl_types[name] = params, args
        acode = :(abstract type $name end)
    else
        _hl_types[name] = vcat(_hl_types[parentname][1], params), :( $(_hl_types[parentname][2].args...) ; $(args.args...))
        acode = :(abstract type $name <: $parentname end)
    end

    tname = Symbol(:_, name)
    tname = :($tname{})
    append!(tname.args, _hl_types[name][1])
    constrcode = quote
        $(Symbol(name,"Builder"))(args...) = args
        $name(args...) = $(Symbol("_",name))($(Symbol(name,"Builder"))(args...)...)
    end

    if mutable
        code = quote
            $acode
            mutable struct $tname <: $name
                $(_hl_types[name][2].args...)
            end
            $constrcode
        end
    else
        code = quote
            $acode
            struct $tname <: $name
                $(_hl_types[name][2].args...)
            end
            $constrcode
        end
    end
    esc(code)
end

end # module


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

V,EV = Lar.randomcuboids(100, .2)
V = Plasm.normalize(V,flag=true)
model2d = V,EV

Sigma =  Lar.spaceindex(model2d);
for k=1:length(Sigma) println(k,Sigma[k]) end

Plasm.view(Plasm.numbering(.15)((V,[[[k] for k=1:size(V,2)], EV])))
Plasm.viewexploded(V,EV)(1.2,1.2,1.2) 	# no numerical errors

W,EW = Lar.fragmentlines((V,EV))
Plasm.viewexploded(W,EW)(1.2,1.2,1.2)
Plasm.view(Plasm.numbering(.015)((W,[[[k] for k=1:size(W,2)], EW])))

V,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells)) # 2-connected components (H & T)
hpcs = [ Plasm.lar2hpc(V,EVs[i]) for i=1:length(EVs) ]
Plasm.view([ Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k]) for k=1:(length(hpcs)) ])

W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EW::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)

triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
W = convert(Lar.Points, V')
Plasm.viewcolor(W::Lar.Points, FVs::Array{Lar.Cells})


text:
include(joinpath("src", "simulation.jl"))
include(joinpath("src", "03runeval.jl"))

####    Transfering Simulation Structs into new format  ####

# If simulation struct changes, this procedure helps to convert old simulation
# objects into new ones.

for file in readdir("results")
    if occursin("transferformat", file) || !occursin(".jld2", file)
        continue
    end
    raw = load(joinpath("results", file))
    current_run = raw[first(keys(raw))]

    current_run_transfer = Tuple{
        String,
        Int64,
        Int64,
        MersenneTwister,
        Config,
        Any,
        Any,
        DataFrame,
        Any,
        Array{AbstractGraph}
    }[]

    for i in 1:length(current_run)
        push!(
            current_run_transfer,
            (
                current_run[i].name,
                current_run[i].runnr,
                current_run[i].repnr,
                current_run[i].rng,
                current_run[i].config,
                current_run[i].init_state,
                current_run[i].final_state,
                current_run[i].agent_log,
                current_run[i].post_log,
                current_run[i].graph_list
            )
        )
    end

    save(
        joinpath("results", current_run[1].name * "_transferformat.jld2"),
        current_run[1].name, current_run_transfer)
end

# Transfer Step 2

for file in readdir("results")
    if !occursin("transferformat", file) || !occursin(".jld2", file)
        continue
    end
    raw = load(joinpath("results", file))
    current_run = raw[first(keys(raw))]

    current_run_transfered = [Simulation(current_run[i]) for i in 1:length(current_run)]

    save(joinpath("results", current_run_transfered[1].name * "_transfered.jld2"), current_run_transfered[1].name, current_run_transfered)
end

####    Appending replications of equal batch simulation runs  ####

# If a batch run shall be continued, this function helps with appending the
# new results to the existing ones. When continuing a batch run on the same
# machine, it will automatically continue with next replication

# To append new results, the existing results should be inside the results
# folder and the new ones in the resultst subfolder "append"

for file in readdir("results")
    if occursin(".jld2", file) && !occursin("appended", file)
        raw = load(joinpath("results", file))
        current_run = raw[first(keys(raw))]
        name = file[1:first(findfirst(".jld2", file))-1]

        raw2 = load(joinpath("results", "append", file))
        appending_run = raw2[first(keys(raw2))]

        if last(current_run).repnr + 1 == first(appending_run).repnr
            append!(current_run, appending_run)
            save(joinpath("results", name * ".jld2"), name, current_run)
        else
            println("Appending runs for $name do not match")
        end

    end
end


text:
function to_rgba(x::UInt32)
    a = ((x & 0xff000000)>>24)/255
    b = ((x & 0x00ff0000)>>16)/255
    g = ((x & 0x0000ff00)>>8)/255
    r = (x & 0x000000ff)/255
    RGBA(r, g, b, a)
end

# plotting
function triplot(x,y,z,t)
    img = to_rgba.(Triplot.rasterize(x,y,z,t)')
    plot([extrema(x)...], [extrema(y)...], img, yflip=false)
end

function uniform_tri_nodes(N)
    r1D = LinRange(-1,1,N+1)
    Np = (N+1)*(N+2)÷2
    r,s = zeros(Np),zeros(Np)
    sk = 1
    for i = 0:N
        for j = 0:N-i
            r[sk] = r1D[i+1]
            s[sk] = r1D[j+1]
            sk += 1
        end
    end
    return r,s
end

function plotTriMesh(mesh)
    x,y = mesh.point[1,:], mesh.point[2,:]
    xmesh = Float64[]
    ymesh = Float64[]
    for vertex_ids in eachcol(mesh.cell)
        ids = vcat(vertex_ids, vertex_ids[1])
        append!(xmesh,[x[ids];NaN])
        append!(ymesh,[y[ids];NaN])
    end
    display(plot(xmesh,ymesh,linecolor=:black,legend=false,ratio=1))
end

# returns vertex coordinates and element-to-vertex connectivities
function unpack_mesh_info(mesh)
    VX,VY = mesh.point[1,:],mesh.point[2,:]
    EToV  = mesh.cell  # element-to-vertex mapping: each column contains vertex ids for a different element
    return VX,VY,EToV
end

function map_triangle_pts(r,s,x,y)
    return λ(r,s)*x, λ(r,s)*y
end

# get list of boundary indices (vector of indices) and boundary faces (stored as list of tuples (f=face,e=element))
function get_boundary_info(reference_face_indices,mesh)
    is_point_on_boundary = mesh.point_marker # 1 if on boundary, 0 otherwise
    is_boundary_face = find_boundary_faces(reference_face_indices,mesh)
    boundary_indices = findall(vec(is_point_on_boundary) .== 1)
    boundary_faces = Tuple.(findall(is_boundary_face .== 1))
    return boundary_indices,boundary_faces
end

# loop through faces and look for matches
function find_boundary_faces(reference_face_indices,mesh)
    num_elements = size(mesh.cell,2)
    is_boundary_face = zeros(Int,3,num_elements)
    for e = 1:num_elements
        for f = 1:3
            face_vertices = mesh.cell[reference_face_indices[f],e]
            if sort(face_vertices) in sort.(eachcol(mesh.segment))
                 is_boundary_face[f,e] = 1
            end
        end
    end
    return is_boundary_face
end

"""
meshgrid(vx), meshgrid(vx,vy)
Computes an (x,y)-grid from the vectors (vx,vx). For more information, see the MATLAB documentation.

Copied and pasted directly from [VectorizedRoutines.jl](https://github.com/ChrisRackauckas/VectorizedRoutines.jl/blob/master/src/matlab.jl).
Using VectorizedRoutines.jl directly causes Pkg versioning issues with SpecialFunctions.jl
"""
meshgrid(v::AbstractVector) = meshgrid(v, v)
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repeat(vx, m, 1), repeat(vy, 1, n))
end

function uniform_tri_mesh(Kx,Ky)
        (VY, VX) = meshgrid(LinRange(-1,1,Ky+1),LinRange(-1,1,Kx+1))
        sk = 1
        EToV = zeros(Int,2*Kx*Ky,3)
        for ey = 1:Ky
                for ex = 1:Kx
                        id(ex,ey) = ex + (ey-1)*(Kx+1) # index function
                        id1 = id(ex,ey)
                        id2 = id(ex+1,ey)
                        id3 = id(ex+1,ey+1)
                        id4 = id(ex,ey+1)
                        EToV[2*sk-1,:] = [id1 id3 id2]
                        EToV[2*sk,:]   = [id3 id1 id4]
                        sk += 1
                end
        end
        return VX[:],VY[:],Matrix(EToV')
end


text:
using DataFrames

function file_data_transform(in_file::AbstractString, out_file::AbstractString, transform::Function, pars::Dict{Any, Any}; element_transform::Function=x->x, header=true)
  in_data = load_data(in_file; header=header)
  out_data = transform(in_data, pars; transform=element_transform)
  save_data(out_data, out_file; header=header)
end

function load_matrix(file_name::AbstractString, N::Int, L::Int)
  if !isfile(file_name)
    return Nullable{Array{Float64, 2}}()
  end
  M = if N > 0 && L >0
      readdlm(file_name, ' ';
        header=false, skipstart=0, skipblanks=true,
        use_mmap=false, quotes=false, dims=(N, L),
        comments=true, comment_char='#')
    else
      readdlm(file_name, ' ';
        header=false, skipstart=0, skipblanks=true,
        use_mmap=false, quotes=false,
        comments=true, comment_char='#')
    end
  return Nullable{Array{Float64, 2}}(M)
end

function load_matrix(file_name::AbstractString, N::Int)
  return load_matrix(file_name, N, N)
end

function load_matrix(file_name::AbstractString)
  return load_matrix(file_name, -1, -1)
end


text:
include(joinpath(@__DIR__, "..", "share", "LogSystem.jl"))
include(joinpath(@__DIR__, "..", "models", "EMOM", "ENGINE_EMOM.jl"))
include(joinpath(@__DIR__, "..", "driver", "driver_working.jl"))
include(joinpath(@__DIR__, "ProgramTunnel", "src", "julia", "BinaryIO.jl"))
include(joinpath(@__DIR__, "ProgramTunnel", "src", "julia", "ProgramTunnel_fs_new.jl"))

using MPI
using CFTime, Dates
using ArgParse
using JSON
using TOML

using .ProgramTunnel_fs
using .BinaryIO
using .LogSystem

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin
        "--config-file"
            help = "Configuration file in TOML format."
            arg_type = String
            required = true
    end

    return parse_args(s)
end

function parseCESMTIME(
        t_str    :: AbstractString,
        timetype :: DataType,
)

    yyyy = parse(Int64, t_str[1:4])
    mm   = parse(Int64, t_str[5:6])
    dd   = parse(Int64, t_str[7:8])
    HH   = parse(Int64, t_str[10:11])
    MM   = parse(Int64, t_str[13:14])
    SS   = parse(Int64, t_str[16:17])

    return timetype(yyyy,mm,dd,HH,MM,SS)

end

function parseMsg(msg::AbstractString)
    pairs = split(msg, ";")
    d = Dict{AbstractString, Any}()
    for i = 1:length(pairs)

        if strip(pairs[i]) == ""
            continue
        end

        key, val = split(pairs[i], ":")
        key = String(key)
        val = String(val)
        d[key] = val
    end
    return d
end

parsed = parse_commandline()

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
is_master = rank == 0

msg = nothing
config = nothing
if is_master

    config = TOML.parsefile(parsed["config-file"])

    timetype = getproperty(CFTime, Symbol(config["MODEL_MISC"]["timetype"]))

    PTI = ProgramTunnelInfo(
        path = joinpath(config["DRIVER"]["caserun"], "x_tmp"),
        reverse_role  = true,
        recv_channels = 2,
        chk_freq = 0.2,
    )
    nullbin  = [zeros(Float64, 1)]

    function recvMsg(;verbose::Bool = true)
        global msg = parseMsg( recvData!(PTI, nullbin, which=1) )
        if verbose
            writeLog("Receive message: ")
            JSON.print(msg, 4)
        end
    end

end

coupler_funcs = (

    master_before_model_init = function()

        global msg

        writeLog("[Coupler] Before model init. My rank = {:d}", rank)

        recvMsg()

        if msg["MSG"] != "INIT"
            throw(ErrorException("Unexpected `MSG` : " * string(msg["MSG"])))
        end

        read_restart = (msg["READ_RESTART"] == "TRUE") ? true : false
        cesm_coupler_time = parseCESMTIME(msg["CESMTIME"], timetype)

        return read_restart, cesm_coupler_time

    end,

    master_after_model_init! = function(OMMODULE, OMDATA)

        global msg

        writeLog("[Coupler] After model init. My rank = {:d}", rank)

        global lsize = parse(Int64, msg["LSIZE"])

        global send_data_list = [OMDATA.o2x["SST"], OMDATA.o2x["Q_FRZMLTPOT"]]
        global recv_data_list = []

        global x2o_available_varnames = split(msg["VAR2D"], ",")
        global x2o_wanted_varnames = keys(OMDATA.x2o)
        global x2o_wanted_flag     = [(x2o_available_varnames[i] in x2o_wanted_varnames) for i = 1:length(x2o_available_varnames)]

        writeLog("# Check if all the requested variables are provided by cesm end...")
        local pass = true
        for varname in x2o_wanted_varnames
            if ! ( varname in x2o_available_varnames )
                writeLog("Error: $(varname) is requested by ocean model but not provided on cesm side.")
                pass = false
            end
        end

        if pass
            writeLog("All variables requested are provided.")
        else
            throw(ErrorException("Some variable are not provided. Please check."))
        end

        writeLog("# List of provided x2o variables:")
        for (i, varname) in enumerate(x2o_available_varnames)
            push!(recv_data_list , ( x2o_wanted_flag[i] ) ? OMDATA.x2o[varname] :  zeros(Float64, lsize))
            println(format(" ({:d}) {:s} => {:s}", i, varname, ( x2o_wanted_flag[i] ) ? "Wanted" : "Dropped" ))
        end

        sendData(PTI, "OK", send_data_list)
        sendData(PTI, "mask",  [OMDATA.o2x["mask"],])

    end,

    master_before_model_run! = function(OMMODULE, OMDATA)
        #writeLog("[Coupler] Before model run")
        #writeLog("[Coupler] This is where flux exchange happens.")
        recvMsg()

        return_values = nothing
        if msg["MSG"] == "RUN"
            Δt = Dates.Second(parse(Float64, msg["DT"]))
            recvData!(
                PTI,
                recv_data_list,
                which=2,
            )

            cesm_coupler_time = parseCESMTIME(msg["CESMTIME"], timetype)
            if OMDATA.clock.time != cesm_coupler_time
                writeLog("Warning: time inconsistent. `cesm_coupler_time` is $(string(cesm_coupler_time)), but ocean model's time is $(string(OMDATA.clock.time)). This can happen if this is a startup run or hybrid run. Because my implementation cannot tell, I have to just forward the model time. Please be extra cautious about this.")

                error_t = Second(cesm_coupler_time - OMDATA.clock.time)
                if error_t.value < 0
                    throw(ErrorException("Error: ocean model time is ahead of `cesm_coupler_time`. This is absolutely an error."))
                end
                advanceClock!(OMDATA.clock, error_t)

            end

            write_restart = msg["WRITE_RESTART"] == "TRUE"

            return_values = ( :RUN,  Δt, write_restart )

        elseif msg["MSG"] == "END"
            return_values = ( :END, 0.0, false  )
        else
            throw(ErrorException("Unexpected `MSG` : " * string(msg["MSG"])))
        end

        return return_values

    end,

    master_after_model_run! = function(OMMODULE, OMDATA)
        #writeLog("[Coupler] After model run")
        global send_data_list = [OMDATA.o2x["SST"], OMDATA.o2x["Q_FRZMLTPOT"]]
        sendData(PTI, "OK", send_data_list)
    end,

    master_finalize! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] Finalize")
        writeLog("Sleep for 30 seconds before archiving to avoid conflicting with CESM archiving process.")
        sleep(30.0)
    end
)

runModel(
    ENGINE_EMOM,
    coupler_funcs,
    config,
)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 896f6942-76a2-11eb-01e9-1156629a4a14
begin
	using Pkg;using Plots; gr()
	Pkg.add("LaTeXStrings")
	using LaTeXStrings
end

# ╔═╡ 0eea2678-7641-11eb-1c89-b53ca4514973
function circleShape(h,k,r)
	θ = LinRange(0,2π,500)
	h .+ r*sin.(θ), k .+ r*cos.(θ)
end

# ╔═╡ 58cc6b66-7641-11eb-2aa0-3b5f7a52eac9
function mohr2D(ϵx,ϵy,ϵxy)
	center = (ϵx+ϵy)/2.
	R = sqrt(((ϵx-ϵy)/2.)^2+ϵxy^2)
	return center,R
end

# ╔═╡ a8d106fa-7641-11eb-0248-1ddac2014484
function mohr2D(ϵ1,ϵ2)
	center = (ϵ1+ϵ2)/2
	R = abs((ϵ1-ϵ2)/2)
	return center,R
end

# ╔═╡ a62626fc-7643-11eb-14ef-c110a134ce50
function plotCircle(ϵ1,ϵ2)
	c,R =  mohr2D(ϵ1,ϵ2)
	p = plot(circleShape(c,0,R),seriestype = [:shape], 		 lw=0.5, c = :blue, linecolor =:black, legend = false, fillalpha = 0.05, aspect_ratio = 1)
	scatter!(p,[c],[0],markercolor = :black)
	scatter!(p,[ϵ1],[0],markercolor = :black)
	scatter!(p,[ϵ2],[0],markercolor = :black)
	xlabel!(L"\varepsilon_n")

	ylabel!(L"\tau_n")

	if R≉0
		annotate!(p, ϵ1 + 0.1*R, 0.1*R, L"\epsilon_1", :color)
		annotate!(p, ϵ2 + 0.1*R, 0.1*R, L"\epsilon_2", :color)
	end

	if R ≈ 0
		ylims!(-maximum([ϵ1,ϵ2]),maximum([ϵ1,ϵ2]))
	end

	hline!(p,[0], linecolor = :black)
	return p
end

# ╔═╡ eaf4f830-7643-11eb-09e3-e1109e9cdb65
function plotCircle(ϵx,ϵy,ϵxy)
	c,R =  mohr2D(ϵx,ϵy,ϵxy)
	ϵ1 = c + R
	ϵ2 = c - R
	p = plot(circleShape(c,0,R),seriestype = [:shape], 		 lw=0.5, c = :blue, linecolor =:black, legend = false, fillalpha = 0.05, aspect_ratio = 1)
	scatter!(p,[c],[0],markercolor = :black)
	scatter!(p,[ϵ1],[0],markercolor = :black)
	scatter!(p,[ϵ2],[0],markercolor = :black)
	scatter!(p,[ϵx],[ϵxy],markercolor = :blue)
	scatter!(p,[ϵy],[-ϵxy],markercolor = :red)

	plot!(p,[ϵx,ϵy],[ϵxy, -ϵxy],linecolor=:black)

	if R≉0
		annotate!(p, ϵ1 + 0.1*R, 0.1*R, L"\epsilon_1", :color)
		annotate!(p, ϵ2 + 0.1*R, 0.1*R, L"\epsilon_2", :color)
		annotate!(p, ϵx+0.06*R, ϵxy+0.1*R, L"X", :color)
		annotate!(p, ϵy+0.06*R, -ϵxy+0.1*R, L"Y", :color)
	end

	xlabel!(L"\varepsilon_n")

	ylabel!(L"\tau_n")
	if R ≈ 0
		ylims!(-maximum([ϵ1,ϵ2]),maximum([ϵ1,ϵ2]))
		annotate!(p, ϵ1 + 0.1*ϵ1, 0.1*ϵ1, L"\epsilon_1 = \epsilon_2 = C", :color)
	end

	hline!(p,[0], linecolor = :black)
	return p
end

# ╔═╡ d56c2fda-7641-11eb-29e8-579e1f95cada
begin
	ϵ1 = 1.
	ϵ2 = 5.
	p = plotCircle(3,5,2)
end

# ╔═╡ Cell order:
# ╠═d56c2fda-7641-11eb-29e8-579e1f95cada
# ╠═896f6942-76a2-11eb-01e9-1156629a4a14
# ╠═0eea2678-7641-11eb-1c89-b53ca4514973
# ╠═58cc6b66-7641-11eb-2aa0-3b5f7a52eac9
# ╠═a8d106fa-7641-11eb-0248-1ddac2014484
# ╠═a62626fc-7643-11eb-14ef-c110a134ce50
# ╠═eaf4f830-7643-11eb-09e3-e1109e9cdb65

max_stars_repo_name: victorjesusamoresmedianero/Prueba

max_issues_repo_name: victorjesusamoresmedianero/Prueba

max_forks_repo_name: victorjesusamoresmedianero/Prueba


text:
module StaticTypeAnalysis

macro lintpragma(ex); end

EQ_METHOD_FALSE = which(==, Tuple{Nothing, Int})
#CONSTRUCTOR_FALLBACK = which(Nothing, Tuple{Nothing})

"""
    StaticTypeAnalysis.canequal(S::Type, T::Type) :: Union{Bool, Nothing}

Given types `S` and `T`, return `false` if it is not possible for
`s::S == t::T`. Return `true` if it is possible, and
`nothing` if it cannot be determined.

```jldoctest
julia> StaticTypeAnalysis.canequal(Int, Float64)
true

julia> StaticTypeAnalysis.canequal(Int, String)
false
```
"""
function canequal(S::Type, T::Type)
    if S == Union{} || T == Union{}
        false
    elseif typeintersect(S, T) ≠ Union{}
        # TODO: this is not fully correct; some types are not Union{} but still
        # not instantiated
        true
    elseif isconcretetype(S) && isconcretetype(T) &&
           EQ_METHOD_FALSE == which(==, Tuple{S, T})
        # == falls back to === here, but we saw earlier that the intersection
        # is empty
        false
    else
        try zero(S) == zero(T)
            true
        catch
            nothing
        end
    end
end

"""
    StaticTypeAnalysis.isknownerror(f, argtypes)

Return `true` if it is known, without inference, that calling `f` with
signature `argtypes` will result in an error.

Consumers of this package are advised to use `infertype` and check the result
against `Union{}`, which covers more cases.
"""
isknownerror(_f, _argtypes) = false
function isknownerror(::typeof(Base.getindex), argtypes::Tuple)
    # correctly forward arguments from `infertype` below
    isknownerror(typeof(Base.getindex), argtypes[1])
end
function isknownerror(::typeof(Base.getindex), argtypes)
    if isempty(argtypes)
        true
    elseif argtypes[1] <: AbstractDict
        if Base.length(argtypes) ≠ 2
            true
        else
            try
                K = keytype(argtypes[1])
                ce = canequal(K, argtypes[2])
                ce !== nothing && !ce
            catch
                false
            end
        end
    else
        false
    end
end

"""
    StaticTypeAnalysis.infertype(f, argtypes)

Given a function `f` and a list of types `argtypes`, use inference and other
static type checking techniques to figure out a type `S` such that the result
of applying `f` to `argtypes` is always of type `S`.
"""
infertype(f, argtypes...) = infertype(f, argtypes)
function infertype(f, argtypes::Tuple)
    if isknownerror(f, argtypes)
        Union{}
    elseif f === Base.getindex && Base.length(argtypes) == 2 &&
           argtypes[1] == UnitRange && argtypes[2] <: UnitRange &&
           eltype(argtypes[2]) <: Integer
        # TODO: would be nice to get rid of this odd special case
        UnitRange
    # elseif isa(f, Type) && Base.length(argtypes) == 1 &&
    #        isconcretetype(argtypes[1]) &&
    #        which(f, Tuple{argtypes[1]}) === CONSTRUCTOR_FALLBACK
    #     # we can infer better code for the constructor `convert` fallback by
    #     # inferring the convert itself
    #     Core.Inference.return_type(convert, Tuple{Type{f}, argtypes[1]})
    else
        try
            typejoin(Base.return_types(f, Tuple{argtypes...})...)
        catch  # error might be thrown if generic function, try using inference
            if all(isconcretetype, argtypes)
                Core.Inference.return_type(f, Tuple{argtypes...})
            else
                Any
            end
        end
    end
end

"""
    StaticTypeAnalysis.getindexable(T::Type)

Return `true` if all objects of type `T` support `getindex`, and the `getindex`
operation on numbers is consistent with iteration order.

Note that, in particular, this is not true for `String` and `Dict`.
"""
getindexable(::Type{T}) where {T <: Union{Tuple,Pair,Array,Number}} = true
getindexable(::Type) = false

"""
    StaticTypeAnalysis.length(T::Type) :: Union{Int, Nothing}

If it can be determined that all objects of type `T` have length `n`, then
return `n`. Otherwise, return `nothing`.
"""
length(::Type{Union{}}) = 0
length(::Type) = nothing
length(::Type{T}) where {T <: Pair} = 2
length(::Type{T}) where T <: NTuple{N, Any} where N = N

"""
    StaticTypeAnalysis.eltype(T::Type)

Return `S` as specific as possible such that all objects of type `T` have
element type `S`.
"""
eltype(::Type{Union{}}) = Union{}
eltype(T::Type) = Base.eltype(T)

_getindex_nth(xs::Any, ::Type{Val{n}}) where {n} = xs[n]
_typeof_nth_getindex(::Type{T}, n::Integer) where {T} =
    infertype(_getindex_nth, Any[T, Type{Val{Int(n)}}])

"""
    StaticTypeAnalysis.typeof_nth(T::Type)

Return `S` as specific as possible such that all objects of type `T`, when
iterated over, have `n`th element type `S`.
"""
typeof_nth(T::Type, n::Integer) =
    if getindexable(T) && 0 < Base.length(T.types)
        if n ≤ Base.length(T.types)
            typeintersect(eltype(T), T.types[n])
        else
            Union{}
        end
    else
        eltype(T)
    end
typeof_nth(::Type{Pair{K,V}}, n::Integer) where {K, V} =
    n == 1 ? K : n == 2 ? V : Union{}
typeof_nth(::Type{Union{}}, ::Integer) = Union{}
typeof_nth(::Type{Tuple{Vararg{Any}}}, ::Integer) = Any

end


text:
"TODO top level doc"
module Multilevel

include("./interface.jl")
include("./cycles.jl")
include("./vcycle.jl")
include("./fcycle.jl")
include("./wcycle.jl")

end # module


text:
function ndgrid{T}(v1::AbstractVector{T}, v2::AbstractVector{T})
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)
    (repmat(v1, 1, n), repmat(v2, m, 1))
end

function pe1()
    N = 1e6
    sum = 0
    for i=1:N
        if i % 3 == 0 || i % 5 == 0
            sum += i
        end
    end
end

function pe10()
    N = 2e6
    primes = [1:N]
    for i in 2:sqrt(N)
        j = (2:N/i)*i
        primes[j] = 0
    end
end

function pe73()
    N = 1e3
    di = [1:N]
    ni = [1:N]
    (n, d) = ndgrid(ni, di)

    f = n ./ d

    f = unique(f)
    f = sort(f)
    f = f[f .< 1/2]
    f = f[f .> 1/3]
end

function soft_threshold()
    N = 1e2
    j = linspace(-1.0, 1.0, convert(Int64, N))
    x, y = ndgrid(j, j)
    z = x.^2 + y.^2
    i = abs(z) .< 0.5
    z[i] = 0
    z[!i] -= 0.5
end

function pi_approx()
    N = 1e6
    pie = 1.0
    for k in 1:N
        pie += (-1)^k / (2*k + 1)
    end
end

# timing
pe1()
pe10()
pe73()
soft_threshold()
pi_approx()

print("pe1: ")
println(@elapsed pe1())
print("pe10: ")
println(@elapsed pe10())
print("pe73: ")
println(@elapsed pe73())
print("st_time: ")
println(@elapsed soft_threshold())
print("pi_approx: ")
println(@elapsed pi_approx())

max_stars_repo_name: stsievert/swix

max_issues_repo_name: stsievert/swix

max_forks_repo_name: scottsievert/swix


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

function get_test_2d_model()
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 1.0],
        2 => [3/4, 1.0],
        3 => [2.0, 1.0],
        4 => [0.0, 1.0],
        5 => [5/4, 1.0],
        6 => [2.0, 1.0])
    sel1 = Element(Seg2, [1, 2])
    sel2 = Element(Seg2, [2, 3])
    mel1 = Element(Seg2, [4, 5])
    mel2 = Element(Seg2, [5, 6])
    update!([mel1, mel2, sel1, sel2], "geometry", X)
    return [sel1, sel2], [mel1, mel2]
end

@testset "calculate flat 2d assembly" begin
    (sel1, sel2), (mel1, mel2) = get_test_2d_model()

    bc = Problem(Mortar, "test interface", 1, "temperature")
    update!([sel1, sel2], "master elements", [mel1, mel2])
    bc.elements = [sel1, sel2, mel1, mel2]

    B_expected = zeros(3, 6)

    S1 = get_gdofs(bc, sel1)
    M1 = get_gdofs(bc, mel1)
    B_expected[S1,S1] += [1/4 1/8; 1/8 1/4]
    B_expected[S1,M1] -= [3/10 3/40; 9/40 3/20]

    S2 = get_gdofs(bc, sel2)
    M2 = get_gdofs(bc, mel1)
    B_expected[S2,S2] += [49/150 11/150; 11/150 2/75]
    B_expected[S2,M2] -= [13/150 47/150; 1/75 13/150]

    S3 = get_gdofs(bc, sel2)
    M3 = get_gdofs(bc, mel2)
    B_expected[S3,S3] += [9/100 27/200; 27/200 39/100]
    B_expected[S3,M3] -= [3/20 3/40; 9/40 3/10]

    assemble!(bc, 0.0)
    B = full(bc.assembly.C1, 3, 6)
    # dump(round(B, 6))
    # dump(B_expected)
    @test isapprox(B, B_expected; rtol=1.0e-9)

end

@testset "solve mortar tie contact with multiple dirichlet boundary conditions and multiple bodies" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0], 2 => [1.0, 0.0],
        3 => [0.0, 0.5], 4 => [1.0, 0.5],
        5 => [0.0, 0.5], 6 => [1.0, 0.5],
        7 => [0.0, 1.0], 8 => [1.0, 1.0])
    T = Dict{Int64, Vector{Float64}}(
        7 => [0.0, 288.0], 8 => [0.0, 288.0]
    )
    e1 = Element(Quad4, [1, 2, 4, 3])
    e2 = Element(Quad4, [5, 6, 8, 7])
    t1 = Element(Seg2, [7, 8])
    update!([e1, e2, t1], "geometry", X)
    update!([e1, e2], "youngs modulus", 288.0)
    update!([e1, e2], "poissons ratio", 1/3)
    update!(t1, "displacement traction force", T)

    body1 = Problem(Elasticity, "block 1", 2)
    body1.properties.formulation = :plane_stress
    push!(body1, e1)
    body2 = Problem(Elasticity, "block 2", 2)
    body2.properties.formulation = :plane_stress
    push!(body2, e2, t1)

    # boundary elements for dirichlet dx=0
    dx1 = Element(Seg2, [1, 3])
    dx2 = Element(Seg2, [5, 7])
    update!([dx1, dx2], "geometry", X)
    update!([dx1, dx2], "displacement 1", 0.0)
    bc1 = Problem(Dirichlet, "symmetry dx=0", 2, "displacement")
    push!(bc1, dx1, dx2)

    # boundary elements for dirichlet dy=0
    dy1 = Element(Seg2, [1, 2])
    update!(dy1, "geometry", X)
    update!(dy1, "displacement 2", 0.0)
    bc2 = Problem(Dirichlet, "symmetry dy=0", 2, "displacement")
    push!(bc2, dy1)

    # mortar boundary between two bodies
    mel1 = Element(Seg2, [3, 4])
    sel1 = Element(Seg2, [5, 6])
    update!([mel1, sel1], "geometry", X)
    update!(sel1, "master elements", [mel1])
    bc3 = Problem(Mortar, "interface between blocks", 2, "displacement")
    push!(bc3, sel1, mel1)

    solver = LinearSolver(body1, body2, bc1, bc2, bc3)
    solver()

    u = e2("displacement", [1.0, 1.0], 0.0)
    u_expected = [-1/3, 1.0]
    info("displacement at tip: $u")
    @test isapprox(u, u_expected)
end


text:
# ---
# title: 193. Valid Phone Numbers
# id: problem193
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/valid-phone-numbers/description/>
# hidden: true
# ---
#
# Given a text file `file.txt` that contains list of phone numbers (one per
# line), write a one liner bash script to print all valid phone numbers.
#
# You may assume that a valid phone number must appear in one of the following
# two formats: (xxx) xxx-xxxx or xxx-xxx-xxxx. (x means a digit)
#
# You may also assume each line in the text file must not contain leading or
# trailing white spaces.
#
# **Example:**
#
# Assume that `file.txt` has the following content:
#
#
#
#     987-123-4567
#     123 456 7890
#     (123) 456-7890
#
#
# Your script should output the following valid phone numbers:
#
#
#
#     987-123-4567
#     (123) 456-7890
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end