text:
@testset "Basic tests" begin
    vs = 1:100
    n = 10
    ch = ConcurrentCircularBuffer{Int}(n)

    # pop
    append!(ch, vs)
    rv = zeros(Int, n)
    for i in 1:n
        v = pop!(ch)
        rv[i] = v
    end
    @test rv == reverse(vs[end-n+1:end])

    # pop several
    append!(ch, vs)
    rv = zeros(Int, n)
    i = pop!(rv, ch)
    @test i == n
    @test rv == reverse(vs[end-n+1:end])

    # popfirst
    append!(ch, vs)
    rv = zeros(Int, n)
    for i in 1:n
        v = popfirst!(ch)
        rv[i] = v
    end
    @test rv == vs[end-n+1:end]

    # popfirst several
    append!(ch, vs)
    rv = zeros(Int, n)
    i = popfirst!(rv, ch)
    @test i == n
    @test rv == vs[end-n+1:end]
end

@testset "Threading" begin
    n = 1000

    # one producer and one consumer
    ch = ConcurrentCircularBuffer{Int}(n)
    rv = zeros(Int, n)
    consumer = Threads.@spawn begin
        try
            i = 1
            while true
                rv[i] = popfirst!(ch)
                i += 1
            end
        catch e
            if e isa InvalidStateException
            else
                rethrow()
            end
        end
    end
    producer = Threads.@spawn for i in 1:n
        push!(ch, i)
    end
    wait(producer)
    close(ch)
    wait(consumer)
    @test rv == 1:n

    # one producer and two consumers
    ch = ConcurrentCircularBuffer{Int}(n)
    rv1 = zeros(Int, n)
    ts1 = zeros(n)
    consumer1 = Threads.@spawn begin
        try
            for i in 1:n
                ts1[i] = @elapsed rv1[i] = popfirst!(ch)
            end
        catch e
            if e isa InvalidStateException
            else
                rethrow()
            end
        end
    end
    rv2 = zeros(Int, n)
    ts2 = zeros(n)
    consumer2 = Threads.@spawn begin
        try
            for i in 1:n
                ts2[i] = @elapsed rv2[i] = popfirst!(ch)
            end
        catch e
            if e isa InvalidStateException
            else
                rethrow()
            end
        end
    end
    tsc = zeros(4n)
    producer = Threads.@spawn for i in 1:4n
        tsc[i] = @elapsed push!(ch, i)
    end
    wait(producer)
    close(ch)
    wait(consumer1)
    wait(consumer2)
    @test issorted(rv1)
    @test issorted(rv2)
    @test length(intersect(Set(rv1), Set(rv2))) == 0
    @info "consumer1 throughput: $(1 / (sum(ts1) / n)), consumer2 throughput: $(1 / (sum(ts2) / n))"

    # create many threads and make sure all exit cleanly
    ch = ConcurrentCircularBuffer{Int}(n)
    ntasks = 100
    tasks = Vector{Task}(undef, ntasks)
    for i in 1:ntasks
        tasks[i] = Threads.@spawn begin
            try
                sleep(0)
                popfirst!(ch)
            catch e
                if e isa InvalidStateException
                else
                    rethrow()
                end
            end
        end
    end
    sleep(0)
    close(ch)
    for i in 1:ntasks
        wait(tasks[i])
    end
end


text:
__precompile__(true)

module SequenceAlignment
export levenshteinAlign
export levenshteinDistance
export ldn

include("levenshtein.jl")

end

ext: jl


text:
using NNlib: scatter, scatter!

