text:
"""
    SumBasis(b1, b2...)

Similar to [`CompositeBasis`](@ref) but for the [`directsum`](@ref) (⊕)
"""
struct SumBasis{S,B} <: Basis
    shape::S
    bases::B
end
SumBasis(bases) = SumBasis(Int[length(b) for b in bases], bases)
SumBasis(shape, bases::Vector) = (tmp = (bases...,); SumBasis(shape, tmp))
SumBasis(bases::Vector) = SumBasis((bases...,))
SumBasis(bases::Basis...) = SumBasis((bases...,))

==(b1::T, b2::T) where T<:SumBasis = equal_shape(b1.shape, b2.shape)
==(b1::SumBasis, b2::SumBasis) = false
length(b::SumBasis) = sum(b.shape)

"""
    directsum(b1::Basis, b2::Basis)

Construct the [`SumBasis`](@ref) out of two sub-bases.
"""
directsum(b1::Basis, b2::Basis) = SumBasis(Int[length(b1); length(b2)], Basis[b1, b2])
directsum(b::Basis) = b
directsum(b::Basis...) = reduce(directsum, b)
function directsum(b1::SumBasis, b2::Basis)
    shape = [b1.shape;length(b2)]
    bases = [b1.bases...;b2]
    return SumBasis(shape, (bases...,))
end
function directsum(b1::Basis, b2::SumBasis)
    shape = [length(b1);b2.shape]
    bases = [b1;b2.bases...]
    return SumBasis(shape, (bases...,))
end
function directsum(b1::SumBasis, b2::SumBasis)
    shape = [b1.shape;b2.shape]
    bases = [b1.bases...;b2.bases...]
    return SumBasis(shape, (bases...,))
end

const ⊕ = directsum

"""
    directsum(x::Ket, y::Ket)

Construct a spinor via the [`directsum`](@ref) of two [`Ket`](@ref)s.
The result is a [`Ket`](@ref) with data given by `[x.data;y.data]` and its
basis given by the corresponding [`SumBasis`](@ref).
**NOTE**: The resulting state is not normalized!
"""
directsum(x::Ket, y::Ket) = Ket(directsum(x.basis, y.basis), [x.data; y.data])
directsum(x::StateVector...) = reduce(directsum, x)

"""
    getblock(x::Ket{<:SumBasis}, i)

For a [`Ket`](@ref) defined on a [`SumBasis`](@ref), get the state as it is defined
on the ith sub-basis.
"""
function getblock(x::Ket{B}, i) where B<:SumBasis
    b_i = x.basis.bases[i]
    inds = cumsum([0;length.(x.basis.bases[1:i])...])
    return Ket(b_i, x.data[inds[i]+1:inds[i+1]])
end

"""
    setblock!(x::Ket{<:SumBasis}, val::Ket, i)

Set the data of `x` on the ith sub-basis equal to the data of `val`.
"""
function setblock!(x::Ket{B}, val::Ket, i) where B<:SumBasis
    check_samebases(x.basis.bases[i], val)
    inds = cumsum([0;length.(x.basis.bases[1:i])...])
    x.data[inds[i]+1:inds[i+1]].data .= val.data
    return x
end

"""
    directsum(x::DataOperator, y::DataOperator)

Compute the direct sum of two operators. The result is an operator on the
corresponding [`SumBasis`](@ref).
"""
function directsum(a::DataOperator, b::DataOperator)
    dType = promote_type(eltype(a),eltype(b))
    data = zeros(dType, size(a,1)+size(b,1), size(a,2)+size(b,2))
    data[1:size(a,1),1:size(a,2)] = a.data
    data[size(a,1)+1:end, size(a,2)+1:end] = b.data
    return Operator(directsum(a.basis_l, b.basis_l), directsum(a.basis_r, b.basis_r), data)
end
function directsum(a::SparseOpType, b::SparseOpType)
    dType = promote_type(eltype(a),eltype(b))
    data = spzeros(dType, size(a,1)+size(b,1), size(a,2)+size(b,2))
    data[1:size(a,1),1:size(a,2)] = a.data
    data[size(a,1)+1:end, size(a,2)+1:end] = b.data
    return Operator(directsum(a.basis_l, b.basis_l), directsum(a.basis_r, b.basis_r), data)
end
directsum(a::AbstractOperator...) = reduce(directsum, a)

"""
    setblock!(op::DataOperator{<:SumBasis,<:SumBasis}, val::DataOperator, i, j)

Set the data of `op` corresponding to the block `(i,j)` equal to the data of `val`.
"""
function setblock!(op::DataOperator{<:SumBasis,<:SumBasis}, val::DataOperator, i, j)
    (bases_l,bases_r) = op.basis_l.bases, op.basis_r.bases
    check_samebases(bases_l[i], val.basis_l)
    check_samebases(bases_r[j], val.basis_r)
    inds_i = cumsum([0;length.(bases_l[1:i])...])
    inds_j = cumsum([0;length.(bases_r[1:j])...])
    op.data[inds_i[i]+1:inds_i[i+1],inds_j[j]+1:inds_j[j+1]] = val.data
    return op
end

"""
    getblock(op::Operator{<:SumBasis,<:SumBasis}, i, j)

Get the sub-basis operator corresponding to the block `(i,j)` of `op`.
"""
function getblock(op::DataOperator{BL,BR}, i, j) where {BL<:SumBasis,BR<:SumBasis}
    (bases_l,bases_r) = op.basis_l.bases, op.basis_r.bases
    inds_i = cumsum([0;length.(bases_l[1:i])...])
    inds_j = cumsum([0;length.(bases_r[1:j])...])
    data = op.data[inds_i[i]+1:inds_i[i+1],inds_j[j]+1:inds_j[j+1]]
    return Operator(bases_l[i],bases_r[j], data)
end

"""
    embed(basis_l::SumBasis, basis_r::SumBasis,
               index::Integer, operator)

Embed an operator defined on a single subspace specified by the `index` into
a [`SumBasis`](@ref).
"""
function embed(basis_l::SumBasis, basis_r::SumBasis,
               index::Integer, op::T) where T<:DataOperator
    @assert length(basis_r.bases) == length(basis_l.bases)

    basis_l.bases[index] == op.basis_l || throw(IncompatibleBases())
    basis_r.bases[index] == op.basis_r || throw(IncompatibleBases())

    embedded_op = SparseOperator(eltype(op), basis_l, basis_r)
    setblock!(embedded_op, op, index, index)
    return embedded_op
end

"""
    embed(basis_l::SumBasis, basis_r::SumBasis,
                indices, operator)

Embed an operator defined on multiple subspaces specified by the `indices` into
a [`SumBasis`](@ref).
"""
function embed(basis_l::SumBasis, basis_r::SumBasis,
               indices, op::T) where T<:DataOperator
    @assert length(basis_r.bases) == length(basis_l.bases)

    embedded_op = SparseOperator(eltype(op), basis_l, basis_r)
    for i=1:length(indices), j=1:length(indices)
        op_ = getblock(op, i, j)
        setblock!(embedded_op, op_, indices[i], indices[j])
    end
    return embedded_op
end

"""
    embed(basis_l::SumBasis, basis_r::SumBasis,
               indices, operators)

Embed a list of operators on subspaces specified by the `indices` into a
[`SumBasis`](@ref).
"""
function embed(basis_l::SumBasis, basis_r::SumBasis,
               indices, ops::Union{Tuple{Vararg{<:DataOperator}},Vector{<:DataOperator}})
    @assert length(basis_r.bases) == length(basis_l.bases)

    T = mapreduce(eltype, promote_type, ops)
    embedded_op = SparseOperator(T, basis_l, basis_r)
    for k=1:length(ops)
        op = ops[k]
        idx = indices[k]
        if length(idx)==1
            setblock!(embedded_op, op, idx[1], idx[1])
        else
            for i=1:length(idx), j=1:length(idx)
                op_ = getblock(op, i, j)
                setblock!(embedded_op, op_, idx[i], idx[j])
            end
        end
    end
    return embedded_op
end

embed(b::SumBasis, indices, ops) = embed(b, b, indices, ops)

"""
    LazyDirectSum <: AbstractOperator

Lazy implementation of `directsum`
"""
mutable struct LazyDirectSum{BL,BR,T} <: AbstractOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    operators::T
end

# Methods
LazyDirectSum(op1::AbstractOperator, op2::AbstractOperator) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1,op2))
LazyDirectSum(op1::LazyDirectSum, op2::AbstractOperator) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1.operators...,op2))
LazyDirectSum(op1::AbstractOperator, op2::LazyDirectSum) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1,op2.operators...))
LazyDirectSum(op1::LazyDirectSum, op2::LazyDirectSum) = LazyDirectSum(directsum(op1.basis_l,op2.basis_l),directsum(op1.basis_r,op2.basis_r),(op1.operators...,op2.operators...))
LazyDirectSum(op::AbstractOperator...) = reduce(LazyDirectSum, op)

function Base.:(==)(op1::LazyDirectSum,op2::LazyDirectSum)
    (op1.basis_l == op2.basis_l && op1.basis_r == op2.basis_r) || return false
    length(op1.operators)==length(op2.operators) || return false
    for (o1,o2) = zip(op1.operators,op2.operators)
        (o1 == o2) || return false
    end
    return true
end

# Algebra
dense(x::LazyDirectSum) = directsum(dense.(x.operators)...)
*(op1::LazyDirectSum{B1,B2},op2::LazyDirectSum{B2,B3}) where {B1<:Basis,B2<:Basis,B3<:Basis}= LazyDirectSum(op1.basis_l,op2.basis_r,Tuple(op1.operators[i]*op2.operators[i] for i=1:length(op1.operators)))
-(op::LazyDirectSum) = LazyDirectSum([-op.operators[1];op.operators[2:end]]...)
+(op1::LazyDirectSum{B1,B2},op2::LazyDirectSum{B1,B2}) where {B1<:Basis,B2<:Basis} = LazyDirectSum((op1.operators .+ op2.operators)...)
dagger(op::LazyDirectSum) = LazyDirectSum(op.basis_r, op.basis_l, dagger.(op.operators))
Base.eltype(x::LazyDirectSum) = mapreduce(eltype, promote_type, x.operators)

directsum(op1::AbstractOperator,op2::LazyDirectSum) = LazyDirectSum(op1,op2)
directsum(op1::LazyDirectSum,op2::AbstractOperator) = LazyDirectSum(op1,op2)
directsum(op1::LazyDirectSum,op2::LazyDirectSum) = LazyDirectSum(op1,op2)

# Use lazy sum for FFTOperator
transform(b1::SumBasis,b2::SumBasis; kwargs...) = LazyDirectSum([transform(b1.bases[i],b2.bases[i];kwargs...) for i=1:length(b1.bases)]...)
directsum(op1::FFTOperator, op2::FFTOperator) = LazyDirectSum(op1,op2)

# Lazy embed
function embed(basis_l::SumBasis, basis_r::SumBasis, indices, op::LazyDirectSum)
    bases_l = basis_l.bases
    bases_r = basis_r.bases
    N = length(bases_r)
    @assert length(bases_l)==N

    T = eltype(op)

    function _get_op(i)
        idx = findfirst(isequal(i), indices)
        if idx === nothing
            return SparseOperator(T, bases_l[i], bases_r[i])
        else
            return op.operators[idx]
        end
    end

    ops = Tuple(_get_op(i) for i=1:N)
    return LazyDirectSum(basis_l,basis_r,ops)
end
# TODO: embed for multiple LazyDirectums?

# Fast in-place multiplication
function mul!(result::Ket{B1},M::LazyDirectSum{B1,B2},b::Ket{B2},alpha_,beta_) where {B1,B2}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    bases_r = b.basis.bases
    bases_l = M.basis_l.bases
    index = cumsum([0;length.(bases_r)...])
    for i=1:length(index)-1
        tmpket = Ket(bases_r[i], b.data[index[i]+1:index[i+1]])
        tmpresult = Ket(bases_l[i], result.data[index[i]+1:index[i+1]])
        mul!(tmpresult,M.operators[i],tmpket,alpha,beta)
        result.data[index[i]+1:index[i+1]] = tmpresult.data
    end
    return result
end
function mul!(result::Bra{B2},b::Bra{B1},M::LazyDirectSum{B1,B2},alpha_,beta_) where {B1,B2}
    alpha = convert(ComplexF64, alpha_)
    beta = convert(ComplexF64, beta_)
    bases_l = b.basis.bases
    bases_r = M.basis_r.bases
    index = cumsum([0;length.(bases_r)...])
    for i=1:length(index)-1
        tmpbra = Bra(bases_l[i], b.data[index[i]+1:index[i+1]])
        tmpresult = Bra(bases_r[i], result.data[index[i]+1:index[i+1]])
        mul!(tmpresult,tmpbra,M.operators[i],alpha,beta)
        result.data[index[i]+1:index[i+1]] = tmpresult.data
    end
    return result
end


text:
###############################################################
# Part of BDisposal.jl package
###############################################################

