text:
function MakieCore.convert_arguments(P::DiscreteSurface, f::ScalarField)
    @debug "ScalarField SurfaceLike"
    grid, data = unwrap(f)

    convert_arguments(P, grid..., data)
end

function MakieCore.convert_arguments(P::VolumeLike, f::ScalarField{3})
    @debug "ScalarField VolumeLike"
    grid, data = unwrap(f)

    convert_arguments(P, grid..., data)
end

function MakieCore.convert_arguments(P::Type{<:Arrows}, f::VectorField{N}) where N
    @debug "VectorField Arrows"
    _f = hasunit(f) ? ustrip(f) : f

    grid = getdomain(_f)
    data = unwrapdata(_f)

    origins = vec(mapgrid(Point{N,Float32}, grid))
    arrowheads = vec(Vec{N,Float32}.(components(data)...))

    convert_arguments(P, origins, arrowheads)
end

function MakieCore.convert_arguments(P::Type{<:Arrows}, v::VectorVariable)
    @debug "VectorVariable Arrows"
    _v = hasunit(v) ? ustrip(v) : v
    N = dimensionality(v)

    grid = getdomain(_v)
    data = unwrapdata(_v)
    origins = mapgrid(Point{N,Float32}, grid)
    arrowheads = Vec{N,Float32}.(components(data)...)

    convert_arguments(P, origins, arrowheads)
end

function MakieCore.convert_arguments(P::Type{<:Contour}, f::ScalarField{3})
    @debug "ScalarField{3} with plot type $P"
    _f = hasunit(f) ? ustrip(f) : f
    grid, data = unwrap(_f)

    # 3D contour plots need the VolumeLike trait
    # plotting from the recipe creates a Contour{T} where T which defaults to SurfaceLike
    # In order to fix the problem, we create the correct type here
    P_fixed = Contour{Tuple{typeof.(grid)..., typeof(data)}}
    convert_arguments(P_fixed, grid..., data)
end

function MakieCore.convert_arguments(P::PointBased, g::ParticlePositions)
    @debug "PointBased ParticlePositions"
    convert_arguments(P, g...)
end

ext: jl


text:
# expcone.jl -- part of the Information Decomposition project
module ExpCone

#incude("infdecomp_base.jl")
using InfDecomp_Base
using MathProgBase
#using CPUTime

type Solution_Stats
    optimum     :: Float64
    primal_val  :: Float64
#    primal_feas :: Float64
#    dual_val    :: Float64
#    dual_feas   :: Float64
    q           :: Vector{Float64}
    time        :: Float64
    model       :: MathProgBase.AbstractMathProgModel
end

function model_L(D::InfDecomp_Base.My_Eval, solver::MathProgBase.AbstractMathProgSolver) :: Solution_Stats
    # min c^T v     #  min (-1,0,0)^T (r^T,q^T,p^T)
    # s.t.          #  s.t.
    # b - Av = 0    #  q_{*xy} = p_{xyz} ;  + marginal eqns
    # v ∈ EXP       #  (r,q,p) ∈ EXP

    local n_vars::Int = D.n*3            # r_{xyz}, q_{xyz}, p_{xyz} ;  ∀ xyz ∈ Variables
                                         # r=3*i-2 ; q=3*i-1 ; p=3*i
    local n_cons::Int = D.n + D.m        #    0          q_{*yz} - p_{xyz}  = 0     ; ∀ xyz
                                         # +  b^y_{xy} - q_{xy*}            = 0     ; ∀ xy
                                         # +  b^z_{xz} - q_{x*z}            = 0     ; ∀ xz

    function make_c_A_b() :: Tuple{Vector{Float64},SparseMatrixCSC{Float64},Vector{Float64}}
        local c::Vector{Float64} = zeros(n_vars)
        local b::Vector{Float64} = zeros(n_cons)

        # fill c: -1 on the r_{xyz} vars:
        for xyz = 1:D.n
            c[3*xyz-2] = -1.
        end

        # fill b:  (0, b^y, b^z)
        for i = 1:D.m
            b[ D.n + i ]  = D.rhs[i]
        end

        #
        # fill A
        # ------

        Eqn   = Array{Int,1}()
        Var   = Array{Int,1}()
        Coeff = Array{Float64,1}()

        # q=p equations:
        for xyz = 1:D.n
            eqn = xyz

            p_xyz = 3*xyz
            append!(Eqn,   eqn)
            append!(Var,   p_xyz)
            append!(Coeff, -1.)

            (x,y,z) = D.xyz[xyz]
            for uvw = 1:D.n
                (u,v,w) = D.xyz[uvw]
                if (v,w)==(y,z)
                    q_uvw = 3*uvw-1
                    append!(Eqn,   eqn)
                    append!(Var,   q_uvw)
                    append!(Coeff, 1.)
                end
            end
        end

        # q_{xy*} = b^y_{xy}
        (K,L,V) = findnz(D.Gt)
        for eq = 1:length(K)
            append!(Eqn,   D.n+L[eq])
            append!(Var,   3*K[eq]-1)
            append!(Coeff, V[eq])
        end
        A = sparse(Eqn,Var,Coeff)
        @show A
        return c,A,b
    end #^ make_c_b_A()

    c,A,b = make_c_A_b()

    constr_cones = [ ( :Zero, 1:n_cons ) ]
    var_cones = []
    for xyz = 1:D.n
        rqp_xyz = [3*xyz-2,3*xyz-1,3*xyz]
        append!(var_cones, [(:ExpPrimal, rqp_xyz)] )
    end

    m = MathProgBase.ConicModel(solver)
    MathProgBase.loadproblem!(m, c,A,b, constr_cones, var_cones)

    # CPUtic()                        # TIC

    MathProgBase.optimize!(m)

    # stats_time = CPUtoq()           # TOQ

    @show MathProgBase.status(m)
    @show stats_time = MathProgBase.getsolvetime(m)

    if MathProgBase.status(m) == :Optimal || MathProgBase.status(m) == :Suboptimal || MathProgBase.status(m) == :Error
        q = Vector{Float64}(D.n)
        for xyz in 1:D.n
            q[xyz] = MathProgBase.getsolution(m)[ 3*xyz-1 ]
        end

        eval_val::Float64 = InfDecomp_Base.condEntropy(D,q,big(0.))

        @show eval_val, MathProgBase.getobjval(m)
        if abs(eval_val + MathProgBase.getobjval(m)) > 1.e-3
            println("BIG DIFFERENCE IN SOLUTION VALUE!! :-( :-( :-( :-( :-( :-(\nProbably a BUG!!!!!!!!!!!!!!!!!!!!!")
        end

        stats = Solution_Stats(-eval_val,MathProgBase.getobjval(m),q,stats_time,m)
        return stats
    else
        print("Fuck!!")
        stats = Solution_Stats(999,999,[],stats_time,m)
        return stats
    end
    ; #
end #^ model_L

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

function model_S(D::InfDecomp_Base.My_Eval, solver::MathProgBase.AbstractMathProgSolver) :: Solution_Stats
    #    min (-1,0,0)(r,q,p)
    # s.t.
    #    q_{*yz} - p_{yz}         = 0     p/q coupling
    #    b - Gq                   = 0     marginal eqns
    #    (r_{xyz},q_{xyz},p_{yz}) ∈ EXP   linear EXPCONE_ieq
    #
    #    rqp free variables

    # count yz's
    local idx_of_yz = Dict{Tuple{Int64,Int64}, Int64}()

    local n_yz = 0
    for xyz = 1:D.n
        (x,y,z) = D.xyz[xyz]
        if (y,z) ∉ keys(idx_of_yz)
            n_yz += 1
            idx_of_yz[(y,z)] = n_yz
        end
    end

    local yz_of_idx = Vector{Tuple{Int64,Int64}}(n_yz)

    for (yz,i) in idx_of_yz
        yz_of_idx[i] = yz
    end

    local n_vars::Int = D.n + D.n + n_yz      # r_{xyz} ∀xyz;  q_{xyz} ∀xyz;  p_{yz} ∀yz
                                              # r=1...n; q=n+1...2n; p=2n+1...2n+n_yz
    local n_cons::Int = n_yz + D.m + 3*D.n    #   0           - ( q_{*yz}  -p_{yz}        ) ∈ {0}   ; ∀ yz
                                              #   / b^y_{xy}  - ( q_{xy*}                 ) ∈ {0}   ; ∀ xy  \  | D.m
                                              #   \ b^z_{xz}  - ( q_{x*z}                 ) ∈ {0}   ; ∀ xz  /  |
                                              #   0           - (-r_{xyz},-q_{xyz},-p_{yz}) ∈ EXPC  ; ∀ xyz

    function make_c_A_b() :: Tuple{Vector{Float64},SparseMatrixCSC{Float64},Vector{Float64}}
        local c::Vector{Float64} = zeros(n_vars)
        local b::Vector{Float64} = zeros(n_cons)

        # fill c: -1 on the r_{xyz} vars:
        for xyz = 1:D.n
            c[xyz] = -1.
        end

        # fill b:  (0, b^y b^z, 0)
        for i = 1:D.m
            b[ n_yz + i ]  = D.rhs[i]
        end

        #
        # fill A
        # ------

        Constr = Array{Int,1}()
        Var    = Array{Int,1}()
        Coeff  = Array{Float64,1}()

        # q=p equations:
        for yz = 1:n_yz
            eqn = yz
            p_yz = 2*D.n + yz

            append!(Constr, eqn)
            append!(Var,    p_yz)
            append!(Coeff,  -1.)

            (y,z) = yz_of_idx[yz]

            for uvw = 1:D.n
                (u,v,w) = D.xyz[uvw]
                if (v,w) == (y,z)
                    q_uvw = D.n + uvw
                    append!(Constr, eqn)
                    append!(Var,    q_uvw)
                    append!(Coeff,  1.)
                end
            end
        end

        # q_{xy*} = b^y_{xy}
        (K,L,V) = findnz(D.Gt)
        for eq = 1:D.m
            append!(Constr, n_yz + L[eq])
            append!(Var,    D.n + K[eq])
            append!(Coeff,  V[eq])
        end

        # EXPC ieqs:
        for xyz = 1:D.n
            cns_base = n_yz + D.m + 3*xyz-3
            r_xyz = xyz
            q_xyz = D.n + xyz
            (x,y,z) = D.xyz[xyz]
            p_yz = 2*D.n + idx_of_yz[ (y,z) ]

            append!(Constr, cns_base +1)
            append!(Var,    r_xyz)
            append!(Coeff,  -1.)

            append!(Constr, cns_base +2)
            append!(Var,    q_xyz)
            append!(Coeff,  -1.)

            append!(Constr, cns_base +3)
            append!(Var,    p_yz)
            append!(Coeff,  -1.)
        end

        A = sparse(Constr,Var,Coeff)

        return c,A,b
    end #^ make_c_b_A()

    (c,A,b) = make_c_A_b()

    var_cones = [ ( :Free, 1:n_vars) ]

    constr_cones = []
    append!(constr_cones, [ ( :Zero, 1:n_yz             ) ] )
    append!(constr_cones, [ ( :NonNeg, (n_yz+1):(n_yz+D.m)) ] )
    for xyz = 1:D.n
        rqp = [ n_yz+D.m+3*xyz-2, n_yz+D.m+3*xyz-1, n_yz+D.m+3*xyz ]
        append!(constr_cones, [(:ExpPrimal, rqp)] )
    end

    m = MathProgBase.ConicModel(solver)
    MathProgBase.loadproblem!(m, c,A,b, constr_cones, var_cones)

    MathProgBase.optimize!(m)

    @show MathProgBase.status(m)

    if MathProgBase.status(m) == :Optimal
        q = Vector{Float64}(D.n)
        for xyz in 1:D.n
            q[xyz] = MathProgBase.getsolution(m)[ D.n + xyz ]
        end
        eval_val::Float64 = InfDecomp_Base.condEntropy(D,q,big(0.))

        @show eval_val, MathProgBase.getobjval(m)
        if abs(eval_val + MathProgBase.getobjval(m)) > 1.e-3
            println("BIG DIFFERENCE IN SOLUTION VALUE!! :-( :-( :-( :-( :-( :-(\nProbably a BUG!!!!!!!!!!!!!!!!!!!!!")
        end

        stats = Solution_Stats(-eval_val,MathProgBase.getobjval(m),q,m)

        return stats
    else
        print("Fuck!!")
        stats = Solution_Stats(0,0,[],m)
        return stats
    end
    ; #
