text:
export getEdgeIntegralOfPolygonalChain

function getEdgeIntegralOfPolygonalChain(mesh::TetraMesh, polygon::Array{Float64,2}; normalize=false)

X     = mesh.Points
T     = mesh.Tetras
edges = mesh.edges
nn    = mesh.nn
ne    = mesh.ne
nc    = mesh.nc

# create connectivity matrix
E = sparse([edges[:,1]; edges[:,2]], [edges[:,2]; edges[:,1]], [edges[:,1]; edges[:,2]], nn, nn)

# number of polygon segments
np = size(polygon,1) - 1

# collect all vertices that discretize the polygon
v  = Int64[]

# locate mesh vertex closest to first point of polygon
X1 = polygon[1,:]
k1 = 0
dmin = Inf
for ii = 1:nn
  d = norm(X[ii,:] - X1)
  if d < dmin
    k1   = ii
    dmin = d
  end
end
X1 = X[k1,:]
push!(v, k1)

# loop over all polygon segments
for ip = 1:np

  # locate mesh vertex closest to end point of current polygon segment
  X2 = polygon[ip+1,:]
  k2 = 0
  dmin = Inf
  for ii = 1:nn
    d = norm(X[ii,:] - X2)
    if d < dmin
      k2   = ii
      dmin = d
    end
  end
  X2 = X[k2,:]

  # find edges connecting k1 (X1) and k2 (X2) by a straight line
  while k1 != k2

    X1 = X[k1,:]
    K  = nonzeros(E[:,k1])  # all vertices connected to k1 by an edge

    # find vertex closest to straight line between k1 (X1) and k2 (X2)
    dmin = Inf
    ll   = norm(X2 - X1)
    for k in K
      Xk = X[k,:]
      l1 = norm(Xk - X1)
      l2 = norm(Xk - X2)
      d  = (l1 + l2) / ll # d should be equal 1 if k (Xi) is on a straight line
      if d < dmin
        k1   = k
        dmin = d
      end
    end
    push!(v, k1)

  end

  # end point of current segment is start point of next segment
  X1 = X2

end

# edge lengths
L = getEdgeLength(mesh)
F = getFaceArea(mesh)
V = getCellVolume(mesh)

# find edge numbers for edges v[1:2], v[2:3], ...
n = length(v) - 1
e = sort(Int64[v[1:n] v[2:n+1]], 2) # sort because edges[:,1] < edges[:,2]
j = findRow(edges, e)

# assemble sparse source vector
s = spzeros(ne,1)
for ii = 1:n
  k = j[ii]
  if v[ii+1] - v[ii] > 0   # line current flows in direction of edge
    s[k] =  L[k]
  else                     # line current flows opposite to direction of edge
    s[k] = -L[k]
  end
end

# normalize
if normalize

	if all(polygon[1,:] .== polygon[np+1,:])

		# closed polygon: divide by enclosed area
		a  = 0.0
		px = polygon[2:np,1] .- polygon[1,1]
		py = polygon[2:np,2] .- polygon[1,2]
		pz = polygon[2:np,3] .- polygon[1,3]
		for ip = 1:np-2
			cx = py[ip] * pz[ip+1] - pz[ip] * py[ip+1]
			cy = pz[ip] * px[ip+1] - px[ip] * pz[ip+1]
			cz = px[ip] * py[ip+1] - py[ip] * px[ip+1]
			a += sqrt(cx * cx + cy * cy + cz * cz)
		end
		a *= 0.5

	else

		# open polygon: divide by length
		a = 0.0
		for ip = 1:np
			dx = polygon[ip+1,1] - polygon[ip,1]
			dy = polygon[ip+1,2] - polygon[ip,2]
			dz = polygon[ip+1,3] - polygon[ip,3]
			a += sqrt(dx * dx + dy * dy + dz * dz)
		end

	end

	s /= a

end

return s

end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
]add DiffEqParamEstim
using DiffEqParamEstim|
]add DiffEqBayes
using DiffEqBayes
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")

ext: jl

max_stars_repo_name: onionpork/sparsemyxo

max_issues_repo_name: onionpork/sparsemyxo

max_forks_repo_name: onionpork/sparsemyxo


text:
module NaturallyUnitful

using Reexport
@reexport using Unitful

const c  = 1u"c"
const ħ  = 1.0545718e-34u"kg*m^2/s"
const ħc = ħ*c
const kB = 1.38064852e-23u"kg*m^2/s^2/K"
const ϵ0 = 8.85418782e-12u"m^-3*kg^-1*s^4*A^2"

setDimPow(D::Dict, x::Unitful.Dimension{T}) where {T} = D[T] = x.power

function dimDict(x::Unitful.Dimensions{T}) where {T}
    D = Dict(:Length => 0,
             :Mass => 0,
             :Temperature => 0,
             :Time => 0,
             :Current => 0)
    for t in T
        setDimPow(D,t)
    end
    return D
end

dimDict(x::Unitful.Quantity{T}) where {T} = dimDict(dimension(x))

dimDict(x) = Dict(:Length => 0, :Mass => 0, :Temperature => 0, :Time => 0, :Current => 0)

gettypeparams(::Unitful.FreeUnits{T, U, V}) where {T, U, V} = T, U, V
const energydimension = gettypeparams(u"GeV")[2]

"""
    natural(q; base=u"eV")

Convert `q` to natural units based on the units specified by `base`. If `base` is
unspecified, `natural` will default to `eV`. Currently, `natural` only supports
`base`s with dimensions of energy. For all other `base` you will need to use `unnatrual` '
to convert.

Examples:

    julia> natural(1u"kg")
    5.609588650020686e35 eV

    julia> natural(1u"kg", base=u"GeV")
    5.609588650020685e26 GeV

    julia> natural(1u"m")
    5.067730759202785e6 eV^-1

    julia> natural(1u"m", base=u"GeV")
    5.067730759202785e15 GeV^-1
"""
natural(q; base=u"eV") = _natural(base, q)

function _natural(base::Unitful.FreeUnits{T, energydimension, U}, q) where {T, U}
    D = dimDict(q)
    (α,β,γ,δ,ϕ) = (D[:Length], D[:Mass], D[:Temperature], D[:Time], D[:Current])
    uconvert(base^(-α-δ+β+γ+ϕ),
             q*ħc^(-α+ϕ/2)*ħ^(-δ)*c^(2(β-ϕ/2))*kB^(γ)*(4π*ϵ0)^(-ϕ/2))
end

function _natural(base::Unitful.FreeUnits{T, U, V}, q) where {T, U, V}
    throw("""natural(q; base)` where `base` has dimensions `$U` has not yet been implemented. Please use a base with dimensions of `energy` and then use `unnatural` to convert to your desired units.""")
end

"""
    unnatural(targetUnit, q)

Convert a quantity `q` to units specified by `targetUnits` automatically inserting whatever
factors  of `ħ`, `c`, `ϵ₀` and `kb` are required to make the conversion work.

Examples:

    julia> unnatural(u"m", 5.067730759202785e6u"eV^-1")
    1.0 m

    julia> unnatural(u"m/s", 1)
    2.99792458e8 m s^-1

    julia> unnatural(u"K", 1u"eV")
    11604.522060401006 K
"""
function unnatural(targetUnit, q)
    natTarget = natural(1targetUnit)
    natQ = natural(q)
    ratio = natQ/natTarget
    if typeof(ratio |> dimension) == Unitful.Dimensions{()}
        (ratio)targetUnit
    else
        throw(Unitful.DimensionError)
    end
end

export natural, unnatural

end # module


text:
function get_raw_comment(cursor::CLCursor)
    stringptr = clang_Cursor_getRawCommentText(cursor.cursor)
    if stringptr.data == C_NULL
        return nothing
    else
        cstring = clang_getCString(stringptr)
        return unsafe_string(cstring)
    end
end

function get_comment(obj::CLCursor)
    clang_Cursor_getParsedComment(obj.cursor)
end

function Clang.children(c::CXComment)

    num_children = Int(clang_Comment_getNumChildren(c))
    list = CXComment[]
    for i in 1:num_children
        child = clang_Comment_getChild(c, i-1)
        clang_Comment_isWhitespace(child) == UInt32(0) && push!(list, child)
    end
    return list
end

function Clang.name(c::CXComment)
end

function Clang.kind(c::CXComment)
    clang_Comment_getKind(c)
end


text:
module DrWatsonSim

using DrWatson
using BSON
using Reexport
@reexport using Dates

include("Locking.jl")
include("Metadata.jl")
include("Metadata_Search.jl")
include("Simulation.jl")
end

max_stars_repo_name: sebastianpech/DrWatsonSim

max_issues_repo_name: sebastianpech/DrWatsonSim

max_forks_repo_name: sebastianpech/DrWatsonSim


text:
using JGraph
using Test
using Statistics

@testset "JGraph.core" begin
    g = Graph();
    @test ne(g) == 0
    @test nv(g) == 0

    # add vertices
    g = Graph();
    add_vertex!(g, 10)
    @test ne(g) == 0
    @test nv(g) == 10

    # add a edge
    g = Graph();
    add_edge!(g, 3, 4);
    @test ne(g) == 1
    @test nv(g) == 4
    @test has_edge(g, 3, 4) == true
    @test has_edge(g, 4, 3) == true

    # test if the edge is already there
    @test add_edge!(g, 3, 4) == false

    # test the handshake lemma
    @test sum(degree(g)) == 2*ne(g)

    # test the neighbors function
    g = Graph();
    @test add_edge!(g, 1, 2) == true
    @test add_edge!(g, 2, 3) == true
    @test add_edge!(g, 3, 1) == true
    @test neighbors(g, 1) == [2, 3]

    # test itertor
    g = Graph();
    add_vertex!(g, 10);
    @test collect(vertices(g)) == collect(1:10)

    # test remove edge
    g = Graph()
    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 1)

    remove_edge!(g, 3, 1)

    @test !has_edge(g, 3, 1)
    @test !has_edge(g, 1, 3)

    e = Edge((1, 2))
    remove_edge!(g, e)

    @test !has_edge(g, 2, 1)
    @test !has_edge(g, 1, 2)

    #  test remove all_edges incident to a node
    g = Graph()
    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 1)

    remove_vertex!(g, 2)
    @test !has_edge(g, 2, 1)
    @test !has_edge(g, 1, 2)
    @test !has_edge(g, 2, 3)
    @test !has_edge(g, 3, 2)
end

@testset "JGraph.edge" begin
    e1 = Edge((1, 2))
    e2 = Edge(Pair(2, 1))
    e3 = Edge((2, 3))
    e4 = Edge((3, 2))
    e5 = Edge((1, 2))
    results = [e4, e3, e2, e1];

    @test src(e1) == 1
    @test dst(e1) == 2

    @test src(e2) == 2
    @test dst(e2) == 1
    @test e1 == e5

    ## test iterator
    g = Graph();
    add_edge!(g, 1, 2);
    add_edge!(g, 2, 3);

    @test length(edges(g)) == 2

    for e in edges(g)
        @test e == pop!(results)
    end
end

@testset "JGraph.generator" begin
    g = erdos_renyi(10, 0.2)
    @test nv(g) == 10

    # test the handshake lemma
    @test sum(degree(g)) == 2*ne(g)

    # test the mean degree
    n = 10000;
    p = 0.0003;
    kmean = (n-1) * p
    g = erdos_renyi(n, p);
    @test mean(degree(g)) ≈ kmean atol=1e-1

    g1 = erdos_renyi(n, p, seed=5);
    g2 = erdos_renyi(n, p, seed=5);
    @test nv(g1) == nv(g2)
    @test ne(g1) == ne(g2)
end

@testset "JGraph.traversal" begin
    g = Graph();
    @test add_edge!(g, 1, 2) == true
    @test add_edge!(g, 2, 3) == true
    @test add_edge!(g, 3, 1) == true
    @test bfs(g, 1) == Set([2, 1, 3])

    g = Graph();
    add_vertex!(g, 6);
    add_edge!(g, 1, 2);
    add_edge!(g, 2, 3);
    add_edge!(g, 1, 3);

    add_edge!(g, 4, 5);
    add_edge!(g, 5, 6);
    add_edge!(g, 4, 6);

    cc = connected_components(g);
    for elem in cc
        @test length(elem) == 3
    end

    @test bfs(g, 1) == Set([2, 3, 1])
    @test bfs(g, 4) == Set([4, 5, 6])
end


text:
"""Julia program to implement Exponential Search algorithm.
Exponential Search Algorithm is an optimized Binary Search to search an element in sorted array.
It is specifically used when the size of array is infinite.
"""

function binary_search(arr, low, n, ele)
    high = n
    while high >= low
        mid = low + (high - low) ÷ 2
        # If the mid element is the required element, return that index
        if(arr[mid] == ele)
            return mid
        # ele is greater than mid element then ele would be present in first half of the array
        elseif (arr[mid] > ele)
            high = mid - 1
        #Else if  ele is smaller than mid element then ele would be present in last half of the array
        else
            low = mid + 1
        end
    end
    # If the element is not found return 0
    return 0
end