dsts = Dict(
    0 => [3, 4, 5, 6, 7],
    1 => [3 3 4 4 5;
          5 5 6 6 7],
)
srcs = Dict(
    (0, true) => ones(Int, 3, 4),
    (0, false) => ones(Int, 3) * collect(1:4)',
    (1, true) => ones(Int, 2, 3, 4),
    (1, false) => [1, 2] .* reshape(ones(Int, 3) * collect(1:4)', 1,3,4),
)
idxs = Dict(
    :int => [1 2 3 4;
             4 2 1 3;
             3 5 5 3],
    :tup => [(1,) (2,) (3,) (4,);
             (4,) (2,) (1,) (3,);
             (3,) (5,) (5,) (3,)],
    :car => CartesianIndex.(
            [(1,) (2,) (3,) (4,);
             (4,) (2,) (1,) (3,);
             (3,) (5,) (5,) (3,)]),
)
res = Dict(
    (+, 0, true) => [5, 6, 9, 8, 9],
    (+, 1, true) => [5 5 8 6 7;
                     7 7 10 8 9],
    (+, 0, false) => [4, 4, 12, 5, 5],
    (+, 1, false) => [4 4 12 5 5;
                      8 8 24 10 10],
    (-, 0, true) => [1, 2, 1, 4, 5],
    (-, 1, true) => [1 1 0 2 3;
                     3 3 2 4 5],
    (-, 0, false) => [-4, -4, -12, -5, -5],
    (-, 1, false) => [-4 -4 -12 -5 -5;
                      -8 -8 -24 -10 -10],
    (max, 0, true) => [3, 4, 5, 6, 7],
    (max, 1, true) => [3 3 4 4 5;
                       5 5 6 6 7],
    (max, 0, false) => [3, 2, 4, 4, 3],
    (max, 1, false) => [3 2 4 4 3;
                        6 4 8 8 6],
    (min, 0, true) => [1, 1, 1, 1, 1],
    (min, 1, true) => [1 1 1 1 1;
                 1 1 1 1 1],
    (min, 0, false) => [1, 2, 1, 1, 2],
    (min, 1, false) => [1 2 1 1 2;
                        2 4 2 2 4],
    (*, 0, true) => [3, 4, 5, 6, 7],
    (*, 1, true) => [3 3 4 4 5;
               5 5 6 6 7],
    (*, 0, false) => [3, 4, 48, 4, 6],
    (*, 1, false) => [3 4 48 4 6;
                      12 16 768 16 24],
    (/, 0, true) => [0.75, 1., 0.3125, 1.5, 1.75],
    (/, 1, true) => [0.75 0.75 0.25 1. 1.25;
               1.25 1.25 0.375 1.5 1.75],
    (/, 0, false) => [1//3, 1//4, 1//48, 1//4, 1//6],
    (/, 1, false) => [1//3 1//4 1//48 1//4 1//6;
                      1//12 1//16 1//768 1//16 1//24],
    (mean, 0, true) => [4., 5., 6., 7., 8.],
    (mean, 1, true) => [4. 4. 5. 5. 6.;
                  6. 6. 7. 7. 8.],
    (mean, 0, false) => [2, 2, 3, 2.5, 2.5],
    (mean, 1, false) => [2. 2. 3. 2.5 2.5;
                         4. 4. 6. 5. 5.],
)

types = [UInt8, UInt16, UInt32, UInt64, UInt128,
         Int8, Int16, Int32, Int64, Int128, BigInt,
         Float16, Float32, Float64, BigFloat, Rational]

@testset "scatter" begin
    for T = types
        @testset "$T" begin
            PT = promote_type(T, Int)
            @testset "+" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(+, T.(copy(dsts[dims])), T.(srcs[(dims, mutated)]), idx) == T.(res[(+, dims, mutated)])
                    @test scatter!(+, T.(copy(dsts[dims])), srcs[(dims, mutated)], idx) == PT.(res[(+, dims, mutated)])
                    @test scatter!(+, copy(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == PT.(res[(+, dims, mutated)])

                    mutated = false
                    @test scatter(+, T.(srcs[(dims, mutated)]), idx) == T.(res[(+, dims, mutated)])
                end
            end

            @testset "-" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(-, T.(copy(dsts[dims])), T.(srcs[(dims, mutated)]), idx) == T.(res[(-, dims, mutated)])
                    @test scatter!(-, T.(copy(dsts[dims])), srcs[(dims, mutated)], idx) == PT.(res[(-, dims, mutated)])
                    @test scatter!(-, copy(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == PT.(res[(-, dims, mutated)])

                    mutated = false
                    if !(T in [UInt8, UInt16, UInt32, UInt64, UInt128])
                        @test scatter(-, T.(srcs[(dims, mutated)]), idx) == T.(res[(-, dims, mutated)])
                    end
                end
            end

            @testset "max" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(max, T.(copy(dsts[dims])), T.(srcs[(dims, mutated)]), idx) == T.(res[(max, dims, mutated)])
                    @test scatter!(max, T.(copy(dsts[dims])), srcs[(dims, mutated)], idx) == PT.(res[(max, dims, mutated)])
                    @test scatter!(max, copy(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == PT.(res[(max, dims, mutated)])

                    mutated = false
                    if !(T in [BigInt])
                        @test scatter(max, T.(srcs[(dims, mutated)]), idx) == T.(res[(max, dims, mutated)])
                    end
                end
            end

            @testset "min" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(min, T.(copy(dsts[dims])), T.(srcs[(dims, mutated)]), idx) == T.(res[(min, dims, mutated)])
                    @test scatter!(min, T.(copy(dsts[dims])), srcs[(dims, mutated)], idx) == PT.(res[(min, dims, mutated)])
                    @test scatter!(min, copy(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == PT.(res[(min, dims, mutated)])

                    mutated = false
                    if !(T in [BigInt])
                        @test scatter(min, T.(srcs[(dims, mutated)]), idx) == T.(res[(min, dims, mutated)])
                    end
                end
            end

            @testset "*" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(*, T.(copy(dsts[dims])), T.(srcs[(dims, mutated)]), idx) == T.(res[(*, dims, mutated)])
                    @test scatter!(*, T.(copy(dsts[dims])), srcs[(dims, mutated)], idx) == PT.(res[(*, dims, mutated)])
                    @test scatter!(*, copy(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == PT.(res[(*, dims, mutated)])

                    mutated = false
                    if !(T in [UInt8, Int8])
                        @test scatter(*, T.(srcs[(dims, mutated)]), idx) == T.(res[(*, dims, mutated)])
                    end
                end
            end
        end
    end

    for T = [Float16, Float32, Float64, BigFloat, Rational]
        @testset "$T" begin
            PT = promote_type(T, Float64)
            @testset "/" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(/, T.(dsts[dims]), T.(srcs[(dims, mutated)].*2), idx) == T.(res[(/, dims, mutated)])
                    @test scatter!(/, T.(dsts[dims]), srcs[(dims, mutated)].*2, idx) == PT.(res[(/, dims, mutated)])
                    @test scatter!(/, T.(dsts[dims]), T.(srcs[(dims, mutated)].*2), idx) == PT.(res[(/, dims, mutated)])

                    mutated = false
                    @test scatter(/, T.(srcs[(dims, mutated)]), idx) == T.(res[(/, dims, mutated)])
                end
            end

            @testset "mean" begin
                for idx = values(idxs), dims = [0, 1]
                    mutated = true
                    @test scatter!(mean, T.(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == T.(res[(mean, dims, mutated)])
                    @test scatter!(mean, T.(dsts[dims]), srcs[(dims, mutated)], idx) == PT.(res[(mean, dims, mutated)])
                    @test scatter!(mean, copy(dsts[dims]), T.(srcs[(dims, mutated)]), idx) == PT.(res[(mean, dims, mutated)])

                    mutated = false
                    @test scatter(mean, T.(srcs[(dims, mutated)]), idx) == T.(res[(mean, dims, mutated)])
                end
            end
        end
    end

    @test_throws AssertionError scatter!(+, dsts[0], srcs[(1, true)], idxs[:int])
    idx = [1 2 3 4; 4 2 1 3; 6 7 8 9]
    @test_throws BoundsError scatter!(+, dsts[1], srcs[(1, true)], idx)
end

@testset "∇scatter" begin
    T = Float64
    fdm(op) = op == min ? :backward : :forward
    # fdm(op) = :forward

    @testset "∂dst" begin
        for op in (+, -, *, /, mean, max, min)
            gradtest(xs -> scatter!(op, copy(xs), srcs[(0, true)], idxs[:int]), T.(dsts[0]), fdm=fdm(op))
            gradtest(xs -> scatter!(op, copy(xs), srcs[(1, true)], idxs[:int]), T.(dsts[1]), fdm=fdm(op))
        end
    end

    @testset "∂src" begin
        for op in (+, -, *, /, mean, max, min)
            gradtest(xs -> scatter!(op, T.(dsts[0]), xs, idxs[:int]), T.(srcs[(0, true)]), fdm=fdm(op))
            gradtest(xs -> scatter!(op, T.(dsts[1]), xs, idxs[:int]), T.(srcs[(1, true)]), fdm=fdm(op))

            gradtest(xs -> scatter(op, xs, idxs[:int]), T.(srcs[(0, false)]), fdm=fdm(op))
            gradtest(xs -> scatter(op, xs, idxs[:int]), T.(srcs[(1, false)]), fdm=fdm(op))
        end
    end
end


text:
# See https://en.wikipedia.org/wiki/Polynomial

using Polynomials

expressions = [
    :(roots(Polynomial([24 + 7im, -50, 35, -10, 1])))
    :(roots(Polynomial([-26, 24 - 44im, 14 + 33im, -9 - 5im, 1])))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
using EffectiveWaves, Test, LinearAlgebra

# # Had to remove from testset block do to strange scoping error
# function randlm(order)
#     l = rand(0:order)
#     m = rand(-l:l)
#     return l, m
# end

@testset "fourth order matrix test" begin

    L = 2
    L1 = 3

    B = [
        ((l,m),(l2,m2),(dl,dm),(l1,m1))
                for dl = 0:L for dm = -dl:dl
            for l1 = 0:L1 for m1 = -l1:l1
        for l = 0:L for m = -l:l
    for l2 = 0:L1 for m2 = -l2:l2]

    len = (L1+1)^2 * (L+1)^2
    B = reshape(B, (:,len))

    eigF = [
        ((l,m),(l2,m2))
        for l = 0:L for m = -l:l
    for l2 = 0:L1 for m2 = -l2:l2]

    function randlm(order)
        local l = rand(0:order) # need unique names due to scoping
        local m = rand(-l:l)
        return l, m
    end

    l,m = randlm(L)
    l2,m2 = randlm(L1)

    dl,dm = randlm(L)
    l1,m1 = randlm(L1)

    n = lm_to_spherical_harmonic_index(l,m)
    n2 = lm_to_spherical_harmonic_index(l2,m2)
    dn = lm_to_spherical_harmonic_index(dl,dm)
    n1 = lm_to_spherical_harmonic_index(l1,m1)

    @test eigF == [b[3:4] for b in B[(n - 1) * (L1+1)^2 + n2,:]]
    @test eigF == [b[1:2] for b in B[:,(dn - 1) * (L1+1)^2 + n1]]

    @test ((l,m),(l2,m2),(dl,dm),(l1,m1)) == B[(n - 1) * (L1+1)^2 + n2, (dn - 1) * (L1+1)^2 + n1]

end
# @testset "MyTest" begin
#
#     function myfun(m)
#         n = m
#         return n
#     end
#
#     n = myfun(1)
#
#     println("n = ", n)
#
#     myfun(3) # <-- changes n to 3
#
#     println("n = ", n)
#
# end

max_stars_repo_name: arturgower/EffectiveWaves


text:
function frag_edge_channel(in_chan, out_chan, V, EV)
    run_loop = true
    while run_loop
        edgenum = take!(in_chan)
        if edgenum != -1
            put!(out_chan, (edgenum, frag_edge(V, EV, edgenum)))
        else
            run_loop = false
        end
    end
end

function frag_edge(V::LinearAlgebraicRepresentation.Points, EV::LinearAlgebraicRepresentation.ChainOp, edge_idx::Int)
    alphas = Dict{Float64, Int}()
    edge = EV[edge_idx, :]
    verts = V[edge.nzind, :]
    for i in 1:size(EV, 1)
        if i != edge_idx
            intersection = intersect_edges(V, edge, EV[i, :])
            for (point, alpha) in intersection
                verts = [verts; point]
                alphas[alpha] = size(verts, 1)
            end
        end
    end

    alphas[0.0], alphas[1.0] = [1, 2]

    alphas_keys = sort(collect(keys(alphas)))
    edge_num = length(alphas_keys)-1
    verts_num = size(verts, 1)
    ev = spzeros(Int8, edge_num, verts_num)

    for i in 1:edge_num
        ev[i, alphas[alphas_keys[i]]] = 1
        ev[i, alphas[alphas_keys[i+1]]] = 1
    end

    verts, ev
end

function intersect_edges(V::LinearAlgebraicRepresentation.Points, edge1::LinearAlgebraicRepresentation.Cell, edge2::LinearAlgebraicRepresentation.Cell)
    err = 10e-8

    x1, y1, x2, y2 = vcat(map(c->V[c, :], edge1.nzind)...)
    x3, y3, x4, y4 = vcat(map(c->V[c, :], edge2.nzind)...)
    ret = Array{Tuple{LinearAlgebraicRepresentation.Points, Float64}, 1}()

    v1 = [x2-x1, y2-y1];
    v2 = [x4-x3, y4-y3];
    v3 = [x3-x1, y3-y1];

    ang1 = dot(normalize(v1), normalize(v2))
    ang2 = dot(normalize(v1), normalize(v3))

    parallel = 1-err < abs(ang1) < 1+err
    colinear = parallel && (1-err < abs(ang2) < 1+err || -err < norm(v3) < err)

    if colinear
        o = [x1 y1]
        v = [x2 y2] - o
        alpha = 1/dot(v,v')
        ps = [x3 y3; x4 y4]
        for i in 1:2
            a = alpha*dot(v',(reshape(ps[i, :], 1, 2)-o))
            if 0 < a < 1
                push!(ret, (ps[i:i, :], a))
            end
        end

    elseif !parallel
        denom = (v2[2])*(v1[1]) - (v2[1])*(v1[2])
        a = ((v2[1])*(-v3[2]) - (v2[2])*(-v3[1])) / denom
        b = ((v1[1])*(-v3[2]) - (v1[2])*(-v3[1])) / denom

        if -err < a < 1+err && -err <= b <= 1+err
            p = [(x1 + a*(x2-x1))  (y1 + a*(y2-y1))]
            push!(ret, (p, a))
        end
    end

    return ret
end

function merge_vertices!(V::LinearAlgebraicRepresentation.Points, EV::LinearAlgebraicRepresentation.ChainOp, edge_map, err=1e-4)
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    kdtree = KDTree(V')

    todelete = []

    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = LinearAlgebraicRepresentation.inrange(kdtree, V[vi, :], err)

            newverts[nearvs] = i

            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)

            i = i + 1
        end
    end

    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    edges = Array{Tuple{Int, Int}, 1}(edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(edgenum)

    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind

        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))

    end
    nedges = union(edges)
    nedges = filter(t->t[1]!=t[2], nedges)

    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))

    etuple2idx = Dict{Tuple{Int, Int}, Int}()

    for ei in 1:nedgenum
        nEV[ei, collect(nedges[ei])] = 1
        etuple2idx[nedges[ei]] = ei
    end

    for i in 1:length(edge_map)
        row = edge_map[i]
        row = map(x->edges[x], row)
        row = filter(t->t[1]!=t[2], row)
        row = map(x->etuple2idx[x], row)
        edge_map[i] = row
    end

    return LinearAlgebraicRepresentation.Points(nV), nEV
end

function biconnected_components(EV::LinearAlgebraicRepresentation.ChainOp)
    ps = Array{Tuple{Int, Int, Int}, 1}()
    es = Array{Tuple{Int, Int}, 1}()
    todel = Array{Int, 1}()
    visited = Array{Int, 1}()
    bicon_comps = Array{Array{Int, 1}, 1}()
    hivtx = 1

    function an_edge(point)
        edges = setdiff(EV[:, point].nzind, todel)
        if length(edges) == 0
            edges = [false]
        end
        edges[1]
    end

    function get_head(edge, tail)
        setdiff(EV[edge, :].nzind, [tail])[1]
    end

    function v_to_vi(v)
        i = findfirst(t->t[1]==v, ps)
        if i == 0
            return false
        else
            return ps[i][2]
        end
    end

    push!(ps, (1,1,1))
    push!(visited, 1)
    exit = false
    while !exit
        edge = an_edge(ps[end][1])
        if edge != false
            tail = ps[end][2]
            head = get_head(edge, ps[end][1])
            hi = v_to_vi(head)
            if hi == false
                hivtx += 1
                push!(ps, (head, hivtx, ps[end][2]))
                push!(visited, head)
            else
                if hi < ps[end][3]
                    ps[end] = (ps[end][1], ps[end][2], hi)
                end
            end
            push!(es, (edge, tail))
            push!(todel, edge)
        else
            if length(ps) == 1
                found = false
                pop!(ps)
                for i in 1:size(EV,2)
                    if !(i in visited)
                        hivtx = 1
                        push!(ps, (i, hivtx, 1))
                        push!(visited, i)
                        found = true
                        break
                    end
                end
                if !found
                    exit = true
                end

            else
                if ps[end][3] == ps[end-1][2]
                    edges = Array{Int, 1}()
                    while true
                        edge, tail = pop!(es)
                        push!(edges, edge)
                        if tail == ps[end][3]
                            if length(edges) > 1
                                push!(bicon_comps, edges)
                            end
                            break
                        end
                    end

                else
                    if ps[end-1][3] > ps[end][3]
                        ps[end-1] = (ps[end-1][1], ps[end-1][2], ps[end][3])
                    end
                end
                pop!(ps)
            end
        end
    end

    bicon_comps
end

function get_external_cycle(V::LinearAlgebraicRepresentation.Points, EV::LinearAlgebraicRepresentation.ChainOp, FE::LinearAlgebraicRepresentation.ChainOp)
    FV = abs.(FE)*EV
    vs = sparsevec(mapslices(sum, abs.(EV), 1)).nzind
    minv_x1 = maxv_x1 = minv_x2 = maxv_x2 = pop!(vs)
    for i in vs
        if V[i, 1] > V[maxv_x1, 1]
            maxv_x1 = i
        elseif V[i, 1] < V[minv_x1, 1]
            minv_x1 = i
        end
        if V[i, 2] > V[maxv_x2, 2]
            maxv_x2 = i
        elseif V[i, 2] < V[minv_x2, 2]
            minv_x2 = i
        end
    end
    cells = intersect(
        FV[:, minv_x1].nzind,
        FV[:, maxv_x1].nzind,
        FV[:, minv_x2].nzind,
        FV[:, maxv_x2].nzind
    )
    if length(cells) == 1
        return cells[1]
    else
        for c in cells
            if LinearAlgebraicRepresentation.face_area(V, EV, FE[c, :]) < 0
                return c
            end
        end
    end
end
function pre_containment_test(bboxes)
    n = length(bboxes)
    containment_graph = spzeros(Int8, n, n)

    for i in 1:n
        for j in 1:n
            if i != j && LinearAlgebraicRepresentation.bbox_contains(bboxes[j], bboxes[i])
                containment_graph[i, j] = 1
            end
        end
    end

    return containment_graph
end
function prune_containment_graph(n, V, EVs, shells, graph)

    for i in 1:n
        an_edge = shells[i].nzind[1]
        origin_index = EVs[i][an_edge, :].nzind[1]
        origin = V[origin_index, :]

        for j in 1:n
            if i != j
                if graph[i, j] == 1
                    shell_edge_indexes = shells[j].nzind
                    ev = EVs[j][shell_edge_indexes, :]

                    if !LinearAlgebraicRepresentation.point_in_face(origin, V, ev)
                        graph[i, j] = 0
                    end
                end
             end
         end

     end
     return graph
end
function transitive_reduction!(graph)
    n = size(graph, 1)
    for j in 1:n
        for i in 1:n
            if graph[i, j] > 0
                for k in 1:n
                    if graph[j, k] > 0
                        graph[i, k] = 0
                    end
                end
            end
        end
    end
end
function cell_merging(n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)
    function bboxes(V::LinearAlgebraicRepresentation.Points, indexes::LinearAlgebraicRepresentation.ChainOp)
        boxes = Array{Tuple{Any, Any}}(indexes.n)
        for i in 1:indexes.n
            v_inds = indexes[:, i].nzind
            boxes[i] = LinearAlgebraicRepresentation.bbox(V[v_inds, :])
        end
        boxes
    end

    sums = Array{Tuple{Int, Int, Int}}(0);

    for father in 1:n
        if sum(containment_graph[:, father]) > 0
            father_bboxes = bboxes(V, abs.(EVs[father]')*abs.(boundaries[father]'))
            for child in 1:n
                if containment_graph[child, father] > 0
                    child_bbox = shell_bboxes[child]
                    for b in 1:length(father_bboxes)
                        if LinearAlgebraicRepresentation.bbox_contains(father_bboxes[b], child_bbox)
                            push!(sums, (father, b, child))
                            break
                        end
                    end
                end
            end
        end
    end

    EV = vcat(EVs...)
    edgenum = size(EV, 1)
    facenum = sum(map(x->size(x,1), boundaries))
    FE = spzeros(Int8, facenum, edgenum)
    shells2 = spzeros(Int8, length(shells), edgenum)
    r_offsets = [1]
    c_offset = 1
    for i in 1:n
        min_row = r_offsets[end]
        max_row = r_offsets[end] + size(boundaries[i], 1) - 1
        min_col = c_offset
        max_col = c_offset + size(boundaries[i], 2) - 1
        FE[min_row:max_row, min_col:max_col] = boundaries[i]
        shells2[i, min_col:max_col] = shells[i]
        push!(r_offsets, max_row + 1)
        c_offset = max_col + 1
    end

    for (f, r, c) in sums
        FE[r_offsets[f]+r-1, :] += shells2[c, :]
    end

    return EV, FE
end

"""
    planar_arrangement(V::Points, EV::ChainOp, [sigma::Chain], [return_edge_map::Bool], [multiproc::Bool])

Compute the arrangement on the given cellular complex 1-skeleton in 2D.

A cellular complex is arranged when the intersection of every possible pair of cell
of the complex is empty and the union of all the cells is the whole Euclidean space.
The basic method of the function without the `sigma`, `return_edge_map` and `multiproc` arguments
returns the full arranged complex `V`, `EV` and `FE`.

## Additional arguments:
- `sigma::Chain`: if specified, `planar_arrangement` will delete from the output every edge and face outside this cell. Defaults to an empty cell.
- `return_edge_map::Bool`: makes the function return also an `edge_map` which maps the edges of the imput to the one of the output. Defaults to `false`.
- `multiproc::Bool`: Runs the computation in parallel mode. Defaults to `false`.
"""
function planar_arrangement(
        V::LinearAlgebraicRepresentation.Points, copEV::LinearAlgebraicRepresentation.ChainOp,
        sigma::LinearAlgebraicRepresentation.Chain=spzeros(Int8, 0),
        return_edge_map::Bool=false,
        multiproc::Bool=false)

    edgenum = size(copEV, 1)
    edge_map = Array{Array{Int, 1}, 1}(edgenum)
    rV = LinearAlgebraicRepresentation.Points(zeros(0, 2))
    rEV = spzeros(Int8, 0, 0)
    finalcells_num = 0

    if (multiproc == true)
        in_chan = RemoteChannel(()->Channel{Int64}(0))
        out_chan = RemoteChannel(()->Channel{Tuple}(0))

        ordered_dict = SortedDict{Int64,Tuple}()

        @schedule begin
            for i in 1:edgenum
                put!(in_chan,i)
            end
            for p in workers()
                put!(in_chan,-1)
            end
        end

        for p in workers()
            @async Base.remote_do(frag_edge_channel, p, in_chan, out_chan, V, copEV)
        end

        for i in 1:edgenum
            frag_done_job = take!(out_chan)
            ordered_dict[frag_done_job[1]] = frag_done_job[2]
        end

        for (dkey, dval) in ordered_dict
            i = dkey
            v, ev = dval
            newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))

            edge_map[i] = newedges_nums

            finalcells_num += size(ev, 1)

            rV, rEV = LinearAlgebraicRepresentation.skel_merge(rV, rEV, v, ev)
        end

    else
        for i in 1:edgenum
            v, ev = frag_edge(V, copEV, i)

            newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))

            edge_map[i] = newedges_nums

            finalcells_num += size(ev, 1)
            rV, rEV = LinearAlgebraicRepresentation.skel_merge(rV, rEV, v, ev)
        end

    end

    V, copEV = rV, rEV

    V, copEV = merge_vertices!(V, copEV, edge_map)

    # Deletes edges outside sigma area
    if sigma.n > 0
        todel = []

        new_edges = []
        map(i->new_edges=union(new_edges, edge_map[i]), sigma.nzind)
        ev = copEV[new_edges, :]

        for e in 1:copEV.m
            if !(e in new_edges)

                vidxs = copEV[e, :].nzind
                v1, v2 = map(i->V[vidxs[i], :], [1,2])
                centroid = .5*(v1 + v2)

                if !LinearAlgebraicRepresentation.point_in_face(centroid, V, ev)
                    push!(todel, e)
                end
            end
        end

        for i in reverse(todel)
            for row in edge_map

                filter!(x->x!=i, row)

                for j in 1:length(row)
                    if row[j] > i
                        row[j] -= 1
                    end
                end
            end
        end

        V, copEV = LinearAlgebraicRepresentation.delete_edges(todel, V, copEV)
    end

    bicon_comps = biconnected_components(copEV)

    if isempty(bicon_comps)
        println("No biconnected components found.")
        if (return_edge_map)
            return (nothing, nothing, nothing, nothing)
        else
            return (nothing, nothing, nothing)
        end
    end

    edges = sort(union(bicon_comps...))
    todel = sort(setdiff(collect(1:size(copEV,1)), edges))

    for i in reverse(todel)
        for row in edge_map

            filter!(x->x!=i, row)

            for j in 1:length(row)
                if row[j] > i
                    row[j] -= 1
                end
            end
        end
    end

    V, copEV = LinearAlgebraicRepresentation.delete_edges(todel, V, copEV)

    bicon_comps = biconnected_components(copEV)

    n = size(bicon_comps, 1)
    shells = Array{LinearAlgebraicRepresentation.Chain, 1}(n)
    boundaries = Array{LinearAlgebraicRepresentation.ChainOp, 1}(n)
    EVs = Array{LinearAlgebraicRepresentation.ChainOp, 1}(n)
    for p in 1:n
        ev = copEV[sort(bicon_comps[p]), :]
        fe = minimal_2cycles(V, ev)
        shell_num = get_external_cycle(V, ev, fe)

        EVs[p] = ev
        tokeep = setdiff(1:fe.m, shell_num)
        boundaries[p] = fe[tokeep, :]
        shells[p] = fe[shell_num, :]
    end

    shell_bboxes = []
    for i in 1:n
        vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
        push!(shell_bboxes, LinearAlgebraicRepresentation.bbox(V[vs_indexes, :]))
    end

    containment_graph = pre_containment_test(shell_bboxes)
    containment_graph = prune_containment_graph(n, V, EVs, shells, containment_graph)

    transitive_reduction!(containment_graph)

    copEV, FE = cell_merging(n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)

    if (return_edge_map)
        return V, copEV, FE, edge_map
    else
        return V, copEV, FE
    end
end


text:
# ps3000CommonElements.jl

export PS3000TimeUnits, PS3000DriverErrorCodes

baremodule PS3000TimeUnits
  const FS = 0   # femptoseconds
  const PS = 1  # picoseconds
  const NS = 2  # nanoseconds
  const US = 3  # microseconds
  const MS = 4  # miliseconds
  const S = 5   # seconds
end # module

baremodule PS3000DriverErrorCodes
  const OK = 0
  const MAX_UNITS_OPENED = 1
  const MEM_FAIL = 2
  const NOT_FOUND = 3
  const FW_FAIL = 4
  const NOT_RESPONDING = 5
  const CONFIG_FAIL = 6
  const OS_NOT_SUPPORTED =7
end # module


text:
using LowRankModels
using Base.Test

TOL = 1e-5

# QuadConstraint
r = QuadConstraint(7)
@test evaluate(r, ones(7)) == 0
@test evaluate(r, ones(100)) == Inf
@test_approx_eq prox(r, ones(100), 1) ones(100)/sqrt(100)*7


text:
import ..Continuations: find_errors

function find_errors(solver::Codim2Solver)
    if augsys(solver.prob) isa BackReferencingAS
        # TODO: find_errors for BackReferencingAS-problems
        @warn("`find_errors` for problem with `BackReferencingAS`",
             " is not implemented yet.")
        return []
    end
    super = as(solver, ContinuationSolver)
    errors = find_errors(super)
    # TODO: do some Codim2Solver-specific checks
    return errors
end


text:
export PointCloud

struct PointCloud <: AbstractObject
    points::Array{Float32,2}
    npoints::Int
end

function PointCloud(points::Array{T,2} where{T<:Number})
    points =convert(Array{Float32,2},points)
    PointCloud(points, size(points,1))
end

PointCloud(points::Array{Float32,2}) = PointCloud(points, size(points,1))

PointCloud(;points) = PointCloud(points)


text:
function OEMFile()
    return OEMFile(())
end

function add_ephemerides_block(obj::OEMFile)
    return jcall(obj, "addEphemeridesBlock", void, ())
end

function check_time_systems(obj::OEMFile)
    return jcall(obj, "checkTimeSystems", void, ())
end

function get_ephemerides_blocks(obj::OEMFile)
    return jcall(obj, "getEphemeridesBlocks", List, ())
end

function get_satellites(obj::OEMFile)
    return jcall(obj, "getSatellites", Map, ())
end


text:
"""
* Assumes Fx0 is linear polynomial mean basis
* Assume btg testing_buffer has been updated using x0, Fx0, y0, i.e. compute_qmc has been called

NOTE:θ theta needs to be passed as an array
"""
function compute_BO_derivs(train_buffer, train_data, test_buffer, θ, λ, x0, Fx0, y0, m, q, βhat, Σθ_inv_y, cdf_deriv, cdf_second_deriv, cdf_eval)
    x0 = reshape(x0, 1, length(x0))
    #unpack pertinent quantities
    Σθ_inv_X  = btg.train_buffer_dict[θ].Σθ_inv_X
    choleskyΣθ = btg.train_buffer_dict[θ].choleskyΣθ
    choleskyXΣX = btg.train_buffer_dict[θ].choleskyXΣX
    (x, Fx, y, d, n, p) = unpack(btg.trainingData)
    (Eθ, Bθ, ΣθinvBθ, Dθ, Hθ, Cθ) = unpack(btg.test_buffer_dict[θ])
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(compute gradients)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #check if theta is an array of length d
    #println("n (number data points in trainingData): ", n)
    #println("d (number length scales): ", d)
    #println("size of Fx: ", size( Fx))
    #println("size of H: ", size( Hθ))
    #println("size of B: ", size    (Bθ))
    function intermediates(x0)::Tuple{Array{T}, Array{T}, Array{T}, Array{T}} where T<:Real
        @assert typeof(θ)<:Array{T, 1} where T
        @assert length(θ) ==d
        S = repeat(x0, n, 1) .- x # n x d
        #println("size of S: ", size(S))
        jacB =   diagm(Bθ[:]) * S * diagm(- θ[:]) #n x d
        #println("size of jacB: ", size(jacB))
        #println("size of Bθ: ", size(Bθ))
        #println("size of choleskyΣθ: ", size(choleskyΣθ))
        #println("choleskyΣθ inv Bθ': ", size(choleskyΣθ \ Bθ'))
        jacD = -2* jacB' * (choleskyΣθ \ Bθ') #d x 1
        #println("size of jacD: ", size(jacD))
        #assuming linear polynomial basis
        jacFx0 = vcat(zeros(1, d), diagm(ones(d))) #p x d
        #println("size of Fx0: ", size(Fx0))
        jacH = jacFx0' - jacB' * Σθ_inv_X #d x p
        #println("size of jacH: ", size(jacH))
        jacC = jacD + 2 * jacH * (choleskyXΣX \ Hθ') #d x 1
        jacm = jacB' * Σθ_inv_y + jacH * βhat  # d x 1
        return (jacB, jacD, jacFx0, jacC)
    end

    function hessians()
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(compute hessian of D)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ΣinvjacB  = choleskyΣθ\jacB
        W = zeros(d, d)
        Stheta = S * diagm(- θ[:]) # n x d
        for i = 1:d
            for j = 1:d
                arg = Bθ' .* Stheta[:, i:i] .* Stheta[:, j:j]
                if i==j
                    eq = Bθ * (choleskyΣθ \ (-θ[i] .* Bθ'))
                else
                    eq = 0
                end
                W[i, j] = (Bθ * (choleskyΣθ\arg) .+ eq)[1]
            end
        end
        hessD = -2*jacB' * ΣinvjacB - 2*W #d x d
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(compute hessian of C)~~~~~~~~~~~~~~~~~~~~~~~~~~~
        V = zeros(d, d)
        for i = 1:d
            for j = 1:d
                arg = Bθ' .* Stheta[:, i:i] .* Stheta[:, j:j]
                if i==j
                    eq = (-θ[i] .* Bθ')
                else
                    eq = 0
                end
                #println("size arg: ", size(arg))
                #println("size eq: ", size(eq))
                # println(arg .+ eq)
                # println((Fx' * (choleskyΣθ \ (arg .+ eq))))
                # println((choleskyXΣX \ (Fx' * (choleskyΣθ \ (arg .+ eq)))))
                # println(Hθ)
                V[i, j] =  (Hθ * (choleskyXΣX \ ( - Fx' * (choleskyΣθ \ (arg .+ eq)))))[1]
            end
        end
        hessC = hessD + 2*jacH*(choleskyXΣX\jacH') + 2*V #d x d
        #~~~~~~~~~~~~~~~~~~~~~~~~~(compute hessian of m)~~~~~~~~~~~~~~~~~~~~~~~~~~
        hess_m = zeros(d, d)
        for i = 1:d
            for j = 1:d
                arg = Bθ' .* Stheta[:, i:i] .* Stheta[:, j:j]
                if i==j
                    eq = (-θ[i] .* Bθ')
                else
                    eq = 0
                end
                #println()
                #println("sigma inv y: ", size(Σθ_inv_y))
                #println("size arg .+ eq: ", size(arg .+ eq))
                expr1 = dot(Σθ_inv_y,  (arg .+ eq))
                expr2 = dot(( - Fx' * (choleskyΣθ \ (arg .+ eq))), βhat)
                #println("size expr1: ", size(expr1))
                #println("size expr2: ", size(expr2))
                hess_m[i, j] = (expr1 + expr2)[1] #d_i d_j H
            end
        end
    end

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(compute gradient of CDF G(u, s) w.r.t s)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #G(u, s)
    gλy0 = btg.g(y0, λ[1])
    #we need to define covariance (sigma here ...)
    #sigma is covariance sqrt(qC/(n-p))
    qC = (q .* Cθ[1])
    sigma = sqrt(qC/(n-p)) #covariance
    dsigma = sqrt(q/(n-p)) .* 0.5 * Cθ[1]^(-1/2) .* jacC
    hess_sigma = zeros(d, d)
    for i = 1:d
        for j = 1:d
            hess_sigma[i, j] = -0.25 * Cθ[1]^(-3/2) .* jacC[i]jacC[j] + 0.5 * Cθ[1]^(-1/2) * hessC[i, j]
        end
    end

    function Y(i, sigma, dsigma)
        return - jacm[i]/sigma + (gλy0 - m)/sigma^2 * dsigma[i]
    end
    function R(i, j, sigma, dsigma, hess_sigma) #assume jacm  and hess_m defined already
        e1 = -hess_m[i, j]/sigma
        e2 = jacm[j] * dsigma[i] / sigma^2
        e3 = - jacm[i] * dsigma[j] / sigma^2
        e4 = - 2*(gλy0 - m) * dsigma[j]/sigma^3
        e5 = (gλy0 - m)*hess_sigma[i, j] / sigma^2
        return e1 + e2 + e3 + e4 + e5
    end
    G = cdf_eval

    jacG = zeros(1, d) #1 x d
    for i = 1:d
        jacG[i] = (cdf_deriv .* Y(i, sigma, dsigma))[1]
    end

    hess_G = zeros(d, d)
    for i = 1:d
        for j = 1:d
            expr1 = cdf_second_deriv * Y(i, sigma, dsigma) * Y(j, sigma, dsigma)
            expr2 = cdf_deriv * R(i, j, sigma, dsigma, hess_sigma)
            hess_G[i, j] = (expr1 .+ expr2)[1]
        end
    end
    #println("size G: ", size(G))
    #println("size jacG: ", size(jacG))
    #println("size hess_G: ", size(hess_G))

    # println("cdf_eval: ", cdf_eval)
    # println("cdf_deriv: ", cdf_deriv)
    # println("cdf_second_deriv: ", cdf_second_deriv)
    # return (cdf_eval, cdf_deriv, cdf_second_deriv)

    return (G, jacG, hess_G)
    #return (m, jacm', hess_m)
    #return (Bθ*(choleskyΣθ\Bθ'), 2*Bθ*(choleskyΣθ\jacB), 2 * jacB' * (choleskyΣθ\ jacB))
    #return (Dθ, jacD', hessD)
    #return (Cθ, jacC', hessC) #dimensions of jacobian needs to be correct with respect to input and output dimensinos
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module Tests
using StandardizedMatrices, StatsBase, Statistics, LinearAlgebra, Test

n, p, k = 100, 5, 7
x = randn(n, p)
x2 = zscore(x, 1)
z = StandardizedMatrix(x)
b = randn(p)
b2 = randn(n)
b3 = randn(p, k)
b4 = randn(n, k)

@testset "Multiplication with Vector" begin
	@test isapprox(z * b, x2 * b)
	storage = zeros(n)
	mul!(storage, z, b)
	@test isapprox(storage, x2 * b)
	storage = zeros(p)
	mul!(storage, z', b2)
	@test isapprox(storage, x2' * b2)
end

@testset "Multiplication with Matrix" begin
	@test isapprox(z * b3, x2 * b3)
	storage = zeros(n, k)
	mul!(storage, z, b3)
	@test isapprox(storage, x2 * b3)
	storage = zeros(p, k)
	mul!(storage, z', b4)
	@test isapprox(storage, x2' * b4)
end

@testset "Indexing" begin
	for i in eachindex(z)
		@test z[i] == x2[i]
	end
	for j in 1:size(z, 2), i in 1:size(z, 1)
		@test z[i, j] == x2[i, j]
	end
end

end #module


text:
using Pkg
pkg"activate ."
pkg"add JuMP"
pkg"add Cbc"
pkg"add Juniper"
pkg"add Ipopt"

using JuMP, Cbc, Juniper

# Escolhendo o modelo
model = Model(Cbc.Optimizer)

# Definindo os parâmetros que, por sinal, coloquei aleatórios, pois não sei as entradas
m = 4 # linha
n = 4 # coluna
c = rand(1:4, n,m) #matriz
# l,c = size(c)
a = rand(0:0, m)
b = rand(1000:1000,n)
lb = [5,4] #limites inferiores
ob = [6,5] # limites superiores
# Implementando realmente o modelo
@variable(model,w[i = 1:n-1, g = 1:m], Bin) #Variável
@variable(model,x[i = 1:n, g = 1:m], Bin) #Variável
@objective(model, Max, sum(w[i,g] for g = 1:m, i = 1:n-1)) #Func. objetivo

#Restrições
@constraint(model, [i = 1:n, g = 1:m], sum(x[i,g] for i = 1:n-1, g = 1:m) == 1) # 1º restrição
@constraint(model, [i = 1:n, g = 1:m], sum(x[i,g] for i = 1:n-1, g = 1:m ) ≤ b[g]) # 2º restrição
@constraint(model, [i = 1:n, g = 1:m], sum(x[i,g] for i = 1:n-1, g = 1:m ) ≥ a[g]) # 3º restrição
@constraint(model, [i = 1:n - 1, g = 1:m], - sum(max(0,(c[i,j])) * x[i,g] + w[i,g] for j = i+1:n) <= 0) # 4º restrição
@constraint(model, [i = 1:n - 1, g = 1:m], - sum(c[i,j]*x[j,g] for j = i+1:n) +  sum(min(0,c[i,j]) * (1 - x[i,g]) for j = i+1:n)+ w[i,g] <= 0) # 5º restrição

optimize!(model) # Otimizando o modelo

println("Valor da função objetivo = ", objective_value(model))
value.(x)


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

using Markdown
using InteractiveUtils

# ╔═╡ 4723a50e-f6d1-11ea-0b1c-3d33a9b92f87
md"# Defining new types"

# ╔═╡ b3bb0212-f6da-11ea-2bea-21a84ea741a9
md"In this notebook we will see how to define our own types in Julia, and why it's useful to do so."

# ╔═╡ bed2a154-f6d1-11ea-0812-1f5c628a9785
md"## Why define new types?"

# ╔═╡ c826737a-f6d1-11ea-188c-8b6435fcd351
md"Many mathematical and other objects can be represented by a pair of numbers, for example:

  - a rectangle has a width and a height;
  - a complex number has a real and imaginary part;
  - a position vector in 2D space has 2 components.

Each of these could naturally be represented by a pair $(x, y)$ of two numbers, or in Julia as a tuple:
"

# ╔═╡ 81113974-f6d2-11ea-2ef8-fb2930402a74
begin
	rectangle = (3, 4)   # (width, height)
	c = (3, 4)   # 3 + 4im
	x = (3, 4)   # position vector
end

# ╔═╡ ac321d80-f6d2-11ea-2951-676e6e1aef56
md"But from the fact that we have to remind ourselves using comments what each oof these numbers represents, and the fact that they all look the same but should behave very differently, that there is a problem here.

For example, we would like to have a function `width` that returns the width of a rectangle, but it makes no sense to apply that to a complex number.
"

# ╔═╡ 645503ce-f6d7-11ea-37a1-915ed24aefb8
md"In other words, we need a way to be able to distinguish between different *types of* objects with different *behaviours*."

# ╔═╡ 9951a870-f6d7-11ea-1d94-ad757a8b196f
md"## Defining new types"

# ╔═╡ 9f38efc6-f6d7-11ea-088c-c5934f8ed969
md"""Julia allows us to do this by **defining new types**. For example, we would like to define a new **type** of, or *kind* of, object, called a `Rectangle`. Any object of type `Rectangle` should have its own `width` and `height` **fields** that are separate from any other `Rectangle` object, and indeed from any variables lying around that happen to be called that.

We can think of a `Rectangle` object as a bag or box that contains two pieces of data, namely the values of that `Rectangle`'s `width` and `height` fields.

When we specify what a `Rectangle` is, we are actually specifying *what all `Rectangle` objects should look like*, i.e. a template or "cookie-cutter" shape for making new objects of that type.

The syntax in Julia to do this is as follows:
"""

# ╔═╡ 66abc04e-f6d8-11ea-27d8-9f8b14659755
struct Rectangle
	width::Float64
	height::Float64
end

# ╔═╡ 7571be3a-f6d8-11ea-174c-9d65d5185153
md"""The keyword `struct` (short for "structure") tells Julia that we are defining a new type. We list the field names together with **type annotations** using `::` that specify which type each field can contain."""

# ╔═╡ 9f384ac2-f6d8-11ea-297e-4bf09acf9fe7
md"Once Julia has this template, we can create objects which have that type as follows:"

# ╔═╡ b0dac516-f6d8-11ea-1bdb-b59723107206
r = Rectangle(1, 2.5)

# ╔═╡ b98b9faa-f6d8-11ea-3610-bf8a84af2b5a
md"We can check that `r` is now a variable whose type is `Rectangle`, in other words `r` *is a* `Rectangle`:"

# ╔═╡ cf1c4aae-f6d8-11ea-3200-c5fb458c7c09
typeof(r)

# ╔═╡ d0749974-f6d8-11ea-2f41-074b6744f3d5
r isa Rectangle

# ╔═╡ d372342c-f6d8-11ea-10cd-573cf7eab992
md"""We can  extract from `r` the information about the values of the fields that it contains using "`.`": """

# ╔═╡ dd8f9e88-f6d8-11ea-15e3-f17f4af0d81b
r.width

# ╔═╡ e3e70064-f6d8-11ea-22fd-892bbc567ed4
r.height

# ╔═╡ e582eb02-f6d8-11ea-1fcc-89bbc9dfbb07
md"We can create a new `Rectangle` with its own width and height:"

# ╔═╡ f2ed18b4-f6d8-11ea-3bc7-0b82eb5e8dc0
r2 = Rectangle(3.0, 4.0)

# ╔═╡ f9d192fe-f6d8-11ea-138d-3dcdff33c034
md"You should check that this does *not* affect the `width` and `height` variables belonging to `r`."

# ╔═╡ 6085144c-f6db-11ea-19fe-ed46dafb4562
md"Types like this are often called **composite types**; they consist of aggregating, or collecting together, different pieces of information that belong to a given object."

# ╔═╡ 1840898e-f6d9-11ea-3035-bb4dac496834
md"## Mutable vs immutable"

# ╔═╡ 299442a2-f6d9-11ea-1ed0-5fb7ccff2b35
md"Now suppose we want to change the width of `r`. We would naturally try the following:"

# ╔═╡ 63f76d28-f6d9-11ea-071c-458528c36008
r.width = 10

# ╔═╡ 68934a2a-f6d9-11ea-37ea-850304f6d3d6
md"But Julia complains that fields of objects of type `Rectangle` *cannot* be modified. This is because `struct` makes **immutable** (i.e. unchangeable) objects. The reason for this is that these usually lead to *faster code*."

# ╔═╡ a38a0164-f6d9-11ea-2e1d-a7a6e2106d0b
md"If we really want to have a **mutable** (modifiable) object, we can declare it using `mutable struct` instead of `struct` -- try it!"

# ╔═╡ 3fa7d8e2-f6d9-11ea-2e82-9f59b5cb9424
md"## Functions on types"

# ╔═╡ 87eeb3fe-f6da-11ea-2fa3-d95454cd50de
md"We can now define functions that act only on given types. To restrict a given function argument to only accept objects of a certain type, we add a type annotation:"

# ╔═╡ 8c4beb6e-f6db-11ea-12d1-cf450181363b
width(r::Rectangle) = r.width

# ╔═╡ 91e21d28-f6db-11ea-1b0f-336719682f28
width(r)

# ╔═╡ 2ef7392e-f6dc-11ea-00e4-770cdf9a102e
md"Applying the function to objects of other types gives an error:"

# ╔═╡ 9371209e-f6db-11ea-3ba2-c3597d42d8ed
width(3)   # throws an error

# ╔═╡ b916acb4-f6dc-11ea-3cdf-2b8ab3c34e03
md"""It is common in Julia to have "generic" versions of functions that apply to any object, and then specialised versions that apply only to certain types.

For example, we might want to define an `area` function with the correct definition for `Rectangle`s, but it might be convenient to fall back to a version that just returns the value itself for any other type:"""

# ╔═╡ 6fe1b332-f6dd-11ea-39d4-1954aeda6f08
begin
	area(r::Rectangle) = r.width * r.height
	area(x) = x
end

# ╔═╡ c34355b8-f6dd-11ea-1089-cf5be2117ba8
area(r)

# ╔═╡ c4a7fae6-f6dd-11ea-1851-3bd445ebf677
area(17)

# ╔═╡ c83180d8-f6dd-11ea-32f7-634b781070f1
md"But since we didn't restrict the type in the second method, we also have"

# ╔═╡ d3dc64b6-f6dd-11ea-1273-7fb3957e4964
area("hello")

# ╔═╡ 6d7ca590-f6f5-11ea-0a00-6128f971b546
area

# ╔═╡ d6e5a276-f6dd-11ea-34aa-b9e2d3805364
md"which does not make much sense."

# ╔═╡ 79fa562e-f6dd-11ea-2e97-df3c62c83685
md"Note that these are different versions of the function with the *same* function name; each version acting on different (combinations of) types is called a **method** of the (generic) function."

# ╔═╡ a1d2ae6c-f6dd-11ea-0216-ef9db5d9e29b
md"Suppose that later we create a `Circle` type. We can then just add a new method `area(c::Circle)` with the corresponding definition, and Julia will continue to choose the correct version (method) when we call `area(x)`, depending on the type of `x`.

[Note that at the time of writing, Pluto requires all methods of a function to be defined in the same cell.]
"

# ╔═╡ 41e21994-f6de-11ea-0e5c-0515a3a52f6f
md"## Multiple dispatch"

# ╔═╡ 46aed906-f6de-11ea-17d8-cb032f59240e
md"""The act of choosing which function to call based on the type of the arguments that are passed to the function is called **dispatch**. A central feature of Julia is **multiple dispatch**: this choice is made based on the types of *all* the arguments to a function.

For example, the following three calls to the `+` function each call *different* methods, defined in different locations. Click on the links to see the definition of each method in the Julia source code on GitHub!
"""

# ╔═╡ 814e328c-f6de-11ea-13c0-d1b97714c4f3
cc = 3 + 4im

# ╔═╡ cdd3e14e-f6f5-11ea-15e2-bd309e658823
cc + cc

# ╔═╡ e01e26f2-f6f5-11ea-13b0-95413a6f7290
+

# ╔═╡ 84cae75c-f6de-11ea-3cd4-1b263e34771f
@which cc + cc

# ╔═╡ 8ac4904a-f6de-11ea-105b-8925016ca6d5
@which cc + 3

# ╔═╡ 8cd9f438-f6de-11ea-2b58-93bbb860a005
@which 3 + cc

# ╔═╡ Cell order:
# ╟─4723a50e-f6d1-11ea-0b1c-3d33a9b92f87
# ╟─b3bb0212-f6da-11ea-2bea-21a84ea741a9
# ╟─bed2a154-f6d1-11ea-0812-1f5c628a9785
# ╟─c826737a-f6d1-11ea-188c-8b6435fcd351
# ╠═81113974-f6d2-11ea-2ef8-fb2930402a74
# ╟─ac321d80-f6d2-11ea-2951-676e6e1aef56
# ╟─645503ce-f6d7-11ea-37a1-915ed24aefb8
# ╟─9951a870-f6d7-11ea-1d94-ad757a8b196f
# ╟─9f38efc6-f6d7-11ea-088c-c5934f8ed969
# ╠═66abc04e-f6d8-11ea-27d8-9f8b14659755
# ╟─7571be3a-f6d8-11ea-174c-9d65d5185153
# ╟─9f384ac2-f6d8-11ea-297e-4bf09acf9fe7
# ╠═b0dac516-f6d8-11ea-1bdb-b59723107206
# ╟─b98b9faa-f6d8-11ea-3610-bf8a84af2b5a
# ╠═cf1c4aae-f6d8-11ea-3200-c5fb458c7c09
# ╠═d0749974-f6d8-11ea-2f41-074b6744f3d5
# ╟─d372342c-f6d8-11ea-10cd-573cf7eab992
# ╠═dd8f9e88-f6d8-11ea-15e3-f17f4af0d81b
# ╠═e3e70064-f6d8-11ea-22fd-892bbc567ed4
# ╟─e582eb02-f6d8-11ea-1fcc-89bbc9dfbb07
# ╠═f2ed18b4-f6d8-11ea-3bc7-0b82eb5e8dc0
# ╟─f9d192fe-f6d8-11ea-138d-3dcdff33c034
# ╟─6085144c-f6db-11ea-19fe-ed46dafb4562
# ╟─1840898e-f6d9-11ea-3035-bb4dac496834
# ╟─299442a2-f6d9-11ea-1ed0-5fb7ccff2b35
# ╠═63f76d28-f6d9-11ea-071c-458528c36008
# ╟─68934a2a-f6d9-11ea-37ea-850304f6d3d6
# ╟─a38a0164-f6d9-11ea-2e1d-a7a6e2106d0b
# ╟─3fa7d8e2-f6d9-11ea-2e82-9f59b5cb9424
# ╟─87eeb3fe-f6da-11ea-2fa3-d95454cd50de
# ╠═8c4beb6e-f6db-11ea-12d1-cf450181363b
# ╠═91e21d28-f6db-11ea-1b0f-336719682f28
# ╟─2ef7392e-f6dc-11ea-00e4-770cdf9a102e
# ╠═9371209e-f6db-11ea-3ba2-c3597d42d8ed
# ╟─b916acb4-f6dc-11ea-3cdf-2b8ab3c34e03
# ╠═6fe1b332-f6dd-11ea-39d4-1954aeda6f08
# ╠═c34355b8-f6dd-11ea-1089-cf5be2117ba8
# ╠═c4a7fae6-f6dd-11ea-1851-3bd445ebf677
# ╟─c83180d8-f6dd-11ea-32f7-634b781070f1
# ╠═d3dc64b6-f6dd-11ea-1273-7fb3957e4964
# ╠═6d7ca590-f6f5-11ea-0a00-6128f971b546
# ╟─d6e5a276-f6dd-11ea-34aa-b9e2d3805364
# ╟─79fa562e-f6dd-11ea-2e97-df3c62c83685
# ╟─a1d2ae6c-f6dd-11ea-0216-ef9db5d9e29b
# ╟─41e21994-f6de-11ea-0e5c-0515a3a52f6f
# ╟─46aed906-f6de-11ea-17d8-cb032f59240e
# ╠═814e328c-f6de-11ea-13c0-d1b97714c4f3
# ╠═cdd3e14e-f6f5-11ea-15e2-bd309e658823
# ╠═e01e26f2-f6f5-11ea-13b0-95413a6f7290
# ╠═84cae75c-f6de-11ea-3cd4-1b263e34771f
# ╠═8ac4904a-f6de-11ea-105b-8925016ca6d5
# ╠═8cd9f438-f6de-11ea-2b58-93bbb860a005

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
a = 0; b = 1
nₓ = 100
m(x) = x
σ = 1
s = 2
V(x) = 0.
left_bc = DirichletBC(0); right_bc = DirichletBC(0)
f(x) = 1.

problem = RandomODEProblem(a, b, nₓ, m, σ, s, left_bc, right_bc, f, potential = V);
assemble_system!(problem)
solve!(problem)

true_solution(x) = -x/2*(x - 1)

error = abs.(true_solution.(problem.x) - problem.u)
norm(error, Inf) < 1e-14


text:
export DummyTemperatureSensor, DummyTemperatureSensorParams, numChannels, getTemperature

Base.@kwdef struct DummyTemperatureSensorParams <: DeviceParams

end
DummyTemperatureSensorParams(dict::Dict) = params_from_dict(DummyTemperatureSensorParams, dict)

Base.@kwdef mutable struct DummyTemperatureSensor <: TemperatureSensor
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::DummyTemperatureSensorParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}
end

init(sensor::DummyTemperatureSensor) = sensor.present = true

neededDependencies(::DummyTemperatureSensor) = []
optionalDependencies(::DummyTemperatureSensor) = []

Base.close(sensor::DummyTemperatureSensor) = nothing

numChannels(sensor::DummyTemperatureSensor) = 1
getTemperatures(sensor::DummyTemperatureSensor) = [42u"°C"]
getTemperature(sensor::DummyTemperatureSensor, channel::Int)::typeof(1u"°C") = 42u"°C"


text:
export NodeHandle, advertise, advertiseService, subscribe, shutdown, serviceClient, hasParam, deleteParam, setParam, param

# -------------------------------------------------------------------
# Hacky helpers, not intended to be accessed by users
# (Ensuring safety from gc)
# -------------------------------------------------------------------
__subs = []
__srvs = []

# -------------------------------------------------------------------
# roscpp functions
# -------------------------------------------------------------------

# TODO support TransportHints in subscribe function

"""
    NodeHandle()

Create a NodeHandle.
"""
NodeHandle() = @cxxnew ros::NodeHandle()

"""
    advertise(nodehandle, topic_name::String, topic_type, queue_size::Int)

Advertise a topic.
"""
advertise(nodehandle, topic_name::String, topic_type, queue_size::Int) = icxx"$(nodehandle)->advertise<$(topic_type)>($(pointer(topic_name)), $(queue_size));"

"""
    serviceClient(nodehandle, topic_name::String, topic_type)

Create a client for a service.
"""
serviceClient(nodehandle, topic_name::String, topic_type) = icxx"$(nodehandle)->serviceClient<$topic_type>($(pointer(topic_name)));"

"""
    advertiseService(nodehandle, service_name::String, service_request_type, service_response_type, callback)

Advertise a service.
"""
function advertiseService(nodehandle, service_name::String, service_request_type, service_response_type, callback)
    push!(__srvs,
        icxx"""
        boost::function<bool ($(service_request_type)&, $(service_response_type)&)> cpp_srv_callback = [&]($(service_request_type) &req, $(service_response_type) &res) {
            return $:(callback(icxx"return &req;", icxx"return &res;")::Bool);
        };
        return $(nodehandle)->advertiseService($(pointer(service_name)), cpp_srv_callback);
        """)
    return last(__srvs)
end

"""
    subscribe(nodehandle, topic_name::String, queue_size::Int, topic_type, callback)

Subscribe to a topic.
"""
function subscribe(nodehandle, topic_name::String, queue_size::Int, topic_type, callback)
    push!(__subs,
        icxx"""
        boost::function<void ($(topic_type))> cpp_callback = [&]($(topic_type) msg) {
            $:(callback(icxx"return &msg;");nothing);
        };
        return $(nodehandle)->subscribe<$(topic_type)>($(pointer(topic_name)), $(queue_size), cpp_callback);
        """)
    return last(__subs)
end

"""
    hasParam(n, key::String)

Check whether a parameter exists on the parameter server.
"""
hasParam(n, key::String) = @cxx n->hasParam(pointer(key))

"""
    deleteParam(n, key::String)

Delete a parameter from the parameter server.
"""
deleteParam(n, key::String) = @cxx n->deleteParam(pointer(key))

"""
    setParam(n, key::String, value::Int32)

Set a value ont the parameter server.
"""
setParam(n, key::String, value::Int32) = @cxx n->setParam(pointer(key), value)

# The cast below will introduce bugs in VERY large numbers.
"""
    setParam(n, key::String, value::Int64)

Set a value ont the parameter server.
"""
setParam(n, key::String, value::Int64) = @cxx n->setParam(pointer(key), Int32(value))

"""
    setParam(n, key::String, value::Float64)

Set a value ont the parameter server.
"""
setParam(n, key::String, value::Float64) = @cxx n->setParam(pointer(key), value)

"""
    setParam(n, key::String, value::Bool)

Set a value ont the parameter server.
"""
setParam(n, key::String, value::Bool) = @cxx n->setParam(pointer(key), value)

"""
    setParam(n, key::String, value::String)

Set a value ont the parameter server.
"""
setParam(n, key::String, value::String) = @cxx n->setParam(pointer(key), pointer(value))

"""
    param(n, key::String, default::Int32)

Assign value from parameter server, with default.
"""
param(n, key::String, default::Int32) = @cxx n->param(pointer(key), default)

# The cast below will introduce bugs in VERY large numbers.
"""
    param(n, key::String, default::Int64)

Assign value from parameter server, with default.
"""
param(n, key::String, default::Int64) = @cxx n->param(pointer(key), Int32(default))

"""
    param(n, key::String, default::Float64)

Assign value from parameter server, with default.
"""
param(n, key::String, default::Float64) = @cxx n->param(pointer(key), default)

"""
    param(n, key::String, default::Bool)

Assign value from parameter server, with default.
"""
param(n, key::String, default::Bool) = @cxx n->param(pointer(key), default)

"""
    param(n, key::String, default::String)

Assign value from parameter server, with default.
"""
param(n, key::String, default::String) = unsafe_string(icxx"return $(n)->param($(pointer(key)), std::string($(default)));")

"""
    shutdown(n)

Shutdown every handle created through `this` NodeHandle.
"""
shutdown(n::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::NodeHandle")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::NodeHandle")},(false, false, false)},N} where N}) = @cxx n->shutdown()


text:
using KitBase, LinearAlgebra
using KitBase.FastGaussQuadrature, KitBase.JLD2, KitBase.Plots
using ProgressMeter: @showprogress

function flux_wall!(
    ff::T1,
    f::T2,
    u::T3,
    dt,
    rot = 1,
) where {
    T1<:AbstractVector{<:AbstractFloat},
    T2<:AbstractVector{<:AbstractFloat},
    T3<:AbstractVector{<:AbstractFloat},
}
    δ = heaviside.(u .* rot)
    fWall = 0.5 .* δ .+ f .* (1.0 .- δ)
    @. ff = u * fWall * dt

    return nothing
end

begin
    # setup
    set = Setup("radiation", "linesource", "1d1f1v", "kfvs", "bgk", 1, 2, "vanleer", "extra", 0.5, 0.3)

    # physical space
    x0 = 0
    x1 = 1
    nx = 100
    nxg = 0
    ps = PSpace1D(x0, x1, nx, nxg)

    # velocity space
    nu = 100
    points, weights = gausslegendre(nu)
    vs = VSpace1D(points[1], points[end], nu, points, ones(nu) .* (points[end] - points[1]) / (nu - 1), weights)

    # material
    σs = ones(Float32, nx)
    σa = zeros(Float32, nx)
    σt = σs + σa
    σq = zeros(Float32, nx)

    # time
    dt = set.cfl * ps.dx[1]
    nt = set.maxTime / dt |> floor |> Int
    global t = 0.0

    # solution
    f0 = 0.0001 * ones(nu)
    phi = zeros(nu, nx)
    for i = 1:nx
        phi[:, i] .= f0
    end

    flux = zeros(nu, nx+1)

    cd(@__DIR__)
end

for iter = 1:nt
    println("iteration $iter of $nt")

    ff = @view flux[:, 1]
    flux_wall!(ff, phi[:, 1], points, dt, 1.0)

    @inbounds for i = 2:nx
        ff = @view flux[:, i]
        KitBase.flux_kfvs!(ff, phi[:, i-1], phi[:, i], points, dt)
    end

    @inbounds for i = 1:nx-1
        for q = 1:nu
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / ps.dx[i] +
                (σs[i] * phi[q, i] - σt[i] * phi[q, i]) * dt +
                σq[i] * dt
        end
    end
    phi[:, nx] .=  phi[:, nx-1]

    global t += dt
end

ρ = [sum(phi[:, i] .* weights) for i = 1:nx]
plot(ps.x[1:nx], ρ)

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
"""
A basis of the classicale Laguerre polynomials. These polynomials are orthogonal
on the positive halfline `[0,∞)` with respect to the weight function
`w(x)=exp(-x)`.
"""
struct Laguerre{T} <: OPS{T}
    n       ::  Int
    α       ::  T

	Laguerre{T}(n::Int, α = 0) where {T} = new{T}(n, α)
end

Laguerre(n::Int; α = 0) = Laguerre(n, α)
Laguerre(n::Int, α::T) where {T <: AbstractFloat} = Laguerre{T}(n, α)
Laguerre(n::Int, α::S) where {S} = Laguerre(n, float(α))

similar(b::Laguerre, ::Type{T}, n::Int) where {T} = Laguerre{T}(n, b.α)

support(b::Laguerre{T}) where {T} = HalfLine{T}()

first_moment(b::Laguerre{T}) where {T} = gamma(b.α+1)

laguerre_α(b::Laguerre) = b.α

measure(b::Laguerre) = LaguerreWeight(b.α)

iscompatible(d1::Laguerre, d2::Laguerre) = d1.α == d2.α
interpolation_grid(dict::Laguerre) = LaguerreNodes(length(dict), dict.α)
iscompatible(dict::Laguerre, grid::LaguerreNodes) = length(dict) == length(grid) && dict.α ≈ grid.α
isorthogonal(dict::Laguerre, measure::GaussLaguerre) = laguerre_α(dict) ≈ laguerre_α(measure) && opsorthogonal(dict, measure)

isorthonormal(dict::Laguerre, measure::LaguerreWeight) = isorthogonal(dict, measure) && laguerre_α(dict) == 0
isorthonormal(dict::Laguerre, measure::GaussLaguerre) where T = isorthogonal(dict, measure) && laguerre_α(dict) == 0
issymmetric(::Laguerre) = false

isorthogonal(dict::Laguerre, measure::LaguerreWeight) = dict.α == measure.α

gauss_rule(dict::Laguerre) = GaussLaguerre(length(dict), dict.α)

function innerproduct_native(d1::Laguerre, i::PolynomialDegree, d2::Laguerre, j::PolynomialDegree, measure::LaguerreWeight; options...)
	T = coefficienttype(d1)
	if iscompatible(d1, d2) && isorthogonal(d1, measure)
		if i == j
			gamma(convert(T, value(i)+1+laguerre_α(d1))) / convert(T, factorial(value(i)))
		else
			zero(T)
		end
	else
		innerproduct1(d1, i, d2, j, measure; options...)
	end
end

# See DLMF, Table 18.9.1
# http://dlmf.nist.gov/18.9#i
rec_An(b::Laguerre{T}, n::Int) where {T} = -T(1) / T(n+1)

rec_Bn(b::Laguerre{T}, n::Int) where {T} = T(2*n + b.α + 1) / T(n+1)

rec_Cn(b::Laguerre{T}, n::Int) where {T} = T(n + b.α) / T(n+1)

## Printing

function show(io::IO, b::Laguerre{Float64})
	if laguerre_α(b) == 0
		print(io, "Laguerre($(length(b)))")
	else
		print(io, "Laguerre($(length(b)); α = $(laguerre_α(b)))")
	end
end

function show(io::IO, b::Laguerre{T}) where {T}
	if laguerre_α(b) == 0
		print(io, "Laguerre{$(T)}($(length(b)))")
	else
		print(io, "Laguerre{$(T)}($(length(b)); α = $(laguerre_α(b)))")
	end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Direct Inversion of the Iterative Subspace
#
# When solving systems of linear (or nonlinear) equations, iterative methods are often employed.  Unfortunately, such methods often suffer from convergence issues such as numerical instability, slow convergence, and significant computational expense when applied to difficult problems.  In these cases, convergence accelleration methods may be applied to both speed up, stabilize and/or reduce the cost for the convergence patterns of these methods, so that solving such problems become computationally tractable.  One such method is known as the direct inversion of the iterative subspace (DIIS) method, which is commonly applied to address convergence issues within self consistent field computations in Hartree-Fock theory (and other iterative electronic structure methods).  In this tutorial, we'll introduce the theory of DIIS for a general iterative procedure, before integrating DIIS into our previous implementation of RHF.
#
# ## I. Theory
#
# DIIS is a widely applicable convergence acceleration method, which is applicable to numerous problems in linear algebra and the computational sciences, as well as quantum chemistry in particular.  Therefore, we will introduce the theory of this method in the general sense, before seeking to apply it to SCF.
#
# Suppose that for a given problem, there exist a set of trial vectors $\{\mid{\bf p}_i\,\rangle\}$ which have been generated iteratively, converging toward the true solution, $\mid{\bf p}^f\,\rangle$.  Then the true solution can be approximately constructed as a linear combination of the trial vectors,
# $$\mid{\bf p}\,\rangle = \sum_ic_i\mid{\bf p}_i\,\rangle,$$
# where we require that the residual vector
# $$\mid{\bf r}\,\rangle = \sum_ic_i\mid{\bf r}_i\,\rangle\,;\;\;\; \mid{\bf r}_i\,\rangle
# =\, \mid{\bf p}_{i+1}\,\rangle - \mid{\bf p}_i\,\rangle$$
# is a least-squares approximate to the zero vector, according to the constraint
# $$\sum_i c_i = 1.$$
# This constraint on the expansion coefficients can be seen by noting that each trial function ${\bf p}_i$ may be represented as an error vector applied to the true solution, $\mid{\bf p}^f\,\rangle + \mid{\bf e}_i\,\rangle$.  Then
# \begin{align}
# \mid{\bf p}\,\rangle &= \sum_ic_i\mid{\bf p}_i\,\rangle\\
# &= \sum_i c_i(\mid{\bf p}^f\,\rangle + \mid{\bf e}_i\,\rangle)\\
# &= \mid{\bf p}^f\,\rangle\sum_i c_i + \sum_i c_i\mid{\bf e}_i\,\rangle
# \end{align}
# Convergence results in a minimization of the error (causing the second term to vanish); for the DIIS solution vector $\mid{\bf p}\,\rangle$ and the true solution vector $\mid{\bf p}^f\,\rangle$ to be equal, it must be that $\sum_i c_i = 1$.  We satisfy our condition for the residual vector by minimizing its norm,
# $$\langle\,{\bf r}\mid{\bf r}\,\rangle = \sum_{ij} c_i^* c_j \langle\,{\bf r}_i\mid{\bf r}_j\,\rangle,$$
# using Lagrange's method of undetermined coefficients subject to the constraint on $\{c_i\}$:
# $${\cal L} = {\bf c}^{\dagger}{\bf Bc} - \lambda\left(1 - \sum_i c_i\right)$$
# where $B_{ij} = \langle {\bf r}_i\mid {\bf r}_j\rangle$ is the matrix of residual vector overlaps.  Minimization of the Lagrangian with respect to the coefficient $c_k$ yields (for real values)
# \begin{align}
# \frac{\partial{\cal L}}{\partial c_k} = 0 &= \sum_j c_jB_{jk} + \sum_i c_iB_{ik} - \lambda\\
# &= 2\sum_ic_iB_{ik} - \lambda
# \end{align}
# which has matrix representation
# \begin{equation}
# \begin{pmatrix}
#   B_{11} & B_{12} & \cdots & B_{1m} & -1 \\
#   B_{21} & B_{22} & \cdots & B_{2m} & -1 \\
#   \vdots  & \vdots  & \ddots & \vdots  & \vdots \\
#   B_{n1} & B_{n2} & \cdots & B_{nm} & -1 \\
#   -1 & -1 & \cdots & -1 & 0
# \end{pmatrix}
# \begin{pmatrix}
# c_1\\
# c_2\\
# \vdots \\
# c_n\\
# \lambda
# \end{pmatrix}
# =
# \begin{pmatrix}
# 0\\
# 0\\
# \vdots\\
# 0\\
# -1
# \end{pmatrix},
# \end{equation}
#
# which we will refer to as the Pulay equation, named after the inventor of DIIS.  It is worth noting at this point that our trial vectors, residual vectors, and solution vector may in fact be tensors of arbitrary rank; it is for this reason that we have used the generic notation of Dirac in the above discussion to denote the inner product between such objects.
#
# ## II. Algorithms for DIIS
# The general DIIS procedure, as described above, has the following structure during each iteration:
# #### Algorithm 1: Generic DIIS procedure
# 1. Compute new trial vector, $\mid{\bf p}_{i+1}\,\rangle$, append to list of trial vectors
# 2. Compute new residual vector, $\mid{\bf r}_{i+1}\,\rangle$, append to list of trial vectors
# 3. Check convergence criteria
#     - If RMSD of $\mid{\bf r}_{i+1}\,\rangle$ sufficiently small, and
#     - If change in DIIS solution vector $\mid{\bf p}\,\rangle$ sufficiently small, break
# 4. Build **B** matrix from previous residual vectors
# 5. Solve Pulay equation for coefficients $\{c_i\}$
# 6. Compute DIIS solution vector $\mid{\bf p}\,\rangle$
#
# For SCF iteration, the most common choice of trial vector is the Fock matrix **F**; this choice has the advantage over other potential choices (e.g., the density matrix **D**) of **F** not being idempotent, so that it may benefit from extrapolation.  The residual vector is commonly chosen to be the orbital gradient in the AO basis,
# $$g_{\mu\nu} = ({\bf FDS} - {\bf SDF})_{\mu\nu},$$
# however the better choice (which we will make in our implementation!) is to orthogonormalize the basis of the gradient with the inverse overlap metric ${\bf A} = {\bf S}^{-1/2}$:
# $$r_{\mu\nu} = ({\bf A}^{\rm T}({\bf FDS} - {\bf SDF}){\bf A})_{\mu\nu}.$$
# Therefore, the SCF-specific DIIS procedure (integrated into the SCF iteration algorithm) will be:
# #### Algorithm 2: DIIS within an SCF Iteration
# 1. Compute **F**, append to list of previous trial vectors
# 2. Compute AO orbital gradient **r**, append to list of previous residual vectors
# 3. Compute RHF energy
# 3. Check convergence criteria
#     - If RMSD of **r** sufficiently small, and
#     - If change in SCF energy sufficiently small, break
# 4. Build **B** matrix from previous AO gradient vectors
# 5. Solve Pulay equation for coefficients $\{c_i\}$
# 6. Compute DIIS solution vector **F_DIIS** from $\{c_i\}$ and previous trial vectors
# 7. Compute new orbital guess with **F_DIIS**

# ## III. Implementation
#
# In order to implement DIIS, we're going to integrate it into an existing RHF program.  Since we just-so-happened to write such a program in the last tutorial, let's re-use the part of the code before the SCF integration which won't change when we include DIIS:

# +
# ==> Basic Setup <==
# Import statements
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using TensorOperations: @tensor
using LinearAlgebra: Diagonal, Hermitian, eigen, tr, norm, dot
using Printf: @printf

# Memory specification
psi4.set_memory(Int(5e8))
numpy_memory = 2

# Set output file
psi4.core.set_output_file("output.dat", false)

# Define Physicist's water -- don't forget C1 symmetry!
mol = psi4.geometry("""
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

# Set computation options
psi4.set_options(Dict("basis" => "cc-pvdz",
                      "scf_type" => "pk",
                      "e_convergence" => 1e-8))

# Maximum SCF iterations
MAXITER = 40
# Energy convergence criterion
E_conv = 1.0e-6
D_conv = 1.0e-3

# +
# ==> Static 1e- & 2e- Properties <==
# Class instantiation
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option("basis"))
mints = psi4.core.MintsHelper(wfn.basisset())

# Overlap matrix
S = np.asarray(mints.ao_overlap()) # we only need a copy

# Number of basis Functions & doubly occupied orbitals
nbf = size(S)[1]
ndocc = wfn.nalpha()

println("Number of occupied orbitals: ", ndocc)
println("Number of basis functions: ", nbf)

# Memory check for ERI tensor
I_size = nbf^4 * 8.e-9
println("\nSize of the ERI tensor will be $I_size GB.")
memory_footprint = I_size * 1.5
if I_size > numpy_memory
    psi4.core.clean()
    throw(OutOfMemoryError("Estimated memory utilization ($memory_footprint GB) exceeds " *
                           "allotted memory limit of $numpy_memory GB."))
end

# Build ERI Tensor
I = np.asarray(mints.ao_eri()) # we only need a copy

# Build core Hamiltonian
T = np.asarray(mints.ao_kinetic()) # we only need a copy
V = np.asarray(mints.ao_potential()) # we only need a copy
H = T + V;

# +
# ==> CORE Guess <==
# AO Orthogonalization Matrix
A = mints.ao_overlap()
A.power(-0.5, 1.e-16) # ≈ Julia's A^(-0.5) after psi4view()
A = np.asarray(A)

# Transformed Fock matrix
F_p = A * H * A

# Diagonalize F_p for eigenvalues & eigenvectors with Julia
e, C_p = eigen(Hermitian(F_p))

# Transform C_p back into AO basis
C = A * C_p

# Grab occupied orbitals
C_occ = C[:, 1:ndocc]

# Build density matrix from occupied orbitals
D = C_occ * C_occ'

# Nuclear Repulsion Energy
E_nuc = mol.nuclear_repulsion_energy()
# -

# Now let's put DIIS into action.  Before our iterations begin, we'll need to create empty lists to hold our previous residual vectors (AO orbital gradients) and trial vectors (previous Fock matrices), along with setting starting values for our SCF energy and previous energy:

# ==> Pre-Iteration Setup <==
# SCF & Previous Energy
SCF_E = 0.0
E_old = 0.0;

# Now we're ready to write our SCF iterations according to Algorithm 2.  Here are some hints which may help you along the way:
#
# #### Starting DIIS
# Since DIIS builds the approximate solution vector $\mid{\bf p}\,\rangle$ as a linear combination of the previous trial vectors $\{\mid{\bf p}_i\,\rangle\}$, there's no need to perform DIIS on the first SCF iteration, since there's only one trial vector for DIIS to use!
#
# #### Building **B**
# 1. The **B** matrix in the Lagrange equation is really $\tilde{\bf B} = \begin{pmatrix} {\bf B} & -1\\ -1 & 0\end{pmatrix}$.
# 2. Since **B** is the matrix of residual overlaps, it will be a square matrix of dimension equal to the number of residual vectors.  If **B** is an $N\times N$ matrix, how big is $\tilde{\bf B}$?
# 3. Since our residuals are real, **B** will be a symmetric matrix.
# 4. To build $\tilde{\bf B}$, make an empty array of the appropriate dimension, then use array indexing to set the values of the elements.
#
# #### Solving the Pulay equation
# 1. Use built-in Julia functionality to make your life easier.
# 2. The solution vector for the Pulay equation is $\tilde{\bf c} = \begin{pmatrix} {\bf c}\\ \lambda\end{pmatrix}$, where $\lambda$ is the Lagrange multiplier, and the right hand side is $\begin{pmatrix} {\bf 0}\\ -1\end{pmatrix}$.

# +
# Start from fresh orbitals
F_p = A * H * A
e, C_p = eigen(Hermitian(F_p))
C = A * C_p
C_occ = C[:, 1:ndocc]
D = C_occ * C_occ' ;

# Trial & Residual Vector Lists
F_list = []
DIIS_RESID = []

# ==> SCF Iterations w/ DIIS <==
println("==> Starting SCF Iterations <==")
SCF_E = let SCF_E = SCF_E, E_old = E_old, D = D

   # Begin Iterations
   for scf_iter in 1:MAXITER
      # Build Fock matrix
      @tensor G[p,q] := (2I[p,q,r,s] - I[p,r,q,s]) * D[r,s]
      F = H + G

      # Build DIIS Residual
      diis_r = A * (F * D * S - S * D * F) * A

      # Append trial & residual vectors to lists
      push!(F_list, F)
      push!(DIIS_RESID, diis_r)

      # Compute RHF energy
      SCF_E = tr((H + F) * D) + E_nuc
      dE = SCF_E - E_old
      dRMS = norm(diis_r)
      @printf("SCF Iteration %3d: Energy = %4.16f dE = %1.5e dRMS = %1.5e \n",
                          scf_iter,        SCF_E, SCF_E - E_old,     dRMS)

      # SCF Converged?
      if abs(SCF_E - E_old) < E_conv && dRMS < D_conv
          break
      end
      E_old = SCF_E

      if scf_iter >= 2
         # Build B matrix
         B_dim = length(F_list) + 1
         B = zeros(B_dim, B_dim)
         B[end,   :] .= -1
         B[:  , end] .= -1
         B[end, end]  =  0
         for i in eachindex(F_list), j in eachindex(F_list)
            B[i, j] = dot(DIIS_RESID[i],  DIIS_RESID[j])
         end

         # Build RHS of Pulay equation
         rhs = zeros(B_dim)
         rhs[end] = -1

         # Solve Pulay equation for c_i's with Julia
         coeff = B \ rhs

         # Build DIIS Fock matrix
         F = zeros(size(F))
         for i in 1:length(coeff) - 1
            F += coeff[i] * F_list[i]
         end
      end

      # Compute new orbital guess with DIIS Fock matrix
      F_p =  A * F * A
      e, C_p = eigen(Hermitian(F_p))
      C = A * C_p
      C_occ = C[:, 1:ndocc]
      D = C_occ * C_occ'

      # MAXITER exceeded?
      if scf_iter == MAXITER
          psi4.core.clean()
          throw(MethodError("Maximum number of SCF iterations exceeded."))
      end
   end
   SCF_E
end

# Post iterations
println("\nSCF converged.")
println("Final RHF Energy: $SCF_E [Eh]")
# -

# Congratulations! You've written your very own Restricted Hartree-Fock program with DIIS convergence accelleration!  Finally, let's check your final RHF energy against  Psi4:

# Compare to Psi4
SCF_E_psi = psi4.energy("SCF")
psi4.compare_values(SCF_E_psi, SCF_E, 6, "SCF Energy")

# ## References
# 1. P. Pulay. *Chem. Phys. Lett.* **73**, 393-398 (1980)
# 2. C. David Sherrill. *"Some comments on accellerating convergence of iterative sequences using direct inversion of the iterative subspace (DIIS)".* Available at: vergil.chemistry.gatech.edu/notes/diis/diis.pdf. (1998)


text:
function print_trial(ex)
    data = ex.trial_data
    println("Trial info: ")
    println("\t  target.................... ", data.target_present)
    println("\t  stimulus.................. ", Crayon(foreground=data.target_color), data.target_shape)
    print(Crayon(foreground=:default))
end

function print_visual_buffer(actr)
    object = actr.visual.buffer[1]
    println("Visual Buffer")
    println("\t  stimulus.................. ", Crayon(foreground=object.features.color.value),
        object.features.shape.value)
    print(Crayon(foreground=:default))
    println("\t  location.................. ", "x: ", round(Int, object.location[1]),
        " y: ", round(Int, object.location[2]))
    println("\t  activation................ ", round(object.activation, digits=2))
    println("\t  attend time .............. ", round.(object.attend_time, digits=3))

end

function get_time(actr)
    return round(actr.time, digits=3)
end

function print_response(actr, response)
    println(get_time(actr), " Motor..................... respond $response")
end

function print_abstract_location(actr, status)
    println(get_time(actr), " Abstract-Location......... $status")
    println("\t  iconic memory size........ ", iconic_memory_size(actr), " out of ", visicon_size(actr))
    println("\t  termination probability... ", "≈ ", round(termination_prob(actr), digits=2))
    status != "object found" ? (return) : nothing
    result = actr.visual_location.buffer[1]
    angular_distance = compute_angular_distance(actr, result)
    angular_size = result.angular_size
    parms = actr.parms.acuity[:color]
    color_threshold = compute_acuity_threshold(parms, angular_distance)
    parms = actr.parms.acuity[:shape]
    shape_threshold = compute_acuity_threshold(parms, angular_distance)
    println("\t  model focus............... ", "x: ", round(Int, actr.visual.focus[1]),
        " y: ", round(Int, actr.visual.focus[2]))
    println("\t  color threshold........... ", round.(color_threshold, digits=2), "°")
    println("\t  shape threshold........... ", round.(shape_threshold, digits=2), "°")
    println("\t  angular size.............. ", round.(angular_size, digits=2), "°")
    println("\t  angular distance.......... ", round.(angular_distance, digits=2), "°")
end

function print_vision(actr)
    println(get_time(actr), " Vision","."^20, " object attended.")
end

print_check(actr, v) = println(get_time(actr), " Procedural","."^16, " target $v")

iconic_memory_size(actr) = sum(x->x.visible, get_iconic_memory(actr))

visicon_size(actr) = length(get_iconic_memory(actr))

function termination_prob(actr)
    σ = actr.parms.noise
    τ = actr.parms.τₐ
    vos = filter(x->relevant_object(actr, x), get_iconic_memory(actr))
    α = map(x->weighted_activation(actr, x), vos)
    push!(α, τ)
    return exp(τ/σ)/sum(exp.(α/σ))
end


text:
import Pkg
Pkg.update()
Pkg.add("CSV")
Pkg.add("ColorBrewer")
Pkg.add("Colors")
Pkg.add("DSP")
Pkg.add("DataFrames")
Pkg.add("EEG")
Pkg.add("FileIO")
Pkg.add("Flux")
Pkg.add("GitHub")
Pkg.add("Gumbo")
Pkg.add("IJulia")
Pkg.add("ImageMagick")
Pkg.add("ImageView")
Pkg.add("Images")
Pkg.add("LaTeX")
Pkg.add("LaTeXStrings")
Pkg.add("LightGraphs")
Pkg.add("MAT")
Pkg.add("Plots")
Pkg.add("Pluto")
Pkg.add("PyCall")
Pkg.add("PyPlot")
# Pkg.add("PySide")
# Pkg.add("RCall")
# Pkg.add("RDatasets")
# Pkg.add("ScikitLearn")
# Pkg.add("ScikitLearnBase")
# Pkg.add("StatsBase")
# Pkg.add("StatsFuns")
Pkg.add("SymPy")
# Pkg.add("Taro")
# Pkg.add("TextAnalysis")

max_stars_repo_name: abcsds/dotfiles

max_issues_repo_name: abcsds/dotfiles

max_forks_repo_name: abcsds/dotfiles


text:
export CVNormal, BoundedWienerProcess

struct CVNormal{CV,T} <: ContinuousUnivariateDistribution
    normal::Normal{T}
    CVNormal{CV}(n::Normal{T}) where {CV, T} = new{CV,T}(n)
end
CVNormal{CV}(μ::Number) where {CV} = CVNormal{CV}(Normal(μ,CV*μ))
CVNormal{CV}() where {CV} = CVNormal{CV}(0)

MacroTools.@forward CVNormal.normal Distributions.mean, Base.rand, Distributions.sampler, Distributions.pdf, Distributions.logpdf, Distributions.cdf, Distributions.quantile, Distributions.minimum,
    Distributions.maximum, Distributions.insupport, Distributions.var, Distributions.std, Distributions.modes, Distributions.mode, Distributions.skewness, Distributions.kurtosis,
    Distributions.entropy, Distributions.mgf, Distributions.cf, Base.eltype
Distributions.params(d::CVNormal) = mean(d)
Base.rand(i::InventoryModels.CVNormal, v::Vararg{Int64, N} where N) = rand(i.normal, v...)

cv(d::CVNormal{CV}) where CV = CV
cv(d::Type{CVNormal{CV}}) where CV = CV

#Base.show(io::IO, n::CVNormal) = show(n.normal)

####BoundedWienerProcess
mutable struct BoundedWienerProcess{T, D <: ContinuousUnivariateDistribution} <: ContinuousUnivariateDistribution
    x::T
    d::D
    lb::T
    ub::T
end

BoundedWienerProcess(d::Distribution, lb, ub) = BoundedWienerProcess(rand(Uniform(lb,ub)), d, eltype(d)(lb), eltype(d)(ub))

function Base.rand(bw::BoundedWienerProcess)
    delta = rand(bw.d)
    bw.x = max(bw.lb, min(bw.ub, bw.x + delta))
    bw.x
end

reset!(b::BoundedWienerProcess) = b.x = rand(Uniform(b.lb,b.ub))


text:
A=[1, 2, 3]
S = Set([1, 2])
D = Dict("a"=>1, "b"=>2, "c"=>3);

# pop!
println(pop!(A))
println(A)

println(pop!(S))
println(S)

println(pop!(D))
println(D)

# popfirst!
println(popfirst!(A))

# pushfirst!
println(pushfirst!([1, 2, 3, 4], 5, 6))

# insert!(a::Vector, index::Integer, item)
println(insert!([6, 5, 2, 1], 3, 3))

# deleteat!(a::Vector, i::Integer)
println(deleteat!([6, 5, 4, 3, 2, 1], 2))

# splice!
# Remove the item at the given index, and return the removed item.
println(splice!([6, 5, 4, 3, 2, 1], 5))

# resize!
println(resize!([6, 5, 4, 3, 2, 1], 3))

# append!
println(append!([1, 2], [5, 6]))

# preappend!
println(prepend!([3],[1,2]))

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia


text:
module TestCollections

using Base.Test
using RLEVectors

@testset begin

x = RLEVector([4,5,6],[3,6,9])
# setdiff, symdiff, union, endof, maxabs, minabs, any, all, in

# push!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
push!(x,9)
@test x == RLEVector([1,1,2,2,3,3,4,4,5,5,9])
x = RLEVector(Int64[], Int64[])
push!(x, 12)
@test x == RLEVector([12])

# pop!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test pop!(x) == 5
@test x == RLEVector([1,1,2,2,3,3,4,4,5])

# shove!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
shove!(x,4)
@test x == RLEVector([4,1,1,2,2,3,3,4,4,5,5])

# shift!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test shift!(x) == 1
@test x == RLEVector([1,2,2,3,3,4,4,5,5])

# vcat
@test RLEVector([1,1,2,2,3,3,4,4]) == vcat( RLEVector([1,1,2,2]), RLEVector([3,3,4,4]) )
@test RLEVector([1,1,2,2,2,2,4,4]) == vcat( RLEVector([1,1,2,2]), RLEVector([2,2,4,4]) )

# insert!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test_throws BoundsError insert!(x,0,5)
@test_throws BoundsError insert!(x,length(x) + 2,5)
@test insert!(x,1,9) == RLEVector([9,1,1,2,2,3,3,4,4,5,5])
@test insert!(x,3,9) == RLEVector([9,1,9,1,2,2,3,3,4,4,5,5])
@test insert!(x,1,7) == RLEVector([7,9,1,9,1,2,2,3,3,4,4,5,5])
@test insert!(x,length(x) + 1,100) == RLEVector([7,9,1,9,1,2,2,3,3,4,4,5,5,100])

# deleterun!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test deleterun!(x,1) == RLEVector([2,2,3,3,4,4,5,5])
@test deleterun!(x,4) == RLEVector([2,2,3,3,4,4])
@test deleterun!(x,2) == RLEVector([2,2,4,4])
x = RLEVector([1,1,2,2,1,1,4,4,5,5])
@test deleterun!(x,2) == RLEVector([1,4,5],[4,6,8])

# deleteat! and decrement_run!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test deleteat!(x,1) == RLEVector([1,2,2,3,3,4,4,5,5])
@test deleteat!(x,1) == RLEVector([2,2,3,3,4,4,5,5])
@test deleteat!(x,4) == RLEVector([2,2,3,4,4,5,5])
@test deleteat!(x,5) == RLEVector([2,2,3,4,5,5])
@test deleteat!(x,6) == RLEVector([2,2,3,4,5])

# splice!
# splice! removing
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,1) == 1
@test x == RLEVector([1,2,2,3,3,4,4,5,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,10) == 5
@test x == RLEVector([1,1,2,2,3,3,4,4,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,3) == 2
@test x == RLEVector([1,1,2,3,3,4,4,5,5])

# splice! replacing
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,1,[9]) == 1
@test x == RLEVector([9,1,2,2,3,3,4,4,5,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,1,[9,11]) == 1
@test x == RLEVector([9,11,1,2,2,3,3,4,4,5,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,5,[10,11,12]) == 3
@test x == RLEVector([1,1,2,2,10,11,12,3,4,4,5,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,5,RLEVector([10,11,12])) == 3
@test x == RLEVector([1,1,2,2,10,11,12,3,4,4,5,5])
x = RLEVector(collect(1:10))
@test splice!(x,10,[100]) == 10
@test x == RLEVector([collect(1:9);100])

# splice! adding
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,1:0,[9,11]) == similar(x,0)
@test x == RLEVector([9,11,1,1,2,2,3,3,4,4,5,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,5:4,[10,11,12]) == similar(x,0)
@test x == RLEVector([1,1,2,2,10,11,12,3,3,4,4,5,5])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,10,[9]) == 5
@test x == RLEVector([1,1,2,2,3,3,4,4,5,9])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,10,[9,11]) == 5
@test x == RLEVector([1,1,2,2,3,3,4,4,5,9,11])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test splice!(x,10:9,[9,11]) == similar(x,0)
@test x == RLEVector([1,1,2,2,3,3,4,4,5,9,11,5])

# resize!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test resize!(x,5) == RLEVector([1,1,2,2,3])
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test resize!(x,12) == RLEVector([1,1,2,2,3,3,4,4,5,5,0,0])

# empty!
x = RLEVector([1,1,2,2,3,3,4,4,5,5])
@test empty!(x) == RLEVector(Int[],Int[])

# intersect
x = RLEVector([9, 9, 12, 12, 13, 13, 9, 9, 14, 14])
@test intersect(x, [12, 9]) == RLEVector([9, 9, 12, 12, 9, 9])
@test intersect(x, [9, 12], [12]) == RLEVector([12, 12])
@test intersect(x, Set([12, 14])) == RLEVector([12, 12, 14, 14])

end # testset

end # module


text:
"""
# Method p41

One dimensional analysis of an axially loaded elastic Rod using 2-node
Line elements.

### Constructors
```julia
p41(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}`  : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Type of  structural fin_el
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Float64}                        : x-coordinate vector
```

### Optional additional data dictionary keys
```julia
* penalty = 1e20               : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}         : Element material vector if np_types > 1
* eq_nodal_forces_and_moments  : Contribution of distributed loads to loaded_nodes
```

### Return values
```julia
* (jfem, dis_df, fm_df)        : Tuple of jFem, dis_df and fm_df
                                 where:
                                    jfem::jFem    : Computational result type
                                    dis_df        : Displacement data table
                                    fm_df         : Forces and moments data table
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Rod                           : Help on a Rod structural element
?FiniteElement                 : List finite element types
?Line                          : Help on Line finite element
```
"""
function p41(data::Dict{Symbol, Any})

  # Parse & check FE problem data input dict

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    throw("No structural element type specified.")
  end

  if typeof(struc_el) == PtFEM.Rod
    ndim = 1
    nst = struc_el.np_types
  else
    throw("p41 expects a Rod structural element.")
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = PtFEM.mesh_size(fin_el, struc_el.nxe)
  else
    throw("p41 expects a Line (Interval) finite element.")
  end

  nodof = fin_el.nodof         # Degrees of freedom per node
  ndof = fin_el.nod * nodof    # Degrees of freedom per fin_el

  # Set penalty

  penalty = :penalty in keys(data) ? data[:penalty] : 1e20

  # All dynamic arrays

  points = zeros(struc_el.nip, ndim)   #
  g = zeros(Int, ndof)                   # Element steering vector
  g_coord = zeros(ndim,nn)                 #
  fun = zeros(fin_el.nod)                 #
  coord = zeros(fin_el.nod, ndim)         #
  gamma = zeros(nels)                      #
  jac = zeros(ndim, ndim)                  #
  g_num = zeros(Int, fin_el.nod, nels)  #
  der = zeros(ndim, fin_el.nod)           #
  deriv = zeros(ndim, fin_el.nod)         #
  bee = zeros(nst,ndof)                    #
  km = zeros(ndof, ndof)                   #
  mm = zeros(ndof, ndof)                   #
  gm = zeros(ndof, ndof)                   #
  kg = zeros(ndof, ndof)                   #
  eld = zeros(ndof)                        #
  weights = zeros(struc_el.nip)        #
  g_g = zeros(Int, ndof, nels)           #
  num = zeros(Int, fin_el.nod)          #
  actions = zeros(nels, ndof)              #
  nf = ones(Int, nodof, nn)              #
  displacements = zeros(size(nf, 1), ndim) #
  gc = ones(ndim, ndim)                    #
  dee = zeros(nst,nst)                     #
  sigma = zeros(nst)                       #
  axial = zeros(nels)                      #
  x_coords = zeros(nn)                     #
  y_coords = zeros(nn)                     # Not used, needed for FEM constructor
  z_coords = zeros(nn)                     #
  etype = ones(Int, nels)                #
  ell = zeros(nels)

  # Start with arrays to be initialized from input dict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    throw("No :properties key found in input dict to FE4_1.")
  end

  # Update nodal freedom array from input dict :support entry
  # nf has dimensions [nodof, nn] and all entries are set to 1 (free)

  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  else
    throw("No :support key found in input dict to FE4_1.")
  end

  if :x_coords in keys(data)
    x_coords = data[:x_coords]
    @assert length(x_coords) == nels + 1
  else
    dl = 1.0/nels
    x_coords = 0.0:dl:1.0
  end

  # Set lengths of elements
  for i in 1:length(x_coords)-1
    ell[i] = x_coords[i+1] - x_coords[i]
  end

  # In multiple np_types, get np_type for each structural element
  if :etype in keys(data)
    etype = data[:etype]
  end

  # Done with input dict and allocations, re-do nf (assign number to each dof)

  PtFEM.formnf!(nodof, nn, nf)
  neq = maximum(nf)
  #kdiag = zeros(Int, neq)

  # Set global numbering, coordinates and array sizes

  for i in 1:nels
    num = [i; i+1]
    PtFEM.num_to_g!(num, nf, g)
    g_g[:, i] = g
  end

  println("There are $(neq) equations.")

  loads = OffsetArray(zeros(neq+1), 0:neq)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end

  gsm = spzeros(neq, neq)
  for i in 1:nels
    km = PtFEM.rod_km!(km, prop[etype[i], 1], ell[i])
    g = g_g[:, i]
    PtFEM.fsparm!(gsm, g, km)
  end

  fixed_freedoms = 0
  if :fixed_freedoms in keys(data)
    fixed_freedoms = size(data[:fixed_freedoms], 1)
  end
  no = zeros(Int, fixed_freedoms)
  node = zeros(Int, fixed_freedoms)
  sense = zeros(Int, fixed_freedoms)
  value = zeros(Float64, fixed_freedoms)
  if :fixed_freedoms in keys(data) && fixed_freedoms > 0
    for i in 1:fixed_freedoms
      node[i] = data[:fixed_freedoms][i][1]
      sense[i] = data[:fixed_freedoms][i][2]
      no[i] = nf[sense[i], node[i]]
      value[i] = data[:fixed_freedoms][i][3]
      gsm[no[i], no[i]] += penalty
      loads[no[i]] = gsm[no[i], no[i]] .* value[i]
    end
  end

  cfgsm = cholesky(gsm)

  loads[1:neq] = cfgsm \ loads[1:neq]
  println()

  displacements = zeros(size(nf))
  for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if nf[i, j] > 0
        displacements[i,j] = loads[nf[i, j]]
      end
    end
  end
  displacements = displacements'

  loads[0] = 0.0
  for i in 1:nels
    km = PtFEM.rod_km!(km, prop[etype[i], 1], ell[i])
    g = g_g[:, i]
    eld = loads[g]
    actions[i, :] = km * eld
  end

  dis_df = DataFrame(
    x_translation = displacements[:, 1],
  )

  if :eq_nodal_forces_and_moments in keys(data)
    fm_df = DataFrame(
      normal_force_1 = actions[:, 1],
      normal_force_2 = actions[:, 2],
      nf_1_uncorrected = actions[:, 1],
      nf_2_uncorrected = actions[:, 2]
    )
  else
    fm_df = DataFrame(
      normal_force_1 = actions[:, 1],
      normal_force_2 = actions[:, 2]
    )
  end

  # Correct element forces and moments for equivalent nodal
  # forces and moments introduced for loading between nodes
  if :eq_nodal_forces_and_moments in keys(data)
    eqfm = data[:eq_nodal_forces_and_moments]
    k = data[:struc_el].fin_el.nod * data[:struc_el].fin_el.nodof
    for t in eqfm
      vals = convert(Array, fm_df[t[1], :])
      for i in 1:k
        fm_df[t[1], i] = round.(vals[i] - t[2][i], digits=2)
      end
    end
  end

  fem = PtFEM.jFEM(struc_el, fin_el, ndim, nels, nst, ndof, nn, nodof,
    neq, penalty, etype, g, g_g, g_num, nf, no,
    node, num, sense, actions, bee, coord, gamma, dee,
    der, deriv, displacements, eld, fun, gc, g_coord, jac,
    km, mm, kg, cfgsm, loads, points, prop, sigma, value,
    weights, x_coords, y_coords, z_coords, axial)

  (fem, dis_df, fm_df)
end

"""
# Method p41

One dimensional analysis of an axially loaded elastic Rod using 2-node
Line elements.

### Constructors
```julia
p41(m, data) # Re-use factored global stiffness matrix
```
### Arguments
```julia
* `m::jFEM`                  : Previously created jFEM model
* `data::Dict{Symbol, Any}`  : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Type of  structural fin_el
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* loaded_nodes::Array{Tuple{Int,Array{Float64,1}},1} : Node load vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Float64}                        : x-coordinate vector
```

### Optional additional data dictionary keys
```julia
* penalty = 1e20               : Penalty used for fixed degrees of freedoms
* etype::Vector{Int}         : Element material vector if np_types > 1
* eq_nodal_forces_and_moments  : Contribution of distributed loads to loaded_nodes
```

### Return values
```julia
* (jfem, dis_df, fm_df)        : Tuple of jFem, dis_df and fm_df
                                 where:
                                    jfem::jFem    : Computational result type
                                    dis_df        : Displacement data table
                                    fm_df         : Forces and moments data table
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Rod                           : Help on a Rod structural element
?FiniteElement                 : List finite element types
?Line                          : Help on Line finite element
```
"""
function p41(m::PtFEM.jFEM, data::Dict)
  loads = OffsetArray(zeros(m.neq+1), 0:m.neq)
  if :loaded_nodes in keys(data)
    for i in 1:size(data[:loaded_nodes], 1)
      loads[m.nf[:, data[:loaded_nodes][i][1]]] = data[:loaded_nodes][i][2]
    end
  end

  loads[1:m.neq] = m.cfgsm \ loads[1:m.neq]
  println()

  displacements = zeros(size(m.nf))
  for i in 1:size(displacements, 1)
    for j in 1:size(displacements, 2)
      if m.nf[i, j] > 0
        displacements[i,j] = loads[m.nf[i, j]]
      end
    end
  end
  displacements = displacements'

  # Set lengths of elements
  ell = zeros(length(m.x_coords)-1)
  for i in 1:length(m.x_coords)-1
    ell[i] = m.x_coords[i+1] - m.x_coords[i]
  end

  km = zeros(m.ndof, m.ndof)
  g = zeros(Int, m.ndof)
  eld = zeros(m.ndof)
  actions = zeros(m.nels, m.ndof)
  loads[0] = 0.0
  for i in 1:m.nels
    km = PtFEM.rod_km!(m.km, m.prop[m.etype[i], 1], ell[i])
    g = m.g_g[:, i]
    eld = loads[g]
    actions[i, :] = km * eld
  end

  dis_df = DataFrame(
    x_translation = displacements[:, 1],
  )

  if :eq_nodal_forces_and_moments in keys(data)
    fm_df = DataFrame(
    normal_force_1 = actions[:, 1],
    normal_force_2 = actions[:, 2],
    nf_1_uncorrected = actions[:, 1],
    nf_2_uncorrected = actions[:, 2]
    )
  else
    fm_df = DataFrame(
      normal_force_1 = actions[:, 1],
      normal_force_2 = actions[:, 2]
    )
  end

  # Correct element forces and moments for equivalent nodal
  # forces and moments introduced for loading between nodes
  if :eq_nodal_forces_and_moments in keys(data)
    eqfm = data[:eq_nodal_forces_and_moments]
    k = data[:struc_el].fin_el.nod * data[:struc_el].fin_el.nodof
    for t in eqfm
      vals = convert(Array, fm_df[t[1], :])
      for i in 1:k
        fm_df[t[1], i] = round.(vals[i] - t[2][i], digits=2)
      end
    end
  end

  fem = PtFEM.jFEM(m.struc_el, m.fin_el, m.ndim, m.nels, m.nst, m.ndof, m.nn, m.nodof,
    m.neq, m.penalty, m.etype, g, m.g_g, m.g_num, m.nf, m.no,
    m.node, m.num, m.sense, actions, m.bee, m.coord, m.gamma, m.dee,
    m.der, m.deriv, displacements, eld, m.fun, m.gc, m.g_coord, m.jac,
    km, m.mm, m.kg, m.cfgsm, loads, m.points, m.prop, m.sigma, m.value,
    m.weights, m.x_coords, m.y_coords, m.z_coords, m.axial)

  (fem, dis_df, fm_df)
end

max_stars_repo_path: src/4 Static Equilibrium/p41.jl

max_issues_repo_path: src/4 Static Equilibrium/p41.jl

max_forks_repo_path: src/4 Static Equilibrium/p41.jl