"""
    efficiencyScores(gI,gO,bO,bI;retToScale,prodStructure,dirGI,dirBI,dirGO,dirBO,startθ,startμ,startλ)

Compute efficiency indicators and convexity test results considering the joined production of good and bad outputs.

Given a set of measures of inputs, "good" ("desiderable") and "bad" ("undesiderable") outputs for different decision making units, compute their
distance to the production frontier, i.e. their degree of efficiency.

## Parameters:
- Positional
  *  `gI`: "Good" inputs (3D matrix by DMUs, input items and periods)
  *  `gO`: "Good" outputs (3D matrix by DMUs, input items and periods)
  *  `bO`: "Bad" outputs (3D matrix by DMUs, input items and periods)
  *  `bI`: "Bad" inputs (optional 3D matrix by DMUs, input items and periods) [default: empty array]
- Keyword
  *  `retToScale`: Wheter the returns to scale should be assumed "constant" or "variable" (default). Non-convex distance and test is computed
      only under the "variable" assumption, except for the `(0,0,1,-1)` distance under multiplicative production function
  *  `prodStructure`: Wheter the production structure should be assumed "additive" (default) or "multiplicative"
  *  `dirGI`,`dirBI`,`dirGO`,`dirBO`: The directions toward where to measure the efficiency (see notes) [default: `(0,0,1,0)`]
  *  `startθ`,`startμ`,`startλ`: Initial values in the convex optimisation problem [default: `(0,0,1.1)`]

## Returns:
- A named tuple with the following items:
  * `λs`: The efficiency indicators coherent with the convexity test (2D matrix by DMUs and periods)
  * `λs_convex`: The efficiency indicators assuming a convex production frontier (2D matrix by DMUs and periods)
  * `λs_nonconvex`:  The efficiency indicators assuming a non-convex production frontier (2D matrix by DMUs and periods)
  * `nonConvTest`: The result of the non-convexity test (2D matrix of boolean by DMUs and periods). Note that "true" here refers to non-convex.
  * `nonConvTest_value`: The value of the non-convexity test (2D matrix by DMUs and periods)

## Description of the function

`efficiencyScores()` defines environmental efficiency indicators for a set
of Decision Making Units (ie., observations). The nature of the DMUs (eg. firms,
countries etc.) depends on the aim of the environmental efficiency analysis. All the
DMUs are characterized by a set of inputs and outputs separated into undesirable (eg.,
polluting) and desirable (eg., no polluting) ones. Note that, inputs separation into
desirable (eg., non emissions-causing) and undesirable (eg., emissions-causing) compo-
nents is optional for the BDisposal package. The BDisposal DEA modelling defines the
smallest environmental production process that contains all DMUs informations (in-
puts/outputs) and satisfies the B-disposal axiomatic pattern.
The BDisposal algorithms for efficiency assessment compute the distance to the best
environmental production procedures; ie., to the efficient production frontier. These
efficiency indices inherit the structure of additive and multiplicative distance functions.
Notice that the BDisposal package computes general shape of additive and multiplicative distance functions. It follows that,
the Bdisposal package introduces a flexible framework for environmental efficiency as-
sessment, with the possibility to analyse different directions for the distance functions.
In addition, the BDisposal package allows to test the (economists’) properties of con-
vexity and disposability.

## Interpretation of the results

### Efficiency indicators

When the additive shape of the distance function is greater than 0 then, the produc-
tion unit doesn’t operate efficiently. If the additive efficiency score is equal to 0 then,
the DMU is efficient. It follows that, the DMU employs one of the best production
technology to transform inputs into outputs; ie., the producer is a benchmark for the
other observations. A similar reasoning holds for the multiplicative distance functions.
If the multiplicative efficiency index is greater than 1 then, the DMU doesn’t perform
efficiently. When the multiplicative efficiency score is equal to 1 then, the production
unit is a benchmark for the other observations.

### Convexity test

When the multiplicative test of convexity is greater than 1
(respectively, equal to 1) then, the production process frontier is non convex (respec-
tively, convex). If the multiplicative test of disposability is greater than 1 (respectively,
equal to 1) then, the frontier of the production set displays B-disposability; ie., there
exists costs for undesirable component decrease. A similar reasoning applies for the
additive framework. The informations provided by these tests define the shape of the
production process boundaries.

## Notes:
* Directions toward where to measure the efficiency distance can be tuned using the `dirGI`,`dirBI`,`dirGO` and `dirBO` parameters.
* The following directions are supported:
  - multiplicative pr. struct.: (-1,0,0,0), (0,-1,0,0), (0,0,1,0), (0,0,0,-1) or (0,0,1,-1)
  - addittive pr. struct.: (1,0,0,0), (0,1,0,0), (0,0,1,0) or (0,0,0,1)
  Other directions can be used to compute the direction under the convex frontier assumption, but the convexity test is not performed and,
  for the multiplicative case, no guarantee is given on solving the underliying (non-linear) problem (different initial value
  `startθ`,`startμ`,`startλ` can be attempted)

"""
function efficiencyScores(gI::Array{Float64,3},gO::Array{Float64,3},bO::Array{Float64,3},bI::Array{Float64,3}=Array{Float64}(undef, (size(gI,1),0,size(gI,3)));
                          retToScale="variable",prodStructure="additive",
                          dirGI=0,dirBI=0,dirGO=1,dirBO=0,
                          startθ=0,startμ=0,startλ=1.1)

    # Data processing
    nDMUs, nPer, nGI, nBI, nGO, nBO = size(gI,1), size(gI,3), size(gI,2), size(bI,2), size(gO,2), size(bO,2)
    λs_convex    = Array{Union{Missing,Float64,String}}(undef,nDMUs,nPer) # Output matrix hosting the results (nDMUs,nPer)
    λs_nonconvex = Array{Union{Missing,Float64,String}}(undef,nDMUs,nPer) # Output matrix hosting the results (nDMUs,nPer)
    directions   = (dirGI,dirBI,dirGO,dirBO)

    # Tests if doing also convex analysis or not
    doNonConvexAnalysis = true

    if retToScale == "constant" && directions != (0,0,1,-1)
        @warn "Non-convex analysis is supported only with variable returns to scale or with direcitons (0,0,1,-1) and multiplicative production structure."
        doNonConvexAnalysis = false
    end

    if  prodStructure == "multiplicative"
        if !( directions in [(-1,0,0,0),(0,-1,0,0),(0,0,1,0),(0,0,0,-1),(0,0,1,-1)])
            @warn "Non-convex analysis not supported with these directions. For the multiplicative case,
            non-convex analysis is supported only with directions (-1,0,0,0),(0,-1,0,0),(0,0,1,0),(0,0,0,-1) or (0,0,1,-1)"
             doNonConvexAnalysis = false
        end
    else
        if !( directions in [(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)])
            @warn "Non-convex analysis not supported with these directions. For the addittive case,
            non-convex analysis is supported only with directions (1,0,0,0),(0,1,0,0),(0,0,1,0) or (0,0,0,1)"
             doNonConvexAnalysis = false
        end
    end
    # temp todo:
    doNonConvexAnalysis = true

    # Loping over the periods
    for t in 1:nPer
        gIₜ = gI[:,:,t]; gOₜ = gO[:,:,t]; bOₜ = bO[:,:,t]
        #bIₜ = nBI > 0 ? bI[:,:,t] : Array{Float64}(undef, 0,0)
        bIₜ = bI[:,:,t]

        #println("Period: $t")
        # Solving for each dmu
        for z in 1:nDMUs
            gIₜ₀ = gIₜ[z,:]; gOₜ₀ = gOₜ[z,:]; bOₜ₀ = bOₜ[z,:]
            #bIₜ₀ = nBI > 0 ? bIₜ[z,:] : Float64[]
            bIₜ₀ = bIₜ[z,:]

        #    println("DMU: $z")

            # CONVEX analysis....
            linearModelPossibleDirections = [(-1,0,0,0),(0,-1,0,0),(0,0,1,0),(0,0,0,-1)]
            forceLinearModel = (dirGI,dirBI,dirGO,dirBO) in linearModelPossibleDirections ?  true : false

            λs_convex[z,t] = problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ;
                                     retToScale=retToScale,prodStructure=prodStructure,
                                     convexAssumption=true,
                                     directions=(dirGI,dirBI,dirGO,dirBO),
                                     startValues=(startθ,startμ,startλ),
                                     forceLinearModel=forceLinearModel)

            # Non-convex analysis (when implemented)...
            λs_nonconvex[z,t] = doNonConvexAnalysis ? problem(gIₜ₀,bIₜ₀,gOₜ₀,bOₜ₀,gIₜ,bIₜ,gOₜ,bOₜ;
                                        retToScale=retToScale,prodStructure=prodStructure,
                                        directions=(dirGI,dirBI,dirGO,dirBO),
                                        convexAssumption=false) : missing
        end # end of each dmu
    end # endof each period

#    println("Done all DMU problems")
    nonConvTest_value = λs_nonconvex ./ λs_convex
    nonConvTest       = nonConvTest_value .< (1.0 - eps())
    λs                = doNonConvexAnalysis ? min.(λs_nonconvex,λs_convex) : λs_convex

    return λs, λs_convex, λs_nonconvex, nonConvTest, nonConvTest_value
end

ext: jl


text:
# We output a message after loading each package to work around a
# ten-minute timeout limit on Travis. Travis assumes the tests have hung
# if the interval between printing something to stdio exceeds ten minutes.
using Images
@info "Finished loading Images package."
using TestImages
@info "Finished loading TestImages package."
using StaticArrays
@info "Finished loading StaticArrays package."
using OffsetArrays
@info "Finished loading OffsetArrays package."
using Random
@info "Finished loading Random package."
using CoordinateTransformations
@info "Finished loading CoordinateTransformations package."

function evaluate_error(dims, flow::Array{SVector{2, Float64}, 1},  Δ,  tol)
    error_count = 0
    maximum_error = 0.0
    for i in eachindex(flow)
        δ = flow[i]
        if abs(δ[1] - Δ[1]) > tol || abs(δ[2] - Δ[2]) > tol
                error_count += 1
        end
        error = sqrt( (δ[1] - Δ[1])^2 + (δ[2] - Δ[2])^2 )
        if error > maximum_error
            maximum_error = error
        end
    end
    return error_count, maximum_error
end

@testset "Lucas-Kanade" begin
    @testset "Horizontal Motion" begin
        @info "Running Horizontal Motion test."
        maximum_percentage_error = 7.5
        number_test_pts = 500
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (0.0, 3.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        corners = imcorner(img1, method=shi_tomasi)
        I = findall(!iszero, corners)
        r, c = (getindex.(I, 1), getindex.(I, 2))
        points = map((ri, ci) -> SVector{2}(Float64(ri), Float64(ci)), r, c)
        Random.seed!(9876)
        points = rand(points, (number_test_pts,))
        flow, status_array = optical_flow(img1, img2, points, LucasKanade(20, 11, 4, 0.000001))

        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100

        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 1}(undef,length(points))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        flow, status_array = optical_flow(img1, img2, points, displacement, LucasKanade(20, 11, 4, 0.000001))

        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100

        @info "Case: Horizontal Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
    end

    @testset "Vertical Motion" begin
        @info "Running Vertical Motion test."
        maximum_percentage_error = 7.5
        number_test_pts = 500
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (3.0, 0.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        corners = imcorner(img1, method=shi_tomasi)
        I = findall(!iszero, corners)
        r, c = (getindex.(I, 1), getindex.(I, 2))
        points = map((ri, ci) -> SVector{2}(Float64(ri), Float64(ci)), r, c)
        Random.seed!(9876)
        points = rand(points, (number_test_pts,))
        flow, status_array = optical_flow(img1, img2, points, LucasKanade(20, 11, 4, 0.000001))

        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100

        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 1}(undef,length(points))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        flow, status_array = optical_flow(img1, img2, points, displacement, LucasKanade(20, 11, 4, 0.000001))

        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100

        @info "Case: Vertical Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
    end

    @testset "Combined Motion" begin
        @info "Running Combined Motion test."
        maximum_percentage_error = 13
        number_test_pts = 500
        img1 = Gray{Float64}.(testimage("mandrill"))

        tol = 0.3
        Δ = (3.0, -1.0)
        trans = Translation(-Δ[1], -Δ[2])
        img2 = warp(img1, trans, axes(img1))

        corners = imcorner(img1, method=shi_tomasi)
        I = findall(!iszero, corners)
        r, c = (getindex.(I, 1), getindex.(I, 2))
        points = map((ri, ci) -> SVector{2}(Float64(ri), Float64(ci)), r, c)
        Random.seed!(9876)
        points = rand(points, (number_test_pts,))
        flow, status_array = optical_flow(img1, img2, points, LucasKanade(20, 11, 4, 0.000001))

        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100

        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error

        # Same as above, except that we pass in an initial displacement field.
        displacement = Array{SVector{2, Float64}, 1}(undef,length(points))
        for i in eachindex(displacement)
                displacement[i] = SVector{2, Float64}(0.0, 0.0)
        end
        flow, status_array = optical_flow(img1, img2, points, displacement, LucasKanade(20, 11, 4, 0.000001))

        error_count, maximum_error = evaluate_error(size(img1), flow[status_array], Δ, tol)
        percentage_error = (error_count / sum(status_array)) * 100

        @info "Case: Combined Motion", "Percentage Error: $(percentage_error)", "Maximum Error: $(maximum_error)"
        @test percentage_error  < maximum_percentage_error
    end

end;


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export LiblinearConfig, LiblinearConfigSpace

using LIBLINEAR

@with_kw struct LiblinearConfig
    C::Float64 = 1.0
    eps::Float64 = 0.1
end