end #^ model_S

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

function model_D(D::InfDecomp_Base.My_Eval, solver::MathProgBase.AbstractMathProgSolver) :: Solution_Stats
    # min c^T v     #  min (0,b)^T (u^T,l^T)
    # s.t.          #  s.t.
    #               #
    # f(v) ∈ EXP*   #  (-1,f(u,l),-u) ∈ EXP*

    local n_vars::Int = D.n + D.m        # u_{xyz} ∀ xyz; l_{xy} ∀ xy; l_{xz}  ∀ xz;
                                         # u=i                  ;      l=D.n+i
    local n_cons::Int = D.n*3            # u_{*yz} + ln(-u_{xyz}) + l_{xy} + l_{xz} \le -1     ; ∀ xyz

    function make_c_A_b() :: Tuple{Vector{Float64},SparseMatrixCSC{Float64},Vector{Float64}}
        local c::Vector{Float64} = zeros(n_vars)
        local b::Vector{Float64} = zeros(n_cons)

        println("the elements are ", D.xyz)
        # fill c: (0, b^y, b^z):
        for i = 1:D.m
            c[ D.n + i ]  = D.rhs[i]
        end

        # fill b: (-1,0,0):
        for xyz = 1:D.n
            b[xyz] = -1.
        end

        #
        # fill A
        # ------

        Constr   = Array{Int,1}()
        Var   = Array{Int,1}()
        Coeff = Array{Float64,1}()

        # -1  EXP* ieqs:
        for xyz = 1:D.n
            ieqn = xyz
            u_xyz = xyz
            append!(Constr, ieqn)
            append!(Var, u_xyz)
            append!(Coeff, 0.)
        end

        # u_{*yz} + l_xy + l_xz EXP* ieqs:
        for xyz = 1:D.n
            (x,y,z) = D.xyz[xyz]
            ieqn = xyz
            # add u_{*yz}
            for rpq = 1:D.n
                (r,p,q) = D.xyz[rpq]
                if (p,q) == (y,z)
                    u_rpq = rpq
                    println("u_rpq = ($r,$p,$q) with idx= $u_rpq")
                    append!(Constr, D.n + ieqn)
                    append!(Var,    u_rpq)
                    append!(Coeff,  -1.)
                end
            end
            # add l_xy
            if D.marg_xy[x,y] > 0
                append!(Constr, D.n + ieqn)
                println("l_xy = ($x ,$y) with idx= ", D.n + D.eqidx["xy",x,y])
                append!(Var,    D.n + D.eqidx["xy",x,y])
                append!(Coeff,  -1.)
            end
            # add l_xz
            # !! D.eqidx adds xz after xy. So no need
            # to have D.n + mr_xy + D.eqidx["xz",x,z] !!
            if D.marg_xz[x,z] > 0
                append!(Constr, D.n + ieqn)
                println("l_xz = ($x ,$z) with idx= ", D.n + D.eqidx["xz",x,z])
                append!(Var,    D.n + D.eqidx["xz",x,z])
                append!(Coeff,  -1.)
            end
        end

        # -u_{xyz} EXP* ieqs:
        for xyz = 1:D.n
            ieqn = xyz
            u_xyz = xyz
            (x,y,z) = D.xyz[xyz]
            append!(Constr, D.n + D.n + ieqn)
            println("u_rpq = ($x,$y,$z) with idx= $u_xyz")
            append!(Var, u_xyz)
            append!(Coeff, 1.)
        end

        println("co", Constr)
        println("va", Var)
        println("cof", Coeff)
        A = sparse(Constr,Var,Coeff)
        println("dimension of A", size(A))
        @show A
        return c,A,b
    end #^ make_c_b_A()

    c,A,b = make_c_A_b()

    var_cones = [ ( :Free, 1:n_vars) ]
    constr_cones = []

    for xyz = 1:D.n
        rqp = [ xyz, D.n + xyz, D.n + D.n + xyz ]
        append!(constr_cones, [(:ExpDual, rqp)] )
    end
    println("cstr cones ", constr_cones)
    println("var cones ", var_cones)

    m = MathProgBase.ConicModel(solver)

    MathProgBase.loadproblem!(m, c,A,b, constr_cones, var_cones)

    println("!!")
    # CPUtic()                        # TIC

    MathProgBase.optimize!(m)

    println("!!!")
    # stats_time = CPUtoq()           # TOQ

    @show MathProgBase.status(m)
    @show stats_time = MathProgBase.getsolvetime(m)

    if MathProgBase.status(m) == :Optimal || MathProgBase.status(m) == :Suboptimal || MathProgBase.status(m) == :Error
        q = Vector{Float64}(D.n)
        mu = Vector{Float64}( Mathprogbase.getdual(m) )

        for xyz in 1:D.n
            q[xyz] = mu[ D.n + xyz ] #or mu[3*xyz - 1] depending on the order of the reduced costs
        end

        eval_val::Float64 = InfDecomp_Base.condEntropy(D,q,big(0.))

        @show eval_val, MathProgBase.getobjval(m)
        if abs(eval_val + MathProgBase.getobjval(m)) > 1.e-3
            println("BIG DIFFERENCE IN SOLUTION VALUE!! :-( :-( :-( :-( :-( :-(\nProbably a BUG!!!!!!!!!!!!!!!!!!!!!")
        end

        stats = Solution_Stats(-eval_val,MathProgBase.getobjval(m),q,stats_time,m)
        return stats
    else
        print("Fuck!!")
        stats = Solution_Stats(999,999,[],stats_time,m)
        return stats
    end
    ; #
end #^ model_D

end #^ module
; #EOF

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
using RegionTrees
using StaticArrays
using PyPlot
using Statistics
using LinearAlgebra
using SparseArrays
using VectorizedRoutines
using ProgressBars
using MAT

struct segment3d{T}
    origin::T
    dir::T
end

function segment3dpoint(p1::T, p2::T) where {T}
    dir::T = p2 - p1
    return segment3d(p1, dir)
end

struct plane3d
    origin::Union{SVector{3,Float64},Vector{Float64}}
    dir1::Union{SVector{3,Float64},Vector{Float64}}
    dir2::Union{SVector{3,Float64},Vector{Float64}}
    n::Union{SVector{3,Float64},Vector{Float64}}
    d::Float64
end

function plane3ddir(
    origin::T,
    dir1::T,
    dir2::T,
) where {T}
    A = [dir1 dir2]
    q = svd(A)
    dir1 = q.U[:,1]
    dir2 = q.U[:,2]
    n = cross(dir1,dir2)
    d = dot(origin,n)
    return plane3d(origin,dir1,dir2,n,d)
end

function plane3d(
    p1::T,
    p2::T,
    p3::T,
) where {T}
    origin = p1
    dir1 = p2-p1
    dir2 = p3-p1
    A = [dir1 dir2]
    q = svd(A)
    dir1 = q.U[:,1]
    dir2 = q.U[:,2]
    n = cross(dir1,dir2)
    d = dot(origin,n)
    return plane3d(origin,dir1,dir2,n,d)
end

@inline function intersects(
    a::HyperRectangle{N,Float64},
    b::HyperRectangle{N,Float64},
) where {N}
    maxA = a.origin + a.widths
    minA = a.origin
    maxB = b.origin + b.widths
    minB = b.origin
    return (all(maxA .>= minB) && all(maxB .> minA))

end

function inittree3d(n)
    empty::Array{Int64,1} = []
    oct = Cell(SVector(-1.0, -1, -1), SVector(2.0, 2, 2), empty)
    for _ = 1:n
        for leaf in allleaves(oct)
            split!(leaf)
        end
    end
    for node in allcells(oct)
        node.data = copy(node.data)
    end
    return oct
end

function inserttotree!(
    box::HyperRectangle{N,Float64},
    tree::Cell,
    index::Int64,
) where {N}
    if (isleaf(tree))
        if (intersects(box, tree.boundary))
            push!(tree.data, index)
        end
    elseif (intersects(box, tree.boundary))
        for i = 1:length(tree.children)
            inserttotree!(box, tree[i], index)
        end
    end
    return nothing
end