function exponential_search(arr, n, ele)
    if(arr[1] == ele)
        return false
    end
    i = 2
    while( i < n && arr[i] <= ele)
        i = i * 2
    end
    mini = i < (n-1) ? i : (n-1)
    return binary_search(arr, i÷2, mini, ele )
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("Array is Empty!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
print("Which number do you want to search in the array? ")
ele = readline()
ele = parse(Int, ele)
# Sort the array in ascending order
arr = sort(arr)
res = exponential_search(arr, n, ele)
if (res == 0)
    print("The number $ele is not present in the array")
else
    print("The number $ele is present in the array.")
end

"""
Time Complexity - O(log(n)), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
How many numbers are present in the array? 5
Enter the numbers: 1 2 3 4 5
Which number do you want to search in the array? 6
The number 6 is not present in the array

SAMPLE II
How many numbers are present in the array? 3
Enter the numbers: 3 1 2
Which number do you want to search in the array? 2
The number 2 is present in the array.
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
struct Cache{TF<:AbstractFloat}
    Ax::Vector{TF}
    resid::Vector{TF}
    dresid::Vector{TF}
    g::Vector{TF}
    dx::Vector{TF}
end

function Cache(TF::Type, M::Integer, N::Integer)
    Ax = Vector{TF}(undef, M)
    resid = Vector{TF}(undef, M)
    dresid = Vector{TF}(undef, M)
    g = Vector{TF}(undef, N)
    dx = Vector{TF}(undef, N)
    return Cache{TF}(Ax, resid, dresid, g, dx)
end

"""
    Cache(A::AbstractMatrix)

Allocate arrays for the solver based on the size and element type of `A`.
"""
Cache(A::AbstractMatrix{TF}) where TF<:AbstractFloat = Cache(TF, size(A)...)

function updaterss!(A, b, x, ca)
    mul!(ca.Ax, A, x)
    @inbounds for i in eachindex(ca.resid)
        ca.resid[i] = ca.Ax[i] - b[i]
    end
    return sum(abs2, ca.resid)
end

function direction!(A, x, λ, ca)
    copyto!(ca.g, x)
    # Compute the half-gradient g = A'(Ax-b) + λx
    mul!(ca.g, A', ca.resid, true, λ)
    gmin, ix = findmin(ca.g)
    mul!(ca.dx, -one(eltype(A)), x)
    ca.dx[ix] += one(eltype(x))
    return ix
end

function stepsize!(ix, A, λ, ca)
    @inbounds for i in eachindex(ca.dresid)
        ca.dresid[i] = A[i,ix] - ca.Ax[i]
    end
    ss = -dot(ca.g, ca.dx) / (sum(abs2,ca.dresid) + λ*sum(abs2,ca.dx))
    TF = eltype(A)
    return min(one(TF), max(zero(TF), ss))
end

function updatex!(ss, x, ca)
    @inbounds for i in eachindex(x)
        x[i] += ss*ca.dx[i]
    end
end

"""
    convexfit(A::Matrix, b::Vector, λ::Real=0; kwargs...)

Fit vector `b` with a convex combination of the columns in matrix `A`
and an optional regularization parameter `λ`.

The associated constrained optimization problem can be expressed as follows:
```math
\\min_x \\|Ax - b\\|_2^2 + λ\\|x\\|_2^2 \\\\
st. x_i ≥ 0 \\forall i, \\, \\sum_i x_i = 1
```
where `x` is the vector of coefficients to be solved.

# Keywords
- `x0::AbstractVector{<:Real}=fill(1/size(A,2),1/size(A,2))`: initial value of `x`.
- `ftol::Real=1e-6`: convergence criterion based on the change in the objective function.
- `xtol::Real=1e-6`: convergence criterion based on the inf-norm change in `x`.
- `maxiter::Integer=1000`: maximum number of iterations.
- `store_trace::Bool=false`: save the solver state in each iteration.
- `show_trace::Bool=false`: print the solver state in each iteration to `stdout`.
- `show_thread::Bool=false`: show the thread ID when printing solver states.
- `cache::Cache=Cache(A)`: `Cache` that stores intermediate results.
"""
function convexfit(A::Matrix{TF}, b::Vector{TF}, λ::Real=zero(TF);
        x0::AbstractVector{<:Real}=fill(convert(TF,1/size(A,2)), size(A,2)),
        ftol::Real=convert(TF, 1e-6),
        xtol::Real=convert(TF, 1e-6),
        maxiter::Integer=1000,
        store_trace::Bool=false,
        show_trace::Bool=false,
        show_thread::Bool=false,
        cache::Cache=Cache(A)) where TF<:AbstractFloat

    sx = sum(x0)
    if !(sx ≈ 1)
        @warn "sum of elements in x0 ($sx) is normalized to 1"
        x0 = x0./sx
    end
    λ = convert(TF, λ)
    λ < 0 && throw(ArgumentError("regularization parameter cannot be negative"))

    f_converged = false
    x_converged = false
    rss = -one(TF)
    xnorm = -one(TF)
    flast = -one(TF)
    f = -one(TF)
    x = Vector{TF}(undef, length(x0))
    copyto!(x, x0)
    dfnorm = convert(TF, NaN)
    dxnorm = convert(TF, NaN)
    iter = 0
    if store_trace || show_trace
        tr = SolverTrace(TF)
    else
        tr = nothing
    end

    for _ in 1:maxiter
        flast = f
        rss = updaterss!(A, b, x, cache)
        xnorm = sum(abs2, x)
        f = rss + λ * xnorm
        dfnorm = abs(f-flast)
        if dfnorm < ftol
            f_converged = true
            break
        elseif x_converged
            break
        elseif iszero(f)
            f_converged = true
            x_converged = true
            dfnorm = zero(TF)
            dxnorm = zero(TF)
            break
        end
        ix = direction!(A, x, λ, cache)
        if iszero(cache.dx[ix])
            f_converged = true
            x_converged = true
            dfnorm = zero(TF)
            dxnorm = zero(TF)
            break
        end
        ss = stepsize!(ix, A, λ, cache)
        dxnorm = ss*maximum(abs, cache.dx)
        if dxnorm < xtol
            x_converged = true
            # Do not immediately break to allow checking f_converged with the new x
        end
        updatex!(ss, x, cache)
        # Ensure that iter always represents the number of steps from x0r
        iter += 1
        # xnorm should not be updated until the next iteration
        update!(tr, iter, f, rss, xnorm, ix, ss, store_trace, show_trace, show_thread)
    end
    return SolverResult{TF}(x0, x, copy(cache.Ax), iter, f, rss, λ, xnorm,
        dfnorm, ftol, f_converged, dxnorm, xtol, x_converged, tr)
end

function convexfit(A::Matrix{TF}, B::Matrix{TF}, λ=0.0;
        multithreads::Bool=false, show_thread::Bool=multithreads,
        kwargs...) where TF<:AbstractFloat
    N = size(B, 2)
    rs = Vector{Any}(undef, N)
    if multithreads
        haskey(kwargs, :cache) && @warn "keyword cache is ignored with multithreads=true"
        Threads.@threads for i in axes(B, 2)
            # Create new Cache for each call no matter whether cache is specified
            rs[i] = convexfit(A, view(B,:,i), λ;
                kwargs..., cache=Cache(A), show_thread=show_thread)
        end
    else
        ca = get(kwargs, :cache, nothing)
        ca === nothing && (ca = Cache(A))
        for i in axes(B, 2)
            get(kwargs, :show_trace, false) && println("target $i of $N:")
            rs[i] = convexfit(A, view(B,:,i), λ;
                kwargs..., cache=ca, show_thread=show_thread)
        end
    end
    return rs
end

"""
    convexfit(A::AbstractMatrix, B::AbstractVecOrMat, λ=0; multithreads=false, kwargs...)

A wrapper method of `convexfit` that accepts multiple `b`s in a matrix `B` and
fits each column in `B` with a convex combination of the columns in matrix `A`.
The actual problems to be solved depend on the type of `λ`.
Results are collected in an array in the order of the columns of `B`.

Problems across different columns can be solved simultaneously
by setting `multithreads=true`.
See also [`convexfit(::Matrix,::Vector,::Real)`](@ref).

!!! compat "Julia 1.2"
    Setting `multithreads=true` requires at least Julia 1.2.

"""
function convexfit(A::AbstractMatrix, b::AbstractVecOrMat, λ=0.0; kwargs...)
    TF = promote_type(eltype(A), eltype(b))
    TF <: Integer && (TF = Float64)
    A = convert(Matrix{TF}, A)
    b = convert(Array{TF}, b)
    # Convert x0 at this stage to avoid making multiple copies of x0 later
    x0 = get(kwargs, :x0, nothing)
    if x0 !== nothing
        x0 = convert(Vector{TF}, x0)
        return convexfit(A, b, λ; kwargs..., x0=x0)
    else
        return convexfit(A, b, λ; kwargs...)
    end
end

function convexfit(A::Matrix{TF}, λ=0.0;
        loo::Bool=true, multithreads::Bool=false, show_thread::Bool=multithreads,
        kwargs...) where TF<:AbstractFloat
    if loo
        N = size(A, 2)
        N > 1 || throw(ArgumentError("matrix A must contain at least two columns if loo=true"))
        rs = Vector{Any}(undef, N)
        if multithreads
            haskey(kwargs, :cache) && @warn "keyword cache is ignored with multithreads=true"
            Threads.@threads for i in axes(A, 2)
                Aloo = view(A,:,axes(A,2).!=i)
                # Create new Cache for each call no matter whether cache is specified
                rs[i] = convexfit(Aloo, view(A,:,i), λ;
                    kwargs..., cache=Cache(Aloo), show_thread=show_thread)
            end
        else
            ca = get(kwargs, :cache, nothing)
            ca === nothing && (ca = Cache(view(A,:,2:N)))
            for i in axes(A, 2)
                get(kwargs, :show_trace, false) && println("target $i of $N:")
                rs[i] = convexfit(view(A,:,axes(A,2).!=i), view(A,:,i), λ;
                    kwargs..., cache=ca, show_thread=show_thread)
            end
        end
        return rs
    else
        ca = get(kwargs, :cache, nothing)
        ca === nothing && (ca = Cache(A))
        return convexfit(A, A, λ; kwargs..., cache=ca)
    end
end

"""
    convexfit(A::AbstractMatrix, λ=0; loo=true, multithreads=false, kwargs...)

A wrapper method of `convexfit` for fitting
each column in matrix `A` with a convex combination of the other columns in matrix `A`.
The actual problems to be solved depend on the type of `λ`.
Results are collected in an array in the order of the columns of `B`.

By default, with `loo=true` the fitted column is not involved in the convex combination.
This can be altered by setting `loo=false`.
Problems across different columns can be solved simultaneously
by setting `multithreads=true`.
See also [`convexfit(::Matrix,::Vector,::Real)`](@ref).

!!! compat "Julia 1.2"
    Setting `multithreads=true` requires at least Julia 1.2.

"""
function convexfit(A::AbstractMatrix, λ=0.0; kwargs...)
    TF = eltype(A)
    TF <: Integer && (TF = Float64)
    A = convert(Matrix{TF}, A)
    # Convert x0 at this stage to avoid making multiple copies of x0 later
    x0 = get(kwargs, :x0, nothing)
    if x0 !== nothing
        x0 = convert(Vector{TF}, x0)
        return convexfit(A, λ; kwargs..., x0=x0)
    else
        return convexfit(A, λ; kwargs...)
    end
end


text:
#!/bin/env julia

#Base.compilecache(GP)

#include("./GP.jl")
import GP
#using GP

if length(ARGS) < 2
    println("give me input file names")
    exit(0)
end

input = ARGS

@time GP.process(input[1])
@time GP.process(input[2])

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
"""
    xcom(com, vel, marker)

Compute the extrapolated center of mass.

`com`, and `marker` are all assumed to be 3-dimensional variables where the order of
dimensions is `[ML, AP, VT]`. Units are assumed to be standard SI units.

See McAndrew et al. (2012)[https://doi.org/10.1016/j.jbiomech.2011.12.027]
"""
function xcom(com, bos)
    size(com, 1) === size(bos, 1) || throw(ArgumentError("Both signals must have equal length"))

    # tsum = sum((com - bos).^2; dims=2)
    # l = mean(sqrt, tsum)
    d = com - bos
    l = mean(filter(!isnan, mapslices(norm, d; dims=2)))
    ω₀ = sqrt(9.81/l)

    vel = centraldiff(com; padding=ForwardBackwardPad())

    xCOM = com + vel ./ ω₀

    return xCOM
end

function margin_of_stability(com, lbos, rbos, levents, revents; axis=1)
    lmos = lbos - xcom(com, lbos)
    rmos = rbos - xcom(com, rbos)
    sp = sortperm([levents; revents])

    lmos_ev, rmos_ev = lmos[levents, axis], rmos[revents, axis]
    return lmos_ev, rmos_ev, [lmos_ev; rmos_ev][sp]
end

function phase_coordination_index(;lfs, rfs, lfo, rfo)
    lswing = mean(swing(lfs, lfo))
    rswing = mean(swing(rfs, rfo))

    if lswing < rswing
        tₗ, tₛ = beginwithevent(rfs, lfs)
    else
        tₗ, tₛ = beginwithevent(lfs, rfs)
    end
    resize!(tₗ, min(length(tₗ), length(tₛ)))
    resize!(tₛ, min(length(tₗ), length(tₛ)))

    φ = 360 .* (tₛ[1:end-1] .- tₗ[1:end-1])./stridetimes(tₗ)
    φ_ABS = mean(abs, φ .- 180)
    Pφ_ABS = (φ_ABS/180)*100
    φ_CV = variation(φ)*100

    return φ_CV + Pφ_ABS
end


text:
#=
---------------------------------------------------------------------
This function computes the RHS vector for a unified CGDG method for the
2D Wave Equation.

Written by F.X. Giraldo on July 12, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

include("create_rhs_volume.jl")
include("create_rhs_flux.jl")

function create_rhs(q,u,M,De_x,De_y,intma,periodicity,face,mapL,mapR,normals,jac_face,ωq,space_method,DFloat)

    #---------------------------Students Add their Contributions Here----------------------------------------
    #---------------------------Students Add their Contributions Here----------------------------------------
    #Construct Volume integral contribution
    #rhs = create_rhs_volume(q,u,De_x,De_y,intma,periodicity,DFloat)

    #Construct Flux integral contribution (will also work for CG but cancels completely due to periodicity)
    #if (space_method == "DG")
    #    rhs = create_rhs_flux!(rhs,q,u,face,normals,jac_face,ωq,mapL,mapR,intma,periodicity,DFloat)
    #end
    #---------------------------Students Add their Contributions Here----------------------------------------
    #---------------------------Students Add their Contributions Here----------------------------------------

    #Multiply by inverse Mass matrix
    rhs=M\rhs

    #Apply Periodicity
    Npoin=length(q)
    for i=1:Npoin
        rhs[i]=rhs[periodicity[i]]
    end

    #Return arrays
    return(rhs)
end


text:
# Demonstrate NTF synthesis (synthesizeNTF)
using RSDeltaSigmaPort
j=im

#==Baseband modulator
===============================================================================#
println("\n*** 5th order, 2-level, baseband modulator")
OSR = 32

@info("Synthesizing NTF without zero-optimization..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
NTF_noopt = synthesizeNTF(5, OSR, opt=0)
println("\tdone.")

plot = plotNTF(NTF_noopt, OSR, color=:blue)
plot.title = "5th-Order Modulator (No Zero-Optimization)"
saveimage(:png, "dsdemo1_o5_noopt.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Synthesizing NTF with optimized zeros..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
NTF_opt = synthesizeNTF(5, OSR, opt=1)
println("\tdone.")

plot = plotNTF(NTF_opt, OSR, color=:red)
plot.title = "5th-Order Modulator (Optimized Zeros)"
saveimage(:png, "dsdemo1_o5_zopt.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Plotting NTF comparison (overlay results)")
#-------------------------------------------------------------------------------
plot = plotNTF(NTF_noopt, OSR, color=:blue)
plot = plotNTF!(plot, NTF_opt, OSR, color=:red)
plot.title = "5th-Order Modulator (Optimized Zeros - Overlay)"
saveimage(:png, "dsdemo1_o5_cmp.png", plot, AR=2/1, width=900)
displaygui(plot)

#==Bandpass modulator
===============================================================================#
println("\n*** 8th order, 2-level, bandpass modulator")
OSR = 64
order = 8
f0 = 0.125 #fs/8

function calcSTF(order, OSR, NTF, f0)
	G = _zpk(zeros(array_round(order/2)),NTF.p,1,1)
	G.k = 1/abs(evalTF(G,exp(2π*j*f0)))
	return G
end

@info("Synthesizing NTF..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
NTF = synthesizeNTF(order, OSR, opt=2, f0=f0)
println("\tdone.")

@info("Plotting NTF")
#-------------------------------------------------------------------------------
STF = calcSTF(order, OSR, NTF, f0)
plot = plotNTF(NTF, OSR, color=:blue, f0=f0, STF=STF)
plot.title = "8th-Order Bandpass Modulator"
saveimage(:png, "dsdemo1_o8_bp.png", plot, AR=2/1, width=900)
displaygui(plot)

:END_OF_DEMO


text:
# Tests for Radial DEA Models
@testset "RadialDEAModel" begin

    ## Test Radial DEA Models with FLS Book data
    # Test against results in R
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    # Input oriented CRS
    deaio = dea(X, Y, orient = :Input, rts = :CRS)

    @test typeof(deaio) == RadialDEAModel

    @test nobs(deaio) == 11
    @test ninputs(deaio) == 2
    @test noutputs(deaio) == 1
    @test efficiency(deaio) ≈ [1.0000000000;
                               0.6222896791;
                               0.8198562444;
                               1.0000000000;
                               0.3103709311;
                               0.5555555556;
                               1.0000000000;
                               0.7576690896;
                               0.8201058201;
                               0.4905660377;
                               1.0000000000]
    @test convert(Matrix, peers(deaio)) ≈
    [1.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.4249783174  0  0 0.10928013877  0  0   0   0;
     1.134321653  0  0 0.4380053908  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 1.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.2573807721  0  0 0.04844814534  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 0.33333333333  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.00000000000  0  0   0   0;
     0.000000000  0  0 1.0348650979  0  0 0.11457435013  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.14814814815  0  0   0   0;
     0.000000000  0  0 0.4905660377  0  0 0.49056603774  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   1.000000000]
     @test slacks(deaio, :X) ≈ [0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                4.444444444   0;
                                0.000000000   0;
                                0.000000000   0;
                                1.640211640   0;
                                0.000000000   0;
                                0.000000000   4]
    @test slacks(deaio, :Y) ≈ zeros(11)

    @test efficiency(dea(targets(deaio, :X), targets(deaio, :Y), orient = :Input, rts = :CRS, slack = false)) ≈ ones(11)
    @test efficiency(deaadd(targets(deaio, :X), targets(deaio, :Y))) ≈ zeros(11) atol=1e-14

    @test peersmatrix(deaio) == deaio.lambda

    # Otuput oriented CRS
    deaoo = dea(X, Y, orient = :Output, rts = :CRS)

    @test nobs(deaoo) == 11
    @test ninputs(deaoo) == 2
    @test noutputs(deaoo) == 1
    @test efficiency(deaoo) ≈ [1.0000000000;
                               1.606968641;
                               1.219726027;
                               1.0000000000;
                               3.221951220;
                               1.800000000;
                               1.0000000000;
                               1.319837398;
                               1.219354839;
                               2.038461538;
                               1.0000000000]
    @test convert(Matrix, peers(deaoo)) ≈
    [1.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.6829268293  0  0 0.1756097561   0  0   0   0;
     1.383561644  0  0 0.5342465753  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 1.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.8292682927  0  0 0.1560975610   0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 0.6000000000   0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.00000000000  0  0   0   0;
     0.000000000  0  0 1.3658536585  0  0 0.1512195122   0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.4000000000   0  0   0   0;
     0.000000000  0  0 1.0000000000  0  0 1.0000000000   0  0   0   0;
     1.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   0]
     @test slacks(deaoo, :X) ≈ [0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                8   0;
                                0.000000000   0;
                                0.000000000   0;
                                2   0;
                                0.000000000   0;
                                0.000000000   4]
    @test slacks(deaoo, :Y) ≈ zeros(11)

    @test efficiency(dea(targets(deaoo, :X), targets(deaoo, :Y), orient = :Output, rts = :CRS, slack = false)) ≈ ones(11)
    @test efficiency(deaadd(targets(deaoo, :X), targets(deaoo, :Y))) ≈ zeros(11) atol=1e-10

    # Input oriented VRS
    deaiovrs = dea(X, Y, orient = :Input, rts = :VRS)

    @test nobs(deaiovrs) == 11
    @test ninputs(deaiovrs) == 2
    @test noutputs(deaiovrs) == 1
    @test efficiency(deaiovrs) ≈ [1.0000000000;
                                  0.8699861687;
                                  1.0000000000;
                                  1.0000000000;
                                  0.7116402116;
                                  1.0000000000;
                                  1.0000000000;
                                  1.0000000000;
                                  1.0000000000;
                                  0.4931209269;
                                  1.0000000000]
    @test convert(Matrix, peers(deaiovrs)) ≈
    [1.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  0  0   0   0;
     0.52558782849  0  0 0.0000000000  0  0.2842323651  0.1901798064   0  0   0   0;
     0.000000000    0  1 0.0000000000  0  0.00000000000 0.00000000000  0  0   0   0;
     0.000000000    0  0 1.0000000000  0  0.00000000000 0.00000000000  0  0   0   0;
     0.56613756614  0  0 0.0000000000  0  0.4338624339  0.00000000000  0  0   0   0;
     0.000000000    0  0 0.0000000000  0  1.00000000000 0.00000000000  0  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 1.00000000000  0  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  1  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  0  1   0   0;
     0.03711078928  0  0 0.4433381608  0  0.00000000000 0.5195510500   0  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  0  0   0   1.000000000]
     @test slacks(deaiovrs, :X) ≈ [0.000000000   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0   0;
                                   0.000000000   0;
                                   0.000000000   4]
    @test slacks(deaiovrs, :Y) ≈ [0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  2.698412698;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000]

    @test efficiency(dea(targets(deaiovrs, :X), targets(deaiovrs, :Y), orient = :Input, rts = :VRS, slack = false)) ≈ ones(11)
    @test efficiency(deaadd(targets(deaiovrs, :X), targets(deaiovrs, :Y))) ≈ zeros(11) atol=1e-12

    # Output oriented VRS
    deaoovrs = dea(X, Y, orient = :Output, rts = :VRS)

    @test nobs(deaoovrs) == 11
    @test ninputs(deaoovrs) == 2
    @test noutputs(deaoovrs) == 1
    @test efficiency(deaoovrs) ≈ [1.0000000000;
                                  1.507518797;
                                  1.0000000000;
                                  1.0000000000;
                                  3.203947368;
                                  1.000000000;
                                  1.0000000000;
                                  1.000000000;
                                  1.000000000;
                                  1.192307692;
                                  1.0000000000]
    @test convert(Matrix, peers(deaoovrs)) ≈
    [1.000000000   0  0 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.38157894737 0  0 0.1710526316  0  0 0.4473684211   0  0   0   0;
     0.000000000   0  1 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000   0  0 1.0000000000  0  0 0.00000000000  0  0   0   0;
     0.03947368421 0  0 0.7763157895  0  0 0.1842105263   0  0   0   0;
     0.000000000   0  0 0.0000000000  0  1 0.00000000000  0  0   0   0;
     0.000000000   0  0 0.0000000000  0  0 1.00000000000  0  0   0   0;
     0.000000000   0  0 0.0000000000  0  0 0.00000000000  1  0   0   0;
     0.000000000   0  0 0.0000000000  0  0 0.00000000000  0  1   0   0;
     0.000000000   0  0 0.0000000000  0  0 0.00000000000  0  1   0   0;
     1.000000000   0  0 0.0000000000  0  0 0.00000000000  0  0   0   0]
    @test slacks(deaoovrs, :X) ≈ [0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                5   11;
                                0.000000000   4]
    @test slacks(deaoovrs, :Y) ≈ zeros(11) atol=1e-10

    @test efficiency(dea(targets(deaoovrs, :X), targets(deaoovrs, :Y), orient = :Output, rts = :VRS, slack = false)) ≈ ones(11)
    @test efficiency(deaadd(targets(deaoovrs, :X), targets(deaoovrs, :Y))) ≈ zeros(11) atol=1e-12

    # Test no slacks
    deaionoslack = dea(X, Y, slack = false)
    @test efficiency(deaionoslack) == efficiency(deaio)
    @test isempty(slacks(deaionoslack, :X)) == 1
    @test isempty(slacks(deaionoslack, :Y)) == 1

    @test efficiency(dea(targets(deaionoslack, :X), targets(deaionoslack, :Y), slack = false)) ≈ ones(11)
    @test efficiency(deaadd(targets(deaionoslack, :X), targets(deaionoslack, :Y))) != zeros(11) # Different as there is no slacks in first model

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deaio_ref_eff = zeros(size(X, 1))
    deaoo_ref_eff = zeros(size(X, 1))

    deaiovrs_ref_eff = zeros(size(X, 1))
    deaoovrs_ref_eff = zeros(size(X, 1))

    deaiovrs_ref_slackX = zeros(size(X))
    deaiovrs_ref_slackY = zeros(size(Y))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deaio_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Input, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaoo_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Output, rts = :CRS, Xref = Xref, Yref = Yref))[1]

        deaiovrs_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref))[1]
        deaoovrs_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Output, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaiovrs_ref_slackX[i,:] = slacks(dea(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref), :X)
        deaiovrs_ref_slackY[i,:] = slacks(dea(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref), :Y)
    end

    @test deaio_ref_eff ≈ efficiency(deaio)
    @test deaoo_ref_eff ≈ efficiency(deaoo)

    @test deaiovrs_ref_eff ≈ efficiency(deaiovrs)
    @test deaoovrs_ref_eff ≈ efficiency(deaoovrs)

    @test deaiovrs_ref_slackX ≈ slacks(deaiovrs, :X) atol=1e-14
    @test deaiovrs_ref_slackY ≈ slacks(deaiovrs, :Y) atol=1e-15

    # Print
    show(IOBuffer(), deaio)
    show(IOBuffer(), deaionoslack)

    # Test errors
    @test_throws DimensionMismatch dea([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch dea([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch dea([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch dea([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError dea([1; 2; 3], [4; 5; 6], orient = :Error) # Invalid orientation
    @test_throws ArgumentError dea([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale

    @test_throws ArgumentError dea([1; 2 ; 3], [4 ; 5; 6], disposX = :Error)  # Invalid inputs disposability
    @test_throws ArgumentError dea([1; 2 ; 3], [4 ; 5; 6], disposY = :Error)  # Invalid outputs disposability

    @test_throws ArgumentError targets(deaio, :Error) # Invalid target
    @test_throws ArgumentError slacks(deaio, :Error) # Invalid slacks
    @test_throws ArgumentError rts(deaio) # rts only for DEA models in multiplier form

    # ------------------
    # Weak Disposability Tests
    # ------------------

    X = [1; 2; 3; 2; 4]
    Y = [2; 3; 4; 1; 3]

    deaioStrong = dea(X, Y, orient = :Input, rts = :VRS)
    @test efficiency(deaioStrong ) ≈ [1.0; 1.0; 1.0; 0.5; 0.5]
    @test slacks(deaioStrong, :X) ≈ [0; 0; 0; 0; 0] atol=1e-15
    @test slacks(deaioStrong, :Y) ≈ [0; 0; 0; 1; 0] atol=1e-15

    deaioWeak = dea(X, Y, orient = :Input, rts = :VRS, disposY = :Weak)
    @test efficiency(deaioWeak ) ≈ [1.0; 1.0; 1.0; 1.0; 0.5]
    @test slacks(deaioWeak, :X) ≈ [0; 0; 0; 0; 0] atol=1e-15
    @test slacks(deaioWeak, :Y) ≈ [0; 0; 0; 0; 0] atol=1e-15

    deaooStrong = dea(X, Y, orient = :Output, rts = :VRS)
    @test efficiency(deaooStrong ) ≈ [1.0; 1.0; 1.0; 3.0; 1.3333333333333333]
    @test slacks(deaooStrong, :X) ≈ [0; 0; 0; 0; 1] atol=1e-15
    @test slacks(deaooStrong, :Y) ≈ [0; 0; 0; 0; 0] atol=1e-15

    deaooWeak = dea(X, Y, orient = :Output, rts = :VRS, disposX = :Weak)
    @test efficiency(deaooWeak ) ≈ [1.0; 1.0; 1.0; 3.0; 1.0]
    @test slacks(deaooWeak, :X) ≈ [0; 0; 0; 0; 0] atol=1e-14
    @test slacks(deaooWeak, :Y) ≈ [0; 0; 0; 0; 0] atol=1e-14

    # Test if weak disposability in inputs in the Input oriented model works
    # In this example, same result as stron disposability
    @test efficiency(dea(X, Y, orient = :Input, disposX = :Strong)) ==
        efficiency(dea(X, Y, orient = :Input, disposX = :Weak))

    # Test if weak disposability in outputs in the Output oriented model works
    # In this example, same result as strong disposability
    @test efficiency(dea(X, Y, orient = :Output, disposY = :Strong)) ==
        efficiency(dea(X, Y, orient = :Output, disposY = :Weak))

    # ------------------
    # DMU names
    # ------------------

    X = [1; 2; 3; 2; 4]
    Y = [2; 3; 4; 1; 3]

    @test names(dea(X, Y)) == ["1"; "2"; "3"; "4"; "5"]
    @test names(dea(X, Y, names = ["A", "B", "C", "D", "E"])) == ["A", "B", "C", "D", "E"]

    logs, value = Test.collect_test_logs() do
        names(dea(X, Y, names = ["A", "B", "C", "D"]))
    end
    @test occursin("Length of names lower than number of observations", string(logs))
    @test value == ["A", "B", "C", "D", "5"]

    logs, value = Test.collect_test_logs() do
        names(dea(X, Y, names = ["A", "B", "C", "D", "E", "F"]))
    end
    @test occursin("Length of names greater than number of observations", string(logs))
    @test value == ["A", "B", "C", "D", "E"]

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(dea(X, Y, orient = :Input)) ≈ [1; 1; 1; 0.6; 0.4; 1; 0.6666666667; 0.625]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    @test efficiency(dea(X, Y, orient = :Output)) ≈ [1; 1; 1; 1.555555556; 2.333333333; 1; 1.272727273; 1.6]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    @test efficiency(dea(X, Y, orient = :Input)) ≈ [0.4; 1; 0.8; 0.6; 0.4; 0.4; 0.5142857143; 0.2]

    # ------------------
    # Run with different progress meter options
    # ------------------
    # Just run the code for code coverage as the progress meter does not affect the result
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    dea(X, Y, progress = true)
    dea(X, Y, progress = false)

end


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

import Optim
import BlackBoxOptim

# Changes block parameters to the optimized values
function fitParameters!(blocksWithOutputs::AbstractArray{<:AbstractBlock},
    parametersAndBounds::AbstractArray{Tuple{String, Float64, Float64}, 1},
    outputsAndStds::AbstractDict{String, Tuple{S, T}}; # Name to value and std.
    kwargs...) where {S, T}
    result = fitParameters(blocksWithOutputs, parametersAndBounds, outputsAndStds; kwargs...);
    for block in blocksWithOutputs
        setParameters!(block, result.parameterToMinimum);
    end
    return result;
end

function fitParameters(blocksWithOutputs::AbstractArray{<:AbstractBlock},
    parametersAndBounds::AbstractArray{Tuple{String, Float64, Float64}, 1},
    outputsAndStds::AbstractDict{String, Tuple{S, T}}; # Name to value and std.
    MaxTime = 10) where {S, T}
    boundBlocks = [BlockWithBindings(block, [name for (name, lower, upper) in parametersAndBounds]) for block in blocksWithOutputs];
    lower = [lower for (_, lower, _) in parametersAndBounds];
    upper = [upper for (_, _, upper) in parametersAndBounds];

    outputs = [computeOutputs(block) for block in blocksWithOutputs];

    expectedValuesArray = [];
    stdsArray = [];
    for output in outputs
        expectedValues = [];
        stds = [];
        for variable in output.variables
            if (haskey(outputsAndStds, variable.name))
                push!(expectedValues, outputsAndStds[variable.name][1]);
                push!(stds, outputsAndStds[variable.name][2]);
            else
                push!(expectedValues, 0.0);
                push!(stds, Inf64);
            end
        end
        push!(expectedValuesArray, expectedValues);
        push!(stdsArray, stds);
    end

    objective = x -> begin
        obj = 0.;
        for i in 1:length(boundBlocks)
            setParameters!(boundBlocks[i], x);
            outputs = computeOutputs(boundBlocks[i]);
            outputs = unfold(outputs.values - expectedValuesArray[i]) ./ unfold(stdsArray[i]);
            obj += sum(outputs .^ 2)
        end
        return obj
    end

    result = BlackBoxOptim.bboptimize(objective; SearchRange = collect(zip(lower, upper)), MaxTime = MaxTime);
    parameterToMinimum = Dict{String, Any}();
    for i = 1:length(lower)
        parameterToMinimum[parametersAndBounds[i][1]] = BlackBoxOptim.best_candidate(result)[i];
    end
    return (rawResult = result, minimizer = BlackBoxOptim.best_candidate(result), parameterToMinimum = parameterToMinimum);
end

function unfold(A)
    V = []
    for x in A
        if x === A
            push!(V, x)
        else
            append!(V, unfold(x))
        end
    end
    V
end

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
function _get_potential_plot_information(p::ElectricPotential)::Tuple{Symbol, Tuple{Real, Real}, String, Unitful.Units}
    return (:viridis, (minimum(p.data), maximum(p.data)), "Electric Potential", internal_voltage_unit)
end

function _get_potential_plot_information(::WeightingPotential)::Tuple{Symbol, Tuple{Real, Real}, String, Unitful.Units}
    return (:viridis, (0, 1), "Weighting Potential", Unitful.NoUnits)
end

function _get_potential_plot_information(ecd::EffectiveChargeDensity)::Tuple{Symbol, Tuple{Real, Real}, String, Unitful.Units}
    return (:inferno, (minimum(ecd.data), maximum(ecd.data)), "Effective Charge Density", u"V * m")
end

function _get_potential_plot_information(::PointTypes)::Tuple{Symbol, Tuple{Real, Real}, String, Unitful.Units}
    return (:viridis, (0, 7), "Point Type Map", Unitful.NoUnits)
end

function get_crosssection_idx_and_value(grid::Grid, ::Any, ::Any, ::Any)
    error(ArgumentError, ": The cross section your are trying to plot is not given in the right format.\nPlease pass a number with matching or no units.")
end

## Potentials with cylindrical grids

function get_crosssection_idx_and_value(grid::CylindricalGrid{T}, r::Union{Real, LengthQuantity, Missing}, φ::Union{Real, AngleQuantity, Missing}, z::Union{Real, LengthQuantity, Missing})::Tuple{Symbol,Int,T,Unitful.Units} where {T <: SSDFloat}

    cross_section::Symbol, idx::Int, units::Unitful.Units = if ismissing(φ) && ismissing(r) && ismissing(z)
        return get_crosssection_idx_and_value(grid, r, T(0.0), z)
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        φ_rad::T = T(to_internal_units(φ isa Real ? deg2rad(φ) : φ))
        while !(grid.φ.interval.left <= φ_rad <= grid.φ.interval.right) && grid.φ.interval.right != grid.φ.interval.left
            if φ_rad > grid.φ.interval.right
                φ_rad -= width(grid.φ.interval)
            elseif φ_rad < grid.φ.interval.left
                φ_rad += width(grid.φ.interval)
            end
        end
        :φ, searchsortednearest(grid.φ, φ_rad), φ isa Real ? u"°" : unit(φ)
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, searchsortednearest(grid.r, T(to_internal_units(r))), r isa Real ? internal_length_unit : unit(r)
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, searchsortednearest(grid.z, T(to_internal_units(z))), z isa Real ? internal_length_unit : unit(z)
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end

    value::T = if cross_section == :φ
        rad2deg(grid.φ[idx])
    elseif cross_section == :r
        grid.r[idx]
    elseif cross_section == :z
        grid.z[idx]
    end

    cross_section, idx, value, units
end

@recipe function f(sp::ScalarPotential{T,3,Cylindrical}; r = missing, φ = missing, z = missing, contours_equal_potential = false, full_det = false) where {T <: SSDFloat}

    gradient::Symbol, clims::Tuple{Real, Real}, name::String, punit::Unitful.Units = _get_potential_plot_information(sp)

    if !(sp.grid[2][end] - sp.grid[2][1] ≈ 2π) sp = get_2π_potential(sp, n_points_in_φ = 72) end

    grid::CylindricalGrid{T} = sp.grid
    cross_section::Symbol, idx::Int, value::T, units::Unitful.Units = get_crosssection_idx_and_value(grid, r, φ, z)

    title --> name * " @ $(cross_section) = $(round(units, Float64(uconvert(units,value*(cross_section == :φ ? u"°" : internal_length_unit))), sigdigits=3))"
    colorbar_title --> name
    seriescolor --> gradient
    clims --> clims
    unitformat --> :slash

    sp, cross_section, idx, value, punit, contours_equal_potential, full_det
end

@recipe function f(sp::ScalarPotential{T,3,Cylindrical}, cross_section::Symbol, idx::Int, value::T, punit::Unitful.Units, contours_equal_potential::Bool, full_det::Bool) where {T <: SSDFloat}

    grid::CylindricalGrid{T} = sp.grid
    data = sp.data
    if eltype(data) == PointType
        data = data .& 0x07
    end
    @series begin
        seriestype := :heatmap
        foreground_color_border --> nothing
        tick_direction --> :out
        unitformat --> :slash
        if cross_section == :φ
            aspect_ratio --> 1
            xguide --> "r"
            yguide --> "z"
            xlims --> (grid.r[1],grid.r[end])
            ylims --> (grid.z[1],grid.z[end])
            gr_ext::Array{T,1} = midpoints(get_extended_ticks(grid.r))
            gz_ext::Array{T,1} = midpoints(get_extended_ticks(grid.z))
            if full_det
                cross_section_dummy, idx_mirror, value_dummy = get_crosssection_idx_and_value(grid, missing, value+180, missing)
                extended_data =  cat(data[end:-1:2, idx_mirror, :]', data[:, idx, :]', dims = 2)
                xlims := (-1*grid.r[end],grid.r[end])
                vcat(-1 .* grid.r[end:-1:2], grid.r)*internal_length_unit, grid.z*internal_length_unit, extended_data*punit
             else
                midpoints(gr_ext)*internal_length_unit, midpoints(gz_ext)*internal_length_unit, data[:,idx,:]'*punit
            end
        elseif cross_section == :r
            xguide --> "φ"
            yguide --> "z"
            ylims --> (grid.z[1],grid.z[end])
            grid.φ*internal_angle_unit, grid.z*internal_length_unit, data[idx,:,:]'*punit
        elseif cross_section == :z
            projection --> :polar
            xguide --> ""
            yguide --> ""
            grid.φ*internal_angle_unit, grid.r*internal_length_unit, data[:,:,idx].*punit
        end
    end

    if contours_equal_potential && cross_section == :φ
        @series begin
            seriescolor := :thermal
            seriestype := :contours
            unitformat --> :slash
            #if cross_section == :φ
                aspect_ratio --> 1
                xguide := "r"
                yguide := "z"
                xlims --> (grid.r[1],grid.r[end])
                ylims --> (grid.z[1],grid.z[end])
                if full_det
                    cross_section_dummy, idx_mirror, value_dummy = get_crosssection_idx_and_value(grid, missing, value+180, missing)
                    extended_data =  cat(data[end:-1:2, idx_mirror, :]', data[:, idx, :]', dims = 2)
                    xlims := (-1*grid.r[end],grid.r[end])
                    vcat(-1 .* grid.r[end:-1:2], grid.r)*internal_length_unit, grid.z*internal_length_unit, extended_data.*punit
                 else
                    # midpoints(gr_ext), midpoints(gz_ext), data[:,idx,:]'*unit
                    grid.r*internal_length_unit, grid.z*internal_length_unit, data[:,idx,:]'.*punit
                end
                #=
            elseif cross_section == :r
                xguide --> "φ / °"
                yguide --> "z / m"
                ylims --> (grid.z[1],grid.z[end])
                grid.φ, grid.z, data[idx,:,:]'
            elseif cross_section == :z
                projection --> :polar
                grid.φ, grid.r, data[:,:,idx]
            end
            =#
        end
    end
end

@recipe function f(ϵ::DielectricDistribution{T,3,Cylindrical}; φ = 0) where {T <: SSDFloat}

    grid::CylindricalGrid{T} = ϵ.grid
    cross_section::Symbol, idx::Int, value::T, units::Unitful.Units = get_crosssection_idx_and_value(grid, missing, φ, missing)

    seriestype --> :heatmap
    seriescolor --> :inferno
    foreground_color_border --> nothing
    tick_direction --> :out
    title --> "Dielectric Distribution @ $(cross_section) = $(round(units, Float64(uconvert(units,value*(cross_section == :φ ? u"°" : internal_length_unit))), sigdigits=3))"
    colorbar_title --> "Dielectric Distribution"
    unitformat --> :slash

    gr_ext::Array{T,1} = midpoints(get_extended_ticks(grid.r))
    gφ_ext::Array{T,1} = midpoints(get_extended_ticks(grid.φ))
    gz_ext::Array{T,1} = midpoints(get_extended_ticks(grid.z))

    @series begin
        if cross_section == :φ
            aspect_ratio --> 1
            xguide --> "r"
            yguide --> "z"
            xlims --> (grid.r[1],grid.r[end])
            ylims --> (grid.z[1],grid.z[end])
            gr_ext*internal_length_unit, gz_ext*internal_length_unit, ϵ.data[:,idx,:]'
        elseif cross_section == :r
            xguide --> "φ"
            yguide --> "z"
            ylims --> (grid.z[1],grid.z[end])
            gφ_ext*internal_angle_unit, gz_ext*internal_length_unit, ϵ.data[idx,:,:]'
        elseif cross_section == :z
            projection --> :polar
            gφ_ext*internal_angle_unit, gr_ext*internal_length_unit, ϵ.data[:,:,idx]
        end
    end
end

## Potentials with Cartesian grids

function get_crosssection_idx_and_value(grid::CartesianGrid3D{T}, x::Union{Real, LengthQuantity, Missing}, y::Union{Real, LengthQuantity, Missing}, z::Union{Real, LengthQuantity, Missing})::Tuple{Symbol,Int,T,Unitful.Units} where {T <: SSDFloat}

    cross_section::Symbol, idx::Int, units::Unitful.Units = if ismissing(x) && ismissing(y) && ismissing(z)
        return get_crosssection_idx_and_value(grid, T(0.0), y, z)
    elseif !ismissing(x) && ismissing(y) && ismissing(z)
        :x, searchsortednearest(grid.x, T(to_internal_units(x))), x isa Real ? internal_length_unit : unit(x)
    elseif ismissing(x) && !ismissing(y) && ismissing(z)
        :y, searchsortednearest(grid.y, T(to_internal_units(y))), y isa Real ? internal_length_unit : unit(y)
    elseif ismissing(x) && ismissing(y) && !ismissing(z)
        :z, searchsortednearest(grid.z, T(to_internal_units(z))), z isa Real ? internal_length_unit : unit(z)
    else
        error(ArgumentError, ": Only one of the keywords `x, y, z` is allowed.")
    end
    value::T = grid[cross_section][idx]
    cross_section, idx, value, units
end

@recipe function f(sp::ScalarPotential{T,3,Cartesian}; x = missing, y = missing, z = missing, contours_equal_potential = false) where {T <: SSDFloat}

    gradient::Symbol, clims::Tuple{Real, Real}, name::String, punit::Unitful.Units = _get_potential_plot_information(sp)

    grid::CartesianGrid3D{T} = sp.grid
    cross_section::Symbol, idx::Int, value::T, units::Unitful.Units = get_crosssection_idx_and_value(grid, x, y, z)

    title --> name * " @ $(cross_section) = $(round(units, Float64(uconvert(units,value*internal_length_unit)), sigdigits=3))"
    colorbar_title --> name
    seriescolor --> gradient
    clims --> clims

    sp, cross_section, idx, value, punit, contours_equal_potential
end

@recipe function f(sp::ScalarPotential{T,3,Cartesian}, cross_section::Symbol, idx::Int, value::T, punit::Unitful.Units, contours_equal_potential::Bool) where {T <: SSDFloat}

    grid::CartesianGrid3D{T} = sp.grid
    data = sp.data
    if eltype(data) == PointType
        data = data .& 0x07
    end
    @series begin
        seriestype := :heatmap
        foreground_color_border --> nothing
        tick_direction --> :out
        unitformat --> :slash
        if cross_section == :x
            aspect_ratio --> 1
            xguide --> "y"
            yguide --> "z"
            xlims --> (grid.y[1],grid.y[end])
            ylims --> (grid.z[1],grid.z[end])
            gy_ext = midpoints(get_extended_ticks(grid.y))
            gz_ext = midpoints(get_extended_ticks(grid.z))
            midpoints(gy_ext)*internal_length_unit, midpoints(gz_ext)*internal_length_unit, data[idx,:,:]'*punit
        elseif cross_section == :y
            aspect_ratio --> 1
            xguide --> "x"
            yguide --> "z"
            xlims --> (grid.x[1],grid.x[end])
            ylims --> (grid.z[1],grid.z[end])
            gx_ext = midpoints(get_extended_ticks(grid.x))
            gz_ext = midpoints(get_extended_ticks(grid.z))
            midpoints(gx_ext)*internal_length_unit, midpoints(gz_ext)*internal_length_unit, data[:,idx,:]'*punit
        elseif cross_section == :z
            aspect_ratio --> 1
            xguide --> "x"
            yguide --> "y"
            xlims --> (grid.x[1],grid.x[end])
            ylims --> (grid.y[1],grid.y[end])
            gx_ext = midpoints(get_extended_ticks(grid.x))
            gy_ext = midpoints(get_extended_ticks(grid.y))
            midpoints(gx_ext)*internal_length_unit, midpoints(gy_ext)*internal_length_unit, data[:,:,idx]'*punit
        end
    end

    if contours_equal_potential
        @series begin
            seriescolor := :thermal
            seriestype := :contours
            aspect_ratio --> 1
            unitformat --> :slash
            if cross_section == :x
                xguide --> "y"
                yguide --> "z"
                xlims --> (grid.y[1],grid.y[end])
                ylims --> (grid.z[1],grid.z[end])
                gy_ext = midpoints(get_extended_ticks(grid.y))
                gz_ext = midpoints(get_extended_ticks(grid.z))
                midpoints(gy_ext)*internal_length_unit, midpoints(gz_ext)*internal_length_unit, data[idx,:,:]'*punit
            elseif cross_section == :y
                xguide --> "x"
                yguide --> "z"
                xlims --> (grid.x[1],grid.x[end])
                ylims --> (grid.z[1],grid.z[end])
                gx_ext = midpoints(get_extended_ticks(grid.x))
                gz_ext = midpoints(get_extended_ticks(grid.z))
                midpoints(gx_ext)*internal_length_unit, midpoints(gz_ext)*internal_length_unit, data[:,idx,:]'*punit
            elseif cross_section == :z
                xguide --> "x"
                yguide --> "y"
                xlims --> (grid.x[1],grid.x[end])
                ylims --> (grid.y[1],grid.y[end])
                gx_ext = midpoints(get_extended_ticks(grid.x))
                gy_ext = midpoints(get_extended_ticks(grid.y))
                midpoints(gx_ext)*internal_length_unit, midpoints(gy_ext)*internal_length_unit, data[:,:,idx]'*punit
            end
        end
    end
end

@recipe function f(ϵ::DielectricDistribution{T,3,Cartesian}; x = missing, y = missing, z = missing) where {T <: SSDFloat}

    grid::CartesianGrid3D{T} = ϵ.grid
    cross_section::Symbol, idx::Int, value::T, units::Unitful.Units = get_crosssection_idx_and_value(grid, x, y, z)
    seriestype --> :heatmap
    seriescolor --> :inferno
    foreground_color_border --> nothing
    tick_direction --> :out
    title --> "Dielectric Distribution @ $(cross_section) = $(round(units, Float64(uconvert(units,value*internal_length_unit)), sigdigits=3))"
    colorbar_title --> "Dielectric Distribution"

    gx_ext::Array{T,1} = midpoints(get_extended_ticks(grid.x))
    gy_ext::Array{T,1} = midpoints(get_extended_ticks(grid.y))
    gz_ext::Array{T,1} = midpoints(get_extended_ticks(grid.z))

    @series begin
        unitformat --> :slash
        if cross_section == :x
            aspect_ratio --> 1
            xguide --> "y"
            yguide --> "z"
            xlims --> (grid.y[1],grid.y[end])
            ylims --> (grid.z[1],grid.z[end])
            gy_ext*internal_length_unit, gz_ext*internal_length_unit, ϵ.data[idx,:,:]'
        elseif cross_section == :y
            aspect_ratio --> 1
            xguide --> "x"
            yguide --> "z"
            xlims --> (grid.x[1],grid.x[end])
            ylims --> (grid.z[1],grid.z[end])
            gx_ext*internal_length_unit, gz_ext*internal_length_unit, ϵ.data[:,idx,:]'
        elseif cross_section == :z
            aspect_ratio --> 1
            xguide --> "x"
            yguide --> "y"
            xlims --> (grid.x[1],grid.x[end])
            ylims --> (grid.y[1],grid.y[end])
            gx_ext*internal_length_unit, gy_ext*internal_length_unit, ϵ.data[:,:,idx]'
        end
    end
end


text:
####################################
### METHODS FOR INFEASIBLE START ###
####################################

"""
$(SIGNATURES)
    Calculate infeasible controls to produce an infeasible state trajectory
"""
function infeasible_controls(solver::Solver,X0::Array{Float64,2},u::Array{Float64,2})
    ui = zeros(solver.model.n,solver.N-1) # initialize
    m = solver.model.m
    m̄,mm = get_num_controls(solver)
    dt = solver.dt

    x = zeros(solver.model.n,solver.N)
    x[:,1] = solver.obj.x0
    for k = 1:solver.N-1
        solver.state.minimum_time ? dt = u[m̄,k]^2 : nothing

        solver.fd(view(x,:,k+1),x[:,k],u[1:m,k], dt)

        ui[:,k] = X0[:,k+1] - x[:,k+1]
        x[:,k+1] += ui[:,k]
    end
    ui
end

function infeasible_controls(solver::Solver,X0::Array{Float64,2})
    u = zeros(solver.model.m,solver.N-1)
    if solver.state.minimum_time
        dt = get_initial_dt(solver)
        u_dt = ones(1,solver.N-1)
        u = [u; u_dt]
    end
    infeasible_controls(solver,X0,u)
end

"""
$(SIGNATURES)
    Infeasible start solution is run through time varying LQR to track state and control trajectories
"""
function get_feasible_trajectory(results0::SolverIterResults,solver::Solver)::SolverIterResults
    # Need to copy to avoid writing over original results
    results = copy(results0)

    get_feasible_trajectory!(results,solver)
    return results
end

function get_feasible_trajectory!(results::SolverIterResults,solver::Solver)::Nothing
    remove_infeasible_controls!(results,solver)

    if solver.opts.feasible_projection
        # backward pass - project infeasible trajectory into feasible space using time varying lqr
        Δv = backwardpass!(results, solver)

        # rollout
        rollout!(results,solver,0.0)
    end

    # update trajectories
    copyto!(results.X, results.X_)
    copyto!(results.U, results.U_)

    # return constrained results if input was constrained
    if !solver.state.unconstrained_original_problem
        update_constraints!(results,solver,results.X,results.U)
        update_jacobians!(results,solver)
    else
        solver.state.constrained = false
    end

    return nothing
end

"""
$(SIGNATURES)
Linear interpolation trajectory between initial and final state(s)
"""
function line_trajectory(solver::Solver, method=:trapezoid)::Array{Float64,2}
    N, = get_N(solver,method)
    line_trajectory(solver.obj.x0,solver.obj.xf,N)
end

function line_trajectory(x0::Array{Float64,1},xf::Array{Float64,1},N::Int64)::Array{Float64,2}
    x_traj = zeros(size(x0,1),N)
    t = range(0,stop=N,length=N)
    slope = (xf-x0)./N
    for i = 1:size(x0,1)
        x_traj[i,:] = slope[i].*t
    end
    x_traj
end

"""
$(SIGNATURES)
Generate a control trajectory that holds a mechanism in the configuration q
"""
function hold_trajectory(solver, mech::Mechanism, q)
    state = MechanismState(mech)
    nn = num_positions(state)
    set_configuration!(state, q[1:nn])
    vd = zero(state.q)
    u0 = dynamics_bias(state)

    n,m,N = get_sizes(solver)
    if length(q) > m
        throw(ArgumentError("system must be fully actuated to hold an arbitrary position ($(length(q)) should be > $m)"))
    end
    U0 = zeros(m,N)
    for k = 1:N
        U0[:,k] = u0
    end
    return U0
end

hold_trajectory(solver::Solver, q) = hold_trajectory(solver, solver.model.mech, q)
hold_trajectory(solver::Solver) = hold_trajectory(solver, solver.model.mech, solver.obj.x0)


text:
function isarmstrong(n::Int)
    digitbydigit = digits(n)
    len = length(digitbydigit)
    n == sum(digitbydigit.^len)
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
"""
Abstract type for Device Formulations (a.k.a Models)

# Example
```julia
import PowerSimulations
const PSI = PowerSimulations
struct MyCustomFormulation <: PSI.AbstractDeviceFormulation
```
"""
abstract type AbstractDeviceFormulation end

"""Formulation that fixes the injection values of devices"""
struct FixedOutput <: AbstractDeviceFormulation end

function _check_device_formulation(
    ::Type{D},
) where {D <: Union{AbstractDeviceFormulation, PSY.Device}}
    if !isconcretetype(D)
        throw(
            ArgumentError(
                "The device model must contain only concrete types, $(D) is an Abstract Type",
            ),
        )
    end
end

"""
Establishes the model for a particular device specified by type. Uses the keyword argument
feedforward to enable passing values between operation model at simulation time

# Arguments
-`::Type{D}`: Power System Device Type
-`::Type{B}`: Abstract Device Formulation

# Accepted Key Words
- `feedforward::Array{<:AbstractAffectFeedForward}` : use to pass parameters between models

# Example
```julia
thermal_gens = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
```
"""
mutable struct DeviceModel{D <: PSY.Device, B <: AbstractDeviceFormulation}
    component_type::Type{D}
    formulation::Type{B}
    feedforward::Union{Nothing, AbstractAffectFeedForward}
    services::Vector{ServiceModel}

    function DeviceModel(
        ::Type{D},
        ::Type{B},
        feedforward = nothing,
    ) where {D <: PSY.Device, B <: AbstractDeviceFormulation}
        _check_device_formulation(D)
        _check_device_formulation(B)
        new{D, B}(D, B, feedforward, Vector{ServiceModel}())
    end
end

get_component_type(m::DeviceModel) = m.component_type
get_formulation(m::DeviceModel) = m.formulation
get_feedforward(m::DeviceModel) = m.feedforward
get_services(m::DeviceModel) = m.services
get_services(m::Nothing) = nothing

DeviceModelForBranches = DeviceModel{<:PSY.Branch, <:AbstractDeviceFormulation}

function _set_model!(
    dict::Dict,
    model::DeviceModel{D, B},
) where {D <: PSY.Device, B <: AbstractDeviceFormulation}
    key = Symbol(D)
    if haskey(dict, key)
        @info("Overwriting $(D) existing model")
    end
    dict[key] = model
end


text:
import QuantumOps.FermiOps
import JLD2

@testset "Jordan-Wigner" begin
    s1 = PauliTerm("ZZZXIIII", complex(-1/2)) +  PauliTerm("ZZZYIIII", im/2)
    @test jordan_wigner(FermiOps.Raise, 4, 8) == s1

    s2 = PauliTerm("ZZZXIIII", complex(-1/2)) +  PauliTerm("ZZZYIIII", -im/2)
    @test jordan_wigner(FermiOps.Lower, 4, 8) == s2

    s3 = PauliTerm("IIIIIIII", 0.5 + 0.0im) +
        PauliTerm("IIIZIIII", -0.5 + 0.0im)
    @test jordan_wigner(FermiOps.NumberOp, 4, 8) == s3
    @test s3 == s1 * s2

    s4 = PauliTerm("IIIIIIII", 0.5 + 0.0im) +
        PauliTerm("IIIZIIII", 0.5 + 0.0im)
    @test jordan_wigner(FermiOps.Empty, 4, 8) == s4
    @test s4 == s2 * s1

    @test jordan_wigner_fermi(FermiOps.Raise, 4, 8) == FermiTerm("ZZZ+IIII", 1.0 + 0.0im)
    @test jordan_wigner_fermi(FermiOps.Lower, 4, 8) == FermiTerm("ZZZ-IIII", 1.0 + 0.0im)
    @test jordan_wigner_fermi(FermiOps.NumberOp, 4, 8) == FermiTerm("IIINIIII", 1.0 + 0.0im)

    h2_hamiltonian = FermiSum([FermiTerm("IIII", 0.7137539936876182),FermiTerm("IIIN", -0.47594871522096355),
                               FermiTerm("IINI", -0.47594871522096355),FermiTerm("IINN", 0.6973937674230275),
                               FermiTerm("INII", -1.2524635735648981),FermiTerm("ININ", 0.482179288212072),
                               FermiTerm("INNI", 0.663468096423568),FermiTerm("NIII", -1.2524635735648981),
                               FermiTerm("NIIN", 0.663468096423568),FermiTerm("NINI", 0.482179288212072),
                               FermiTerm("NNII", 0.6744887663568375),FermiTerm("++--", -0.181288808211496),
                               FermiTerm("+--+", 0.181288808211496),FermiTerm("-++-", 0.181288808211496),
                               FermiTerm("--++", -0.181288808211496)])

    data = JLD2.load(joinpath(@__DIR__, "iop_h2.jld2"))
    iop = data["iop"]
    @test OpSum{FermiOp}(iop) == h2_hamiltonian

    jw_h2_hamiltonian = PauliSum([PauliTerm("IIII", -0.09886396933545721 + 0.0im),PauliTerm("IIIZ", -0.2227859304041851 + 0.0im),PauliTerm("IIZI", -0.2227859304041851 + 0.0im),PauliTerm("IIZZ", 0.17434844185575687 + 0.0im),PauliTerm("IZII", 0.17119774903432972 + 0.0im),PauliTerm("IZIZ", 0.120544822053018 + 0.0im),PauliTerm("IZZI", 0.165867024105892 + 0.0im),PauliTerm("XXYY", -0.045322202052874 + 0.0im),PauliTerm("XYYX", 0.045322202052874 + 0.0im),PauliTerm("YXXY", 0.045322202052874 + 0.0im),PauliTerm("YYXX", -0.045322202052874 + 0.0im),PauliTerm("ZIII", 0.17119774903432966 + 0.0im),PauliTerm("ZIIZ", 0.165867024105892 + 0.0im),PauliTerm("ZIZI", 0.120544822053018 + 0.0im),PauliTerm("ZZII", 0.16862219158920938 + 0.0im)])

    @test jordan_wigner(h2_hamiltonian) == jw_h2_hamiltonian

    @test jordan_wigner(h2_hamiltonian, Pauli) == jordan_wigner(h2_hamiltonian, PauliI)
end


text:
using LinearForms
using Base.Test

j, = hilbertspace(:j)
m, = hilbertspace(:m)

# The form compiler does not force the
# operators and functionals to derive
# from any type. They can be any object
# as long as they're defined.
T = :T
K = :K
I = :I
e = :e
h = :h
dvg = :dvg
trc = :trc

# Build a bilinear form
L1 = @varform T[trc(m), dvg(j)]
println(L1)

#B Build a functional
R1 = @varform e[dvg(j)]
println(R1)

# Build an equation
EFIE       = @varform T[m, j]              == e[m]
MFIE_inner = @varform K[m, j] + 0.5I[m, j] == h[m]
MFIE_outer = @varform K[m, j] - 0.5I[m, j] == h[m]
println(EFIE)
println(MFIE_inner)
println(MFIE_outer)

# An example with a non trivial opstack
PL = @varform T[div(m), div(j)] == e[m]
println(PL)


text:
"""
    HalfInteger(lnℓ::AbstractFloat, p::Integer) <: RadialKernel{SqEuclidean}

The matern kernel with
``\\exp(\\ln\\ell) = \\ell > 0`` and
``p + \\frac{1}{2} = \\mathscr{v}
\\in \\Big\\lbrace \\frac{1}{2}, \\frac{3}{2}, \\ldots\\Big\\rbrace``.

```math
k_{\\ell,p}(x, y)
= k_{\\ell,p}(\\lVert x - y\\rVert_2) =
\\exp\\bigg(−\\sqrt{2\\mathscr{v}}\\frac{\\lVert x - y\\rVert_2}{\\ell}\\bigg)
\\sum^p_{i = 0}\\binom{p}{i}\\frac{(p + i)!}{(2p)!}
\\bigg(\\sqrt{8\\mathscr{v}}\\frac{\\lVert x - y\\rVert_2}{\\ell}\\bigg)^{p - i}
```

External links
* [Matérn covariance function on Wikipedia](https://en.wikipedia.org/wiki/Matérn_covariance_function)

"""
mutable struct HalfInteger{
    F<:AbstractFloat,
    Z<:Integer
} <: RadialKernel{Euclidean}
    dist::Euclidean
    lnℓ::F
    p::Z

    sqrttwo𝓋_ℓ::F

    function HalfInteger(
        dist::Euclidean,
        lnℓ::AbstractFloat,
        p::Integer
    )
        new{typeof(lnℓ), typeof(p)}(
            dist,
            lnℓ,
            p,
            sqrt(2 * p + oneunit(lnℓ)) * exp(-lnℓ),
        )
    end
end

HalfInteger(lnℓ::AbstractFloat, p::Integer) =
    HalfInteger(Euclidean(), lnℓ, p)

@inline function (k::HalfInteger)(τ::AbstractFloat)
    acc = 1
    coef = 1
    for j=k.p:-1:1
        coef *= (2 * k.sqrttwo𝓋_ℓ * τ) / (k.p + j)
        acc += binomial(k.p, j) * coef
    end
    return exp(-k.sqrttwo𝓋_ℓ * τ) * acc
end

numparams(::HalfInteger) = (1, 1)
paramtypes(::HalfInteger{F, Z}) where F where Z = (F, Z)
params(k::HalfInteger) = (lnℓ = k.lnℓ, p = k.p)
function setparams!(k::HalfInteger{F,Z}, lnℓ::F, p::Z) where F where Z
    k.lnℓ = lnℓ
    k.p = p
    k.sqrttwo𝓋_ℓ = sqrt(2 * p + oneunit(lnℓ)) * exp(-lnℓ)
end


text:
function intersect_QW(ep::Vector{Point{Float64}}, w1, w2)
    # TODO: Need to consider special cases where it is not binding on one side
    # but is on the other. It should be as simple as returning all extrema to
    # the right side of the binding extreme. Right thing to do here is break
    # this function into cases and write the cases out, that will make it
    # look much less messy

    # Number of extreme points
    npts = length(ep)
    o = Point(w1, w2)

    # Counters to determine whether points belong in the
    # set or not. vcount and hcount are for number of intersections
    # the line makes with the convex hull each way. If either are
    # even then origin is not in intersection
    checknext = true
    storecurr = false
    vcount = 0
    hcount = 0
    # Start a convex hull to store outpoints in
    outhull = []

    # Will check between each point individually
    # Points need to be ordered either clockwise
    # or counter-clockwise
    for i=2:npts+1
        # Pull out points
        if i<=npts
            x1 = ep[i-1]
            x2 = ep[i]
        else
            x1 = ep[npts]
            x2 = ep[1]
        end

        # First add point if we are supposed to store
        if storecurr
            push!(outhull, x1)
        end

        # Check whether intersects
        w1intersect = evalliney(x1, x2, o.x)
        w2intersect = evallinex(x1, x2, o.y)
        cond1 = ((x2.x - o.x)*(x1.x - o.x) < 0) && (w1intersect >= o.y)
        cond2 = ((x2.y - o.y)*(x1.y - o.y) < 0) && (w2intersect >= o.x)

        # Check first condition
        if cond1
            storecurr = !storecurr
            push!(outhull, Point(w1, w1intersect))
            vcount += 1
        # Account for if point is right on top of value
        elseif isapprox(x2.x, o.x) && isapprox(x2.y, w1intersect)
            if storecurr
                push!(outhull, x2)
            end
            storecurr = !storecurr
            vcount += 1
        end

        # Check second condition
        if cond2
            storecurr = !storecurr
            # Make sure we don't double count if both conditions hold
            push!(outhull, Point(w2intersect, w2))
            hcount += 1
        elseif isapprox(x2.y, o.y) && isapprox(x2.x, w2intersect)
            if storecurr
                push!(outhull, x2)
            end
            storecurr = !storecurr
            hcount += 1
        end
    end

    if (hcount*vcount)%2!=0 || (o.x<minx(ep)[1].x && o.y<miny(ep)[1].y)
        push!(outhull, o)
    end

    return outhull
end


text:
unique_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) =
  siteinds(uniqueind, A, B, j)

unique_siteinds(A::AbstractMPS, B::AbstractMPS, j::Integer) =
  siteinds(uniqueinds, A, B, j)

unique_siteinds(A::AbstractMPS, B::AbstractMPS) =
  siteinds(uniqueind, A, B)

common_siteind(A::AbstractMPS, B::AbstractMPS, j::Integer) =
  siteinds(commonind, A, B, j)

common_siteinds(A::AbstractMPS, B::AbstractMPS, j::Integer) =
  siteinds(commoninds, A, B, j)

common_siteinds(A::AbstractMPS, B::AbstractMPS) =
  siteinds(commoninds, A, B)

firstsiteind(M::AbstractMPS, j::Integer; kwargs...) =
  siteind(first, M, j; kwargs...)

map_linkinds!(f::Function, M::AbstractMPS) = map!(f, linkinds, M)

map_linkinds(f::Function, M::AbstractMPS) = map(f, linkinds, M)

map_common_siteinds!(f::Function, M1::AbstractMPS, M2::AbstractMPS) =
  map!(f, siteinds, commoninds, M1, M2)

map_common_siteinds(f::Function, M1::AbstractMPS, M2::AbstractMPS) =
  map(f, siteinds, commoninds, M1, M2)

map_unique_siteinds!(f::Function, M1::AbstractMPS, M2::AbstractMPS) =
  map!(f, siteinds, uniqueinds, M1, M2)

map_unique_siteinds(f::Function, M1::AbstractMPS, M2::AbstractMPS) =
  map(f, siteinds, uniqueinds, M1, M2)

for fname in (:sim, :prime, :setprime, :noprime, :addtags, :removetags,
              :replacetags, :settags)
  @eval begin
    $(Symbol(fname, :_linkinds))(M::AbstractMPS, args...; kwargs...) =
      map(i -> $fname(i, args...; kwargs...), linkinds, M)
    $(Symbol(fname, :_linkinds!))(M::AbstractMPS, args...; kwargs...) =
      map!(i -> $fname(i, args...; kwargs...), linkinds, M)
      $(Symbol(fname, :_common_siteinds))(M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...) =
       map(i -> $fname(i, args...; kwargs...), siteinds, commoninds, M1, M2)
     $(Symbol(fname, :_common_siteinds!))(M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...) =
      map!(i -> $fname(i, args...; kwargs...), siteinds, commoninds, M1, M2)
    $(Symbol(fname, :_unique_siteinds))(M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...) =
      map(i -> $fname(i, args...; kwargs...), siteinds, uniqueinds, M1, M2)
    $(Symbol(fname, :_unique_siteinds!))(M1::AbstractMPS, M2::AbstractMPS, args...; kwargs...) =
      map!(i -> $fname(i, args...; kwargs...), siteinds, uniqueinds, M1, M2)
  end
end


text:
using TargetFunction
using DrawFunc
using OptMethods
using PyCall
@pyimport matplotlib.pyplot as plt

plt.figure(TargetFunction.rosenbrock.name)
DrawFunc.draw2D(TargetFunction.rosenbrock)
plt.annotate("optimal", xy = (1.0, 1.0), xytext=(1.0, 1.0))

start_x = [-1.5; -4.0]
track = OptMethods.Fletcher_Reeves(start_x, TargetFunction.rosenbrock, 10)
plt.plot(track[1,:]', track[2,:]', "o--")
println(track[:,end])
println("num_steps: ", size(track, 2))

plt.show()

max_stars_repo_name: wyli/gradient

max_issues_repo_name: wyli/gradient

max_forks_repo_name: wyli/gradient


text:
abstract Regularizer

## l2-regularizer

immutable L2reg <: Regularizer
    w::Vector
    λ::Float64
end

function value(r::L2reg)
    norm(r.w)^2 / (2 * r.λ)
end

function gradient(r::L2reg)
    r.w / r.λ
end


text:
#define patterned stuctures

"""
    Pattern

Abstract type for structured patterns of molecules
"""
abstract type Pattern end

"""
    Nmer2D <: Pattern

N molecules symmetricaly organized around a circle with diameter d

Nmer2D(;n::Int=8, d::AbstractFloat=.1)

# Fields
- 'n': Numbor of Points = 1
- 'd': Diameter
- 'x': X position
- 'y': Y position

"""
mutable struct Nmer2D <: Pattern
    n::Int
    d::AbstractFloat
    x::Vector{AbstractFloat}
    y::Vector{AbstractFloat}
end
function Nmer2D(;n::Int=8, d::AbstractFloat=.1)

    nmer = Nmer2D(n, d, zeros(n), zeros(n))
    for nn = 1:n
        θ = 2 * pi / n * (nn - 1)
        nmer.x[nn] = d / 2 * cos(θ)
        nmer.y[nn] = d / 2 * sin(θ)
    end
    return nmer
end

"""
    Point2D <: Pattern

    A single 2D point.

Point2D() = new(1, [0.0], [0.0])

# Fields
- 'n': Numbor of Points = 1
- 'x': X position
- 'y': Y position

"""
mutable struct Point2D <: Pattern
    n::Int
    x::Vector{AbstractFloat}
    y::Vector{AbstractFloat}
    Point2D() = new(1, [0.0], [0.0])
end

"""
    Line2D <: Pattern

Points with uniform random distribution between 2 endpoints.

Line2D(;λ::AbstractFloat=10.0, endpoints=[(-1.0,0.0),(1.0,0.0)])

# Fields
- `λ`: linear molecule density
- 'endpoints': Vector of Tuple
- 'n': Numbor of Points = 1
- 'x': X position
- 'y': Y position

"""
mutable struct Line2D <: Pattern
    n::Int
    x::Vector{AbstractFloat}
    y::Vector{AbstractFloat}
    λ::AbstractFloat
    endpoints::Vector{Tuple{<:AbstractFloat,<:AbstractFloat}}
end
function Line2D(;λ::AbstractFloat=10.0, endpoints=[(-1.0,0.0),(1.0,0.0)])

    lx=(endpoints[2][1]-endpoints[1][1])
    ly=(endpoints[2][2]-endpoints[2][1])
    l=sqrt( lx^2+ ly^2)

    pois=Poisson(λ*l)
    n=rand(pois)

    line = Line2D(n, zeros(n), zeros(n),λ,endpoints)
    for nn = 1:n
        d=l*rand()
        line.x[nn] = endpoints[1][1] + d/l*lx
        line.y[nn] = endpoints[1][2] + d/l*ly
    end
    return line
end

"""
    function uniformPattern2D(ρ,xsize::AbstractFloat,ysize::AbstractFloat, p::Pattern)

Create positions of molecules from uniformly randomly placed and rotated patterns
"""
function uniform2D(ρ, p::Pattern, xsize::Real, ysize::Real)

    npatterns = rand(Poisson(xsize * ysize * ρ))
    ntotal = npatterns * p.n

    #make smd
    smd = SMLMData.SMLD2D(ntotal)
    smd.datasize = Int.(ceil.([ysize; xsize]))
    for nn = 1:npatterns
        θ = 2 * pi * rand()
        x0=rand()*xsize
        y0=rand()*ysize

        for mm = 1:p.n
            idx = (p.n) * (nn - 1) + mm
            smd.x[idx] = p.x[mm] * cos(θ) - p.y[mm] * sin(θ) + x0
            smd.y[idx] = p.x[mm] * sin(θ) + p.y[mm] * cos(θ) + y0
        end
    end

    return smd
end


text:
"""
Restricted Two-body Model tests.
"""
module R2BPTests

using LinearAlgebra, Unitful, GeneralAstrodynamics, Test

@testset verbose=false "R2BP Determination" begin

    @testset "Unitful" begin
        rᵢ = [0.0, 11681.0, 0.0] * u"km"
        vᵢ = [5.134, 4.226, 2.787] * u"km/s"
        corbit = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s")

        @test all(
            keplerian(corbit) .≈ (
                0.723452708202361,
                24509.265399338536u"km",
                151.50460766373865u"°",
                90.0u"°",
                270.0034742609256u"°",
                89.99652573907436u"°"
            )
        )

        @test CartesianOrbit(KeplerianOrbit(corbit)) ≈ corbit

        korbit = KeplerianOrbit(
            0.723452708202361,
            24509.265399338536u"km",
            151.50460766373865u"°",
            90.0u"°",
            270.0034742609256u"°",
            89.99652573907436u"°",
            Earth,
            0.0u"s"
        )

        @test KeplerianOrbit(CartesianOrbit(korbit)) ≈ KeplerianOrbit(corbit)

    end

end

@testset verbose=false "Kepler's Algorithm" begin

    @testset "Unitful" begin

        rᵢ = [0.0, 11681.0, 0.0] * u"km"
        vᵢ = [5.134, 4.226, 2.787] * u"km/s"
        corbit = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s")

        @test kepler(corbit) ≈ corbit

    end

end

@testset verbose=true "Lambert Solvers" begin

    @testset "Universal" begin

        rᵢ = [0.0, 11681.0, 0.0]u"km"
        vᵢ = [5.134, 4.226, 2.787]u"km/s"
        initial = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s")

        Δt = 1000u"s"
        final = kepler(initial, Δt; tol=1e-12)

        v₁, v₂ = lambert_universal(position_vector(initial), position_vector(final), mass_parameter(Earth), Δt; trajectory=:short, tolerance=1e-6, max_iter=1000)

        @test all(v₁ .≈ velocity_vector(initial))
        @test all(v₂ .≈ velocity_vector(final))

    end

    @testset "Oldenhuis" begin

        rᵢ = [0.0, 11681.0, 0.0]u"km"
        vᵢ = [5.134, 4.226, 2.787]u"km/s"
        initial = CartesianOrbit(rᵢ, vᵢ, Earth, 0.0u"s")

        Δt = 1000u"s"
        final = kepler(initial, Δt; tol=1e-12)

        m = 0
        v₁, v₂ = lambert(position_vector(initial), position_vector(final), Δt, m, mass_parameter(Earth))

        @test_skip all(v₁ .≈ velocity_vector(initial))
        @test_skip all(v₂ .≈ velocity_vector(final))

    end

end

end # module


text:
using Documenter, Pages

Documenter.makedocs(sitename="Pages.jl",build="")


text:
using CVChannel
using Test
using DelimitedFiles
"""
This script looks at the communication value of the
generalized Siddhu channel using SDP relaxations.
Specifically, it shows that 3-sym cv is the same
as cvPPT.
WARNING: This takes hours.
WARNING: Saving data sub-routine only works if length(s_range)>3
"""

@testset "3-symCV of Gen Siddhu not helpful" begin
    println("Here we see the 3-sym cv of the generalized Siddhu channel doesn't")
    println("improve cvPPT though we know cvPPT is loose.")
    s_range = [0:0.1:0.5;]
    μ_range = [0:0.1:1;]
    cv_table = zeros(length(s_range),length(μ_range))
    threesym_cv_table = zeros(length(s_range),length(μ_range))
    s_ctr, μ_ctr = 1,1
    for s_id in s_range
        println("---Now scanning for s=",s_id,"---")
        for μ_id in μ_range
            println("μ=",μ_id)
            genSidChan(X) = generalizedSiddhu(X,s_id,μ_id)
            gensid_chan = Choi(genSidChan,3,3)
            cv_ppt, = pptCV(gensid_chan, :dual)
            cv_3sym, = threeSymCVPrimal(gensid_chan)
            cv_table[s_ctr,μ_ctr] = cv_ppt
            threesym_cv_table[s_ctr,μ_ctr] = cv_3sym
            μ_ctr += 1
        end
        s_ctr += 1
        μ_ctr = 1
    end

    println("Here we verify that cvPPT is approximately 3symCV over the whole range.")
    diff = threesym_cv_table - cv_table
    @test all(diff -> diff < 1e-5, diff[:,:])

    info_vec = Vector{Union{Nothing,String}}(nothing, length(s_range)+1)
    label_vec = hcat("s↓ μ:",μ_range')
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by three-sym-CV-gen-Siddhu-channel.jl"
    info_vec[3] = "s_range = " * string(s_range)
    info_vec[4] = "μ_range = " * string(μ_range)
    data_to_save = hcat(vcat(label_vec,hcat(s_range,threesym_cv_table)),info_vec)

    println("Here is the 3symCV of the channel:")
    show(stdout, "text/plain", data_to_save)
    println("\nPlease name the file you'd like to write these results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
#=
# 3D Poisson, Dirichlet bc
# Uses Dendro imported as a custom gen target
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("dendrotest");

useLog("dendrotestlog", level=3)

# default values (max_depth=6, wavelet_tol = 0.1, partition_tol = 0.3, solve_tol = 1e-6, max_iters = 100)
generateFor(DENDRO, params=(8, 0.001, 0.3, 0.000001, 500))

domain(3)
functionSpace(order=2)

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, "0")

# Write the weak form
coefficient("k1", "1.1")
coefficient("k2", "0.1")
coefficient("k3", "0.1")
coefficient("k4", "0.1")
coefficient("k5", "0.1")
f = coefficient("f", "-1.1*14*pi*pi*sin(3*pi*x)*sin(2*pi*y)*sin(pi*z) + 0.1*3*pi*cos(3*pi*x)*sin(2*pi*y)*sin(pi*z) + 0.1*2*pi*sin(3*pi*x)*cos(2*pi*y)*sin(pi*z) + 0.1*pi*sin(3*pi*x)*sin(2*pi*y)*cos(pi*z) - 0.1*sin(3*pi*x)*sin(2*pi*y)*sin(pi*z)")
weakForm(u, "-k1*dot(grad(u),grad(v)) + k2*deriv(u,1)*v + k3*deriv(u,2)*v + k4*deriv(u,3)*v - k5*u*v - f*v")

build_octree_with(f);
solve(u);

finalize_finch()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# #fucntion to get argument names
# function method_argnames(m::Method)
#     argnames = ccall(:jl_uncompress_argnames, Vector{Symbol}, (Any,), m.slot_syms)
#     isempty(argnames) && return argnames
#     return argnames[1:m.nargs]
# end

# function assert_assembly_callbacks(biomass_func::Function, new_parameters::Function)
#     #assert functions have correct args and return types
#     @assert Base.return_types(biomass_func)[1] == Float64 "`biomass_func` must return biomass of new consumer as a Float64)"
#     @assert method_argnames(collect(methods(biomass_func))[1]) == [Symbol("#self#"), :integrator] "`biomass_func` must take only `integrator` as argument"

#     @assert Base.return_types(new_parameters)[1] == Tuple{Float64,Any} "`new_parameters` must return new uptake matrix and maintenance cost as a `Tuple{Any,Float64}`"
#     @assert method_argnames(collect(methods(new_parameters))[1]) == [Symbol("#self#"), :integrator] "`new_parameters` must take only `integrator` as argument"
# end

# """
#     add_consumer!(integrator, biomass_func::Function, new_parameters::Function)

# Function to add consumer to system within a callback. Takes two functions `biomass_func` and `new_parameters` which give the biomass of the new consumer and the updated parameters respectively. Both functions should take `integrator` as the argument and return either a single `Float64` for biomass or a new `u` array and `Rm` value for new parameters.
# """
# function add_consumer!(integrator, biomass_func::Function, new_parameters::Function)

#     #find extinct species
#     extant = findall(integrator.u[integrator.p.M+1 : end] .> 1e-5)
#     # println(extant)
#     # println(integrator.u)
#     # println(minimum(integrator.u[integrator.p.M+1 : end]))
#     #move biomasses
#     integrator.u[(integrator.p.M+1) : (integrator.p.M+length(extant))] .= integrator.u[extant .+ integrator.p.M]
#     # #resize cache to add new species ....
#     resize!(integrator, length(extant) + integrator.p.M + 1)
#     integrator.u[end] = biomass_func(integrator)
#     # println(integrator.u)
#     #alter parameters
#     new_p = new_parameters(integrator)
#     #N - get from extant species
#     integrator.p.N = length(extant) + 1

#     #Rm - resize to add new and then remove
#     integrator.p.Rm[1:length(extant)] .= integrator.p.Rm[extant]
#     resize!(integrator.p.Rm, integrator.p.N)
#     integrator.p.Rm[end] = new_p[1]

#     # #U
#     new_u = Array{Float64,2}(undef, integrator.p.N, integrator.p.M)
#     new_u[1:(end-1), :] .= integrator.p.u[extant,:]
#     new_u[end,:] .= new_p[2]
#     integrator.p.u = new_u
# end

# """
#      consumer_equilibrium(integrator, tmax::Float64, abstol::Float64, reltol::Float64)

# function to detect equilbirium based on an absolute tolerance `abstol` and relative tolerance `reltol`. Will not detect equilibrium past a given time `tmax`
# """
# function consumer_equilibrium(integrator, tmax::Float64, abstol::Float64, reltol::Float64)
#     #test time - dont return true when over Tmax
#     # if integrator.t > tmax
#     #    return(false)
#     # end

#     #get derivative - from DiffEqCallbacks TerminateSteadyState
#     if DiffEqBase.isinplace(integrator.sol.prob)
#         testval = first(DiffEqBase.get_tmp_cache(integrator))
#         DiffEqBase.get_du!(testval, integrator)
#         if typeof(integrator.sol.prob) <: DiffEqBase.DiscreteProblem
#             @. testval =  testval - integrator.u
#         end
#     else
#         testval = get_du(integrator)
#         if typeof(integrator.sol.prob) <: DiffEqBase.DiscreteProblem
#             testval =  testval - integrator.u
#         end
#     end
#     #set resource du to 0

#     testval[1:integrator.p.M] .= 0.0
#     # println(testval)

#     any(abs(d) > abstol && abs(d) > reltol*abs(u) for (d,abstol, reltol, u) = zip(testval, Iterators.cycle(abstol), Iterators.cycle(reltol), integrator.u)) && (return false)

#     return(true)

# end

# """

# """
# function add_at_equilibrium(biomass_func::Function, new_parameters::Function; tmax::Float64 = Inf, abstol::Float64 = 1e-8, reltol::Float64 = 1e-6)
#     # #warn if tolerances are high
#     # if abstol > 1e-2
#     #     @warn "abstol is set greater than `1e-2` (`$abstol`). Are you sure you dont want tighter tolerances?"
#     # elseif reltol > 1e-2
#     #     @warn "reltol is set greater than `1e-2` (`$reltol`). Are you sure you dont want tighter tolerances?"
#     # end

#     #assert functions have correct args and return types
#     assert_assembly_callbacks(biomass_func, new_parameters)

#     #define condtion and affect!
#     condition = (u, t, integrator) -> consumer_equilibrium(integrator, tmax, abstol, reltol)
#     affect! = (integrator) -> add_consumer!(integrator, biomass_func, new_parameters)
#     return(:equilibrium , DiffEqCallbacks.DiscreteCallback(condition, affect!; save_positions = (true, true)))
# end

# function add_at_t(biomass_func::Function, new_parameters::Function, times::Vector{Float64})
#     #assert functions have correct args and return types
#     assert_assembly_callbacks(biomass_func, new_parameters)

#     affect! = (integrator) -> add_consumer!(integrator, biomass_func, new_parameters)
#     return(:at_time_stop, DiffEqCallbacks.PresetTimeCallback(times, affect!), times)
# end


text:
# This file is a part of BitArraynge.jl, licensed under the MIT License (MIT).
println("building...")


text:
using ANYmalCRobot
using MeshCat
using MeshCatMechanisms
using RigidBodyDynamics

# Initialize anymal mechanism and MeshCat visualizer
anymal = ANYmalCRobot.mechanism()
vis = Visualizer()
mvis = MechanismVisualizer(
    anymal,
    URDFVisuals(ANYmalCRobot.urdfpath(), package_path=[ANYmalCRobot.packagepath()]),
    vis)

# set_nominal!() is a helper function to put the robot into a default config
anymal_state = MechanismState(anymal)
ANYmalCRobot.set_nominal!(anymal_state)
set_configuration!(mvis, configuration(anymal_state))

# Helper functions to visualise joint frames and contacts
# Try running visualize_frames() with optional kwarg "show_all_frames = true" to
# show all frames, by default only revolute are shown!
ANYmalCRobot.visualize_frames(anymal, mvis)
ANYmalCRobot.visualize_contacts(anymal, mvis)


text:
module ProximalAlgorithms

using AbstractOperators
using ProximalOperators
using LinearAlgebra
using Printf

import Base: iterate

include("utilities/identity.jl")
include("utilities/zero.jl")
include("utilities/conjugate.jl")
include("utilities/tomove.jl")

abstract type ProximalAlgorithm{I,T} end

# The following methods give `ProximalAlgorithm` objects the iterable behavior.

function iterate(solver::ProximalAlgorithm{I,T}) where {I,T}
    point = initialize!(solver)
    return (point, one(I))
end

function iterate(solver::ProximalAlgorithm{I,T}, it::I) where {I,T}
    if done(solver, it)
        return nothing
    else
        point = iterate!(solver, it)
        return (point, it + one(I))
    end
end

function done(solver::ProximalAlgorithm{I,T}, it::I)::Bool where {I,T}
    return it >= maxit(solver) || converged(solver, it)
end

# Running a `ProximalAlgorithm` unrolls the iterations

function run!(solver::ProximalAlgorithm{I,T})::Tuple{I,T} where {I,T}
    local it, point
    if verbose(solver) display(solver) end
    # NOTE: the following loop is translated into:
    #   it = start(solver)
    #   while !done(solver, it)
    #       (point, it) = next(solver, it)
    #       [...]
    #   end
    # See: https://docs.julialang.org/en/stable/manual/interfaces
    for (it_, point_) in enumerate(solver)
        if verbose(solver, it_) display(solver, it_) end
        it, point = it_, point_
    end
    if verbose(solver) display(solver, it) end
    return it, point
end

# Functions `verbose` and `display` are used for inspecting the iterations.
# Here we provide their default behavior (no output).

verbose(sol::ProximalAlgorithm) = false
verbose(sol::ProximalAlgorithm, it) = false
function display(sol::ProximalAlgorithm) end
function display(sol::ProximalAlgorithm, it) end

# It remains to define what concrete ProximalAlgorithm types are and how
# `initialize!`, `iterate!`, `maxit`, `converged` work for each specific solver.
# This is done in the following included files.

include("algorithms/ForwardBackward.jl")
include("algorithms/ZeroFPR.jl")
include("algorithms/PANOC.jl")
include("algorithms/DouglasRachford.jl")
include("algorithms/AsymmetricForwardBackwardAdjoint.jl")

# include("algorithms/VuCondat.jl")
# include("algorithms/ChambollePock.jl")
# include("algorithms/DavisYin.jl")

# The following template can be copy-pasted to implement new algorithms.

include("template/Template.jl")

end # module


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Compat, Random
include("convexpenalties.jl")

results_template = DataFrame(n=Int[],k=Int[], theSeed=Int[], gap_persp=Real[], LB_persp=Real[], UB_persp=Real[], UB_BM=Real[], t_persp=Real[], t_BM=Real[],
mse_persp=Real[], mse_BM=Real[], acc_persp=Real[], acc_BM=Real[])

results_fact=similar(results_template, 0)

for ARG in ["5"]

array_num = parse(Int, ARG)
seeds=1:20

n=50
k_true=10
ks=2:2:40
gamma=1e6

theSeed=seeds[(array_num-1)÷20+1]
k=ks[(array_num-1)%20+1]
for theSeedI in (theSeed*10):(theSeed*10+9)
    Random.seed!(theSeedI)
    U=rand(n,k_true)
    E=0.05*sqrt(k)*randn(n,n)
    A=U*U'.+E
    @show t_persp=@elapsed X_persp, LB_persp, UB_persp, U_persp=solveNNMFRelax(A, k)
    @show gap_persp=(UB_persp-LB_persp)/UB_persp
    mse_persp=norm(A.-X_persp)/norm(A)
    acc_persp=norm(U*U'.-X_persp)/norm(U*U')

    # Solve using BM with random initialization

    t_BM=@elapsed X_BM, UB_BM=solveBM(A, k, rand(n,k))
    mse_BM=norm(A.-X_BM)/norm(A)
    acc_BM=norm(U*U'.-X_BM)/norm(U*U')

    # # Solve using BM with rounding initialization
    #
    t_BM=@elapsed X_BM, UB_BM=solveBM(A, k, U_persp)
    mse_BM=norm(A.-X_BM)/norm(A)
    acc_BM=norm(U*U'.-X_BM)/norm(U*U')
    # Write results

    push!(results_fact, [n, k, theSeedI, gap_persp, LB_persp, UB_persp, UB_BM, t_persp, t_BM,
    mse_persp, mse_BM, acc_persp, acc_BM])
end

CSV.write("createFig1_times_acc.csv", results_fact, append=true)

end

max_stars_repo_path: Matrix Factorization/createFigure1.jl

max_issues_repo_path: Matrix Factorization/createFigure1.jl

max_forks_repo_path: Matrix Factorization/createFigure1.jl