struct LibLinearWrapper{LinearModel}
    map::Dict{UInt64,Int}
    cls::LinearModel
end

function create(config::LiblinearConfig, train_X, train_y, dim)
    map = Dict{UInt64,Int}()
    map[0] = 1 # the vectorization procedure use the zero id as an special symbol
    train_X_ = [Dict(get!(map, k, length(map)+1) => v for (k, v) in x) for x in train_X]
    # dim + 1 for out-of-vocabulary tokens
    cls = linear_train(train_y.refs, sparse(train_X_, dim + 1); C=config.C, eps=config.eps, bias=1.0)
    LibLinearWrapper(map, cls)
end

@with_kw struct LiblinearConfigSpace <: AbstractSolutionSpace
    C = [1.0]
    eps = [0.1]
    scale_C = (lower=0.001, s=3.0, upper=1000.0)
    scale_eps = (lower=0.0001, s=3.0, upper=0.99)
end

Base.eltype(::LiblinearConfigSpace) = LiblinearConfig

function Base.rand(space::LiblinearConfigSpace)
    LiblinearConfig(rand(space.C), rand(space.eps))
end

function combine(a::LiblinearConfig, b::LiblinearConfig)
    LiblinearConfig(a.C, b.eps)
end

function mutate(space::LiblinearConfigSpace, a::LiblinearConfig, iter)
    C = space.scale_C === nothing ? a.C : SearchModels.scale(a.C; space.scale_C...)
    eps = space.scale_eps === nothing ? a.eps : SearchModels.scale(a.eps; space.scale_eps...)
    LiblinearConfig(C, eps)
end

function predict(w::LibLinearWrapper, vec::SVEC)
    D = Dict(get(w.map, k, 1) => v for (k, v) in vec)
    ypred = linear_predict(w.cls, sparse([D], w.cls.nr_feature))
    ypred[1][1]
end


text:
using Flux, BSON, Images, DelimitedFiles, Base.Iterators

labels = readdlm("./char_dict.txt",String)[:,2]
BSON.@load "./trained_model.bson" predict

function recognition(imgs,batch_size=128)
    x = (cat(map(x->imresize(1 .- float(Gray.(load(x))),64,64),img)..., dims=4) for img in partition(imgs,batch_size))
    y = labels[vcat([Flux.onecold(predict(i)) for i in x])...]
    for i=1:length(imgs)
        println(imgs[i], '\t', y[i])
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    recognition(ARGS)
end

max_stars_repo_name: afternone/HCCR

max_issues_repo_name: afternone/HCCR

max_forks_repo_name: afternone/HCCR


text:
module ComfyJulia

    # Usage
    # living_repl(filename, time, redef_filter)
    # reads the file at "filename" and checks it for changes
    # if it hasnt changed, sleep for "time" and check again, forever
    # redef_filter filters out the "WARNING: Method definition..." error
    function living_repl(filename, time, redef_filter = true)
        if is_windows()
            # read the file as a string to compare with future reads for file modification
            basefile = String(read(filename))
        elseif is_unix()
            # Unix timestamp of when the file was last modified
            # (dont know if this works with windows)
            basefile = stat(filename).mtime
        else
            warn("what filesystem are you on?")
        end

        while true
            sleep(time)
            #file = String(read(filename))
            file = stat(filename).mtime

            if basefile != file
                err_rd, err_wt = redirect_stderr()
                try
                    include(filename)
                catch error
                    println(error)
                finally
                    basefile = file
                end
                redirect_stderr(STDOUT)

                outerr = String(readavailable(err_rd))
                outerr = split(outerr, "\n")

                for i in outerr
                    if contains(i, "WARNING: Method definition") & redef_filter
                        continue
                    end
                    println(i)
                end
            end
        end
    end

    export living_repl
end

max_stars_repo_name: isebarn/ComfyJulia

max_issues_repo_name: isebarn/ComfyJulia

max_forks_repo_name: isebarn/ComfyJulia


text:
# PORT
struct Port
    node::AbstractNode
    index::Int
    function Port(node::AbstractNode, index)
        ports = freeports(node)
        any(ports) || error("Node $node has no free ports")
        ports[index] || error("Port $index in node $node is already connected")
        new(node, index)
    end
    Port(node::Junction, index) = new(node, index)
end
Port(node::AbstractNode) = Port(node, nextfreeport(node))

# BOND
struct Bond <: g.AbstractSimpleEdge{Int}
    srcport::Port
    dstport::Port
end
function Bond(srcnode::AbstractNode, dstnode::AbstractNode)
    Bond(Port(srcnode), Port(dstnode))
end

# Source and Destination
srcnode(b::Bond) = b.srcport.node
dstnode(b::Bond) = b.dstport.node

# Base functions
in(n::AbstractNode, b::Bond) = n === srcnode(b) || n === dstnode(b)

iterate(b::Bond) = (b.srcport, true)
iterate(b::Bond, state) = state ? (b.dstport, false) : nothing

show(io::IO, port::Port) = print(io, "Port $(port.node) ($(port.index))")
show(io::IO, b::Bond) = print(io, "Bond $(srcnode(b)) ⇀ $(dstnode(b))")

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


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

using UpROOT
using Test

using ArraysOfArrays
using TypedTables
using StatsBase

@testset "testfiles" begin
    @testset "nesteddirs" begin
        file = TFile(UpROOT.testfiles["nesteddirs"])
        @test "one;1" in keys(file)
        @test length(file) == 2
        @test file["one"] isa TDirectory
        @test file["one"]["two"] isa TDirectory
    end

    @testset "leaflist" begin
        file = TFile(UpROOT.testfiles["leaflist"])
        @test "tree;1" in keys(file)
        tree = file["tree"]

        @test propertynames(tree) == [:leaflist]

        @test tree[1] isa NamedTuple
        @test tree[1].leaflist isa NamedTuple
        @test tree[1:5] isa Table
        @test tree[1:5].leaflist isa Table
        @test tree[:] isa Table
        @test tree[:].leaflist isa Table

        @test tree.leaflist[1] isa NamedTuple
        @test tree.leaflist[1:5] isa Table
        @test tree.leaflist[:] isa Table

        @test tree[1].leaflist == tree.leaflist[1]
    end

    @testset "HZZ" begin
        file = TFile(UpROOT.testfiles["HZZ"])
        @test "events;1" in keys(file)
        tree = file["events"]

        @test tree[1] isa NamedTuple
        @test tree[1:5] isa Table
        @test tree[:] isa Table

        @test tree.Jet_E[1] isa AbstractVector
        @test tree.Jet_E[1:5] isa VectorOfVectors
        @test tree.Jet_E[:] isa VectorOfVectors

        @test copy(tree.Jet_E) == tree.Jet_E[:]
        @test Array(tree.Jet_E) == tree.Jet_E[:]
        @test convert(Array, tree.Jet_E) == tree.Jet_E[:]

        @test Table(tree) == tree[:]

        @test tree[1:5, (:Jet_Px, :Jet_Py, :Jet_Pz)] isa Table{<:NamedTuple{(:Jet_Px, :Jet_Py, :Jet_Pz)}}
        @test tree[:, [:Jet_Px, :Jet_Py, :Jet_Pz]] isa Table{<:NamedTuple{(:Jet_Px, :Jet_Py, :Jet_Pz)}}
        @test tree[[2, 4, 7, 11], ["Jet_Px", "Jet_Py", "Jet_Pz"]] isa Table{<:NamedTuple{(:Jet_Px, :Jet_Py, :Jet_Pz)}}
    end

    @testset "histograms" begin
        file = TFile(UpROOT.testfiles["hepdata-example"])
        @test "hpx;1" in keys(file)
        @test "hpxpy;1" in keys(file)
        hpx = file["hpx;1"]
        hpxpy = file["hpxpy;1"]

        @test hpx isa Histogram{<:Real, 1}
        @test hpxpy isa Histogram{<:Real, 2}
    end
end


text:
similar_range(::Union{OneToSRange,UnitSRange,StepSRange,LinSRange,StepSRangeLen}, ::Union{OneToSRange,UnitSRange,StepSRange,LinSRange,StepSRangeLen}) = srange
similar_range(::Union{OneToSRange,UnitSRange,StepSRange,LinSRange,StepSRangeLen}, ::Union{OneTo,UnitRange,StepRange,LinRange,StepRangeLen}) = range
similar_range(::Union{OneTo,UnitRange,StepRange,LinRange,StepRangeLen}, ::Union{OneToSRange,UnitSRange,StepSRange,LinSRange,StepSRangeLen}) = range
similar_range(::Union{OneToMRange,UnitMRange,StepMRange,LinMRange,StepMRangeLen}, ::Union{OneTo,UnitRange,StepRange,LinRange,StepRangeLen}) = mrange
similar_range(::Union{OneTo,UnitRange,StepRange,LinRange,StepRangeLen}, ::Union{OneToMRange,UnitMRange,StepMRange,LinMRange,StepMRangeLen}) = mrange
similar_range(::Union{OneToSRange,UnitSRange,StepSRange,LinSRange,StepSRangeLen}, ::Union{OneToMRange,UnitMRange,StepMRange,LinMRange,StepMRangeLen}) = mrange
similar_range(::Union{OneToMRange,UnitMRange,StepMRange,LinMRange,StepMRangeLen}, ::Union{OneToSRange,UnitSRange,StepSRange,LinSRange,StepSRangeLen}) = mrange

for RANGE_TYPE in (:OneTo,:OneToMRange,:OneToSRange)
    @eval begin
        function StaticArrays.similar_type(
            ::$(RANGE_TYPE){T},
            element_type=T,
           ) where {T}
            return $(RANGE_TYPE){element_type}
        end
    end
end

lower_rangetype(::Type{<:OneToSRange{T,E}}) where {T,E} = OneTo{T}
lower_rangetype(::Type{OneTo{T}}) where {T} = OneToMRange{T}
lower_rangetype(::Type{OneToMRange{T}}) where {T} = OneToMRange{T}

for RANGE_TYPE in (:UnitRange,:UnitMRange,:UnitSRange)
    @eval begin
        function StaticArrays.similar_type(
            ::$(RANGE_TYPE){T},
            element_type=T,
           ) where {T}
            return $(RANGE_TYPE){element_type}
        end
    end
end

lower_rangetype(::Type{<:UnitSRange{T}}) where {T} = UnitRange{T}
lower_rangetype(::Type{UnitRange{T}}) where {T} = UnitMRange{T}
lower_rangetype(::Type{UnitMRange{T}}) where {T} = UnitMRange{T}

for RANGE_TYPE in (:LinRange,:LinMRange,:LinSRange)
    @eval begin
        function StaticArrays.similar_type(
            ::$(RANGE_TYPE){T},
            element_type=T,
           ) where {T}
            return $(RANGE_TYPE){element_type}
        end
    end
end

lower_rangetype(::Type{<:LinSRange{T}}) where {T} = LinRange{T}
lower_rangetype(::Type{LinRange{T}}) where {T} = LinMRange{T}
lower_rangetype(::Type{LinMRange{T}}) where {T} = LinMRange{T}

for RANGE_TYPE in (:StepRangeLen,:StepMRangeLen,:StepSRangeLen)
    @eval begin
        function StaticArrays.similar_type(
            ::$(RANGE_TYPE){T,R,S},
            element_type=T,
            reference_type=R,
            step_type=S
           ) where {T,R,S}
            return $(RANGE_TYPE){element_type,reference_type,step_type}
        end
    end
end

lower_rangetype(::Type{<:StepSRangeLen{T,R,S}}) where {T,R,S} = StepRangeLen{T,R,S}
lower_rangetype(::Type{StepRangeLen{T,R,S}}) where {T,R,S} = StepMRangeLen{T,R,S}
lower_rangetype(::Type{StepMRangeLen{T,R,S}}) where {T,R,S} = StepMRangeLen{T,R,S}

for RANGE_TYPE in (:StepRange,:StepMRange,:StepSRange)
    @eval begin
        function StaticArrays.similar_type(
            ::$(RANGE_TYPE){T,S},
            element_type=T,
            step_type=S
           ) where {T,R,S}
            return $(RANGE_TYPE){element_type,step_type}
        end
    end
end

lower_rangetype(::Type{<:StepSRange{T,S}}) where {T,S} = StepRange{T,S}
lower_rangetype(::Type{StepRange{T,S}}) where {T,S} = StepMRange{T,S}
lower_rangetype(::Type{StepMRange{T,S}}) where {T,S} = StepMRange{T,S}

for S in (:OneToSRange,:UnitSRange,:StepSRange,:LinSRange,:StepSRangeLen)
    for M in (:OneToMRange,:UnitMRange,:StepMRange,:LinMRange,:StepMRangeLen)
        @eval begin
            Base.promote_rule(a::Type{<:$S}, b::Type{<:$M}) = promote_rule(lower_rangetype(a), b)
            Base.promote_rule(b::Type{<:$M}, a::Type{<:$S}) = promote_rule(lower_rangetype(a), b)
        end
    end
end

for S in (:OneToSRange,:UnitSRange,:StepSRange,:LinSRange,:StepSRangeLen)
    for M in (:OneTo,:UnitRange,:StepRange,:LinRange,:StepRangeLen)
        @eval begin
            Base.promote_rule(a::Type{<:$S}, b::Type{<:$M}) = promote_rule(lower_rangetype(a), b)
            Base.promote_rule(b::Type{<:$M}, a::Type{<:$S}) = promote_rule(lower_rangetype(a), b)
        end
    end