function possiblevoxels(r::plane3d, tree::Cell)::Vector{Int64}
    if (isleaf(tree) && intersects(r, tree.boundary) > 0.0)
        return tree.data

    elseif (~(tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblevoxels(r, tree[i]))
        end
        return v

    elseif ((tree.parent === nothing) && intersects(r, tree.boundary) > 0.0)
        N = length(tree.children)
        v = Vector{Int64}()
        for i = 1:N
            append!(v, possiblevoxels(r, tree[i]))
        end
        return collect(Set(v))

    else
        return []
    end

end

function sub2ind(sizes, i::Int64, j::Int64, k::Int64)
    @assert i > 0 && i <= sizes[1]
    @assert j > 0 && j <= sizes[2]
    @assert k > 0 && k <= sizes[3]
    return (k - 1) * sizes[1] * sizes[2] + sizes[1] * (j - 1) + i
end

function intersects(p::plane3d,r::segment3d)
    t = dot(p.n,r.origin-p.origin)/(-dot(r.dir,p.n))
    #println(t)
    if (0<=t<=1)
        return r.origin + t*r.dir
    else
        return nothing
    end
end

function setupvoxels(sizex::Int64,sizey::Int64,sizez::Int64,octsize::Int64)
    oct = inittree3d(octsize)
    l = Array([[1.0, 0.0, 0.0] [0,0,1] [1, 1, 0] [0, 0, 0] [0, 1, 1] [1, 0, 1] [1, 1, 1] [0, 1, 0]]')
    pm = mean(l, dims = 1)
    size = [sizex, sizey, sizez]
    cellsize = 2.0 ./ size
    width = @SVector[cellsize[1], cellsize[2], cellsize[3]]
    for i = 1:3
        l[:, i] =cellsize[i] .* (l[:, i] .- pm[i]) .- cellsize[i] / 2.0 * size[i] .+ 0.5 * cellsize[i]
    end
    voxelvector = Vector{HyperRectangle{3,Float64}}(undef, sizex * sizey * sizez)
    t = time()
    for k = 1:sizez
        for j = 1:sizey
            for i = 1:sizex
                ll = l .+ reshape(cellsize .* [i - 1, j - 1, k - 1], 1, 3)
                origin = @SVector[
                    minimum(ll[:, 1]),
                    minimum(ll[:, 2]),
                    minimum(ll[:, 3]),
                ]
                voxel = HyperRectangle(origin, width)
                index = j + (i - 1) * sizey + (k - 1) * sizey * sizex
                inserttotree!(voxel, oct, index)
                voxelvector[index] = voxel
            end
        end
    end
    println("Voxels initialized in ", time() - t, " seconds.")
    return (oct, voxelvector)
end

#https://gdbooks.gitbooks.io/3dcollisions/content/Chapter2/static_aabb_plane.html
@inline @inbounds function intersects(
    p::plane3d,
    b::HyperRectangle{3,Float64})
    e = b.widths.*0.5
    c = b.origin + e

    r = e[1]*abs(p.n[1]) + e[2]*abs(p.n[2]) +e[3]*abs(p.n[3])
    s = dot(p.n,c) -p.d

    return  -r <= s <r

end

@inbounds @inline function planeintegal(p::plane3d,h::HyperRectangle{3,Float64})::Float64
    e = Vector{segment3d{SArray{Tuple{3},Float64,1,3}}}(undef,12)
    e[1] = segment3d(SVector(h.origin),@SVector[h.widths[1],0,0])
    e[2] = segment3d(SVector(h.origin),@SVector[0,h.widths[2],0])
    e[3] = segment3d(SVector(h.origin + [h.widths[1],0,0]),@SVector[0,h.widths[2],0])
    e[4] = segment3d(SVector(h.origin + [0,h.widths[2],0]),@SVector[h.widths[1],0,0])

    e[5] = segment3d(SVector(h.origin + [0,0,h.widths[3]]),@SVector[h.widths[1],0,0])
    e[6] = segment3d(SVector(h.origin + [0,0,h.widths[3]]),@SVector[0,h.widths[2],0])
    e[7] = segment3d(SVector(h.origin + [h.widths[1],0,h.widths[3]]),@SVector[0,h.widths[2],0])
    e[8] = segment3d(SVector(h.origin + [0,h.widths[2],h.widths[3]]),@SVector[h.widths[1],0,0])

    e[9] = segment3d(SVector(h.origin ),@SVector[0,0,h.widths[3]])
    e[10] = segment3d(SVector(h.origin + [h.widths[1],0,0]),@SVector[0,0,h.widths[3]])
    e[11] = segment3d(SVector(h.origin + [0,h.widths[2],0]),@SVector[0,0,h.widths[3]])
    e[12] = segment3d(SVector(h.origin + [h.widths[1],h.widths[2],0]),@SVector[0,0,h.widths[3]])

    points = Vector{SArray{Tuple{3},Float64,1,3}}()
    for i = 1:12
        q = intersects(p,e[i])
        if(!isnothing(q))
            push!(points,q)
            #scatter3D(q[1],q[2],q[3])
        end
        #plot3D([e[i].origin[1],e[i].origin[1]+ e[i].dir[1]],[e[i].origin[2],e[i].origin[2]+ e[i].dir[2]],[e[i].origin[3],e[i].origin[3]+ e[i].dir[3]])
    end

    N = length(points)

    if(N<=2)
        return 0.0
    end
    pointsplane = Vector{SArray{Tuple{2},Float64,1,2}}(undef,N)
    for i = 1:N
        pointsplane[i] = SVector(dot(p.dir1,points[i]),dot(p.dir2,points[i]))
    end
    c = mean(pointsplane)
    sort!(pointsplane,by=x->atan((x[2]-c[2]),(x[1]-c[1])))
    val = shoelace(pointsplane)
    return val

end

@inbounds @inline function shoelace(v::T)::Float64 where {T}
    N = length(v)
    A = 0.0
    for i = 1:N-1
        A = A+ v[i][1]*v[i+1][2] - v[i+1][1]*v[i][2]
    end
    A = A + v[N][1]*v[1][2] - v[1][1]*v[N][2]
    return 0.5*abs(A)
end

function voxelplane(
    r::plane3d,
    vv::Vector{HyperRectangle{3,Float64}},
    checklist::Vector{Int64},
)
    N = length(checklist)
    indices = Vector{Int64}()
    tlens = Vector{Float64}()
    for i = 1:N
        q = planeintegal(r, vv[checklist[i]])
        if (q>0.0)
            push!(tlens, q)
            push!(indices, checklist[i])
        end

    end
    return (indices, tlens)
end

function onlyintersections(r::plane3d,vv::Vector{HyperRectangle{3,Float64}},checklist::Vector{Int64})
    N = length(checklist)
    indices = Vector{Int64}()
    for i = 1:N
        if( intersects(r, vv[checklist[i]]))
            push!(indices, checklist[i])
        end
    end

    return indices

end

function constructmatrix(
    tree::Cell,
    vv::Vector{HyperRectangle{3,Float64}},
    N::Int64,
    planevector::Vector{plane3d})

    r = 0.5 * sqrt(12)
    Na = length(planevector)
    if (N > 1)
        span = range(r, stop = -r, length = N)
    else
        span = [0.0]
    end
    Ntotal = Na * N

    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()

    Nth = Threads.nthreads()
    rows = Vector{Vector{Int64}}(undef, Nth)
    cols = Vector{Vector{Int64}}(undef, Nth)
    vals = Vector{Vector{Float64}}(undef, Nth)
    for p = 1:Nth
        rows[p] = []
        cols[p] = []
        vals[p] = []
    end

    t = time()
    if (Na > 1)
        pb = ProgressBar(1:Na)
    else
        pb = 1:1
    end
    for a in pb
        pl = planevector[a]
        Threads.@threads for j = 1:N
            #for j = 1:N
                planeindex = (a - 1) * N + j
                newp = plane3ddir(pl.origin + pl.n.*span[j],pl.dir1,pl.dir2)
                checklist = possiblevoxels(newp, tree)
                checklist = onlyintersections(newp,vv,checklist)
                indices, tlens = voxelplane(newp, vv, checklist)
                # append!(rows,planeindex*ones(Int64,length(tlens)))
                # append!(cols,indices)
                # append!(vals,tlens)

                append!(rows[Threads.threadid()],planeindex * ones(Int64, length(tlens)))
                append!(cols[Threads.threadid()], indices)
                append!(vals[Threads.threadid()], tlens)

            end
        end

    rows = vcat(rows...)
    cols = vcat(cols...)
    vals = vcat(vals...)
    M = sparse(rows, cols, vals, Ntotal, length(vv))
    println("Matrix constructed in ", time() - t, " seconds.")
    return M
end

function test()
    Nx = 100; Ny = 100; Nz =100; Os = 6; Ns = 100
    (oct,voxelvector)=setupvoxels(Nx,Ny,Nz,Os)
    p1=plane3ddir([0,0,0.0],[1,1.0,1.0],[0.0,1,1.0])
    p2=plane3ddir([0,0,0.0],[1,0.0,0.0],[0.0,0,1.0])
    planevector=[p1,p2]

    # Plotting a part of the first plane.
    c1 = p1.origin + p1.dir1*3 + p1.dir2*3
    c2 = p1.origin + p1.dir1*3 - p1.dir2*3
    c3 = p1.origin - p1.dir1*3 + p1.dir2*3
    c4 = p1.origin - p1.dir1*3 - p1.dir2*3
    xx = reshape([c1[1],c2[1],c3[1],c4[1]],2,2)
    yy = reshape([c1[2],c2[2],c3[2],c4[2]],2,2)
    z = reshape([c1[3],c2[3],c3[3],c4[3]],2,2)
    plot_surface(xx,yy,z)

    M=constructmatrix(oct,voxelvector,Ns,planevector)

    l = Array([[1.0, 0.0, 0.0] [0,0,1] [1, 1, 0] [0, 0, 0] [0, 1, 1] [1, 0, 1] [1, 1, 1] [0, 1, 0]]')
    pm = mean(l, dims = 1)
    size = [Nx,Ny,Nz]
    cellsize = 2.0 ./ size
    qx=Vector(range(-1 + cellsize[1]/2.0,stop=1- cellsize[1]/2.0,length=Nx))
    qy=Vector(range(-1+ cellsize[2]/2.0,stop=1- cellsize[2]/2.0,length=Ny))
    qz=Vector(range(-1+ cellsize[3]/2.0,stop=1- cellsize[3]/2.0,length=Nz))

    mgrid = Matlab.meshgrid(qx,qy,qz)
    x = mgrid[1]; y = mgrid[2]; z = mgrid[3]

    object = @. 1.0.*(x^2 + y^2 + z^2 < 1) #Testing Radon transform of a unit ball.
    objf = object[:]
    sino = M*objf
    sino1 = sino[1:Ns]
    sino2 = sino[Ns+1:end]
    fig = plt.figure()
    plot(sino1)
    plot(sino2)

    return  nothing
end

test()

max_stars_repo_name: suurj/radon

max_issues_repo_name: suurj/radon

max_forks_repo_name: suurj/radon


text:
using AAindex
using Test

const testdir = dirname(@__FILE__)

tests = [
    "index",
    "parse",
    "search"
]

include("testdata/testvars.jl")

@testset "AAindex" begin
    for t in tests
        @info "Testing $t"
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ClusteringTask(x, y)

A clustering task consists of grouping training examples based on
features `x ∈ Rⁿ`. The resulting clusters are saved in spatial
variable `y`.
"""
struct ClusteringTask{N} <: AbstractLearningTask
  features::NTuple{N,Symbol}
  output::Symbol
end

ClusteringTask(x::Symbol, out::Symbol) = ClusteringTask{1}((x,), out)

ClusteringTask(x::AbstractVector{Symbol}, out::Symbol) =
  ClusteringTask{length(x)}(Tuple(x), out)

"""
    features(task)

Return the features of the clustering `task`.
"""
features(task::ClusteringTask) = task.features

outputvars(task::ClusteringTask) = (task.output,)

# ------------
# IO methods
# ------------
function Base.show(io::IO, task::ClusteringTask)
  x = task.features
  y = task.output
  lhs = length(x) > 1 ? "("*join(x, ", ")*")" : "$(x[1])"
  print(io, "Clustering $lhs → $y")
end


text:
using ..ChemistryFeaturization.AbstractType: AbstractCodec

"""
    SimpleCodec(encode_f, decode_f)

SimpleCodec is a most simple, atomic (as in indivisible), and straightforward Codec.

A SimpleCodec is characterized by only the encoding and decoding functions it has as fields.
This means that the encoding-decoding scheme is entirely dependent only upon these functions, and the fields of the FeatureDescriptor to which this Codec is attached. No external parameters (such as number of bins, etc.) are required
by the encoding-decoding implementations.

This Codec is very useful when the encoding-decoding scheme that you want to define for a `FeatureDescriptor` is very simple.
"""
struct SimpleCodec <: AbstractCodec
    encode_f::Function
    decode_f::Function
end


text:
module MeteomaticsAPI

include("QueryFunctions.jl")

export querytimeseries, querygrid, querytimeranges, queryinitdatetime,
       data_url, meta_url

end

max_stars_repo_name: meteomatics/julia-connector-api

max_issues_repo_name: meteomatics/julia-connector-api

max_forks_repo_name: meteomatics/julia-connector-api


text:
using DataFrames
using NetCDF
using IcanProj
using Statistics
using CSV
using ProgressMeter

function write_netcdfs(path, savefolder, res, frac_landuse, lai_landuse, hcan_landuse)

    filename = joinpath(path, savefolder, "params_$(res).nc")

    isfile(filename) && rm(filename, force = true)

    dim_space = collect(1:size(frac_landuse, 1))
    dim_classes = collect(1:size(frac_landuse, 2))

    nccreate(filename, "frac_landuse", "dim_space", dim_space, "dim_classes", dim_classes)
    nccreate(filename, "lai_landuse", "dim_space", dim_space, "dim_classes", dim_classes)
    nccreate(filename, "hcan_landuse", "dim_space", dim_space, "dim_classes", dim_classes)

    ncwrite(frac_landuse, filename, "frac_landuse")
    ncwrite(lai_landuse, filename, "lai_landuse")
    ncwrite(hcan_landuse, filename, "hcan_landuse")

    ncclose(filename)

end

function aggregate_params(path, res, name)

    @info res

    # Metadata from table

    df_landuse = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "landuse.csv"), rows_for_type_detect = typemax(Int)) |> DataFrame

    # Metadata from input netcdfs

    file_in = joinpath(path, "forcings_st", "tair_$(res).nc")

    # Grid cell indicies in netcdf files

    ind_nc = Int.(ncread(file_in, "id"))

    # Leaf area index and canopy height

    lai = [1.4, 4.3, 6.7, 9.1, 0.9, 2.4, 2.3, 4.4, 0.1, 0.2, 0.3, 0.4, 0]
    hcan = [7.5, 12.3, 16.8, 22, 7.5, 11.6, 17, 17.2, 4.9, 8.4, 12.2, 18.3, 0]

    # Output arrays for forest runs

    landuse_classes = 3

    frac_landuse = Array{Float64}(undef, length(ind_nc), landuse_classes)
    lai_landuse = Array{Float64}(undef, length(ind_nc), landuse_classes)
    hcan_landuse = Array{Float64}(undef, length(ind_nc), landuse_classes)

    @showprogress 1 "Running for ..." for i in eachindex(ind_nc)

        ind_cell = ind_nc[i]

        irows = findall(df_landuse[Symbol(name)] .== ind_cell)

        # No data available - define as open areas

        if isempty(irows)

            frac_landuse[i, 1] .= 1

            frac_landuse[i, 2:end] .= 0

            lai_landuse[i, :] .= 0

            hcan_landuse[i, :] .= 0

            continue

        end

        # Compute total fraction

        tmp = colwise(mean, df_landuse[irows, 9:end])

        isorted = sortperm(tmp, rev=true)

        ifinal = isorted[1:landuse_classes]

        frac_total = sum(tmp[ifinal])

        # Total fraction is too small - define as open areas

        if frac_total < 0.1

            frac_landuse[i, 1] .= 1

            frac_landuse[i, 2:end] .= 0

            lai_landuse[i, :] .= 0

            hcan_landuse[i, :] .= 0

            continue

        end

        # Data is valid - use landuse classes

        frac_landuse[i, :] = tmp[ifinal] ./ frac_total

        lai_landuse[i, :] = lai[ifinal]

        hcan_landuse[i, :] = hcan[ifinal]

    end

    # Save to netcdfs for forest runs

    write_netcdfs(path, "params_forest", res, frac_landuse, lai_landuse, hcan_landuse)

    # Deal with simulations omitting forest processes

    frac_landuse = ones(length(ind_nc), 1)
    lai_landuse = zeros(length(ind_nc), 1)
    hcan_landuse = zeros(length(ind_nc), 1)

    # Save to netcdfs for forest runs

    write_netcdfs(path, "params_open", res, frac_landuse, lai_landuse, hcan_landuse)

end

# Settings

path = "/data04/jmg/fsm_simulations/netcdf"

res_name= Dict("1km" => "ind_senorge",
               "5km" => "ind_5km",
               "10km" => "ind_10km",
               "25km" => "ind_25km",
               "50km" => "ind_50km")

for (res, name) in res_name

    aggregate_params(path, res, name)

end


text:
function solveKnapNaive(data::KnapData)
    profits = zeros(Int64, length(data.items) + 1, data.capacity + 1)

    for i in 2:length(data.items) + 1
        w = data.items[i - 1].weight
        p = data.items[i - 1].profit

        profits[i, :] = profits[i - 1, :]

        for j in w + 1:data.capacity + 1
            profits[i, j] = max(profits[i, j], profits[i - 1, j - w] + p)
        end
    end

    best = 1
    for j in 2:data.capacity + 1
        if profits[end, j] > profits[end, best]
            best = j
        end
    end

    best_val = profits[end, best]
    solution = Int64[]
    for i in length(data.items) + 1:-1:2
        if profits[i, best] != profits[i - 1, best]
            push!(solution, i - 1)
            best -= data.items[i - 1].weight
        end
    end
    return best_val, reverse(solution)
end


text:
@testset "518.coin-change-2.jl" begin
    @test change(5, [1, 2, 5]) == 4
    @test change(3, [2]) == 0
    @test change(10, [10]) == 1
end


text:
__precompile__(false)

module Cassette

using Core: CodeInfo, SlotNumber, NewvarNode, LabelNode, GotoNode, SSAValue, arrayref, arrayset

using Logging

struct Unused end

abstract type Context end
abstract type Tag{C,T} end

@inline tagtype(x) = tagtype(Nothing, x)

@generated function tagtype(::Type{C}, ::T) where {C,T}
    return quote
        $(Expr(:meta, :inline))
        Tag{C,$(objectid(T))}
    end
end

const unused = Unused()
const MAX_ARGS = 20

include("utilities.jl")
include("metadata.jl")
include("reflection.jl")
include("execution.jl")
include("macros.jl")
include("workarounds.jl")

function __init__()
    # FIXME: Base should provide a mechanism for this (eg. Julia/julia#26265)
    DEBUG = parse(Bool, get(ENV, "DEBUG", "false"))
    if DEBUG
        global_logger(ConsoleLogger(global_logger().stream, Logging.Debug))
    end
end

end # module


text:
module SimpleTrees

# using AbstractTrees
using FMMTrees

struct TreeNode{T}
    num_children::Int
    data::T
end

struct SimpleTree{V <: (AbstractVector{N} where {N<:TreeNode})}
    nodes::V
end

struct ChildView{T} tree::T end
FMMTrees.children(tree::SimpleTree, node=root(tree)) = collect(ChildView(node))
FMMTrees.root(tree::SimpleTree) = tree
FMMTrees.data(tree::SimpleTree, node=root(tree)) = first(node.nodes).data
FMMTrees.haschildren(tree::SimpleTree, node=root(tree)) = (first(node.nodes).num_children > 0)

Base.IteratorSize(cv::ChildView) = Base.SizeUnknown()

function Base.iterate(itr::ChildView, state=(0,2))
    state[1] == first(itr.tree.nodes).num_children && return nothing
    child = itr.tree.nodes[state[2]]
    newstate = (state[1] + child.num_children + 1, state[2] + child.num_children + 1)
    return SimpleTree(view(itr.tree.nodes, state[2]:lastindex(itr.tree.nodes))), newstate
end

Base.getindex(tree::SimpleTree, i::Int) = tree.nodes[i]

# AbstractTrees.printnode(io::IO, tree::SimpleTree) = show(io, data(tree))

end # module SimpleTrees


text:
function pull_request_build(::NewVersion,
                            pr::GitHub.PullRequest,
                            current_pr_head_commit_sha::String,
                            registry::GitHub.Repo;
                            auth::GitHub.Authorization,
                            authorized_authors::Vector{String},
                            registry_head::String,
                            registry_master::String,
                            suggest_onepointzero::Bool,
                            whoami::String)
    # first check if authorized author - if not, then quit
    # then check rules 1-3. if fail, post comment.
    # then check rules 4-5. if fail, post comment.
    # if everything passed, merge the pull request now
    # 1. Sequential version number - if the last version was 1.2.3 then the next can be 1.2.4, 1.3.0 or 2.0.0
    # 2. Compat for all dependencies - all [deps] should also have [compat] entries (and Julia itself) - [compat] entries should have upper bounds
    # 3. If it is a patch release, then it does not narrow the Julia compat range
    # 4. Version can be installed - given the proposed changes to the registry, can we resolve and install the new version of the package?
    # 5. Version can be loaded - once it's been installed (and built?), can we load the code?
    pkg, version = parse_pull_request_title(NewVersion(), pr)
    @info("This is a new version pull request", pkg, version)
    pr_author_login = author_login(pr)
    if is_open(pr)
        if pr_author_login in authorized_authors
            my_retry(() -> delete_all_of_my_reviews!(registry, pr; auth = auth, whoami = whoami))
            my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "pending", "context" => "automerge/new-version")))
            g0, m0 = pr_only_changes_allowed_files(NewVersion(), registry, pr, pkg; auth = auth)
            newv_g1, newv_m1, release_type = meets_sequential_version_number(pkg,
                                                                             version;
                                                                             registry_head = registry_head,
                                                                             registry_master = registry_master)
            newv_g2, newv_m2 = meets_compat_for_all_deps(registry_head,
                                                         pkg,
                                                         version)
            if release_type == :patch
                newv_g3, newv_m3 = meets_patch_release_does_not_narrow_julia_compat(pkg,
                                                                                version;
                                                                                registry_head = registry_head,
                                                                                registry_master = registry_master)
            else
                newv_g3 = true
                newv_m3 = ""
            end
            newv_g1through3 = [g0, newv_g1, newv_g2, newv_g3]
            @info("Only modifies the files that it's allowed to modify", meets_this_guideline = g0, message = m0)
            @info("Sequential version number", meets_this_guideline = newv_g1, message = newv_m1)
            @info("Compat (with upper bound) for all dependencies", meets_this_guideline = newv_g2, message = newv_m2)
            @info("If it is a patch release, then it does not narrow the Julia compat range", meets_this_guideline = newv_g3, message = newv_m3)
            if all(newv_g1through3)
                newv_g4and5, newv_m4and5 = meets_version_can_be_loaded(registry_head,
                                                                       pkg,
                                                                       version)
                @info("Version can be installed and loaded", meets_this_guideline = newv_g4and5, message = newv_m4and5)
                if newv_g4and5
                    newv_commenttextpass = comment_text_pass(NewVersion(),
                                                             suggest_onepointzero,
                                                             version)
                    my_retry(() -> delete_all_of_my_reviews!(registry, pr; auth = auth, whoami = whoami))
                    my_retry(() -> approve!(registry, pr, current_pr_head_commit_sha; auth = auth, body = newv_commenttextpass, whoami = whoami))
                    my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "success", "context" => "automerge/new-version", "description" => "$(current_pr_head_commit_sha)")))
                    return nothing
                else
                    newv_commenttext4and5 = comment_text_fail(NewVersion(),
                                                              [newv_m4and5],
                                                              suggest_onepointzero,
                                                              version)
                    my_retry(() -> post_comment!(registry, pr, newv_commenttext4and5; auth = auth))
                    my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "failure", "context" => "automerge/new-version")))
                    return nothing
                end
            else
                newv_allmessages1through3 = [m0, newv_m1, newv_m2, newv_m3]
                newv_failingmessages1through3 = newv_allmessages1through3[.!newv_g1through3]
                newv_commenttext1through3 = comment_text_fail(NewVersion(),
                                                              newv_failingmessages1through3,
                                                              suggest_onepointzero,
                                                              version)
                my_retry(() -> post_comment!(registry, pr, newv_commenttext1through3; auth = auth))
                my_retry(() -> GitHub.create_status(registry, current_pr_head_commit_sha; auth=auth, params=Dict("state" => "failure", "context" => "automerge/new-version")))
                return nothing
            end
        else
            @info("Author $(pr_author_login) is not authorized to automerge. Exiting...")
            return nothing
        end
    else
        @info("The pull request is not open. Exiting...")
        return nothing
    end
    return nothing
end


text:
function constructMIRK_IV(S::BVPSystem{T,U}) where {T,U}
    c = [0, 1, 1//2, 3//4]
    v = [0, 1, 1//2, 27//32]
    b = [1//6, 1//6, 2//3, 0]
    x = [0      0       0 0
         0      0       0 0
         1//8   -1//8   0 0
         3//64  -9//64  0 0]
    MIRKTableau(T.(c),T.(v),T.(b),T.(x))
end

constructMIRK(S::BVPSystem) = MIRK_dispatcher(S, Val{S.order})
MIRK_dispatcher(S::BVPSystem, ::Type{Val{4}}) = constructMIRK_IV(S)


text:
#############################################################
#
# Timelines Functions
#
#############################################################

function get_mentions_timeline(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/statuses/mentions_timeline.json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function get_user_timeline(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/statuses/user_timeline.json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function get_home_timeline(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/statuses/home_timeline.json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function get_retweets_of_me(; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/statuses/retweets_of_me.json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using SuiteSparse.SPQR
using SuiteSparse.CHOLMOD
using LinearAlgebra: rmul!, lmul!, Adjoint, Transpose

@testset "Sparse QR" begin
m, n = 100, 10
nn = 100

@test size(qrfact(sprandn(m, n, 0.1)).Q) == (m, m)

@testset "element type of A: $eltyA" for eltyA in (Float64, Complex{Float64})
    if eltyA <: Real
        A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], randn(nn), m, n)
    else
        A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], complex.(randn(nn), randn(nn)), m, n)
    end

    F = qrfact(A)
    @test size(F) == (m,n)
    @test size(F, 1) == m
    @test size(F, 2) == n
    @test size(F, 3) == 1
    @test_throws ArgumentError size(F, 0)

    @testset "getindex" begin
        @test istriu(F.R)
        @test isperm(F.pcol)
        @test isperm(F.prow)
        @test_throws ErrorException F.T
    end

    @testset "apply Q" begin
        Q = F.Q
        Imm = Matrix{Float64}(I, m, m)
        @test Q' * (Q*Imm) ≈ Imm
        @test (Imm*Q) * Q' ≈ Imm

        # test that Q'Pl*A*Pr = R
        R0 = Q'*Array(A[F.prow, F.pcol])
        @test R0[1:n, :] ≈ F.R
        @test norm(R0[n + 1:end, :], 1) < 1e-12

        offsizeA = Matrix{Float64}(I, m+1, m+1)
        @test_throws DimensionMismatch lmul!(Q, offsizeA)
        @test_throws DimensionMismatch lmul!(adjoint(Q), offsizeA)
        @test_throws DimensionMismatch rmul!(offsizeA, Q)
        @test_throws DimensionMismatch rmul!(offsizeA, adjoint(Q))
    end

    @testset "element type of B: $eltyB" for eltyB in (Int, Float64, Complex{Float64})
        if eltyB == Int
            B = rand(1:10, m, 2)
        elseif eltyB <: Real
            B = randn(m, 2)
        else
            B = complex.(randn(m, 2), randn(m, 2))
        end

        @inferred A\B
        @test A\B[:,1] ≈ Array(A)\B[:,1]
        @test A\B ≈ Array(A)\B
        @test_throws DimensionMismatch A\B[1:m-1,:]
        C, x = A[1:9, :], fill(eltyB(1), 9)
        @test C*(C\x) ≈ x # Underdetermined system
    end

    # Make sure that conversion to Sparse doesn't use SuiteSparse's symmetric flag
    @test qrfact(SparseMatrixCSC{eltyA}(I, 5, 5)) \ fill(eltyA(1), 5) == fill(1, 5)
end

@testset "basic solution of rank deficient ls" begin
    A = sprandn(m, 5, 0.9)*sprandn(5, n, 0.9)
    b = randn(m)
    xs = A\b
    xd = Array(A)\b

    # check that basic solution has more zeros
    @test count(!iszero, xs) < count(!iszero, xd)
    @test A*xs ≈ A*xd
end

end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
"""
    run_sims(simtype::Type{<:SpatQSimSpec},
             repl_range::AbstractRange{<:Integer};
             prep_file = "prep.h5",
             checkpoint = true)

Run simulations of a fishery for replicates in `repl_range`, using `simtype`
spec. If `checkpoint` is `true` any replicates where the simulation result files
are already present will be skipped.
"""
function run_sims(simtype::Type{<:SpatQSimSpec},
                  repl_range::AbstractRange{<:Integer};
                  prep_file = "prep.h5",
                  checkpoint = true,
                  csv = false,
                  feather = true)
    for rlz in repl_range
        run_sims(simtype,
                 rlz;
                 prep_file = prep_file,
                 checkpoint = checkpoint,
                 csv = csv,
                 feather = feather)
    end
    nothing
end

function run_sims(simtype::Type{<:SpatQSimSpec},
                  rlz::Integer;
                  prep_file = "prep.h5",
                  checkpoint = true,
                  csv = false,
                  feather = true)
    study_range = sim_values(simtype)
    for val in study_range
        spec = simtype(rlz, val, prep_file)
        files = file_paths(spec)

        # If checkpointing, skip the rest of the inner for loop body if the
        # results files are present.
        if (checkpoint && all(isfile.(files)))
            continue
        end

        setup = SpatQSimSetup(spec)
        result = simulate(setup)
        save(result; csv = csv, feather = feather)
    end
    nothing
end


text:
module PointPatternStatistics
using LinearAlgebra: norm
using StaticArrays
using RecipesBase
using Requires

struct PointPattern{PT, WT}
    data::Vector{PT}
    window::WT
end
Base.length(pp::PointPattern) = npoints(pp)
npoints(pp) = length(pp.data)
window(pp) = pp.window

Base.iterate(pp::PointPattern) = iterate(pp.data)
Base.iterate(pp::PointPattern, state) = iterate(pp.data, state)
Base.eltype(pp::PointPattern) = eltype(pp.data)

getx(x::NamedTuple) = x.x
gety(x::NamedTuple) = x.y
getx(x) = x[1]
gety(x) = x[2]

include("window.jl")
include("weighted_distance_histogram.jl")
include("pair_correlation_function.jl")
include("K.jl")
include("kaplanmeier.jl")
#include("distance_transform.jl")
#include("bdist.jl")
#include("Fdt.jl")
include("Fnn.jl")
#include("F.jl")
include("G.jl")
include("envelope_erl.jl")

include("simple_sequential_inhibition.jl")

include("plotrecipes.jl")
Fest(pp, r, N=100) = Fkmnn(pp.data, window(pp), r, N)
Gest = Gkm
Kest = Ktrans
globalenvelope = erlenvelope
export PointPattern
export Kest, pcf, Fest, Gest, globalenvelope, Lest, L12
export SimpleSequentialInhibition
export inside, window

function __init__()
    @require RCall="6f49c342-dc21-5d91-9882-a32aef131414" include("convertR.jl")
end

end # module


text:
# # 2D forced-dissipative turbulence
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/twodnavierstokes_stochasticforcing.ipynb).
#
# A simulation of forced-dissipative two-dimensional turbulence. We solve the
# two-dimensional vorticity equation with stochastic excitation and dissipation in
# the form of linear drag and hyperviscosity.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Random, Printf, Plots"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using GeophysicalFlows, Random, Printf, Plots
parsevalsum = FourierFlows.parsevalsum

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical, domain, and simulation parameters
#
# First, we pick some numerical and physical parameters for our model.

 n, L  = 256, 2π             # grid resolution and domain length
 ν, nν = 2e-7, 2             # hyperviscosity coefficient and hyperviscosity order
 μ, nμ = 1e-1, 0             # linear drag coefficient
    dt = 0.005               # timestep
nsteps = 4000                # total number of steps
 nsubs = 20                  # number of steps between each plot
nothing # hide

# ## Forcing

# We force the vorticity equation with stochastic excitation that is delta-correlated in time
# and while spatially homogeneously and isotropically correlated. The forcing has a spectrum
# with power in a ring in wavenumber space of radius ``k_f`` (`forcing_wavenumber`) and width
# ``δ_f`` (`forcing_bandwidth`), and it injects energy per unit area and per unit time
# equal to ``\varepsilon``. That is, the forcing covariance spectrum is proportional to
# ``\exp{[-(|\bm{k}| - k_f)^2 / (2 δ_f^2)]}``.

forcing_wavenumber = 14.0 * 2π/L  # the forcing wavenumber, `k_f`, for a spectrum that is a ring in wavenumber space
forcing_bandwidth  = 1.5  * 2π/L  # the width of the forcing spectrum, `δ_f`
ε = 0.1                           # energy input rate by the forcing

grid = TwoDGrid(dev, n, L)

K = @. sqrt(grid.Krsq)             # a 2D array with the total wavenumber

forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
@CUDA.allowscalar forcing_spectrum[grid.Krsq .== 0] .= 0 # ensure forcing has zero domain-average

ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
@. forcing_spectrum *= ε/ε0        # normalize forcing to inject energy at rate ε
nothing # hide

# We reset of the random number generator for reproducibility
if dev==CPU(); Random.seed!(1234); else; CUDA.seed!(1234); end
nothing # hide

# Next we construct function `calcF!` that computes a forcing realization every timestep.
# First we make sure that if `dev=GPU()`, then `CUDA.rand()` function is called for random
# numbers uniformly distributed between 0 and 1.
random_uniform = dev==CPU() ? rand : CUDA.rand

function calcF!(Fh, sol, t, clock, vars, params, grid)
  Fh .= sqrt.(forcing_spectrum) .* exp.(2π * im * random_uniform(eltype(grid), size(sol))) ./ sqrt(clock.dt)

  return nothing
end
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. The
# `stepper` keyword defines the time-stepper to be used.
prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4",
                                calcF=calcF!, stochastic=true)
nothing # hide

# Define some shortcuts for convenience.
sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

x, y = grid.x, grid.y
nothing # hide

# First let's see how a forcing realization looks like. Function `calcF!()` computes
# the forcing in Fourier space and saves it into variable `vars.Fh`, so we first need to
# go back to physical space.

# Note that when plotting, we decorate the variable to be plotted with `Array()` to make sure
# it is brought back on the CPU when the variable lives on the GPU.
calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid)

heatmap(x, y, Array(irfft(vars.Fh, grid.nx)'),
     aspectratio = 1,
               c = :balance,
            clim = (-200, 200),
           xlims = (-L/2, L/2),
           ylims = (-L/2, L/2),
          xticks = -3:3,
          yticks = -3:3,
          xlabel = "x",
          ylabel = "y",
           title = "a forcing realization",
      framestyle = :box)

# ## Setting initial conditions

# Our initial condition is a fluid at rest.
TwoDNavierStokes.set_ζ!(prob, ArrayType(dev)(zeros(grid.nx, grid.ny)))

# ## Diagnostics

# Create Diagnostics; the diagnostics are aimed to probe the energy budget.
E  = Diagnostic(TwoDNavierStokes.energy,                prob, nsteps=nsteps) # energy
Z  = Diagnostic(TwoDNavierStokes.enstrophy,             prob, nsteps=nsteps) # enstrophy
diags = [E, Z] # a list of Diagnostics passed to `stepforward!` will  be updated every timestep.
nothing # hide

# ## Visualizing the simulation

# We initialize a plot with the vorticity field and the time-series of
# energy and enstrophy diagnostics. To plot energy and enstrophy on the same
# axes we scale enstrophy with ``k_f^2``.

p1 = heatmap(x, y, Array(vars.ζ'),
         aspectratio = 1,
                   c = :balance,
                clim = (-40, 40),
               xlims = (-L/2, L/2),
               ylims = (-L/2, L/2),
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "vorticity, t=" * @sprintf("%.2f", clock.t),
          framestyle = :box)

p2 = plot(2, # this means "a plot with two series"
               label = ["energy E(t)" "enstrophy Z(t) / k_f²"],
              legend = :right,
           linewidth = 2,
               alpha = 0.7,
              xlabel = "μ t",
               xlims = (0, 1.1 * μ * nsteps * dt),
               ylims = (0, 0.55))

l = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout = l, size = (900, 420))

# ## Time-stepping the `Problem` forward

# Finally, we time-step the `Problem` forward in time.

startwalltime = time()

anim = @animate for j = 0:round(Int, nsteps / nsubs)
  if j % (1000/nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Z: %.4f, walltime: %.2f min",
          clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)
    println(log)
  end

  p[1][1][:z] = Array(vars.ζ)
  p[1][:title] = "vorticity, μt = " * @sprintf("%.2f", μ * clock.t)
  push!(p[2][1], μ * E.t[E.i], E.data[E.i])
  push!(p[2][2], μ * Z.t[Z.i], Z.data[Z.i] / forcing_wavenumber^2)

  stepforward!(prob, diags, nsubs)
  TwoDNavierStokes.updatevars!(prob)
end

mp4(anim, "twodturb_forced.mp4", fps=18)


text:
function copy(obj::UnivariateStatistic)
    return jcall(obj, "copy", UnivariateStatistic, ())
end

function evaluate(obj::UnivariateStatistic, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using DataFrames
using SISLES
using RLESUtils, MathUtils

const COLTYPES =
        [
            Bool,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Float64,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Float64,
            Bool ,
            Bool,
            Bool,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Bool,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Float64,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Bool,
            Float64,
            Bool,
            Bool,
            Bool,
            Float64,
            Float64,
            Float64,
            Float64,
            Float64,
            Bool,
            Float64,
            Float64
        ]

const COLNAMES =
        [
            :RA_1,
            :vert_rate_1,
            :alt_diff_1,
            :psi_1,
            :intr_sr_1,
            :intr_chi_1,
            :intr_vrc0_1,
            :intr_vrc1_1,
            :intr_vrc2_1,
            :cc0_1,
            :cc1_1,
            :cc2_1,
            :vc0_1,
            :vc1_1,
            :vc2_1,
            :ua0_1,
            :ua1_1,
            :ua2_1,
            :da0_1,
            :da1_1,
            :da2_1,
            :target_rate_1,
            :crossing_1,
            :alarm_1,
            :alert_1,
            :intr_out_vrc0_1,
            :intr_out_vrc1_1,
            :intr_out_vrc2_1,
            :intr_out_tds_1,
            :response_none_1,
            :response_stay_1,
            :response_follow_1,
            :response_timer_1,
            :response_h_d_1,
            :response_psi_d_1,
            :v_1,
            :h_1,
            :RA_2,
            :vert_rate_2,
            :alt_diff_2,
            :psi_2,
            :intr_sr_2,
            :intr_chi_2,
            :intr_vrc0_2,
            :intr_vrc1_2,
            :intr_vrc2_2,
            :cc0_2,
            :cc1_2,
            :cc2_2,
            :vc0_2,
            :vc1_2,
            :vc2_2,
            :ua0_2,
            :ua1_2,
            :ua2_2,
            :da0_2,
            :da1_2,
            :da2_2,
            :target_rate_2,
            :crossing_2,
            :alarm_2,
            :alert_2,
            :intr_out_vrc0_2,
            :intr_out_vrc1_2,
            :intr_out_vrc2_2,
            :intr_out_tds_2,
            :response_none_2,
            :response_stay_2,
            :response_follow_2,
            :response_timer_2,
            :response_h_d_2,
            :response_psi_d_2,
            :v_2,
            :h_2,
            :converging,
            :abs_alt_diff,
            :horizontal_range
        ]

extract_state(sim::ACASX_GM) = extract_state(Vector, sim::ACASX_GM)

function extract_state(::Type{DataFrame}, sim::ACASX_GM)
    D = DataFrame(COLTYPES, COLNAMES, 0)
    push!(D, extract_state(Vector, sim))
    D
end

function extract_state(::Type{Vector}, sim::ACASX_GM)
    ra_active_1 = (sim.cas[1].output.dh_min > -9999.0 || sim.cas[1].output.dh_max < 9999.0)::Bool
    ra_active_2 = (sim.cas[2].output.dh_min > -9999.0 || sim.cas[2].output.dh_max < 9999.0)::Bool
    converging = is_converging(sim.cas[1].input.ownInput.psi,
        sim.cas[1].input.intruders[1].chi,
        sim.cas[2].input.ownInput.psi,
        sim.cas[2].input.intruders[1].chi)
    abs_alt_diff = abs(sim.dm[2].state.h - sim.dm[1].state.h)
    horizontal_range = sqrt(max(0.0, sim.cas[1].input.intruders[1].sr^2 - abs_alt_diff^2))
    v = [
           ra_active_1,
           sim.cas[1].input.ownInput.dz,
           sim.cas[1].input.ownInput.z,
           sim.cas[1].input.ownInput.psi,
           sim.cas[1].input.intruders[1].sr,
           sim.cas[1].input.intruders[1].chi,
           sim.cas[1].input.intruders[1].vrc == 0,
           sim.cas[1].input.intruders[1].vrc == 1,
           sim.cas[1].input.intruders[1].vrc == 2,
           bin(sim.cas[1].output.cc, 3)[1] == '1',
           bin(sim.cas[1].output.cc, 3)[2] == '1',
           bin(sim.cas[1].output.cc, 3)[3] == '1',
           bin(sim.cas[1].output.vc, 3)[1] == '1',
           bin(sim.cas[1].output.vc, 3)[2] == '1',
           bin(sim.cas[1].output.vc, 3)[3] == '1',
           bin(sim.cas[1].output.ua, 3)[1] == '1',
           bin(sim.cas[1].output.ua, 3)[2] == '1',
           bin(sim.cas[1].output.ua, 3)[3] == '1',
           bin(sim.cas[1].output.da, 3)[1] == '1',
           bin(sim.cas[1].output.da, 3)[2] == '1',
           bin(sim.cas[1].output.da, 3)[3] == '1',
           sim.cas[1].output.target_rate,
           sim.cas[1].output.crossing,
           sim.cas[1].output.alarm,
           sim.cas[1].output.alert,
           sim.cas[1].output.intruders[1].vrc == 0,
           sim.cas[1].output.intruders[1].vrc == 1,
           sim.cas[1].output.intruders[1].vrc == 2,
           sim.cas[1].output.intruders[1].tds,
           sim.pr[1].state == "none",
           sim.pr[1].state == "stay",
           sim.pr[1].state == "follow",
           sim.pr[1].timer,
           sim.pr[1].output.h_d,
           sim.pr[1].output.psi_d,
           sim.dm[1].state.v,
           sim.dm[1].state.h,
           ra_active_2,
           sim.cas[2].input.ownInput.dz,
           sim.cas[2].input.ownInput.z,
           sim.cas[2].input.ownInput.psi,
           sim.cas[2].input.intruders[1].sr,
           sim.cas[2].input.intruders[1].chi,
           sim.cas[2].input.intruders[1].vrc == 0,
           sim.cas[2].input.intruders[1].vrc == 1,
           sim.cas[2].input.intruders[1].vrc == 2,
           bin(sim.cas[2].output.cc, 3)[1] == '1',
           bin(sim.cas[2].output.cc, 3)[2] == '1',
           bin(sim.cas[2].output.cc, 3)[3] == '1',
           bin(sim.cas[2].output.vc, 3)[1] == '1',
           bin(sim.cas[2].output.vc, 3)[2] == '1',
           bin(sim.cas[2].output.vc, 3)[3] == '1',
           bin(sim.cas[2].output.ua, 3)[1] == '1',
           bin(sim.cas[2].output.ua, 3)[2] == '1',
           bin(sim.cas[2].output.ua, 3)[3] == '1',
           bin(sim.cas[2].output.da, 3)[1] == '1',
           bin(sim.cas[2].output.da, 3)[2] == '1',
           bin(sim.cas[2].output.da, 3)[3] == '1',
           sim.cas[2].output.target_rate,
           sim.cas[2].output.crossing,
           sim.cas[2].output.alarm,
           sim.cas[2].output.alert,
           sim.cas[2].output.intruders[1].vrc == 0,
           sim.cas[2].output.intruders[1].vrc == 1,
           sim.cas[2].output.intruders[1].vrc == 2,
           sim.cas[2].output.intruders[1].tds,
           sim.pr[2].state == "none",
           sim.pr[2].state == "stay",
           sim.pr[2].state == "follow",
           sim.pr[2].timer,
           sim.pr[2].output.h_d,
           sim.pr[2].output.psi_d,
           sim.dm[2].state.v,
           sim.dm[2].state.h,
           converging,
           abs_alt_diff,
           horizontal_range
        ]
    v
end

function is_converging(psi1::Float64, chi1::Float64, psi2::Float64, chi2::Float64)
    #println("psi1=$psi1, chi1=$chi1, psi2=$psi2, chi2=$chi2")
    if abs(chi1) > pi/2 && abs(chi2) > pi/2 #flying away from each other
        return false
    end
    z1 = to_plusminus_pi(psi2 - psi1)
    z2 = to_plusminus_pi(psi1 - psi2)
    isconverge = z1 * chi1 <= 0 && z2 * chi2 <= 0
    isconverge
end

function print_state(filename::AbstractString, sim::ACASX_GM)
    if !endswith(filename, ".csv")
        filename *= ".csv"
    end
    state = extract_state(DataFrame, sim)
    writetable(filename, state)
end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using HTTP

function callapi(
    url::String,
    initial_delay::Int64=3,
    delay::Int64=3
)
    sleep(initial_delay)
    try
        response = HTTP.request("GET", url; verbose=1, retries=4)
        if response.status != 200
            sleep(delay)
            error("GoogleApiInvalidRequestError")
        else
            @info "response status: $(response.status)"

            return response
        end
    catch e
        @debug "response status: $(response.status)"
    end
    error("Failed to call $url")
end


text:
# module Lupin

struct Argument
    args
    kwargs
    function Argument(args...; kwargs...)
        new(args, kwargs)
    end
end

function apply(f, a::Argument)
    f(a.args...; a.kwargs...)
end

# module Lupin


text:
using DataDeps

register(DataDep(
  "lpga-2009",
  "http://www.stat.ufl.edu/~winner/data/lpga2009.txt",
  "http://www.stat.ufl.edu/~winner/data/lpga2009.dat",
  "9f578d992e8e60a2727fe738bda266dfb0d518a0f13ba81d08bc79ad0fcd695d",
  post_fetch_method=(path -> begin
    # First char of each row is whitespace, remove
    write(path, join(strip.(readlines(path)), "\n"))

    UCIData.process_dataset(path,
        target_index=8,
        feature_indices=[2:6; 9:14],
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using UncertaintyQuantification

l = Parameter(1.8, :l) # length
b = Parameter(0.12, :b) # width

h = RandomVariable(Normal(0.24, 0.01), :h) # height

μ = log(10e9^2 / sqrt(1.6e9^2 + 10e9^2))
σ = sqrt(log(1.6e9^2 / 10e9^2 + 1))
E = RandomVariable(LogNormal(μ, σ), :E) # young's modulus

μ = log(5000^2 / sqrt(400^2 + 5000^2))
σ = sqrt(log(400^2 / 5000^2 + 1))
P = RandomVariable(LogNormal(μ, σ), :P) # tip load

μ = log(600^2 / sqrt(140^2 + 600^2))
σ = sqrt(log(140^2 / 600^2 + 1))
ρ = RandomVariable(LogNormal(μ, σ), :ρ) # density

c = GaussianCopula([1 0.8; 0.8 1])
jd = JointDistribution([E ρ], c)

inputs = [l, b, h, P, jd]

inertia = Model(df -> df.b .* df.h .^ 3 / 12, :I)

displacement = Model(
    df ->
        (df.ρ .* 9.81 .* df.b .* df.h .* df.l .^ 4) ./ (8 .* df.E .* df.I) .+
        (df.P .* df.l .^ 3) ./ (3 .* df.E .* df.I),
    :w,
)

max_displacement = 0.01

# Compute probability of failure using standard Monte Carlo
mc = MonteCarlo(10^6)

mc_pf, mc_samples = probability_of_failure(
    [inertia, displacement], df -> max_displacement .- df.w, inputs, mc
)

println(
    "Monte Carlo probability of failure $mc_pf ($(size(mc_samples, 1)) model evaluations)"
)

# Compute probability of failure using Line Sampling
ls = LineSampling(50)

ls_pf, ls_samples = probability_of_failure(
    [inertia, displacement], df -> max_displacement .- df.w, inputs, ls
)

println(
    "Line Sampling probability of failure: $ls_pf ($(size(ls_samples, 1)) model evaluations)",
)

subset = UncertaintyQuantification.SubSetSimulation(2000, 0.1, 10, Uniform(-0.5, 0.5))

subset_pf, subset_samples = probability_of_failure(
    [inertia, displacement], df -> max_displacement .- df.w, inputs, subset
)

println(
    "Subset Simulation probability of failure: $subset_pf ($(size(subset_samples, 1)) model evaluations)",
)

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
@testset "Constant" begin
    # Instantiation
    N1 = 10
    A1 = rand(Float64, N1) * 100
    A2 = rand(Float64, N1, N1) * 100
    A3 = rand(Float64, N1, N1, N1) * 100

    for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
        isinplace = constructor == interpolate!
        itp1 = @inferred(constructor(copier(A1), BSpline(Constant())))
        itp2 = @inferred(constructor(copier(A2), BSpline(Constant())))
        itp3 = @inferred(constructor(copier(A3), BSpline(Constant())))

        @test parent(itp1) === itp1.coefs
        @test Interpolations.lbounds(itp1) == (1,)
        @test Interpolations.ubounds(itp1) == (N1,)

        # Evaluation on provided data points
        for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
            check_axes(itp, A, isinplace)
            check_inbounds_values(itp, A)
            check_oob(itp)
            can_eval_near_boundaries(itp)
        end

        # Evaluation between data points (tests constancy)
        for i in 2:N1-1
            @test A1[i] == itp1(i+.3) == itp1(i+.3) == itp1(i-.3) == itp1(i-.3)
        end
        # 2D
        for i in 2:N1-1, j in 2:N1-1
            @test A2[i,j] == itp2(i+.4,j-.3) == itp2(i+.4,j-.3)
        end
        # 3D
        for i in 2:N1-1, j in 2:N1-1, k in 2:N1-1
            @test A3[i,j,k] == itp3(i+.4,j-.3,k+.1) == itp3(i+.4,j-.3,k+.2)
        end
    end
end


text:
function pack_bits(b::BitArray)::Vector{UInt8}
    paddingsize = 8 - (length(b) % 8)
    with_padding = b
    if paddingsize > 0
        padding::BitArray = ones(paddingsize)
        with_padding = vcat(b, padding)
    end
    bytes = []
    for bidx=1:8:length(with_padding)-8
        subarray = b[bidx:bidx+7]
        str = *([x ? "1" : "0" for x in subarray]...)
        push!(bytes, parse(UInt8, str, base=2))
    end
    return bytes
end

function unpack_bits(v::Vector{UInt8})::BitArray
    b::BitArray = []
    for byte in v
        for char in lpad(string(byte, base=2), 8, '0')
            push!(b, char == '1' ? 1 : 0)
        end
    end
    @assert length(v) == length(b)/8
    return b
end

function segments_into_bytes(s::Vector{Vector{UInt8}})
    bytes = []
    for segment in s
        push!(bytes, segment)
    end
    return vcat(bytes...)
end

function split_into_segments(b::Vector{UInt8})::Vector{Vector{UInt8}}
    markers = []
    for i=1:length(b)-1
        if b[i] == 0xFF && b[i+1] != 0x00
            push!(markers, i)
        end
    end
    segments = []
    for m=1:length(markers)-1
        start = markers[m]
        stop = markers[m+1]-1
        push!(segments, b[start:stop])
    end
    push!(segments, b[markers[end]:end])
    return segments
end

function embed_jpeg_scan(scan::BitArray)::Vector{Vector{UInt8}}
    @assert scan[end-3:end] == [1,1,0,0]
    soi = [0xFF, 0xD8]
    eoi = [0xFF, 0xD9]
    scan_seg = [
    0xFF, 0xDA, #Marker
    0x00, 0x0C, #Length
    0x03, #Number of Components
    0x01, 0x00, #Y component
    0x02, 0x11, #Cb component
    0x03, 0x11, #Cr component
    0x00, 0x00, 0x00 # mandatory skip these
    ]
    sos_data = pack_bits(scan)
    fixed_scan_data = []
    for b in sos_data
        push!(fixed_scan_data, b)
        if b == 0xFF
            push!(fixed_scan_data, 0x00) # Sincd 0xFF is used for markers, add 0x00 after a 0xFF in scan data
        end
    end
    with_scan_data = vcat(scan_seg, fixed_scan_data)
    return [soi, with_scan_data, eoi]
end

function extract_jpeg_scan(jpeg::Vector{Vector{UInt8}})::BitArray
    with_scan_data = jpeg[2]
    scan_data = with_scan_data[15:end]
    fixed_scan_data::Vector{UInt8} = []
    idx = 1
    while idx <= length(scan_data)
        push!(fixed_scan_data, scan_data[idx])
        if scan_data[idx] == 0xFF && scan_data[idx+1] == 0x00
            idx +=1
        end
        idx += 1
    end
    bits = unpack_bits(fixed_scan_data)
    return bits
end

function serialize_jpeg_scan(scan::BitArray)::Vector{UInt8}
    return segments_into_bytes(embed_jpeg_scan(scan))
end

function deserialize_into_jpeg_scan(v::Vector{UInt8})::BitArray
    jpeg_scan = extract_jpeg_scan(split_into_segments(v))
    return jpeg_scan
end


text:
module DispatchOps

export
    # simplenexus
    Graph, Digraph, Arc, ArcIter,
    add_arc!, rem_arc!, add_node!, rem_node!,
    nn, na, nodes, arcs, fadj, badj, has_arc, has_node,
    src, tgt, key,

    # metanexus
    MetaGraph, MetaDigraph, props, get_prop, has_prop,
    set_props!, set_prop!, rem_prop!, clear_props!,
    filter_arcs, filter_nodes

include("Nexus/Nexus.jl")
include("Model/Model.jl")
using .Nexus
using .Model

end


text:
include("utils.jl")
include("autoutils.jl")

export transferentropy, transferentropy!

function get_marginals(s, t, emb::EmbeddingTE)
    pts, vars, τs, js = te_embed(s, t, emb)

    # Get marginals
    ST = pts[:, [vars.S; vars.T]]
    T𝒯 = pts[:, [vars.𝒯; vars.T]]
    T = pts[:, vars.T]
    joint = pts

    return joint, ST, T𝒯, T
end

function get_marginals(s, t, c, emb::EmbeddingTE)
    pts, vars, τs, js = te_embed(s, t, c, emb)

    # Get marginals
    ST = pts[:, [vars.S; vars.T; vars.C]]
    T𝒯 = pts[:, [vars.𝒯; vars.T; vars.C]]
    T = pts[:, [vars.T; vars.C]]
    joint = pts

    return joint, ST, T𝒯, T
end

# map a set of pre-embedded points to the correct marginals for transfer entropy computation
function get_marginals(pts, emb::TEVars)
    # Get marginals
    ST = pts[:, [vars.S; vars.T]]
    T𝒯 = pts[:, [vars.𝒯; vars.T]]
    T = pts[:, vars.T]
    joint = pts

    return joint, ST, T𝒯, T
end

"""
Abstract types for transfer entropy estimators that are not already implemented as basic
entropy estimators in Entropies.jl, but needs some kind of special treatment to  work
for transfer entropy. """
abstract type TransferEntropyEstimator <: EntropyEstimator end

"""
    transferentropy(s, t, [c,] est; base = 2, q = 1,
        τT = -1, τS = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1, [τC = -1, dC = 1]
    )

Estimate transfer entropy[^Schreiber2000] from source `s` to target `t`, ``TE^{q}(s \\to t)``, using the
provided entropy/probability estimator `est` with logarithms to the given `base`. Optionally, condition
on `c` and estimate the conditional transfer entropy ``TE^{q}(s \\to t | c)``.
The input series `s`, `t`, and `c` must be equal-length real-valued vectors.

Compute either Shannon transfer entropy (`q = 1`, which is the default) or the order-`q`
Rényi transfer entropy[^Jizba2012] by setting `q` different from 1.

All possible estimators that can be used are described in the online documentation.

## Keyword Arguments
Keyword arguments tune the embedding that will be done to each of the timeseries
(with more details following below).
In short, the embedding lags `τT`, `τS`, `τC` must be zero or negative, the
prediction lag `η𝒯` must be positive, and the embedding dimensions `dT`, `dS`, `dC`, `d𝒯`
must be greater than or equal to 1. Thus, the convention is to use negative lags to
indicate embedding delays for past state vectors (for the ``T``, ``S`` and ``C`` marginals,
detailed below), and positive lags to indicate embedding delays for future state vectors
(for the ``\\mathcal T`` marginal, also detailed below).

The default behaviour is to use scalar timeseries for past state vectors
(in that case, the `τT`, `τS` or `τC` does not affect the analysis).

## Description

### Transfer entropy on scalar time series

Transfer entropy[^Schreiber2000] between two simultaneously measured scalar time series ``s(n)`` and ``t(n)``,
``s(n) = \\{ s_1, s_2, \\ldots, s_N \\} `` and ``t(n) = \\{ t_1, t_2, \\ldots, t_N \\} ``, is
is defined as

```math
TE(s \\to t) = \\sum_i p(s_i, t_i, t_{i+\\eta}) \\log \\left( \\dfrac{p(t_{i+\\eta} | t_i, s_i)}{p(t_{i+\\eta} | t_i)} \\right)
```

### Transfer entropy on generalized embeddings

By defining the vector-valued time series, it is possible to include more than one
historical/future value for each marginal (see 'Uniform vs. non-uniform embeddings' below
for embedding details):

- ``\\mathcal{T}^{(d_{\\mathcal T}, \\eta_{\\mathcal T})} = \\{t_i^{(d_{\\mathcal T}, \\eta_{\\mathcal T})} \\}_{i=1}^{N}``,
- ``T^{(d_T, \\tau_T)} = \\{t_i^{(d_T, \\tau_T)} \\}_{i=1}^{N}``,
- ``S^{(d_S, \\tau_S)} = \\{s_i^{(d_T, \\tau_T)} \\}_{i=1}^{N}``,  and
- ``C^{(d_C, \\tau_C)} = \\{s_i^{(d_C, \\tau_C)} \\}_{i=1}^{N}``.

The non-conditioned generalized and conditioned generalized forms of the transfer entropy are then

```math
TE(s \\to t) = \\sum_i p(S,T, \\mathcal{T}) \\log \\left( \\dfrac{p(\\mathcal{T} | T, S)}{p(\\mathcal{T} | T)} \\right)
```

```math
TE(s \\to t | c) = \\sum_i p(S,T, \\mathcal{T}, C) \\log \\left( \\dfrac{p(\\mathcal{T} | T, S, C)}{p(\\mathcal{T} | T, C)} \\right)
```

### Uniform vs. non-uniform embeddings

The `N` state vectors for each marginal are either

- uniform, of the form ``x_{i}^{(d, \\omega)} = (x_i, x_{i+\\omega}, x_{i+2\\omega}, \\ldots x_{i+(d - 1)\\omega})``,
    with equally spaced state vector entries. *Note: When constructing marginals for ``T``, ``S`` and ``C``,
    we need ``\\omega \\leq 0`` to get present/past values, while ``\\omega > 0`` is necessary to get future states
    when constructing ``\\mathcal{T}``.*
- non-uniform, of the form ``x_{i}^{(d, \\omega)} = (x_i, x_{i+\\omega_1}, x_{i+\\omega_2}, \\ldots x_{i+\\omega_{d}})``,
    with non-equally spaced state vector entries ``\\omega_1, \\omega_2, \\ldots, \\omega_{d}``,
    which can be freely chosen. *Note: When constructing marginals for ``T``, ``S`` and ``C``,
    we need ``\\omega_i \\leq 0`` for all ``\\omega_i`` to get present/past values, while ``\\omega_i > 0`` for all ``\\omega_i``
    is necessary to get future states when constructing ``\\mathcal{T}``.*

In practice, the `dT`-dimensional, `dS`-dimensional and `dC`-dimensional state vectors
comprising ``T``, ``S`` and ``C`` are constructed with embedding lags `τT`,
`τS`, and `τC`, respectively. The `d𝒯`-dimensional future states ``\\mathcal{T}^{(d_{\\mathcal T}, \\eta_{\\mathcal T})}``
are constructed with prediction lag `η𝒯` (i.e. predictions go from present/past states to
future states spanning a maximum of `d𝒯*η𝒯` time steps).
*Note: in Schreiber's paper, only the historical states are defined as
potentially higher-dimensional, while the future states are always scalar.*

### Estimation

Transfer entropy is here estimated by rewriting the above expressions as a sum of marginal
entropies, and extending the definitions above to use Rényi generalized entropies of order
`q` as

```math
TE^{q}(s \\to t) = H^{q}(\\mathcal T, T) + H^{q}(T, S) - H^{q}(T) - H^{q}(\\mathcal T, T, S),
```

```math
TE^{q}(s \\to t | c) = H^{q}(\\mathcal T, T, C) + H^{q}(T, S, C) - H^{q}(T, C) - H^{q}(\\mathcal T, T, S, C),
```

where ``H^{q}(\\cdot)`` is the generalized Rényi entropy of order ``q``. This is equivalent
to the Rényi transfer entropy implementation in Jizba et al. (2012)[^Jizba2012].

## Examples

Default estimation (scalar marginals):

```julia
# Symbolic estimator, motifs of length 4, uniform delay vectors with lag 1
est = SymbolicPermutation(m = 4, τ = 1)

x, y = rand(100), rand(100)
transferentropy(x, y, est)
```

Increasing the dimensionality of the ``T`` marginal (present/past states of the target
variable):

```julia
# Binning-based estimator
est = VisitationFrequency(RectangularBinning(4))
x, y = rand(100), rand(100)

# Uniform delay vectors when `τT` is an integer (see explanation above)
# Here t_{i}^{(dT, τT)} = (t_i, t_{i+τ}, t_{i+2τ}, \\ldots t_{i+(dT-1)τ})
# = (t_i, t_{i-2}, t_{i-4}, \\ldots t_{i-6τ}), so we need zero/negative values for `τT`.
transferentropy(x, y, est, dT = 4, τT = -2)

# Non-uniform delay vectors when `τT` is a vector of integers
# Here t_{i}^{(dT, τT)} = (t_i, t_{i+τ_{1}}, t_{i+τ_{2}}, \\ldots t_{i+τ_{dT}})
# = (t_i, t_{i-7}, t_{i-25}), so we need zero/negative values for `τT`.
transferentropy(x, y, est, dT = 3, τT = [0, -7, -25])
```

Logarithm bases and the order of the Rényi entropy can also be tuned:

```julia
x, y = rand(100), rand(100)
est = NaiveKernel(0.3)
transferentropy(x, y, est, base = MathConstants.e, q = 2) # TE in nats, order-2 Rényi entropy
```

[^Schreiber2000]: Schreiber, T. (2000). Measuring information transfer. Physical review letters, 85(2), 461.
[^Jizba2012]: Jizba, P., Kleinert, H., & Shefaat, M. (2012). Rényi’s information transfer between financial time series. Physica A: Statistical Mechanics and its Applications, 391(10), 2971-2989.
"""
function transferentropy end
function transferentropy! end

# Estimate transfer entropy from time series by first embedding them and getting required
# marginals.
function transferentropy(s, t, est::Est; base = 2, q = 1,
        τT = -1, τS = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1)

    emb = EmbeddingTE(τT = τT, τS = τS, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯)
    joint, ST, T𝒯, T = get_marginals(s, t, emb)

    _transferentropy(joint, ST, T𝒯, T, est; base = base, q = q)
end

function transferentropy(s, t, c, est::Est; base = 2, q = 1,
        τT = -1, τS = -1, τC = -1, η𝒯 = 1, dT = 1, dS = 1, dC = 1, d𝒯 = 1)

    emb = EmbeddingTE(τT = τT, τS = τS, τC = τC, η𝒯 = η𝒯, dT = dT, dS = dS, dC = dC, d𝒯 = d𝒯)
    joint, ST, T𝒯, T = get_marginals(s, t, c, emb)

    _transferentropy(joint, ST, T𝒯, T, est; base = base, q = q)
end

transferentropy(s::AbstractVector{<:Real}, t::AbstractVector{<:Real}) =
    error("Estimator missing. Please provide a valid estimator as the third argument.")

transferentropy(s::AbstractVector{<:Real}, t::AbstractVector{<:Real}, c::AbstractVector{<:Real}) =
    error("Estimator missing. Please provide a valid estimator as the fourth argument.")

# estimate transfer entropy from marginal entropies, as described in docstring
function _transferentropy(joint, ST, T𝒯, T, est::Est; base = 2, q = 1)
    te = genentropy(T𝒯, est, base = base, q = q) +
        genentropy(ST, est, base = base, q = q) -
        genentropy(T, est, base = base, q = q) -
        genentropy(joint, est, base = base, q = q)
end

# TODO: estimate using mutual information decomposition,
# function transferentropy(marginal1, marginal2, est; base = 2, q = 1)

include("symbolic.jl")
include("hilbert.jl")
include("nearestneighbor.jl")
include("transferoperator.jl")

# automated approaches
include("bbnue.jl")


text:
# This file implements:
#   - LinkedList: a linked list datastructure used by DependencyGraph.
#   - DependencyGraph: a leightweight adjacency-list based directed graph datastructure.
#   - children: a DFS-like algorithm to find all children of a set of vertices in a DependencyGraph.

# DependencyGraph is used as an intermediate datastructure by ForneyLab's message scheduling algorithms.
# The general recipe for building a message scheduling involves at least two steps:
#   1. Build a DependencyGraph (DG) based on a FactorGraph.
#      The vertices in the DG are (for example) node interfaces.
#      An edge v -> w in the DG denotes that the message out of interface v depends on the message out of w.
#   2. Invoke the children algorithm on the DG to find a valid sequential message passing schedule.

###############################
# Linked list datastructure
###############################

mutable struct LinkedListElement{T}
    payload::T
    next::LinkedListElement{T}

    LinkedListElement{T}(payload::T) where T = new(payload)
end

LinkedListElement(payload::T) where T = LinkedListElement{T}(payload)

mutable struct LinkedList{T}
    first::LinkedListElement{T}
    last::LinkedListElement{T}

    LinkedList{T}() where T = new{T}()
end

function Base.push!(list::LinkedList, payload)
    if isdefined(list, :first)
        list.last.next = LinkedListElement(payload)
        list.last = list.last.next
    else
        list.first = list.last = LinkedListElement(payload)
    end

    return list
end

Base.isempty(list::LinkedList) = !isdefined(list, :first)

function Base.collect(list::LinkedList{T}) where T
    items = T[]
    if isdefined(list, :first)
        entry = list.first
        while isdefined(entry, :next)
            push!(items, entry.payload)
            entry = entry.next
        end
        push!(items, entry.payload) # last entry
    end

    return items
end

###############################
# DependencyGraph
###############################

"""
A `DependencyGraph` is a directed graph in which an edge `v -> w`
represents a dependency of vertex `v` on vertex `w`.
Dependency graphs are used for example by message scheduling algorithms.
"""
mutable struct DependencyGraph{VT}
    vertices::Vector{VT}
    neighbors::Vector{LinkedList{Int}}

    DependencyGraph{VT}() where VT = new(Vector{VT}(), Vector{LinkedList{Int}}())
end

function addVertex!(dg::DependencyGraph{VT}, v::VT) where VT
    push!(dg.vertices, v)
    push!(dg.neighbors, LinkedList{Int}())
    return dg
end

function addEdge!(dg::DependencyGraph{VT}, from::VT, to::VT) where VT
    v = something(findfirst(isequal(from), dg.vertices), 0)
    w = something(findfirst(isequal(to), dg.vertices), 0)
    push!(dg.neighbors[v], w)
    return dg
end

neighbors(vertex_idx::Int, dg::DependencyGraph) = collect(dg.neighbors[vertex_idx])

function neighbors(vertex::VT, dg::DependencyGraph{VT}) where VT
    vertex_idx = something(findfirst(isequal(vertex), dg.vertices), 0)
    return dg.vertices[neighbors(vertex_idx, dg)]
end

#########################################
# Graph algorithms for DependencyGraph
#########################################

"""
children(vertices, graph; allow_cycles=false, breaker_sites=[], restrict_to=[])

Return a vector consisting of `vertices` and all their children in `graph`.
`v` is a child of `u` iff there exists a path from `u` to `v`.
The resulting array is sorted in reverse topological order,
i.e. for each directed edge `u -> v`, `v` (child of `u`) appears before `u`.

Optional keyword arguments:

- `allow_cycles`: set to true to accept cycles.
- `breaker_sites`: a Set of vertices on which the search will terminate.
- `restrict_to`: a Set of vertices to restrict the search to.

This function can be used to generate message passing schedules
if `graph` is a dependency graph.
"""
function find_vertex_indexes(vertex::V, graph::DependencyGraph{V}) where V
    return something(findfirst(isequal(vertex), graph.vertices), 0)
end

function children(  vertices::Vector{V},
                    graph::DependencyGraph{V};
                    allow_cycles::Bool=false,
                    breaker_sites::Set{V}=Set{V}(),
                    restrict_to::Set{V}=Set{V}()) where V

    # Find vertex indexes of breaker_sites
    breaker_vertices = Set{Int}((find_vertex_indexes(v, graph) for v=breaker_sites))
    restrict_to_vertices = Set{Int}((find_vertex_indexes(v, graph) for v=restrict_to))

    visited = Int[] # Hold topologically sorted list of indices of child vertices

    for root in vertices
        # Iterative depth-first search
        stack = Int[something(findfirst(isequal(root), graph.vertices), 0)]
        visited_from_root = Int[]
        while !isempty(stack)
            v = pop!(stack)

            # Check search restrictions
            (v in breaker_vertices) && continue
            (isempty(restrict_to_vertices) || (v in restrict_to_vertices)) || continue

            # Check if vertex v has already been visited
            (v in visited) && continue
            if v in visited_from_root
                if allow_cycles
                    continue
                else
                    throw(ArgumentError("The input graph contains a loop around $(graph.vertices[v])."))
                end
            end

            push!(visited_from_root, v) # Mark v as visited

            # Add all adjacent vertices to the stack
            for w in neighbors(v, graph)
                push!(stack, w)
            end
        end

        # Add the newly visited children to the global visited list
        visited = vcat(visited, reverse(visited_from_root))
    end

    return graph.vertices[visited] # Return actual vertices instead of their indexes
end

children(vertex::V, graph::DependencyGraph{V}; kwargs...) where V = children(V[vertex], graph; kwargs...)


text:
#!/usr/bin/env python
# vim: set fileencoding=utf-8
""" This file is a Julia translation of the MATLAB file prony.m

 Julia version by Adrien Pajon 28 Sept 2020
 Copyright notice from prony.m:
 copyright 1996, by M.H. Hayes.  For use with the book
 "Statistical Digital Signal Processing and Modeling"
 (John Wiley & Sons, 1996).
"""

include("convm.jl")

function prony(x, p, q)
    # Model a signal using Prony's method
    #
    # Usage: [b,a,err] = prony(x,p,q)
    #
    # The input sequence x is modeled as the unit sample response of
    # a filter having a system function of the form
    #     H(z) = B(z)/A(z)
    # The polynomials B(z) and A(z) are formed from the vectors
    #     b=[b(0), b(1), ... b(q)]
    #     a=[1   , a(1), ... a(p)]
    # The input q defines the number of zeros in the model
    # and p defines the number of poles. The modeling error is
    # returned in err.
    #
    # This comes from Hayes, p. 149, 153, etc

    # x = x[:]
    # N = length(x)
    # if p + q >= length(x)
    #     print("ERROR: model order too large")
    #     print(
    #         "p q length(x) " +
    #         convert(String, p) +
    #         " " +
    #         convert(String, q) +
    #         " " +
    #         convert(String, length(x)),
    #     )
    #     return false
    # end
    #
    # # This formulation uses eq. 4.50, p. 153
    # # Set up the convolution matrices
    # X = convm(x, p + 1)
    # Xq = X[q+1:N+p-1, 1:p]
    # xq1 = -X[q+2:N+p, 1]
    # # Solve for denominator coefficients
    # a = []
    # if p > 0
    #     a = (Xq \ xq1)
    #     a = [1; a] # a(1) is 1
    # end
    #
    # # Solve for the model error
    # err = Array(conj(x[q+2:N])') * X[q+2:N, 1:p+1]
    # err = err * a
    #
    # # Solve for numerator coefficients
    # if q > 0
    #     # (This is the same as for Pad?)
    #     b = X[1:q+1, 1:p+1] * a
    # else
    #     # all-pole model
    #     # b(0) is x(0), but a better solution is to match energy
    #     b = sqrt(err)
    # end

    x   = x[:];

    N   = length(x);

    if p+q>=length(x)
        error("Model order too large")
        print(
            "p q length(x) " +
            convert(String, p) +
            " " +
            convert(String, q) +
            " " +
            convert(String, length(x))
        )
    end

    X   = convm(x,p+1);

    Xq  = X[q+1:N+p-1,1:p];

    a   = [1;-Xq\X[q+2:N+p,1]];

    b   = X[1:q+1,1:p+1]*a;

    err = x[q+2:N]'*X[q+2:N,1:p+1]*a;

    return (b, a, err)
end

####################################################################
# Example of code to use prony()

# function main()
#     # Test driver
#     # From pp. 149-150
#     x = ones(21)
#     p = q = 1
#     b, a, err = prony(x, p, q)
#     print("a: ")
#     println(a)
#     print("b: ")
#     println(b)
#     print("err: ")
#     println(err)
#
#     # From pp. 152-153
#     # Note that these results don't match the book, but they do match the
#     # MATLAB version. So I'm either setting things up wrong or this is an
#     # errata in the book.
#     p = q = 5
#     nd = 5
#     n = 0:11
#     i = sinc.((n .- nd) ./ 2) ./ 2
#     b, a, err = prony(i, p, q)
#     print("a: ")
#     println(a)
#     print("b: ")
#     println(b)
#     print("err: ")
#     println(err)
# end
# main()