end

for S in (:OneTo,:UnitRange,:StepRange,:LinRange,:StepRangeLen)
    for M in (:OneToMRange,:UnitMRange,:StepMRange,:LinMRange,:StepMRangeLen)
        @eval begin
            Base.promote_rule(a::Type{<:$S}, b::Type{<:$M}) = promote_rule(lower_rangetype(a), b)
            Base.promote_rule(b::Type{<:$M}, a::Type{<:$S}) = promote_rule(lower_rangetype(a), b)
        end
    end
end

# Necessary to avoid ambiguities
Base.promote_rule(a::Type{StepRangeLen{T,R,S}}, b::Type{<:StepMRangeLen}) where {T,R,S} = promote_rule(lower_rangetype(a), b)

function Base.promote_rule(a::Type{StepRangeLen{T,R,S}}, b::Type{StepMRangeLen{T,R,S}}) where {T,R,S}
    return StepMRangeLen{T,R,S}
end

###
### AbstractStepRangeLen
###
function Base.promote_rule(
    ::Type{StepMRangeLen{T1,R1,S1}},
    ::Type{StepMRangeLen{T2,R2,S2}}
   ) where {T1,T2,R1,R2,S1,S2}
    return el_same(
        promote_type(T1,T2),
        StepMRangeLen{T1,promote_type(R1,R2),promote_type(S1,S2)},
        StepMRangeLen{T2,promote_type(R1,R2),promote_type(S1,S2)}
       )
end

function Base.promote_rule(
    ::Type{<:StepSRangeLen{T1,Tr1,Ts1}},
    ::Type{<:StepSRangeLen{T2,Tr2,Ts2}},
   ) where {T1,T2,Tr1,Tr2,Ts1,Ts2,R1,R2,S1,S2,L1,L2,F1,F2}
    return el_same(
        promote_type(T1,T2),
        StepSRangeLen{T1,promote_type(Tr1,Tr2),promote_type(Ts1,Ts2)},
        StepSRangeLen{T2,promote_type(Tr1,Tr2),promote_type(Ts1,Ts2)}
       )
end

# has to be included to avoid ambiguities
Base.promote_rule(a::Type{StepRangeLen{T,R,S}}, b::Type{A}) where {A<:StepSRangeLen,T,R,S} = promote_rule(a, lower_rangetype(b))

Base.promote_rule(a::Type{StepMRangeLen{T1,R,S}}, ::Type{LinMRange{T2}}) where {T1,R,S,T2} = promote_rule(a, StepMRangeLen{T2,T2,T2})

Base.promote_rule(a::Type{<:StepSRangeLen{T1,R,S}}, ::Type{<:LinSRange{T2}}) where {T1,R,S,T2} = promote_rule(StepSRangeLen{T2,T2,T2}, a)

###
### AbstractLinRange
###
Base.promote_rule(a::Type{LinMRange{T1}}, b::Type{LinMRange{T2}}) where {T1,T2} = LinMRange{promote_type(T1,T2)}
Base.promote_rule(a::Type{<:LinSRange{T1}}, b::Type{<:LinSRange{T2}}) where {T1,T2} = LinSRange{promote_type(T1,T2)}

###
### AbstractStepRange
###
function Base.promote_rule(
    ::Type{StepMRange{T1a,T1b}},
    ::Type{StepMRange{T2a,T2b}}
   ) where {T1a,T1b,T2a,T2b}
    return Base.el_same(
        promote_type(T1a,T2a),
        # el_same only operates on array element type, so just promote
        # second type parameter
        StepMRange{T1a, promote_type(T1b,T2b)},
        StepMRange{T2a, promote_type(T1b,T2b)}
       )
end

function Base.promote_rule(
    ::Type{<:StepSRange{T1a,T1b}},
    ::Type{<:StepSRange{T2a,T2b}}
   ) where {T1a,T1b,T2a,T2b}
    return Base.el_same(
        promote_type(T1a,T2a),
        # el_same only operates on array element type, so just promote
        # second type parameter
        StepSRange{T1a, promote_type(T1b,T2b)},
        StepSRange{T2a, promote_type(T1b,T2b)}
       )
end

Base.promote_rule(a::Type{LinMRange{T1}}, ::Type{StepMRange{T2,Ts2}}) where {T1,T2,Ts2} = promote_rule(a, LinMRange{T2})

Base.promote_rule(a::Type{StepMRangeLen{T1,R,S1}}, ::Type{StepMRange{T2,S2}}) where {T1,T2,R,S1,S2} = promote_rule(a, StepMRangeLen{T2,T2,T2})
Base.promote_rule(a::Type{<:StepSRangeLen{T1,R,S1}}, ::Type{<:StepSRange{T2,S2}}) where {T1,T2,R,S1,S2} = promote_rule(a, StepMRangeLen{T2,T2,T2})
Base.promote_rule(a::Type{<:LinSRange{T1}}, ::Type{<:StepSRange{T2,Ts2}}) where {T1,T2,Ts2} = promote_rule(a, LinSRange{T2})

###
### UnitRange
###

Base.promote_rule(a::Type{UnitMRange{T1}}, b::Type{UnitMRange{T2}}) where {T1,T2} = el_same(promote_type(T1,T2), a, b)
Base.promote_rule(a::Type{<:UnitSRange{T1}}, b::Type{<:UnitSRange{T2}}) where {T1,T2} = el_same(promote_type(T1,T2), a, b)

Base.promote_rule(::Type{LinMRange{T2}}, ::Type{UnitMRange{T1}}) where {T1,T2} = LinMRange{promote_type(T1,T2)}
Base.promote_rule(::Type{UnitMRange{T1}}, ::Type{LinMRange{T2}}) where {T1,T2} = LinMRange{promote_type(T1,T2)}

Base.promote_rule(::Type{StepMRange{T1,S1}}, ::Type{UnitMRange{T2}}) where {T1,S1,T2} = promote_rule(StepMRange{T1,S1}, StepMRange{T2,T2})
Base.promote_rule(::Type{<:StepSRange{T1,S1}}, ::Type{<:UnitSRange{T2}}) where {T1,S1,T2} = promote_rule(StepSRange{T1,S1}, StepSRange{T2,T2})

Base.promote_rule(::Type{<:LinSRange{T2}}, ::Type{<:UnitSRange{T1}}) where {T1,T2} = LinSRange{promote_type(T1,T2)}
Base.promote_rule(::Type{<:UnitSRange{T1}}, ::Type{<:LinSRange{T2}}) where {T1,T2} = LinSRange{promote_type(T1,T2)}

###
### OneToRange
###

Base.promote_rule(a::Type{LinRange{T}}, ::Type{OR}) where {T,OR<:OneToMRange} = promote_rule(LinMRange{T},LinMRange{eltype(OR)})
Base.promote_rule(a::Type{<:OneToSRange{T1,Any}}, b::Type{<:OneToSRange{T2,Any}}) where {T1,T2} = OneToSRange{promote_type(T1,T2)}
Base.promote_rule(a::Type{OneToMRange{T1}}, b::Type{OneToMRange{T2}}) where {T1,T2} = OneToMRange{promote_type(T1,T2)}

Base.promote_rule(a::Type{<:UnitSRange{T1}}, b::Type{<:OneToSRange{T2}}) where {T1,T2} = UnitSRange{promote_type(T1,T2)}
Base.promote_rule(a::Type{UnitMRange{T1}}, b::Type{OneToMRange{T2}}) where {T1,T2} = UnitMRange{promote_type(T1,T2)}

Base.promote_rule(::Type{OneToMRange{T1}}, ::Type{OneTo{T2}}) where {T1,T2} = promote_rule(OneToMRange{T1},OneToMRange{T2})
Base.promote_rule(::Type{OneTo{T2}}, ::Type{OneToMRange{T1}}) where {T1,T2} = promote_rule(OneToMRange{T1},OneToMRange{T2})

Base.promote_rule(::Type{OneToMRange{T1}}, ::Type{LinMRange{T2}}) where {T1,T2} = promote_rule(LinMRange{T1},LinMRange{T2})
Base.promote_rule(::Type{LinMRange{T1}}, ::Type{OneToMRange{T2}}) where {T1,T2} = promote_rule(LinMRange{T1},LinMRange{T2})

Base.promote_rule(::Type{<:OneToSRange{T1}}, ::Type{<:LinSRange{T2}}) where {T1,T2} = promote_rule(LinSRange{T1},LinSRange{T2})
Base.promote_rule(::Type{<:LinSRange{T1}}, ::Type{<:OneToSRange{T2}}) where {T1,T2} = promote_rule(LinSRange{T1},LinSRange{T2})


text:
#File: fStd.jl
#Author:
#Date: 29-June-2020

#STANDARD DEVIATION
function fStd(x,flag,f)
#--------------------------------------------------------------------------
#x:     returns
#flag:  0 = sample, 1 = population
#f:     reporting frequency
#m2:    standard deviation (sample or population)
#--------------------------------------------------------------------------

    m2 = std(x,flag).*sqrt(f)
    return m2
end


text:
#=
The following file contains various function to help smoke test and benchmark
arbitrary AbstractRBMs and configurations.
=#

using Base.Test

"""
The TestReporter uses the ratio test or D'Alembert's criterion
to monitor convergence. This is helpful in automating tests that
confirm that a given RBM is learning something.

NOTE: current we just store the calculated ratios from each epoch,
but we could probably just do an online mean calculation.
"""
type TestReporter
    ratios::Vec{Float64}
    prev::Float64
    log::Bool

    function TestReporter(log=false)
        new(Array{Float64}(0), NaN, log)
    end
end

"""
Computes and adds a new ratio to the reporter. If log is set to
true in the reporter a message containing the score, time and mean ratio
will be printed.
"""
function report(r::TestReporter, rbm::AbstractRBM,
                epoch::Int, epoch_time::Float64, score::Float64)
    if epoch == 1
        r.ratios = Array{Float64}(0)
        r.prev = NaN
    end

    if !isnan(r.prev)
        # We're ratio test, also known as D'Alembert's criterion
        ratio = r.prev == 0.0 ? abs(score / nextfloat(r.prev)) : abs(score / r.prev)

        push!(r.ratios, ratio)
    end
    r.prev = score

    if r.log
        println("[Epoch $epoch] Score: $score [$(epoch_time)s] (Mean Ratio: $(mean(r.ratios)) )")
    end
end

"""
Simply return true or false as to whether the mean ratio is less than 1
"""
converged(r::TestReporter) = mean(r.ratios) < 1.0

"""
Simply returns a DefaultContext Dict.
"""
function DefaultContext()
    return Dict(
        :weight_decay_kind => :l2,
        :weight_decay_rate => 0.001,
        :sparsity_cost => 0.001,
        :sparsity_target => 0.02,
        :lr => 0.1,
        :momentum => 0.9,
        :batch_size => 100,
        :n_epochs => 5,
        :n_gibbs => 1,
        :reporter => TestReporter()
    )
end

"""
Generates synthetic random datasets with several modifiable properties.

Args:
* T::Type - the type (or precision) of the resulting matrix
* n_features::Int - the number of features in each observation

Optional:
* n_classes::Int - the number of unique classes or categories in the resulting dataset. (default=10)
* n_obs::Int - total number of observations in resulting dataset. (default=1000)
* sparsity::Float64 - specifies the density if a sparse matrix is desired. (default=-1.0)
                      if less than 0.0 a dense matrix is created.
* binary::Bool - whether or not to round the result to 0.0 or 1.0

Returns:
* Mat{T}(n_features, n_obs) - with the various properties specified.
"""
function generate_dataset(T::Type, n_features; n_classes=10, n_obs=1000, sparsity=-1.0, binary=true)
    function rand_fill!(X::Mat, prototypes::Mat)
        for i in 1:size(X, 2)
            j = rand(1:size(prototypes, 2))
            X[:, i] = prototypes[:, j]
        end

        if binary
            return round(X)
        else
            return X
        end
    end

    if sparsity > 0
        prototypes = map(T, sprand(n_features, n_classes, sparsity))
        X = map(T, spzeros(n_features, n_obs))
        return rand_fill!(X, prototypes)
    else
        prototypes = rand(T, n_features, n_classes)
        X = zeros(n_features, n_obs)
        return rand_fill!(X, prototypes)
    end
end

"""
Runs simple smoke tests on the provided RBM.

Args:
* rbm::AbstractRBM - the rbm to test

Optional:
* opts::Dict - an alternate context to use when calling `fit`. (default=DEFAULT_CONTEXT)
* n_obs::Int - the number of observations to generate for the synthetic datasets. (default=1000)
* debug::Bool - whether or not to print each epoch. (default=false)
                Only applies if the reporter in `ctx` is `TestReporter`

NOTE: Only using dense arrays for the dataset cause the conditional rbm doesn't support
sparse ones yet.
"""
function test{T,V,H}(rbm::AbstractRBM{T,V,H}; opts::Dict=DefaultContext(), input_size=-1, n_obs=1000, debug=false)
    n_hid, n_vis = size(rbm.W)

    if input_size > 0
        n_vis = input_size
    end

    ctx = deepcopy(opts)
    DX = generate_dataset(T, n_vis; n_obs=n_obs)
    SX = generate_dataset(T, n_vis; n_obs=n_obs, sparsity=0.3)

    if debug && isa(ctx[:reporter], TestReporter)
        ctx[:reporter].log = true
    end

    for X in (DX, SX)
        info("Testing against $(typeof(X).name)")
        info("Testing fit")
        fit(rbm, X, ctx)

        info("Testing for convergence")
        @test converged(ctx[:reporter])

        # Check that none of the array fields contain Infs or NaNs
        for name in fieldnames(rbm)
            if fieldtype(typeof(rbm), name) <: AbstractArray
                field = getfield(rbm, name)

                @test !any(isnan, field) && !any(isinf, field)
            end
        end

        info("Testing transform")
        Y = transform(rbm, X)
        @test !any(isnan, Y) && !any(isinf, Y)

        info("Testing generate")
        Y = generate(rbm, X)
        @test !any(isnan, Y) && !any(isinf, Y)
    end
end

"""
Runs a benchmark comparison between various rbms on a given function.

Args:
* func::Function - the function to compare
* rbms::Array{R<:AbstractRBM,1} - the RBMs to test
* args::Array{T<:Tuple,1} - the arguments to pass to each RBM's function

Returns: a DataFrame of the results
"""
function compare{R<:AbstractRBM, T<:Tuple}(func::Function, rbms::Array{R,1}, args::Array{T,1})
    @eval import Benchmark
    @assert length(rbms) == length(args)

    funcs = Array{Function}(length(rbms))
    names = Array{TypeName}(length(rbms))

    for i in eachindex(rbms)
        funcs[i] = () -> func(rbms[i], args[i]...)
        names[i] = typeof(rbms[i]).name
    end

    df = Benchmark.compare(funcs, 1)
    df[:Function] = fill(last(split(string(func), ".")), length(funcs))
    df[:RBM] = names

    return df
end

"""
Runs a full benchmark comparison on the public interface of an RBM
against those in Boltzmann.jl

Args:
* rbm::AbstractRBM - the rbm to test

Optional:
* n_obs::Int - number of observation in the synthetic dataset. (default=1000)

Returns: a DataFrame of the results

NOTE: For now we're only using dense matrices cause the conditional rbm doesn't
work with sparse ones.
"""
function compare{T,V,H}(rbm::AbstractRBM{T,V,H}; opts::Dict=DefaultContext(), input_size=-1, debug=false, n_obs=1000)
    @eval using DataFrames

    ctx = deepcopy(opts)
    models = AbstractRBM[rbm]
    n_hid, n_vis = size(rbm.W)

    if input_size > 0
        n_vis = input_size
    end

    DX = generate_dataset(T, n_vis; n_obs=n_obs)
    SX = generate_dataset(T, n_vis; n_obs=n_obs, sparsity=0.3)

    if !isa(rbm, RBM)
        push!(models, RBM(T, V, H, n_vis, n_hid))
    end

    if !isa(rbm, ConditionalRBM)
        n_cond = round(Int, n_vis / 2)
        push!(models, ConditionalRBM(T, V, H, n_vis - n_cond, n_hid, n_cond))
    end

    n = length(models)
    results = DataFrame()

    if debug && isa(ctx[:reporter], TestReporter)
        ctx[:reporter].log = true
    end

    for X in (DX, SX)
        fit_result = compare(fit, models, fill((X, deepcopy(ctx)), n))
        trans_result = compare(transform, models, fill((X,), n))

        tmp_result = vcat(
            fit_result,
            trans_result
        )

        tmp_result[:Input] = fill(typeof(X).name, n * 2)
        tmp_result[:Type] = fill(T, n * 2)
        results = vcat(results, tmp_result)
    end

    return results
end

"""
Runs a benchmark comparison on an RBM's fit method with different contexts.
This can be useful if you want benchmark various update or gradient methods
specified in the context.

Args:
* rbm::AbstractRBM - the rbm to test
* options::Dict... - the contexts to compare

Optional:
* n_obs::Int - the number of observations to generate in the synthetic dataset.

Returns: a DataFrame of the results
"""
function compare{T,V,H}(rbm::AbstractRBM{T,V,H}, options::Dict...; input_size=-1, n_obs=1000)
    n_hid, n_vis = size(rbm.W)
    if input_size > 0
        n_vis = input_size
    end

    X = generate_dataset(T, n_vis; n_obs=n_obs)
    args = [map(opt -> (X, opt), options)...]
    rbms = map(i -> deepcopy(rbm), 1:length(args))
    return compare(fit, rbms, args)
end

"""
Runs a series of basic smoke and benchmark tests on the interface of any RBM.

Default methods tested:

* `fit(model, X)`
* `transform(model, X)`
* `generate(model, X)`

Args:
* model::Union{AbstractRBM, Net} - the model to test

Optional:
* opts::Dict - the context passed when running fit. (default=DEFAULT_CONTEXT)
* n_obs::Int - total number of observations in generated dataset (default=1000)

Returns: a DataFrame of the results
"""
function benchmark{T,V,H}(rbm::AbstractRBM{T,V,H}; opts::Dict=DefaultContext(), input_size=-1, debug=false, n_obs=1000)
    @eval begin
        using DataFrames
        import Benchmark
    end

    ctx = deepcopy(opts)
    n_hid, n_vis = size(rbm.W)

    if input_size > 0
        n_vis = input_size
    end

    DX = generate_dataset(T, n_vis; n_obs=n_obs)
    SX = generate_dataset(T, n_vis; n_obs=n_obs, sparsity=0.3)

    if debug && isa(ctx[:reporter], TestReporter)
        ctx[:reporter].log = true
    end

    results = DataFrame()

    for X in (DX, SX)
        # For now we're just benchmarking the fit method
        # but we could also get transform, generate, etc
        fit_func() = fit(rbm, X, ctx)
        trans_func() = transform(rbm, X)
        gen_func() = generate(rbm, X)

        results = vcat(results, Benchmark.benchmark(fit_func, "fit", "$(typeof(rbm).name):$(typeof(X).name)", 10))
        results = vcat(results, Benchmark.benchmark(trans_func, "transform", "$(typeof(rbm).name):$(typeof(X).name)", 10))
        results = vcat(results, Benchmark.benchmark(gen_func, "generate", "$(typeof(rbm).name):$(typeof(X).name)", 10))
    end

    return results
end


text:
module MarkovBounds

using Reexport
using SumOfSquares

# SumOfSquares and DynamicPolynomials are the main modeling tools to set up
# moment bounding problems --> Reexport those
@reexport using SumOfSquares: @set, SemialgebraicSets, FullSpace

@reexport using DynamicPolynomials: @polyvar, PolyVar, AbstractPolynomialLike,
                                    Polynomial, Term, subs, polynomial, differentiate,
                                    MonomialVector, maxdegree

# MTK/Symbolics/SymbolicUtils/Catalyst can also be used to specify moment problem
# data but not as its main approach
using Catalyst: @reaction_network, @parameters, species, speciesmap,
                reactions, paramsmap, prodstoichmat, substoichmat,
                ReactionSystem, Reaction

using Symbolics: Num, expand
using SymbolicUtils: Add, Mul, Pow, Sym, arguments
using LightGraphs: edges, vertices, SimpleGraph, add_edge!, add_vertex!, Edge
using MetaGraphs: MetaGraph, MetaDiGraph, props, set_prop!

import LinearAlgebra: qr, nullspace, diag, Diagonal, tr, transpose
import Base: show
import Parameters: @unpack
import SumOfSquares.MultivariatePolynomials.polynomial
import SumOfSquares.SemialgebraicSets.inequalities

const PV = PolyVar
const APL = AbstractPolynomialLike

include("distributed.jl")
include("processes.jl")
include("symbolics.jl")
include("print_wrapper.jl")
include("utils.jl")
include("problem_setup.jl")
include("sos_programs_stationary.jl")
include("sos_programs_transient.jl")
include("sos_programs_control.jl")

end


text:
choleskyvectens(a,b) = permutedims(sum(a.*b,1), [3 2 1])

macro setupQTIC()
    quote
        m     = size(u,1)
        n,_,N = size(fx)

        @assert size(cx) == (n, N) "size(cx) should be (n, N)"
        @assert size(cu) == (m, N) "size(cu) should be (m, N)"
        @assert size(cxx) == (n, n) "size(cxx) should be(n, n) "
        @assert size(cxu) == (n, m) "size(cxu) should be(n, m) "

        k          = zeros(m,N)
        K          = zeros(m,n,N)
        Vx         = zeros(n,N)
        Vxx        = zeros(n,n,N)
        Quu        = Array{T}(undef,m,m,N)
        Quui       = Array{T}(undef,m,m,N)
        dV         = [0., 0.]

        @views Vx[:,N] .= cx[:,N]
        Vxx[:,:,N] = cxx
        Quu[:,:,N] = cuu
        diverge    = 0
    end |> esc
end

macro end_backward_pass()
    @views quote
        QuF = Qu
        if isempty(lims) || lims[1,1] > lims[1,2]
            # debug("#  no control limits: Cholesky decomposition, check for non-PD")
            local R
            try
                R = cholesky(Hermitian(QuuF))
            catch
                diverge  = i
                return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx, dV
            end
            # debug("#  find control law")
            k_i = -(R\QuF)
            K_i = -(R\Qux_reg)
        else
            # debug("#  solve Quadratic Program")
            lower = lims[:,1]-u[:,i]
            upper = lims[:,2]-u[:,i]
            local k_i,result,free
            try
                k_i,result,R,free = boxQP(QuuF,QuF,lower,upper,k[:,min(i+1,N-1)])
            catch
                result = 0
            end
            if result < 1
                diverge  = i
                return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx, dV
            end
            K_i  = zeros(m,n)
            if any(free)
                Lfree         = -R\(R'\Qux_reg[free,:])
                K_i[free,:]   = Lfree
            end
        end
        # debug("#  update cost-to-go approximation")
        Quuki       = Quu[:,:,i]*k_i
        kiQuuki     = dot(k_i, Quuki)
        KiQuuki     = (K_i'Quuki)
        KiQuuKi     = K_i'Quu[:,:,i]*K_i
        dV         .= dV .+ [dot(k_i,Qu); .5*kiQuuki]
        Vx[:,i:i]  .= Qx .+ KiQuuki .+ K_i'Qu .+ Qux'k_i
        Vxx[:,:,i] .= Qxx .+ KiQuuKi .+ K_i'Qux .+ Qux'K_i
        Vxx[:,:,i] .+= Vxx[:,:,i]'
        Vxx[:,:,i] ./= 2

        # debug("# save controls/gains")
        k[:,i]   .= vec(k_i)
        K[:,:,i] .= K_i

    end |> esc
end

function back_pass(cx,cu,cxx::AbstractArray{T,3},cxu,cuu,fx::AbstractArray{T,3},fu,fxx,fxu,fuu,λ,regType,lims,x,u) where T # nonlinear time variant

    m,N        = size(u)
    n          = size(cx,1)
    @assert size(cx) == (n, N)
    @assert size(cu) == (m, N)
    @assert size(cxx) == (n, n, N)
    @assert size(cxu) == (n, m, N)
    @assert size(cuu) == (m, m, N)
    k          = zeros(m,N)
    K          = zeros(m,n,N)
    Vx         = zeros(n,N)
    Vxx        = zeros(n,n,N)
    Quu        = Array{T}(undef,m,m,N)
    Quui       = Array{T}(undef,m,m,N)
    dV         = [0., 0.]
    Vx[:,N]    = cx[:,N]
    Vxx[:,:,N] = cxx[:,:,N]
    Quu[:,:,N] = cuu[:,:,N]

    diverge  = 0
    for i = N-1:-1:1
        Qu  = cu[:,i]      + fu[:,:,i]'Vx[:,i+1]
        Qx  = cx[:,i]      + fx[:,:,i]'Vx[:,i+1]
        Qux = cxu[:,:,i]'  + fu[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        if !isempty(fxu)
            fxuVx = vectens(Vx[:,i+1],fxu[:,:,:,i])
            Qux   = Qux + fxuVx
        end

        Quu[:,:,i] = cuu[:,:,i]   + fu[:,:,i]'Vxx[:,:,i+1]*fu[:,:,i]
        if !isempty(fuu)
            fuuVx = vectens(Vx[:,i+1],fuu[:,:,:,i])
            Quu[:,:,i] .+= fuuVx
        end

        Qxx = cxx[:,:,i] + fx[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        isempty(fxx) || (Qxx .+= vectens(Vx[:,i+1],fxx[:,:,:,i]))
        Vxx_reg = Vxx[:,:,i+1] .+ (regType == 2 ? λ*eye(n) : 0)
        Qux_reg = cxu[:,:,i]'  + fu[:,:,i]'Vxx_reg*fx[:,:,i]
        isempty(fxu) || (Qux_reg .+= fxuVx)
        QuuF = cuu[:,:,i]  + fu[:,:,i]'Vxx_reg*fu[:,:,i] .+ (regType == 1 ? λ*eye(m) : 0)
        isempty(fuu) || (QuuF .+= fuuVx)

        @end_backward_pass
    end

    return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx,dV
end

function back_pass(cx,cu,cxx::AbstractArray{T,2},cxu,cuu,fx::AbstractArray{T,3},fu,fxx,fxu,fuu,λ,regType,lims,x,u) where T # quadratic timeinvariant cost, dynamics nonlinear time variant

    @setupQTIC

    for i = N-1:-1:1
        Qu  = cu[:,i]   + fu[:,:,i]'Vx[:,i+1]
        Qx  = cx[:,i]   + fx[:,:,i]'Vx[:,i+1]
        Qux = cxu' + fu[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        if !isempty(fxu)
            fxuVx = vectens(Vx[:,i+1],fxu[:,:,:,i])
            Qux   = Qux + fxuVx
        end

        Quu[:,:,i] = cuu + fu[:,:,i]'Vxx[:,:,i+1]*fu[:,:,i]
        if !isempty(fuu)
            fuuVx = vectens(Vx[:,i+1],fuu[:,:,:,i])
            Quu[:,:,i]   = Quu[:,:,i] + fuuVx
        end
        Qxx = cxx  + fx[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        isempty(fxx) || (Qxx .+= vectens(Vx[:,i+1],fxx[:,:,:,i]))
        Vxx_reg = Vxx[:,:,i+1] .+ (regType == 2 ? λ*eye(n) : 0)
        Qux_reg = cxu'   + fu[:,:,i]'Vxx_reg*fx[:,:,i]
        isempty(fxu) || (Qux_reg .+= fxuVx)
        QuuF = cuu  + fu[:,:,i]'Vxx_reg*fu[:,:,i] .+ (regType == 1 ? λ*eye(m) : 0)
        isempty(fuu) || (QuuF .+= fuuVx)
        @end_backward_pass
    end
    return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx,dV
end

function back_pass(cx,cu,cxx::AbstractArray{T,2},cxu,cuu,fx::AbstractArray{T,3},fu,λ,regType,lims,x,u) where T # quadratic timeinvariant cost, linear time variant dynamics
    @setupQTIC
    for i = N-1:-1:1
        Qu         = cu[:,i] + fu[:,:,i]'Vx[:,i+1]
        Qx         = cx[:,i] + fx[:,:,i]'Vx[:,i+1]
        Qux        = cxu' + fu[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        Quu[:,:,i] = cuu + fu[:,:,i]'Vxx[:,:,i+1]*fu[:,:,i]
        Qxx        = cxx + fx[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        Vxx_reg    = Vxx[:,:,i+1] .+ (regType == 2 ? λ*eye(n) : 0)
        Qux_reg    = cxu' + fu[:,:,i]'Vxx_reg*fx[:,:,i]
        QuuF       = cuu + fu[:,:,i]'Vxx_reg*fu[:,:,i] .+ (regType == 1 ? λ*eye(m) : 0)
        @end_backward_pass
    end

    return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx,dV
end

function back_pass(cx,cu,cxx::AbstractArray{T,3},cxu,cuu,fx::AbstractArray{T,3},fu,λ,regType,lims,x,u) where T # quadratic timeVariant cost, linear time variant dynamics
    m          = size(u,1)
    n,N        = size(fx, 1), size(fx, 3)

    @assert size(cx) == (n, N)
    @assert size(cu) == (m, N)
    @assert size(cxx) == (n, n, N)
    @assert size(cxu) == (n, m, N)
    @assert size(cuu) == (m, m, N)

    k          = zeros(m,N)
    K          = zeros(m,n,N)
    Vx         = zeros(n,N)
    Vxx        = zeros(n,n,N)
    Quu        = Array{T}(undef,m,m,N)
    Quui       = Array{T}(undef,m,m,N)
    dV         = [0., 0.]

    Vx[:,N]    = cx[:,N]
    Vxx[:,:,N] = cxx[:,:,end]
    Quu[:,:,N] = cuu[:,:,N]

    diverge    = 0
    for i = N-1:-1:1
        Qu          = cu[:,i] + fu[:,:,i]'Vx[:,i+1]
        Qx          = cx[:,i] + fx[:,:,i]'Vx[:,i+1]
        Vxx_reg     = Vxx[:,:,i+1] .+ (regType == 2 ? λ*eye(n) : 0)
        Qux_reg     = cxu[:,:,i]' + fu[:,:,i]'Vxx_reg*fx[:,:,i]
        QuuF        = cuu[:,:,i]  + fu[:,:,i]'Vxx_reg*fu[:,:,i] .+ (regType == 1 ? λ*eye(m) : 0)
        Qux         = cxu[:,:,i]' + fu[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        Quu[:,:,i] .= cuu[:,:,i] .+ fu[:,:,i]'Vxx[:,:,i+1]*fu[:,:,i]
        Qxx         = cxx[:,:,i]  + fx[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        @end_backward_pass
    end

    return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx,dV
end

@views function back_pass(cx,cu,cxx::AbstractArray{T,2},cxu,cuu,fx::AbstractMatrix{T},fu,λ,regType,lims,x,u) where T # cost quadratic and cost and LTI dynamics

    m,N = size(u)
    n   = size(fx,1)
    @assert size(cx) == (n, N)
    @assert size(cu) == (m, N)
    @assert size(cxx) == (n, n)
    @assert size(cxu) == (n, m)
    @assert size(cuu) == (m, m)
    k   = zeros(m,N)
    K   = zeros(m,n,N)
    Vx  = zeros(n,N)
    Vxx = zeros(n,n,N)
    Quu = Array{T}(undef,m,m,N)
    Quui = Array{T}(undef,m,m,N)
    dV  = [0., 0.]

    Vx[:,N]    = cx[:,end]
    Vxx[:,:,N] = cxx
    Quu[:,:,N] = cuu

    diverge    = 0
    for i = N-1:-1:1
        Qu         = cu[:,i] + fu'Vx[:,i+1]
        Qx         = cx[:,i] + fx'Vx[:,i+1]
        Qux        = cxu' + fu'Vxx[:,:,i+1]*fx
        Quu[:,:,i] = cuu + fu'Vxx[:,:,i+1]*fu
        Qxx        = cxx + fx'Vxx[:,:,i+1]*fx
        Vxx_reg    = Vxx[:,:,i+1] .+ (regType == 2 ? λ*eye(n) : 0)
        Qux_reg    = cxu' + fu'Vxx_reg*fx
        QuuF       = cuu + fu'Vxx_reg*fu .+ (regType == 1 ? λ*eye(m) : 0)
        @end_backward_pass
    end

    return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx,dV
end

function graphics(x...)
    return 0
end

function back_pass_gps(cx,cu,cxx::AbstractArray{T,3},cxu,cuu, fx::AbstractArray{T,3},fu,lims,x,u,kl_cost_terms) where T # quadratic timeVariant cost, linear time variant dynamics
    m          = size(u,1)
    n,_,N      = size(fx)
    ηbracket = kl_cost_terms[2]
    η = isa(ηbracket,AbstractMatrix) ? ηbracket[2,N] : ηbracket[2]
    cxkl,cukl,cxxkl,cxukl,cuukl = kl_cost_terms[1]

    @assert size(cx) == (n, N)
    @assert size(cu) == (m, N)
    @assert size(cxx) == (n, n, N)
    @assert size(cxu) == (n, m, N)
    @assert size(cuu) == (m, m, N)

    k          = zeros(m,N)
    K          = zeros(m,n,N)
    Vx         = zeros(n,N)
    Vxx        = zeros(n,n,N)
    Quu        = Array{T}(undef,m,m,N)
    Quui       = Array{T}(undef,m,m,N)
    dV         = [0., 0.]

    Vx[:,N]    = cx[:,N]
    Vxx[:,:,N] = cxx[:,:,end]
    Quu[:,:,N] = cuu[:,:,N]./ η .+ cuukl[:,:,N]
    Quui[:,:,N] = inv(Quu[:,:,N])

    diverge    = 0
    for i = N-1:-1:1
        Qu          = cu[:,i] + fu[:,:,i]'Vx[:,i+1]
        Qx          = cx[:,i] + fx[:,:,i]'Vx[:,i+1]
        Qux         = cxu[:,:,i]' + fu[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]
        Quu[:,:,i] .= cuu[:,:,i] .+ fu[:,:,i]'Vxx[:,:,i+1]*fu[:,:,i]
        Qxx         = cxx[:,:,i]  + fx[:,:,i]'Vxx[:,:,i+1]*fx[:,:,i]

        ηbracket = kl_cost_terms[2]
        η = isa(ηbracket,AbstractMatrix) ? ηbracket[2,i] : ηbracket[2]
        Qu         = Qu ./ η + cukl[:,i]
        Qux        = Qux ./ η + cxukl[:,:,i]
        Quu[:,:,i] = Quu[:,:,i] ./ η + cuukl[:,:,i]
        Qx         = Qx ./ η + cxkl[:,i]
        Qxx        = Qxx ./ η + cxxkl[:,:,i]

        Quu[:,:,i] = 0.5*(Quu[:,:,i] + Quu[:,:,i]')

        if isempty(lims) || lims[1,1] > lims[1,2]
            # debug("#  no control limits: Cholesky decomposition, check for non-PD")
            local R
            try
                R = cholesky(Hermitian(Quu[:,:,i]))
            catch
                diverge  = i
                return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx, dV
            end

            # debug("#  find control law")
            k_i = -(R\Qu)
            K_i = -(R\Qux)
        else
            # debug("#  solve Quadratic Program")
            lower = lims[:,1]-u[:,i]
            upper = lims[:,2]-u[:,i]
            local k_i,result,free
            try
                k_i,result,R,free = boxQP(Quu[:,:,i],Qu,lower,upper,k[:,min(i+1,N-1)])
            catch
                result = 0
            end
            if result < 1
                diverge  = i
                return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx, dV
            end
            K_i  = zeros(m,n)
            if any(free)
                Lfree         = -R\(R'\Qux[free,:])
                K_i[free,:]   = Lfree
            end
        end
        # debug("#  update cost-to-go approximation")

        dV         = dV  + [k_i'Qu; .5*k_i'Quu[:,:,i]*k_i]
        Vx[:,i]    = Qx  + K_i'Quu[:,:,i]*k_i + K_i'Qu + Qux'k_i
        Vxx[:,:,i] = Qxx + K_i'Quu[:,:,i]*K_i + K_i'Qux + Qux'K_i
        Vxx[:,:,i] = .5*(Vxx[:,:,i] + Vxx[:,:,i]')

        # debug("# save controls/gains")
        k[:,i]   = k_i
        K[:,:,i] = K_i
        Quui[:,:,i] = inv(Quu[:,:,i])
    end

    return diverge, GaussianPolicy(N,n,m,K,k,Quui,Quu), Vx, Vxx,dV
end


text:
using Documenter
using QuickPOMDPs

makedocs(
    sitename = "QuickPOMDPs",
    format = Documenter.HTML(),
    modules = [QuickPOMDPs],
    pages = ["index.md",
             "discrete_explicit.md",
             "quick.md"
            ]
)

deploydocs(
    repo = "github.com/JuliaPOMDP/QuickPOMDPs.jl.git"
)


text:
using SpinDoctor
using Documenter

DocMeta.setdocmeta!(SpinDoctor, :DocTestSetup, :(using SpinDoctor); recursive = true)

makedocs(;
    modules = [SpinDoctor],
    authors = "Syver Døving Agdestein <syverda@icloud.com> and contributors",
    repo = "https://github.com/agdestein/SpinDoctor.jl/blob/{commit}{path}#{line}",
    sitename = "SpinDoctor.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://agdestein.github.io/SpinDoctor.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Examples" => [
            "Solve BTPDE" => "examples/solve_btpde.md",
            "Custom gradients" => "examples/custom_gradient.md",
            "Compare ADCs" => "examples/compute_adc.md",
            "High Angular Resolution" => "examples/hardi.md",
        ],
        "Theory" => [
            "Definitions" => "theory/definitions.md",
            "BTPDE" => "theory/btpde.md",
            "ADC" => "theory/adc.md",
            "Matrix Formalism" => "theory/matrix_formalism.md",
            "Discretization" => "theory/discretization.md",
        ],
        "Neurons" => "neurons.md",
        "API Reference" => "api.md",
    ],
)

deploydocs(; repo = "github.com/agdestein/SpinDoctor.jl", devbranch = "main")


text:
using LinearAlgebra, Jets, JetPack, Test

n1,n2 = 33,44

@testset "JopLog, correctness T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    F = JopLog(JetSpace(T,n1,n2))
    x1 = rand(domain(F)) .+ T(0.0001)
    @test F*x1 ≈ log.(x1)
end

@testset "JopLog, linearity test, T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    F = JopLog(JetSpace(T,n1,n2))
    J = jacobian(F, rand(domain(F)) .+ T(0.0001))
    lhs,rhs = linearity_test(J)
    @test lhs ≈ rhs
end

@testset "JopLog, dot product test, T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    F = JopLog(JetSpace(T,n1,n2))
    J  = jacobian!(F, rand(domain(F)) .+ T(0.0001))
    lhs, rhs = dot_product_test(J, -1 .+ 2 .* rand(domain(J)), -1 .+ 2 .* rand(range(J)))
    @test lhs ≈ rhs
end

# note the key here is to increase the size of the nonlinear vector
@testset "JopLog, linearization test, T=$(T)" for T in (Float64,Float32,Complex{Float64},Complex{Float32})
    F = JopLog(JetSpace(T,n1,n2))
    m0 = 1000 .* rand(domain(F)) .+ T(1)
    μ  = sqrt.([1/1,1/2,1/4,1/8,1/16,1/32,1/64,1/128,1/256,1/512,1/1024,1/2048])
    δm = rand(domain(F)) .+ T(0.0001)
    observed, expected = linearization_test(F, m0, μ = μ, δm = δm)
    δ = minimum(abs, observed - expected)
    #write(stdout, @sprintf("\nLinearization test -- type(%s)\n", T))
    #for i = 1:length(observed)
        #write(stdout, @sprintf("mu,observed,expected,diff; %12.6f %12.6f %12.6f %12.6f\n", μ[i], observed[i], expected[i], abs(observed[i] - expected[i])))
    #end
    #write(stdout, @sprintf("minimum difference %12.6f\n", minimum(abs,observed .- expected)))
    @test δ < 0.1
end


text:
mutable struct MovingWindowSolution{windows_T} <: InitialValueSolution
    windows::windows_T
end

function MovingWindowSolution(problem, solver::MovingWindowSolver)
    @↓ (t0, tN) ← tspan = problem
    @↓ τ, Δτ = solver
    M = 1 # if solver.Δτ == 0 then T = τ and M = 1
    if Δτ > 0
        M = trunc(Int, (tN - τ) / Δτ) + 1
    end
    windows = Vector{TimeParallelSolution}(undef, M)
    return MovingWindowSolution(windows)
end

Base.length(solution::MovingWindowSolution) = length(solution.windows)
Base.getindex(solution::MovingWindowSolution, m::Int) = solution.windows[m] # solution[m] ≡ solution.windows[m]
Base.setindex!(solution::MovingWindowSolution, value::TimeParallelSolution, m::Int) = solution.windows[m] = value
Base.lastindex(solution::MovingWindowSolution) = lastindex(solution.windows)


text:
"""
Documentation goes here
"""
function new_slic(img, k)

    no_rows, no_columns = size(img)
    n = length(img)
    s = round(Int,sqrt(n/k))

    img_height, img_width = size(img)
    error = 0.0

    # Initalise arrays for the first time
    counts = zeros(k)
    totals = zeros(1,k,5)
    centers = Array{CartesianIndex}(undef,k)
    center_colors = Array{Lab{Float64}}(undef,k)
    pixel_associactions = zeros(Int,img_height,img_width)

    distances = fill(Inf,size(img))

    # Initial stuff
    centers=initialise_centers(centers,s,k,no_rows, no_columns)
    centers=seed_locations(centers,img)

    for i=1:9
        pixel_associactions,counts,distances,totals=calculate_distances(centers,img,pixel_associactions,counts,distances,totals,s)
        centers,center_colors,img = calc_new_centers(centers,center_colors,counts,totals,img)
        fill!(counts,0)
        fill!(totals,0.0)
        fill!(pixel_associactions,0)
        fill!(distances,Inf)
    end

    pixel_associactions,counts,distances,totals=calculate_distances(centers,img,pixel_associactions,counts,distances,totals,s)
    centers,center_colors,img=calc_new_centers(centers,center_colors,counts,totals,img)

    # Create the superpixels
    img2 = similar(img)
    #img2 = copy(img)
    for i in CartesianIndices(img2)
        if pixel_associactions[i] !=0
            img2[i] = center_colors[pixel_associactions[i]]
        end
    end
    for i in centers
        if i != CartesianIndex(0,0)
            img2[i] = Lab{Float64}(0,128,-128)
        end
    end
    return img2

end

function calc_gradient(img)::Array{Float64}
    rgb_img = RGB.(img)
    grad_y, grad_x = imgradients(rgb_img,KernelFactors.sobel, "replicate")
    magnitudes = magnitude(grad_x, grad_y)
    return magnitudes
end

# Calclates colour distance between 2 pixels.
function d_lab(𝑙ₖ::Float64, 𝑙ᵢ::Float64, 𝑎ₖ::Float64, 𝑎ᵢ::Float64, 𝑏ₖ::Float64, 𝑏ᵢ::Float64)
    return sqrt( (𝑙ₖ - 𝑙ᵢ)^2 + (𝑎ₖ - 𝑎ᵢ)^2 + (𝑏ₖ - 𝑏ᵢ)^2 )
end
# Calculates xy distance between 2 pixels.
function d_xy(𝑥ₖ::Real, 𝑥ᵢ::Real, 𝑦ₖ::Real, 𝑦ᵢ::Real)
    return sqrt( (𝑥ₖ - 𝑥ᵢ)^2 + (𝑦ₖ - 𝑦ᵢ)^2 )
end

# Calculate the distance measure between two CartesianIndexes.
function calc_dist(img,center::CartesianIndex, px::CartesianIndex,s)
    return d_lab(img[center].l, img[px].l,img[center].a, img[px].a,img[center].b, img[px].b) + (10 * d_xy(center[2], px[2], center[1], px[1]) ) / s
end

#=Calculates the top-left and bottom-right bounds of a neighbourhood
of width w x w around the pixel passed in.
=#
function calc_neighbourhood(pixel::CartesianIndex, w,img)
    return CartesianIndices((max(pixel[1] - w, 1):min(pixel[1] + w, size(img)[1]), max(pixel[2] - w, 1): min(pixel[2] + w, size(img)[2])))
end

function calculate_distances(centers,img,pixel_associactions,counts,distances,totals,s)
    temp_dist = 0.0
    for i = 1:size(centers)[1]
        if centers[i] != CartesianIndex(0,0)
            a = img[centers[i]]
            c1 = centers[i][2]
            c2 = centers[i][1]
            search_area = calc_neighbourhood(centers[i],s,img)
            # Initialise shortest distance
            for pixel in search_area
                b = img[pixel]
                dlab=sqrt((a.l-b.l)^2+(a.a-b.a)^2+(a.b-b.b)^2)
                dxy=sqrt((c2-pixel[2])^2+(c1-pixel[1])^2)
                temp_dist =  dlab + (1/s) * dxy
                if temp_dist < distances[pixel]
                    pixel_data = img[pixel]
                    # IF the pixel was previously associated with another center
                    if pixel_associactions[pixel] != 0
                        # Store the previously associated center
                        prev_center = pixel_associactions[pixel]
                        # Take away the count
                        counts[prev_center] -= 1
                        # How to do this in a more efficient way?
                        # Remove pixel's data
                        totals[1,prev_center,1] -= pixel_data.l
                        totals[1,prev_center,2] -= pixel_data.a
                        totals[1,prev_center,3] -= pixel_data.b
                        totals[1,prev_center,4] -= pixel[1]
                        totals[1,prev_center,5] -= pixel[2]
                    end
                    # Associate pixel with new center
                    pixel_associactions[pixel] = i
                    # Update shortest distance
                    distances[pixel] = temp_dist
                    # Update counts
                    counts[i] +=1
                    # Update average calculations
                    totals[1,i,1] += pixel_data.l
                    totals[1,i,2] += pixel_data.a
                    totals[1,i,3] += pixel_data.b
                    totals[1,i,4] += pixel[1]
                    totals[1,i,5] += pixel[2]
                end
            end
        end
    end
    for i in eachindex(counts)
        if counts[i]<=0
            centers[i]=CartesianIndex(0,0)
        end
    end
    return pixel_associactions,counts,distances,totals
end

function calc_new_centers(centers,center_colors,counts,totals,img)
    # Compute new centers
    for i in eachindex(centers)
        if centers[i] != CartesianIndex(0,0)
            # Assign new centers to their respective arrays.
            center_colors[i] = Lab{Float64}(totals[1,i,1]/counts[i],totals[1,i,2]/counts[i],totals[1,i,3]/counts[i])
            centers[i] = CartesianIndex(round(Int,totals[1,i,4]/counts[i]),round(Int,totals[1,i,5]/counts[i]))
        end
    end
    return centers,center_colors,img
end

#= Generate starting locations for centers, at regularly
spaced intervals.
=#
function initialise_centers(centers,s,k,no_rows, no_columns)
    row_pos = col_pos = 1
    for i = 1:k
        col_pos += s
        if col_pos > no_columns
            row_pos = min(row_pos +=s, no_rows)
            col_pos = 1
        end
        centers[i] = CartesianIndex(row_pos,col_pos)
    end
    return centers
end

function seed_locations(centers,img)
    # Calculate gradients
    gradients = calc_gradient(img)
    min_gradient = CartesianIndex(0,0)
    for i in eachindex(centers)
        indexes = CartesianIndices((max(centers[i][1] - 1, 1):min(centers[i][1] + 1, size(img)[1]), max(centers[i][2] - 1, 1): min(centers[i][2] + 1, size(img)[2])))
        min_gradient = findmin(gradients[indexes])[2]
        centers[i] = CartesianIndex(min_gradient[1]+max(centers[i][1] - 1, 0),min_gradient[2]+max(centers[i][2] - 1, 0))
        #img[centers[i]]=Lab{Float64}(100,-128,-128)
    end
    return centers
end


text:
module AllocatorModel

using IterTools

export  allocate,
        free,
        isfree,
        MemoryAllocator,
        MemoryNode,
        showfree

# moDNN Memory allocater
mutable struct MemoryNode
    isfree::Bool
    size::Int64
end
isfree(m::MemoryNode) = m.isfree
Base.sizeof(m::MemoryNode) = m.size

struct MemoryAllocator
    node_list::Vector{MemoryNode}
    alignment::Int64
end

MemoryAllocator(limit::Number, alignment::Number) =
    MemoryAllocator(
        [MemoryNode(true, convert(Int, limit))],
        convert(Int, alignment)
    )

function showfree(M::MemoryAllocator)
    offset = 0
    for node in M.node_list
        if isfree(node)
            @show (offset, sizeof(node))
        end
        offset += sizeof(node)
    end
    return nothing
end

function allocate(M::MemoryAllocator, sz)
    # Make the allocation size match the alignment
    sz = align(sz, M.alignment)
    offset = 0
    min_delta = typemax(Int64)
    best_offset = offset

    best_node = nothing
    best_index = zero(Int64)

    # Go through the list see if we can find the an exact fit
    for (index, node) in enumerate(M.node_list)
        if isfree(node) && sizeof(node) >= sz
            delta = sizeof(node) - sz
            if delta < min_delta
                min_delta   = delta
                best_node   = node
                best_offset = offset
                best_index  = index
            end
        end
        offset += sizeof(node)
    end

    # Check if we were able to allocate.
    isnothing(best_node) && return nothing

    # Do the allocation
    if iszero(min_delta)
        best_node.isfree = false
    else
        insert!(M.node_list, best_index, MemoryNode(false, sz))
        best_node.size -= sz
    end
    return best_offset
end

free(M::MemoryAllocator, ::Nothing) = nothing
function free(M::MemoryAllocator, offset)
    search_offset = 0
    found = false
    for node in M.node_list
        if (offset == search_offset)
            # Mark this node as free
            node.isfree = true
        end
        search_offset += sizeof(node)
    end

    # Coalesce free nodes together
    while true
        changed = false
        for (index, (a, b)) in enumerate(IterTools.partition(M.node_list, 2, 1))
            # If both nodes are free, merge the size into b and delete a
            #
            # `index` is the index of `a`.
            if isfree(a) && isfree(b)
                b.size += sizeof(a)
                deleteat!(M.node_list, index)
                changed = true
                break
            end
        end
        changed || break
    end
    return nothing
end

function align(sz, alignment)
    if iszero(sz)
        sz = alignment
    else
        remainder = rem(sz, alignment)
        if remainder > zero(remainder)
            sz += (alignment - remainder)
        end
    end
    return sz
end

end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
list_test_str = split(
    "00100
11110
10110
10111
10101
01111
00111
11100
10000
11001
00010
01010"
)

list_str = readlines("day3_input.txt")

function part1(list)
    mat = reshape([c=='1' for l in list for c in l], length(list[1]), :)
    bin = sum(mat, dims=2) .> div(size(mat, 2), 2)
    gamma = sum(bin[end-j]*2^j for j in 0:length(bin)-1)
    eps = 2^size(mat, 1) - 1 - gamma
    return gamma * eps
end

@info "Part 1"
@info "Test"
@show part1(list_test_str) == 198
@info "Challenge"
@show part1(list_str) == 1131506

# Part 2

function part2(list)
    mat = reshape([c=='1' for l in list for c in l], length(list[1]), :)

    j = 1
    cmat = copy(mat)
    while j ≤ size(mat,1) && size(cmat,2) > 1
        num_ones = sum(cmat[j,:])
        bit = sum(cmat[j,:]) ≥ size(cmat,2) - num_ones
        inds = [i for i in 1:size(cmat,2) if cmat[j,i] == bit]
        cmat = cmat[:,inds]
        j += 1
    end

    oxygen = sum([cmat[end+1-j]*2^(j-1) for j in eachindex(cmat)])

    j = 1
    cmat = mat[:,:]
    while j ≤ size(mat,1) && size(cmat,2) > 1
        num_ones = sum(cmat[j,:])
        bit = sum(cmat[j,:]) ≥ size(cmat,2) - num_ones
        inds = [i for i in 1:size(cmat,2) if cmat[j,i] != bit]
        cmat = cmat[:,inds]
        j += 1
    end

    scrubber = sum([cmat[end+1-j]*2^(j-1) for j in eachindex(cmat)])

    return oxygen * scrubber
end

@info "Part 2"
@info "Test"
@show part2(list_test_str) == 230
@info "Challenge"
@show part2(list_str) == 7863147

nothing


text:
abstract type RequestTrigger end
abstract type RegisterTrigger <: RequestTrigger end

struct PullRequestTrigger <: RegisterTrigger
    prid::Int
end

struct IssueTrigger <: RegisterTrigger
    branch::String
end

struct CommitCommentTrigger <: RegisterTrigger
end

struct ApprovalTrigger <: RequestTrigger
    prid::Int
end

struct EmptyTrigger <: RequestTrigger
end


text:
# introduce embedding between matrices (vectors) of different dimensions

"""
    mat_embed(X, dim; zero_pos = -1, tpin = "her", tpout = "her")

This function embeds X into new dimensions dim, with the Cartesian
order of the entries are preserved. Fill in zeros at zero_pos if necessary.

# Arguments
-   `X`: a Hermitian or rectangular matrice
-   `dim`: target dimension, array-like [Nnew,Mnew]
-   `zero_pos`: an array of user specify Cartesian position for zeros. if
               it is -1, then use randomly generated position.
-   `tpin`: the type of input matrix, valid options are "her" and "rec".
-   `tpout`: the type of ouput matrix, valid options are "her" and "rec"

# Outputs
-   `Xnew`: a matrix of dimension dim and the type tpin. The entries
          Xnew are from X, with possible extra zeros at zero_pos,
"""
function mat_embed(X, dim; zero_pos = -1, tpin = "her", tpout = "her")

    # get the original dimension
    X = X[:,:]
    N, M = size(X)

    # get the target dimension
    Nnew, Mnew = dim

    # check the input dimension if mat is Hermitian
    if (tpin == "her") & (N != M)
        error("Input dimension is not for a Hermitian matrix!")
    end

    # reshape X into vector according to type in
    if tpin == "her"
        in_size = N*(N + 1)/2
        Xvec = [X[i,j] for i=1:N for j=i:N]
    elseif tpin == "rec"
        in_size = N*M
        Xvec = vec(X)
    end

    # check whether there is enough room in the target shape
    if tpout == "her"
        out_size = Int((Nnew - 1)*Nnew/2)
    elseif tpout == "rec"
        out_size = Nnew*Mnew
    end

    if out_size < in_size
        error("The target dimension is too small!")
    end

    # if the position of zeros are not preset, generate one
    if zero_pos == -1
        zero_pos = sample(1:out_size, out_size - in_size, replace=false)
    end

    # Construct Xnewvec
    Xnewvec = zeros(out_size)
    Xnewvec[[i for i in 1: out_size if ~(i in zero_pos)]] = Xvec

    # reshape Xnewvec to Xnew
    if tpout == "her"
        k = 0
        Xnew = [i < j ? (k += 1; Xnewvec[k]) : 0 for i = 1: Nnew, j = 1: Mnew]
        Xnew = Xnew + Xnew'
    elseif tpout == "rec"
        Xnew = reshape(Xnewvec, (Nnew, Mnew))
    end

    return convert(Array{typeof(X[1,1]), 2}, Xnew)
end


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

mutable struct IoK8sApiAutoscalingV2beta2ResourceMetricSource <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    target::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricTarget } # spec name: target

    function IoK8sApiAutoscalingV2beta2ResourceMetricSource(;name=nothing, target=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2ResourceMetricSource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAutoscalingV2beta2ResourceMetricSource, Symbol("target"), target)
        setfield!(o, Symbol("target"), target)
        o
    end
end # type IoK8sApiAutoscalingV2beta2ResourceMetricSource

const _property_map_IoK8sApiAutoscalingV2beta2ResourceMetricSource = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("target")=>Symbol("target"))
const _property_types_IoK8sApiAutoscalingV2beta2ResourceMetricSource = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("target")=>"IoK8sApiAutoscalingV2beta2MetricTarget")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2ResourceMetricSource }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2ResourceMetricSource))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2ResourceMetricSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2ResourceMetricSource[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2ResourceMetricSource }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2ResourceMetricSource[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2ResourceMetricSource)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("target")) === nothing) && (return false)
    true
end

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


text:
import Base.+
Base.:+(s::AbstractString, n::Real) = string((x = tryparse(Int, s)) isa Int ? x + 1 : parse(Float64, s) + 1)
@show "125" + 1
@show "125.15" + 1
@show "1234567890987654321" + 1

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Tullio, LoopVectorization

I, J, K = 400, 200, 100
f1 = rand(I, 1, K)
f2 = rand(I, J, 1)
dev = rand(I, 1, 1)

function make_R(f1, f2, dev)
    I, J, K = axes(dev, 1), axes(f2, 2), axes(f1, 3)
    R = similar(dev, I, J, K)
    for i in I, j in J, k in K
        @inbounds R[i, j, k] = f1[i, 1, k] * f2[i, j, 1] / dev[i, 1, 1]
    end
    R
end

R1 = make_R(f1, f2, dev)
@tullio R2[i, j, k] := f1[i, 1, k] * f2[i, j, 1] / dev[i, 1, 1]
R3 = @. f1 * f2 / dev
R4 = @turbo @. f1 * f2 / dev
R5 = @tturbo @. f1 * f2 / dev

@show R1 ≈ R2 ≈ R3 ≈ R4 ≈ R5
@show typeof(R1) size(R1);

# %%
using BenchmarkTools
@btime R1 = make_R($f1, $f2, $dev)
@btime @tullio R2[i, j, k] := $f1[i, 1, k] * $f2[i, j, 1] / $dev[i, 1, 1]
@btime R3 = @. $f1 * $f2 / $dev
@btime R4 = @turbo @. $f1 * $f2 / $dev
@btime R5 = @tturbo @. $f1 * $f2 / $dev
;

# %%

max_stars_repo_path: 0009/Tullio.jl and LoopVectorization.jl example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using RvSpectML
 using Statistics
 using Dates

make_plots = false
include("expres_1_read.jl")
order_list_timeseries = RvSpectML.make_order_list_timeseries(expres_data)
order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
lambda_range_with_data = (min = maximum(d->minimum(d.λ[expres_data[1].metadata[:excalibur_mask]]),expres_data), max = minimum(d->maximum(d.λ[expres_data[1].metadata[:excalibur_mask]]),expres_data) )
RvSpectML.normalize_spectra!(order_list_timeseries,expres_data);

espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks","G2.espresso.mas")
 espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
 line_list_df = espresso_df |>
     @filter(lambda_range_with_data.min <= _.lambda ) |>
    @filter( _.lambda < lambda_range_with_data.max) |>
#    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
#    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
    DataFrame

# Setup to run CCF
mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=1.6)
 line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
 ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=-5e3)
 v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)

# Compute CCF's & measure RVs
tstart = now()
 @time ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
 println("# CCF runtime: ", now()-tstart)
 # Write CCFs to file
 using CSV
 CSV.write("101501_ccfs.csv",Tables.table(ccfs',header=Symbol.(v_grid)))

make_plots = true
# Plot CCFs (normalized by their maximimum)
if make_plots
  using Plots
 plot(v_grid,ccfs./maximum(ccfs,dims=1),label=:none)
 xlabel!("v (m/s)")
 ylabel!("CCF")
end

# Compute RVs
make_plots = true
 rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "gaussian") for i in 1:length(order_list_timeseries) ]
 rvs_ccf_quad  = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i], fit_type = "quadratic") for i in 1:length(order_list_timeseries) ]
 println("RMS of RVs (ESPRESSO lines):  Gaussian: ",std(rvs_ccf_gauss), "   Quadratic: ", std(rvs_ccf_quad))
 # Write to RVs to CSV file
 CSV.write("101501_rvs_ccf_std.csv",DataFrame("Time [MJD]"=>order_list_timeseries.times,"CCF RV [m/s]"=>rvs_ccf_gauss))
 if make_plots
   t_idx = 1
   using Plots
   plot(v_grid,ccfs[:,t_idx],label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
 end

rvs_yale = CSV.read(joinpath(expres_data_path,target_subdir,"101501_activity.csv"), DataFrame)

plt1 = plot()
 rv_ref = mean(rvs_yale[!,"CBC RV [m/s]"])
 scatter!(plt1,rvs_yale[!,"Time [MJD]"],rvs_yale[!,"CBC RV [m/s]"].-rv_ref,yerr=rvs_yale[!,"CBC RV Error [m/s]"], color=1, label="Yale CBC")
 rv_ref = mean(rvs_yale[!,"CCF RV [m/s]"])
 scatter!(plt1,rvs_yale[!,"Time [MJD]"],rvs_yale[!,"CCF RV [m/s]"].-rv_ref,yerr=rvs_yale[!,"CCF RV Error [m/s]"], color=2, label="Yale CCF")
 times_plt = order_list_timeseries.times .- 40000
 rv_ref = mean(rvs_ccf_gauss)
 scatter!(plt1,times_plt,rvs_ccf_gauss.-mean(rvs_ccf_gauss),color=3,label="Gaussian")
 plt2 = plot( legend=:topleft)
 rv_ref = rvs_yale[!,"CCF RV [m/s]"].-mean(rvs_yale[!,"CCF RV [m/s]"]).+mean(rvs_yale[!,"CBC RV [m/s]"])
 scatter!(plt2,rvs_yale[!,"Time [MJD]"],rvs_yale[!,"CBC RV [m/s]"].-rv_ref,yerr=rvs_yale[!,"CBC RV Error [m/s]"], color=1, label="Yale CBC-Yale CCF")
 #rv_ref = rvs_yale[!,"CCF RV [m/s]"].-mean(rvs_yale[!,"CCF RV [m/s]"]).+mean(rvs_yale[!,"CCF RV [m/s]"])
 #scatter!(plt2,rvs_yale[!,"Time [MJD]"],rvs_yale[!,"CCF RV [m/s]"].-rv_ref,yerr=rvs_yale[!,"CCF RV Error [m/s]"], color=2, label="Yale CCF-Yale CCF")
 times_plt = order_list_timeseries.times .- 40000   # Why do I need to subtract this off?
 rv_ref = rvs_yale[!,"CCF RV [m/s]"].-mean(rvs_yale[!,"CCF RV [m/s]"]).+mean(rvs_ccf_gauss)
 scatter!(plt2,times_plt,rvs_ccf_gauss.-rv_ref,color=3,label="RvSpecML CCF - Yale CCF")
 plot(plt1,plt2,layout=(2, 1))

 plt1 = plot()# legend=:topleft)
  plt2 = plot()# legend=:topleft)
  rv_ref1 = mean(rvs_yale[!,"CCF RV [m/s]"])
  rv_ref2 = mean(rvs_ccf_gauss)
  scatter!(plt1,rvs_yale[!,"CCF RV [m/s]"].-rv_ref1,rvs_ccf_gauss.-rv_ref2,xerr=rvs_yale[!,"CCF RV Error [m/s]"],yerr=rvs_yale[!,"CCF RV Error [m/s]"],color=1, label=:none)
  xlabel!(plt1,"RV (m/s) from Yale's CCF")
  ylabel!(plt1,"RV (m/s) from RvSpectML's CCF")
  x_line = range(minimum(rvs_yale[!,"CCF RV [m/s]"]),stop=maximum(rvs_yale[!,"CCF RV [m/s]"]),length=20).-mean(rvs_yale[!,"CCF RV [m/s]"])
  plot!(plt1,x_line,x_line, color=6,label=:none)
  scatter!(plt2,rvs_yale[!,"CCF RV [m/s]"].-rv_ref1,rvs_ccf_gauss.-rv_ref2.-(rvs_yale[!,"CCF RV [m/s]"].-rv_ref1),xerr=rvs_yale[!,"CCF RV Error [m/s]"],yerr=rvs_yale[!,"CCF RV Error [m/s]"],color=1, label=:none)
  xlabel!(plt2,"RV (m/s) from Yale's CCF")
  ylabel!(plt2,"ΔRV (m/s)")
  plot(plt1,plt2,layout=(2, 1))

# Plot Nightly mean RVs
if make_plots
   nbin = 4
   plt_t = (order_list_timeseries.times .- minimum(order_list_timeseries.times) )
   plt = scatter(plt_t,rvs_ccf_gauss.-mean(rvs_ccf_gauss),label="Gaussian")
   scatter!(plt,plt_t,rvs_ccf_quad.-mean(rvs_ccf_quad),label="Quadratic")
   times_binned = RvSpectML.bin_times_nightly(plt_t)
   rvs_ccf_gauss_binned = RvSpectML.bin_rvs_nightly(times=plt_t,rvs=rvs_ccf_gauss)
   rvs_ccf_quad_binned = RvSpectML.bin_rvs_nightly(times=plt_t,rvs=rvs_ccf_quad)
   println("# RMS of nightly binned RVs: ", std(rvs_ccf_gauss_binned), "  ", std(rvs_ccf_quad_binned) ) #, "  ", std(rvs_ccf_cent_binned))
   scatter!(times_binned,rvs_ccf_gauss_binned.-mean(rvs_ccf_gauss),label="Binned Gauss")
   scatter!(times_binned,rvs_ccf_quad_binned.-mean(rvs_ccf_quad),label="Binned Quad")
   xlabel!(plt,"Time (days)")
   ylabel!(plt,"RV (m/s)")
   display(plt)
 end


text:
#!/bin/env julia
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.test("NodeJS") # This forces .julia/artifacts to gain a node installation

max_stars_repo_name: PolyCement/gcvt

max_issues_repo_name: PolyCement/gcvt

max_forks_repo_name: PolyCement/gcvt


text:
using FlightGNC
using Test

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