text:
# Additions to Base's searchsorted functionality
"""
    searchsortednearest(vec::AbstractVector, x)

Like `searchsortedfirst` or `searchsortedlast`, this returns the the index of
the element in the sorted vector `vec` whose value is closest to `x`, rounding
up. If there are multiple elements that are equally close to `x`, this will
return the first index if `x` is less than or equal to those in the vector or
the last index if `x` is greater.
"""
function searchsortednearest(vec::AbstractVector, x)
    idx = searchsortedfirst(vec, x) # Returns the first idx | vec[idx] >= x
    if idx > 1 && (idx > length(vec) || (vec[idx] - x) > (x - vec[idx-1]))
        idx -= 1 # The previous element is closer
    end
    return idx
end
# Base only specializes searching ranges by Numbers; so optimize for Intervals
function Base.searchsorted(a::AbstractRange, I::ClosedInterval)
    searchsortedfirst(a, I.left):searchsortedlast(a, I.right)
end

"""
The internal `Extrapolated` module contains implementations for indexing and
searching into ranges beyond their bounds. The `@inbounds` macro is not
sufficient since it can be turned off by `--check-bounds=yes`.
"""
module Extrapolated
using IntervalSets: ClosedInterval

function searchsortednearest(vec::AbstractRange, x)
    idx = searchsortedfirst(vec, x) # Returns the first idx | vec[idx] >= x
    if (getindex(vec, idx) - x) > (x - getindex(vec, idx-1))
        idx -= 1 # The previous element is closer
    end
    return idx
end

"""
    searchsorted(a::AbstractRange, I::ClosedInterval)

Return the indices of the range that fall within an interval without checking
bounds, possibly extrapolating outside the range if needed.
"""
function searchsorted(a::AbstractRange, I::ClosedInterval)
    searchsortedfirst(a, I.left):searchsortedlast(a, I.right)
end

# When running with `--check-bounds=yes` (like on Travis), the bounds-check isn't elided
@inline function getindex(v::AbstractRange{T}, i::Integer) where T
    convert(T, first(v) + (i-1)*step(v))
end
@inline function getindex(r::AbstractRange, s::AbstractRange{<:Integer})
    f = first(r)
    st = oftype(f, f + (first(s)-1)*step(r))
    range(st, step=step(r)*step(s), length=length(s))
end
getindex(r::AbstractRange, I::Array) = [getindex(r, i) for i in I]
@inline getindex(r::StepRangeLen, i::Integer) = Base.unsafe_getindex(r, i)
@inline function getindex(r::StepRangeLen, s::AbstractUnitRange)
    soffset = 1 + (r.offset - first(s))
    soffset = clamp(soffset, 1, length(s))
    ioffset = first(s) + (soffset-1)
    if ioffset == r.offset
        StepRangeLen(r.ref, r.step, length(s), max(1,soffset))
    else
        StepRangeLen(r.ref + (ioffset-r.offset)*r.step, r.step, length(s), max(1,soffset))
    end
end

function searchsortedlast(a::AbstractRange, x)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    n = round(Integer,(x-first(a))/step(a))+1
    isless(x, getindex(a, n)) ? n-1 : n
end
function searchsortedfirst(a::AbstractRange, x)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    n = round(Integer,(x-first(a))/step(a))+1
    isless(getindex(a, n), x) ? n+1 : n
end
function searchsortedlast(a::AbstractRange{<:Integer}, x)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    fld(floor(Integer,x)-first(a),step(a))+1
end
function searchsortedfirst(a::AbstractRange{<:Integer}, x)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    -fld(floor(Integer,-x)+first(a),step(a))+1
end
function searchsortedfirst(a::AbstractRange{<:Integer}, x::Unsigned)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    -fld(first(a)-signed(x),step(a))+1
end
function searchsortedlast(a::AbstractRange{<:Integer}, x::Unsigned)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    fld(signed(x)-first(a),step(a))+1
end
end


text:
# Test script to test solver for a qp

using Test
using COSMO, SparseArrays, LinearAlgebra

# Quadratic program example from OSQP Doc
# min 0.5 * x'Px +  q'x
# s.t. l <= Ax <= u
q = [1; 1.]
P = sparse([4. 1; 1 2])
A = [1. 1; 1 0; 0 1]
l = [1.; 0; 0]
u = [1; 0.7; 0.7]

# Define the constraint l <= Ax <= u with the help of a Nonnegatives set
Aa = [-A; A]
ba = [u; -l]
constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives)

# define example problem
settings = COSMO.Settings(verbose=true, eps_abs = 1e-4, eps_rel = 1e-4)

model = COSMO.Model()
assemble!(model, P, q, constraint1, settings = settings)
res = COSMO.optimize!(model);

# solve again by defining the constraints with the help of a box
constraint1 = COSMO.Constraint(A, zeros(3), COSMO.Box(l, u))

model = COSMO.Model()
assemble!(model, P, q, constraint1, settings = settings)
res_box = COSMO.optimize!(model);

@testset "QP Problem" begin
  @test norm(res.x[1:2] - [0.3; 0.7], Inf) < 1e-3
  @test norm(res_box.x[1:2] - [0.3; 0.7], Inf) < 1e-3
  @test abs(res.obj_val - 1.88) < 1e-3
  @test abs(res_box.obj_val - 1.88) < 1e-3
end
nothing

ext: jl


text:
push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")

println("On startup, we have ", nprocs(), " processes.")
println("On startup, we have ", nworkers(), " workers.")
nchild = nworkers()

@everywhere using parallel

@everywhere function f(dset, key, p)
    println("Likelihood call to $dset, $p, $key on process ", myid(), " Hostname: ", readall(`hostname`))
    return 0.0
end

@everywhere function initfunc(key::Int)
    println("Loading dset $key")
    return key
end

const global keylist = Int[i for i=1:nchild]
pipes = initialize(nchild, keylist, initfunc, f)

for j=1:nworkers()
    distribute!(pipes, [0.0, 0.0])
    result = gather!(pipes)
    println(result)
end

println("completing all results")
quit!(pipes)
println("all results should be completed")

# For some reason the workers() command is slow, and if you call it multiple times
# in quick succession it will give outdated information. Also, if you call it while
# there are no workers, it will return [1]
# for j=1:3
#     println(workers())
#     sleep(2)
# end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using Optim

#function find_energy(V, L, initial)
function find_energy(z0, ini)

#    a = L/2
#    z0 = a * sqrt(2*V)

    function f(z)

        L = sqrt.( abs.((z0./z).^2 .- 1))
        R = tan.(z)
        return sum( (L - R).^2)
        #        alpha = sqrt.( 2.0 / L^2 * abs.(V .+ E))
#        k = sqrt.(abs.(2.0* E))

#        return sum(( k .* tan.(k*L/2.0) .- alpha ).^2)
    end

    ret = optimize(f, ini, BFGS())
    println(ret)
    println()
    println("min val ", f(Optim.minimizer(ret)))
    println(Optim.minimizer(ret))

    return Optim.minimizer(ret)

end


text:
import RDatasets # R datasets from GitHub
import DataFrames
import Plots
import StatsPlots
import KernelDensity
import LinearAlgebra
import HypothesisTests
import PyCall
import MLBase

Plots.plotly()

DF = DataFrames

# Get data from RDatasets on GitHub
D = RDatasets.dataset("datasets", "faithful")
@show DF.names(D)
@show DF.describe(D)

# what does ! mean
eruptions = D[!, :Eruptions]
eruptions = D[:, :Eruptions] # same as above?
waittime = D[!, :Waiting]
pt = Plots.scatter(eruptions, label="Eruptions")
Plots.scatter!(pt, waittime, label="Wait time")
Plots.gui()

pt = StatsPlots.boxplot(["Eruption Length"], eruptions, legend=true,
    size=(200,400), whisker_width=1, ylabel="Time")
Plots.gui(pt)

pt = StatsPlots.histogram(eruptions, label="Eruptions")
Plots.gui(pt)

# Kernel density estimates
p = KernelDensity.kde(eruptions)
pt = StatsPlots.histogram(eruptions, label="Eruptions")
Plots.plot!(pt, p.x, p.density .* length(eruptions), linewidth=3, color=2,
    label="KDE Fit")
Plots.gui(pt)

A = [1 6; 2 5; 3 4]
p = KernelDensity.kde(A)
pt = Plots.plot(A[:,1], A[:,2])
Plots.plot!(pt, p.x, p.density , linewidth=3, color=2,
    label="KDE Fit")
Plots.gui(pt)

# Skip


text:
# Copy and paste these line into a REPL and screen capture to produce the
# sample usage image.
using Memento
logger = Memento.config!("trace"; fmt="[{level} | {name}]: {msg}")
trace(logger, "Some code you're tracing.")
debug(logger, "Something to help you track down a bug.")
info(logger, "Something you might want to know.")
notice(logger, "This is probably pretty important.")
warn(logger, ErrorException("A caught exception that we want to log as a warning."))
error(logger, "Something that should throw an error.")
critical(logger, "The world is exploding")
child_logger = getlogger("Foo.bar")
setlevel!(child_logger, "warn")
push!(child_logger, DefaultHandler(tempname(), DefaultFormatter("[{date} | {level} | {name}]: {msg}")));
debug(child_logger, "Something that should only be printed to stdout on the root_logger.")
warn(child_logger, "Warning to stdout and the log file.")
exit()


text:
# Utility functions for CliffordAlgebras.jl

"""
    nextselector(s::NTuple{N,<:Integer},M) where N

Helper function for enumerating base vectors.
"""
function nextselector(s::NTuple{N,<:Integer}, M) where {N}
    @assert issorted(s)
    @assert all(s .<= M) && all(s .>= 1)
    for n = 1:(N-1)
        if s[n] + 1 < s[n+1]
            return ntuple(i -> i < n ? i : i == n ? s[i] + 1 : s[i], N)
        end
    end
    if s[N] < M
        return ntuple(i -> i < N ? i : i == N ? s[i] + 1 : s[i], N)
    else
        return nothing
    end
end

"""
    enumeratebase(N::Integer)

Generates the basis vector product combinations for the Clifford algebra of an N dimensional vector space.
"""
function enumeratebase(N::Integer)
    result = Tuple[]
    push!(result, ())
    for n = 1:N
        s = ntuple(identity, n)
        while true
            push!(result, s)
            s = nextselector(s, N)
            !isnothing(s) || break
        end
    end
    Tuple(result)
end

"""
    adaptbasefordual(base)

Shuffles the base vector products so that the Poincaré dual of every basis vector is another basis vector.
"""
function adaptbasefordual(base)
    N = length(base)
    Nh = N ÷ 2
    Tuple(
        n > Nh && length(base[n]) >= 2 && levicivita([base[N-n+1]..., base[n]...]) < 0 ?
        Tuple(
            k == length(base[n]) ? base[n][k-1] :
            k + 1 == length(base[n]) ? base[n][k+1] : base[n][k] for k = 1:length(base[n])
        ) : base[n] for n = 1:N
    )
end

"""
    findpermutation(a,b)

Find the permutation that takes the a to b.
"""
function findpermutation(a, b)
    @assert allunique(a) "Permutations can only be found in sequences of unique elements."
    @assert length(a) == length(b) "Permutations can only be found between sequences of identical length."
    @assert length(intersect(a, b)) == length(a) "Permutations can only be found between sequences with the same elements."
    [findfirst(isequal(a[n]), b) for n = 1:length(a)]
end

max_stars_repo_name: ATell-SoundTheory/CliffordAlgebras

max_issues_repo_name: ATell-SoundTheory/CliffordAlgebras


text:
"""
    q1_q0

The manufactured-solution colliding flow example from Elman et al 2014. The
linear quadrilaterals for the velocity and discontinuous
pressure.

The formulation is the one derived in Reddy, Introduction to the finite element
method, 1993. Page 486 ff.
"""
module q1_q0

using LinearAlgebra
using StaticArrays
using MeshCore
using MeshCore: retrieve, nrelations, nentities, ir_identity, attribute
using MeshSteward: Q4block
using MeshSteward: Mesh, attach!, baseincrel, boundary, increl
using MeshSteward: vselect, geometry, summary, transform
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEL2_Q4, FEH1_Q4, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs
using Elfel.FEIterators: jacjac
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # This mesh will be both for the velocities and for the pressure
    mesh = Mesh()
    attach!(mesh, Q4block(2 * A, 2 * A, N, N), "velocity+pressure")
    ir = baseincrel(mesh)
    transform(ir, x -> x .- A)
    eidir = ir_identity(ir)
    attach!(mesh, eidir)
    return mesh
end

function assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    function integrateK!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxuy = LocalMatrixAssembler(uxnedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxuy, eldofs(uxel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(uxel, uxqp) # L2 pressure: must integrate over the velocity mesh
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * (2 * gradNux[i][1] * gradNux[j][1] + gradNux[i][2] * gradNux[j][2])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * (gradNuy[i][1] * gradNuy[j][1] + 2 * gradNuy[i][2] * gradNuy[j][2])
                end
                for j in 1:uynedof, i in 1:uxnedof
                    kuxuy[i, j] += (mu * JxW) * (gradNux[i][1] * gradNuy[j][2])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuxuy)
            assemble!(ass, transpose(kuxuy))
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(uxfesp), FEIterator(uyfesp), FEIterator(pfesp))
    qargs = (kind = :default, order = 2,)
    qpits = (QPIterator(uxfesp, qargs), QPIterator(uyfesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    @time integrateK!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    @time KT = K * U
    @time U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

centroid(p, c) = begin
    (p[c[1]] + p[c[2]] + p[c[3]] + p[c[4]]) / 4
end

function evaluate_error(uxfesp, uyfesp, pfesp)
    geom = geometry(pfesp.mesh)
    bir = baseincrel(pfesp.mesh)
    points = MeshCore.attribute(bir.right, "geom")
    centroids = [centroid(points, retrieve(bir, i)) for i in 1:nrelations(bir)]
    p = attribute(increl(pfesp.mesh, (2, 2)).right, "p")
    pt = [truep(centroids[i]...) for i in 1:length(centroids)]
    return norm(p - pt) / norm(pt)
end

function run(N)
    mesh = genmesh(N)
    # Velocity spaces
    uxfesp = FESpace(Float64, mesh, FEH1_Q4(), 1)
    uyfesp = FESpace(Float64, mesh, FEH1_Q4(), 1)
    locs = geometry(mesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(uxfesp, 0, i, 1, trueux(locs[i]...))
            setebc!(uyfesp, 0, i, 1, trueuy(locs[i]...))
        end
    end
    # Pressure space
    pfesp = FESpace(Float64, mesh, FEL2_Q4(), 1)
    atcenter = vselect(geometry(mesh); nearestto = [0.0, 0.0])
    setebc!(pfesp, 2, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([uxfesp, uyfesp, pfesp])
    @show tndof = ndofs(uxfesp) + ndofs(uyfesp) + ndofs(pfesp)
    @show tnunk = nunknowns(uxfesp) + nunknowns(uyfesp) + nunknowns(pfesp)
    # Assemble the coefficient matrix
    K = assembleK(uxfesp, uyfesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [uxfesp, uyfesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([uxfesp, uyfesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(uxfesp, "ux", 1)
    makeattribute(uyfesp, "uy", 1)
    @show evaluate_error(uxfesp, uyfesp, pfesp)
    vtkwrite("q1_q0-p", baseincrel(mesh), [(name = "p",), ])
    vtkwrite("q1_q0-v", baseincrel(mesh), [(name = "ux",), (name = "uy",)])
    true
end

end

let
    N = 4
    for loop in 1:5
        q1_q0.run(N)
        N = N * 2
    end
end


text:
_nrow(df::AbstractDataFrame) = nrow(df)
_nrow(x::NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}) =
    isempty(x) ? 0 : length(x[1])
_ncol(df::AbstractDataFrame) = ncol(df)
_ncol(x::Union{NamedTuple, DataFrameRow}) = length(x)

function _combine_multicol((firstres,)::Ref{Any}, wfun::Ref{Any}, gd::GroupedDataFrame,
                           wincols::Ref{Any})
    @assert only(wfun) isa Base.Callable
    @assert only(wincols) isa Union{Nothing, AbstractVector, Tuple, NamedTuple}
    firstmulticol = firstres isa MULTI_COLS_TYPE
    if !(firstres isa Union{AbstractVecOrMat, AbstractDataFrame,
                            NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}})
        idx_agg = Vector{Int}(undef, length(gd))
        fillfirst!(nothing, idx_agg, 1:length(gd.groups), gd)
    else
        idx_agg = NOTHING_IDX_AGG
    end
    return _combine_with_first(Ref{Any}(wrap(firstres)), wfun, gd, wincols,
                               Val(firstmulticol), idx_agg)
end

function _combine_with_first((first,)::Ref{Any},
                             (f,)::Ref{Any}, gd::GroupedDataFrame,
                             (incols,)::Ref{Any},
                             firstmulticol::Val, idx_agg::Vector{Int})
    @assert f isa Base.Callable
    @assert incols isa Union{Nothing, AbstractVector, Tuple, NamedTuple}
    @assert first isa Union{NamedTuple, DataFrameRow, AbstractDataFrame}
    extrude = false

    if first isa AbstractDataFrame
        n = 0
        eltys = eltype.(eachcol(first))
    elseif first isa NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}
        n = 0
        eltys = map(eltype, first)
    elseif first isa DataFrameRow
        n = length(gd)
        eltys = [eltype(parent(first)[!, i]) for i in parentcols(index(first))]
    elseif firstmulticol == Val(false) && first[1] isa Union{AbstractArray{<:Any, 0}, Ref}
        extrude = true
        first = wrap_row(first[1], firstmulticol)
        n = length(gd)
        eltys = (typeof(first[1]),)
    else # other NamedTuple giving a single row
        n = length(gd)
        eltys = map(typeof, first)
        if any(x -> x <: AbstractVector, eltys)
            throw(ArgumentError("mixing single values and vectors in a named tuple is not allowed"))
        end
    end
    idx = idx_agg === NOTHING_IDX_AGG ? Vector{Int}(undef, n) : idx_agg
    local initialcols
    let eltys=eltys, n=n # Workaround for julia#15276
        initialcols = ntuple(i -> Tables.allocatecolumn(eltys[i], n), _ncol(first))
    end
    targetcolnames = tuple(propertynames(first)...)
    if !extrude && first isa Union{AbstractDataFrame,
                                   NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}}
        outcols, finalcolnames = _combine_tables_with_first!(first, initialcols, idx, 1, 1,
                                                             f, gd, incols, targetcolnames,
                                                             firstmulticol)
    else
        outcols, finalcolnames = _combine_rows_with_first!(first, initialcols,
                                                           f, gd, incols, targetcolnames,
                                                           firstmulticol)
    end
    return idx, outcols, collect(Symbol, finalcolnames)
end

function _names_match(r1, t::NTuple{N, Symbol}) where N
    n1 = _getnames(r1)
    length(n1) == length(t) || return false
    for (a, b) in zip(n1, t)
        a == b || return false
    end
    return true
end

function fill_row!(row, outcols::NTuple{N, AbstractVector},
                   i::Integer, colstart::Integer,
                   colnames::NTuple{N, Symbol}) where N
    if !_names_match(row, colnames)
        throw(ArgumentError("return value must have the same column names " *
                            "for all groups (got $colnames and $(propertynames(row)))"))
    end
    @inbounds for j in colstart:length(outcols)
        col = outcols[j]
        val = row[j]
        S = typeof(val)
        T = eltype(col)
        if S <: T || promote_type(S, T) <: T
            col[i] = val
        else
            return j
        end
    end
    return nothing
end

function _combine_rows_with_first_task!(tid::Integer,
                                        rowstart::Integer,
                                        rowend::Integer,
                                        rownext::Integer,
                                        outcols::NTuple{<:Any, AbstractVector},
                                        outcolsref::Ref{NTuple{<:Any, AbstractVector}},
                                        type_widened::AbstractVector{Bool},
                                        widen_type_lock::ReentrantLock,
                                        f::Base.Callable,
                                        gd::GroupedDataFrame,
                                        starts::AbstractVector{<:Integer},
                                        ends::AbstractVector{<:Integer},
                                        incols::Union{Nothing, AbstractVector,
                                                      Tuple, NamedTuple},
                                        colnames::NTuple{<:Any, Symbol},
                                        firstmulticol::Val)
    j = nothing
    gdidx = gd.idx
    local newoutcols
    for i in rownext:rowend
        row = wrap_row(do_call(f, gdidx, starts, ends, gd, incols, i), firstmulticol)
        j = fill_row!(row, outcols, i, 1, colnames)
        if j !== nothing # Need to widen column
            # If another thread is already widening outcols, wait until it's done
            lock(widen_type_lock)
            try
                newoutcols = outcolsref[]
                # Workaround for julia#15276
                newoutcols = let i=i, j=j, newoutcols=newoutcols, row=row
                    ntuple(length(newoutcols)) do k
                        S = typeof(row[k])
                        T = eltype(newoutcols[k])
                        U = promote_type(S, T)
                        if S <: T || U <: T
                            newoutcols[k]
                        else
                            type_widened .= true
                            Tables.allocatecolumn(U, length(newoutcols[k]))
                        end
                    end
                end
                for k in 1:length(outcols)
                    if outcols[k] !== newoutcols[k]
                        copyto!(newoutcols[k], rowstart,
                                outcols[k], rowstart, i - rowstart + (k < j))
                    end
                end
                j = fill_row!(row, newoutcols, i, j, colnames)
                @assert j === nothing # eltype is guaranteed to match

                outcolsref[] = newoutcols
                type_widened[tid] = false
            finally
                unlock(widen_type_lock)
            end
            return _combine_rows_with_first_task!(tid, rowstart, rowend, i+1, newoutcols, outcolsref,
                                                  type_widened, widen_type_lock,
                                                  f, gd, starts, ends,
                                                  incols, colnames, firstmulticol)
        end
        # If other thread widened columns, copy already processed data to new vectors
        # This doesn't have to happen immediately (hence type_widened isn't atomic),
        # but the more we wait the more data will have to be copied
        if type_widened[tid]
            lock(widen_type_lock) do
                type_widened[tid] = false
                newoutcols = outcolsref[]
                for k in 1:length(outcols)
                    # Check whether this particular column has been widened
                    if outcols[k] !== newoutcols[k]
                        copyto!(newoutcols[k], rowstart,
                                outcols[k], rowstart, i - rowstart + 1)
                    end
                end
            end
            return _combine_rows_with_first_task!(tid, rowstart, rowend, i+1, newoutcols, outcolsref,
                                                  type_widened, widen_type_lock,
                                                  f, gd, starts, ends,
                                                  incols, colnames, firstmulticol)
        end
    end
    return outcols
end

# CategoricalArray is thread-safe only when input and output levels are equal
# since then all values are added upfront. Otherwise, if the function returns
# CategoricalValues mixed from different pools, one thread may add values,
# which may put the invpool in an invalid state while the other one is reading from it
function isthreadsafe(outcols::NTuple{<:Any, AbstractVector},
                      incols::Union{Tuple, NamedTuple})
    anycat = any(outcols) do col
        T = typeof(col)
        # If the first result is missing, widening can give a CategoricalArray
        # if later results are CategoricalValues
        eltype(col) === Missing ||
            (nameof(T) === :CategoricalArray &&
             nameof(parentmodule(T)) === :CategoricalArrays)
    end
    if anycat
        levs = nothing
        for col in incols
            T = typeof(col)
            if nameof(T) === :CategoricalArray &&
                nameof(parentmodule(T)) === :CategoricalArrays
                if levs !== nothing
                    levs == levels(col) || return false
                else
                    levs = levels(col)
                end
            end
        end
    end
    return true
end
isthreadsafe(outcols::NTuple{<:Any, AbstractVector}, incols::AbstractVector) =
    isthreadsafe(outcols, (incols,))
isthreadsafe(outcols::NTuple{<:Any, AbstractVector}, incols::Nothing) = true

function _combine_rows_with_first!(firstrow::Union{NamedTuple, DataFrameRow},
                                   outcols::NTuple{N, AbstractVector},
                                   f::Base.Callable, gd::GroupedDataFrame,
                                   incols::Union{Nothing, AbstractVector, Tuple, NamedTuple},
                                   colnames::NTuple{N, Symbol},
                                   firstmulticol::Val) where N
    len = length(gd)
    gdidx = gd.idx
    starts = gd.starts
    ends = gd.ends

    # handle empty GroupedDataFrame
    len == 0 && return outcols, colnames

    # Handle first group
    j1 = fill_row!(firstrow, outcols, 1, 1, colnames)
    @assert j1 === nothing # eltype is guaranteed to match

    # Handle groups other than the first one
    # Create up to one task per thread
    # This has lower overhead than creating one task per group,
    # but is optimal only if operations take roughly the same time for all groups
    if VERSION >= v"1.4" && isthreadsafe(outcols, incols)
        basesize = max(1, cld(len - 1, Threads.nthreads()))
        partitions = Iterators.partition(2:len, basesize)
    else
        partitions = (2:len,)
    end
    widen_type_lock = ReentrantLock()
    outcolsref = Ref{NTuple{<:Any, AbstractVector}}(outcols)
    type_widened = fill(false, length(partitions))
    tasks = Vector{Task}(undef, length(partitions))
    for (tid, idx) in enumerate(partitions)
        tasks[tid] =
            @spawn _combine_rows_with_first_task!(tid, first(idx), last(idx), first(idx),
                                                  outcols, outcolsref,
                                                  type_widened, widen_type_lock,
                                                  f, gd, starts, ends, incols, colnames,
                                                  firstmulticol)
    end

    # Workaround JuliaLang/julia#38931:
    # we want to preserve the exception type thrown in user code,
    # and print the backtrace corresponding to it
    for t in tasks
        try
            wait(t)
        catch e
            throw(t.exception)
        end
    end

    # Copy data for any tasks that finished before others widened columns
    oldoutcols = outcols
    outcols = outcolsref[]
    if outcols !== oldoutcols # first group
        for k in 1:length(outcols)
            outcols[k][1] = oldoutcols[k][1]
        end
    end
    for (tid, idx) in enumerate(partitions)
        if type_widened[tid]
            oldoutcols = fetch(tasks[tid])
            for k in 1:length(outcols)
                # Check whether this particular column has been widened
                if oldoutcols[k] !== outcols[k]
                    copyto!(outcols[k], first(idx), oldoutcols[k], first(idx),
                            last(idx) - first(idx) + 1)
                end
            end
        end
    end

    return outcols, colnames
end

# This needs to be in a separate function
# to work around a crash due to JuliaLang/julia#29430
if VERSION >= v"1.1.0-DEV.723"
    @inline function do_append!(do_it, col, vals)
        do_it && append!(col, vals)
        return do_it
    end
else
    @noinline function do_append!(do_it, col, vals)
        do_it && append!(col, vals)
        return do_it
    end
end

_get_col(rows::AbstractDataFrame, j::Int) = rows[!, j]
_get_col(rows::NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}, j::Int) = rows[j]

function append_rows!(rows, outcols::NTuple{N, AbstractVector},
                      colstart::Integer, colnames::NTuple{N, Symbol}) where N
    if !_names_match(rows, colnames)
        throw(ArgumentError("return value must have the same column names " *
                            "for all groups (got $colnames and $(propertynames(rows)))"))
    end
    @inbounds for j in colstart:length(outcols)
        col = outcols[j]
        cn = colnames[j]
        vals = _get_col(rows, j)
        S = eltype(vals)
        T = eltype(col)
        if !do_append!(S <: T || promote_type(S, T) <: T, col, vals)
            return j
        end
    end
    return nothing
end

function _combine_tables_with_first!(first::Union{AbstractDataFrame,
                                     NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}}},
                                     outcols::NTuple{N, AbstractVector},
                                     idx::Vector{Int}, rowstart::Integer, colstart::Integer,
                                     f::Base.Callable, gd::GroupedDataFrame,
                                     incols::Union{Nothing, AbstractVector, Tuple, NamedTuple},
                                     colnames::NTuple{N, Symbol},
                                     firstmulticol::Val) where N
    len = length(gd)
    gdidx = gd.idx
    starts = gd.starts
    ends = gd.ends
    # Handle first group

    @assert _ncol(first) == N
    if !isempty(colnames) && length(gd) > 0
        j = append_rows!(first, outcols, colstart, colnames)
        @assert j === nothing # eltype is guaranteed to match
        append!(idx, Iterators.repeated(gdidx[starts[rowstart]], _nrow(first)))
    end
    # Handle remaining groups
    @inbounds for i in rowstart+1:len
        rows = wrap_table(do_call(f, gdidx, starts, ends, gd, incols, i), firstmulticol)
        _ncol(rows) == 0 && continue
        if isempty(colnames)
            newcolnames = tuple(propertynames(rows)...)
            if rows isa AbstractDataFrame
                eltys = eltype.(eachcol(rows))
            else
                eltys = map(eltype, rows)
            end
            initialcols = ntuple(i -> Tables.allocatecolumn(eltys[i], 0), _ncol(rows))
            return _combine_tables_with_first!(rows, initialcols, idx, i, 1,
                                               f, gd, incols, newcolnames, firstmulticol)
        end
        j = append_rows!(rows, outcols, 1, colnames)
        if j !== nothing # Need to widen column type
            local newcols
            let i = i, j = j, outcols=outcols, rows=rows # Workaround for julia#15276
                newcols = ntuple(length(outcols)) do k
                    S = eltype(rows isa AbstractDataFrame ? rows[!, k] : rows[k])
                    T = eltype(outcols[k])
                    U = promote_type(S, T)
                    if S <: T || U <: T
                        outcols[k]
                    else
                        copyto!(Tables.allocatecolumn(U, length(outcols[k])), outcols[k])
                    end
                end
            end
            return _combine_tables_with_first!(rows, newcols, idx, i, j,
                                               f, gd, incols, colnames, firstmulticol)
        end
        append!(idx, Iterators.repeated(gdidx[starts[i]], _nrow(rows)))
    end
    return outcols, colnames
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module BlsData

using Compat
using DataFrames

import Dates
import JSON
import HttpCommon
import HTTP
import Printf: @printf

export
    # Bls type
    Bls, get_api_url, set_api_url!, get_api_key, get_api_version, requests_made,
        requests_remaining,
    # BlsSeries type
    BlsSeries,
    # get data
    get_data

const DEFAULT_API_URL            = "https://api.bls.gov/publicAPI/v2/timeseries/data/"
const API_KEY_LENGTH             = 32
const BLS_RESPONSE_SUCCESS       = "REQUEST_SUCCEEDED"
const BLS_RESPONSE_CATALOG_FAIL1 = "unable to get catalog data"
const BLS_RESPONSE_CATALOG_FAIL2 = "catalog has been disabled"

# See https://www.bls.gov/developers/api_faqs.htm for status code reasons.
const BLS_STATUS_CODE_REASONS    = Dict(400 => "Your request did not follow the correct syntax.",
                                        401 => "You are not authorized to make this request.",
                                        404 => "Your request was not found and/or does not exist.",
                                        429 => "You have made too many requests.",
                                        500 => "The server has encountered an unexpected condition, and the request cannot be completed.")

const LIMIT_DAILY_QUERY          = [25, 500]
const LIMIT_SERIES_PER_QUERY     = [25, 50]
const LIMIT_YEARS_PER_QUERY      = [10, 20]

# Debugging info
const DEBUG                      = true
function log_error(url, payload, headers, response)
    try
        open(joinpath(homedir(), ".blsdatajl.log"), "a") do f
            println(f, "--------")
            println(f, "Request")
            println(f, url)
            println(f, JSON.json(payload))
            println(f, headers)
            println(f, "Response")
            println(f, String(response))
        end
    catch e
        # pass
    end
end

"""
A connection to the BLS API.

Constructors
------------
* `Bls()`
* `Bls(key::AbstractString)`

Arguments
---------
* `key`: Registration key provided by the BLS.

Methods
-------
* `api_url(b::Bls)`: Get the base URL used to connect to the server
* `set_api_url!(b::Bls, url::AbstractString)`: Set the base URL used to connect to the server
* `api_key(b::Bls)`: Get the API key
* `api_version(b::Bls)`: Get the API version (v1 or v2) used
* `requests_made(b::Bls)`: Get the number of requests made today
* `requests_remaining(b::Bls)`: Get the number of requests remaining today

Notes
-----
* A default API key can be specified in a ~/.blsdatarc file.
* A valid registration key increases the allowable number of requests per day as well making
  catalog metadata available.

"""
mutable struct Bls
    url::AbstractString
    key::AbstractString
    n_requests::Int16
    t_created::Dates.DateTime
end
function Bls(key="")
    if isempty(key)
        try
            key = open(joinpath(homedir(),".blsdatarc"), "r") do f
                read(f, String)
            end
            key = rstrip(key)
            @printf "API key loaded.\n"
            # Key validation
            if length(key) > API_KEY_LENGTH
                key = key[1:API_KEY_LENGTH]
                warn("Key too long. First ", API_KEY_LENGTH, " chars used.")
            end
            if !isxdigit(key)
                error("Invalid BLS API key: ", key)
            end
        catch err
        end
    end

    url = DEFAULT_API_URL
    n_requests = 0
    t_created = Dates.now()
    Bls(url, key, n_requests, t_created)
end

get_api_url(b::Bls) = b.url
set_api_url!(b::Bls, url::AbstractString) = setfield!(b, :url, url)
get_api_key(b::Bls) = b.key
get_api_version(b::Bls) = 1 + !isempty(get_api_key(b))
requests_made(b::Bls) = b.n_requests
requests_remaining(b::Bls) = LIMIT_DAILY_QUERY[get_api_version(b)] - requests_made(b)

function increment_requests!(b::Bls)
    # Reset request if we are in a new day!
    if Dates.day(Dates.now()) ≠ Dates.day(b.t_created)
        b.t_created = Dates.now()
        b.n_requests = 0
    end

    b.n_requests += 1
end

function Base.show(io::IO, b::Bls)
    @printf io "BLS API v%d Connection\n"   get_api_version(b)
    @printf io "\turl: %s\n"                get_api_url(b)
    @printf io "\tkey: %s\n"                get_api_key(b)
    @printf io "\trequests made (this session): %d\n"       requests_made(b)
    @printf io "\trequests remaining (this session): ≤%d\n" requests_remaining(b)
end

"""
A time series with metadata returned from a `get_data` call.

For a series `s`, access fields with
```
s.id
s.data
s.catalog
```
"""
mutable struct BlsSeries
    id::AbstractString
    data::DataFrame
    catalog::AbstractString
end

function Base.show(io::IO, s::BlsSeries)
    @printf io "BlsSeries\n"
    @printf io "\tid: %s\n" s.id
    @printf io "\tdata: %dx%d DataFrame with columns %s\n" size(s.data)...  names(s.data)
    @printf io "\tcatalog: %s\n" s.catalog
end

EMPTY_RESPONSE() = BlsSeries("",DataFrame(),"")
function Base.isempty(s::BlsSeries)
    for name in propertynames(s)
        if !isempty(getproperty(s, name))
            return false
        end
    end
    return true
end

# deprecated
export
    api_url, api_key, api_version,
    id, series, catalog
@deprecate api_url(b::Bls) get_api_url(b)
@deprecate api_key(b::Bls) get_api_key(b)
@deprecate api_version(b::Bls) get_api_version(b)
@deprecate id(s::BlsSeries) getfield(s, :id)
@deprecate series(s::BlsSeries) getfield(s, :data)
@deprecate catalog(s::BlsSeries) getfield(s, :catalog)

include("get_data.jl")
end # module


text:
#Video:
function make_video(all_traj::Traj_Data,params::Params, i_iter= 300::Int)
    i_best = 1
    max_len = 0

    #Preselect good trajectory
    for i_traj = 1: params.n_traj
        if  all_traj.traj[i_traj].len > max_len
            max_len = all_traj.traj[i_traj].len
            i_best =  i_traj
        end
    end

    #Take data for video
    t_plot = 0:params.dt:params.t_horizon;
    t_plot = t_plot[1:max_len]
    x1 = all_traj.traj[i_best].s_t[1:max_len,:]
    u1 = all_traj.traj[i_best].a_t[1:max_len,:]'
    filename = string("Results/Videos/Actor_Critic_",Dates.today(),"_iter_",i_iter,".mp4")

    #plot_all_segways(all_traj, filename, params; options = Dict(:U => u1, :size => (1920, 1080)))

    #Select the Environment to visualise here!
    visualise(:InvPendulum, x1', t_plot, filename; options = Dict(:U => u1, :size => (1920, 1080)))
end

#Visualization:
#-------------------------------------------------------------------------------
@everywhere function plot_stuff(params::Params, traj::Trajectory, a_t_old::Array{Float64,1},
    weights_old::Array{Float64,1}, i_iter::Int,θ::Array{Any,1},
    res::Array{Float64,1},precision::Array{Float64,1})
    #Plotting part
    #--------------------------------------------------
    figure(1)
    clf()
    subplot(311)
    a_t_new = traj.a_t
    plot(0:params.dt:params.t_horizon,a_t_old,label="a_old");
    plot(0:params.dt:params.t_horizon,a_t_new,label="a_new");
    #plot(t,uout,label="perfect");
    xlabel("t")
    ylabel("u_motor")
    legend()
    a_t_old = a_t_new

    subplot(312)
    lab = ("x","v","dφ","φ", )

    for i_pl = 1: length(traj.s_t[1,:])
        plot(0:params.dt:params.t_horizon,traj.s_t[:,i_pl],label=lab[i_pl]);
    end
    #plot(t,x,label="perfect");
    xlabel("t")
    ylabel("s_t")
    legend()

    subplot(313)
    weights = net_to_vec(θ)
    plot(weights_old[:,1],label="θ_old");
    plot(weights,label="θ_new")
    #legend()
    ylabel("Weights")

    figure(2)
    clf()
    subplot(211)
    plot(res[1:i_iter])
    xlabel("Iterations")
    ylabel("Cost")

    subplot(212)
    plot(precision[1:i_iter])
    xlabel("Iterations")
    ylabel("Advantage-Precision")

    sleep(0.5)

    return a_t_new ,weights
    #--------------------------------------------------
end

#-------------------------------------------------------------------------------
@everywhere function plot_stuff(params::Params, all_traj::Traj_Data, a_t_old::Array{Float64,2},
    weights_old::Array{Float64,2}, i_iter::Int,θa,θv,
    res::Array{Float64,1},precision::Array{Float64,1}, ĝ_a, ĝ_v)
    #Plotting part
    #--------------------------------------------------
    #select best trajectory:
    max_len = 0
    i_best = 1
    for i_traj = 1: params.n_traj
        if max_len < all_traj.traj[i_traj].len
            i_best = i_traj
            max_len = all_traj.traj[i_traj].len
        end
    end

    len = all_traj.traj[i_best].len
    ##################################

    figure("Trajectory Infos")
    clf()
    subplot(311)
    #plot(0:params.dt:params.t_horizon,a_t_old,label="a_old");
    a_t_new = all_traj.traj[i_best].a_t
    plot(0:params.dt:params.t_horizon,a_t_new,label="a_new");
    #plot(t,uout,label="perfect");
    xlabel("t")
    ylabel("u_motor")
    #legend()
    a_t_old = copy(a_t_new)

    subplot(312)
    lab = ("v","x","dφ","φ","a","b","c","d","e","f","g","h", "i")

    for i_pl = 1: length(all_traj.traj[i_best].s_t[1,:])
        plot(0:params.dt:params.t_horizon,all_traj.traj[i_best].s_t[:,i_pl])#,label=lab[i_pl]);
    end
    xlabel("t")
    ylabel("s_t")
    #legend()

    subplot(313)
    plot(weights_old[:,1],label="θ_old")
    weights_old[:,1] = net_to_vec(θa)
    plot(weights_old[:,1],label="θ_new")
    #legend()
    ylabel("Weights-Policy")

    ##################################
    figure("Improvement: Actor improvement")
    clf()
    #subplot(211)
    plot(res[1:i_iter])
    xlabel("Iterations")
    ylabel("Cost")

    figure("Improvement: Value improvement")
    clf()
    if i_iter < 11
        start = 1
    else
        start = i_iter - 10
    end
    plot(precision[1:i_iter])
    xlabel("Iterations")
    ylabel("Value-Precision")

    sleep(0.5)

    #Plot Advantage:
    ##################################
    figure("Advantage")
    clf()
    for i_traj = 1: params.n_traj
        len = all_traj.traj[i_traj].len
        plot(all_traj.traj[i_traj].A_t[1:len])
        ylabel("Advantage")
    end

    figure("Reward_sum")
    clf()
    rew_sum = zeros(params.l_traj)
    for i_t = 1:params.l_traj
        active_traj = 0

        for i_traj = 1:params.n_traj
            if all_traj.traj[i_traj].r_t[i_t] != 0
                active_traj += 1
                rew_sum[i_t] += all_traj.traj[i_traj].r_t[i_t]
            end
        end

        if active_traj > 0
            rew_sum[i_t] /= active_traj
        end
    end
    plot(rew_sum)

    #Plot Rewards:
    ##################################
    figure("Rewards")
    clf()
    for i_traj = 1: params.n_traj
        len = all_traj.traj[i_traj].len
        plot(all_traj.traj[i_traj].r_t[1:len])
        ylabel("Costs")
    end

    len = all_traj.traj[i_best].len
    #Plot Trajectory and Advantage, Value comparison
    ##################################
    figure("State Details")
    clf()
    state_len = length(all_traj.traj[i_best].s_t[1,:])
    #single comparison of one trajectory:
    R_t = discount(all_traj.traj[i_best].r_t[1:len],params.γ_disc)

    subplot(2,1,1)
    lab = ("x","y","φz_Ball","ϕ_y","ϕ_x","ϕ_z","δx","δy","δφz_Ball","δϕ_y","δϕ_x","δϕ_z", "i","j","k")

    for i_pl = 1: Int(floor(state_len/2))
        plot(all_traj.traj[i_best].s_t[1:len,i_pl])#,label=lab[i_pl]);
    end

    subplot(2,1,2)
    for i_pl = (Int(floor(state_len/2))+1) : state_len
        plot(all_traj.traj[i_best].s_t[1:len,i_pl])#,label=lab[i_pl]);
    end

    figure("Advantage and Value")
    clf()
    plot(all_traj.traj[i_best].val_t[1:len], label="Qf_1 (half- lr)")
    plot(all_traj.traj[i_best].A_t[1:len], label="Qf_2 (double- lr)")

    plot(R_t, label="R_t")
    legend()
    return a_t_old, weights_old
    #--------------------------------------------------
end

#3D Plot of the Advantage:
function polt_surface_single(params::Params, θv::Array{Any,1}, θa::Array{Any,1})

    n = 100
    x = linspace(-2, 2, n)
    y = linspace(-2, 2,n)

    xgrid = repmat(x',n,1)
    ygrid = repmat(y,1,n)

    values = zeros(n,n)
    out    = zeros(n,n)

    for i in 1:n
        for j in 1:n
            s_t = [x[i];y[j]]
            values[i:i,j:j] =  mlp_nlayer_value(θv,s_t,params)[1]
            out[i:i,j:j]    =  mlp_nlayer_policy(θa,s_t,params)[1]
        end
    end

    fig = figure("Advantage_Visulization_separated")
    clf()
    plot_surface(xgrid,ygrid,values)
    plot_surface(xgrid,ygrid,out)
    xlabel("Velocity")
    ylabel("Position")
    title("Value and Policy Output")
end

function polt_surface_combined(params::Params, θv::Array{Any,1})

    n = 100
    x = linspace(-0.07, 0.07, n)
    y = linspace(-2, 2,n)

    xgrid = repmat(x',n,1)
    ygrid = repmat(y,1,n)

    values = zeros(n,n)
    out    = zeros(n,n)

    for i in 1:n
        for j in 1:n
            s_t = [x[i];y[j]]
            out_s, values[i:i,j:j] =  mlp_nlayer(θv,s_t,params)
            out[i:i,j:j] = out_s[1]
        end
    end

    fig = figure("Advantage_Visulization_combined")
    clf()
    plot_surface(xgrid,ygrid,values,label ="Value")
    plot_surface(xgrid,ygrid,out,label="Action")
    xlabel("Velocity")
    ylabel("Position")
    title("Value Estimation")
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
"""
Reads pcap data from a stream.
"""
mutable struct PcapStreamReader{Src <: IO} <: PcapReader
    src::Src
    raw_header::Vector{UInt8}
    header::PcapHeader
    usec_mul::Int64
    bswapped::Bool
    record_buffer::Vector{UInt8}

    @doc """
        PcapStreamReader(src::IO)

    Create reader over `src`. Will read and process pcap header,
    and yield records through `read(::PcapStreamReader)`.
    """
    function PcapStreamReader(src::Src) where {Src <: IO}
        raw_header = read(src, sizeof(PcapHeader))
        length(raw_header) != sizeof(PcapHeader) && throw(EOFError())
        h = GC.@preserve raw_header unsafe_load(Ptr{PcapHeader}(pointer(raw_header)))
        header, bswapped, nanotime = process_header(h)
        new{Src}(src, raw_header, header, nanotime ? 1 : 1000, bswapped, zeros(UInt8, 9000 + sizeof(RecordHeader)))
    end
end

"""
    PcapStreamReader(path)

Open file at `path` and create PcapStreamReader over its content.
"""
PcapStreamReader(path::AbstractString) = PcapStreamReader(open(path))

Base.close(x::PcapStreamReader) = close(x.src)
Base.position(x::PcapStreamReader) = position(x.src)
Base.seek(x::PcapStreamReader, pos) = seek(x.src, pos)
Base.mark(x::PcapStreamReader) = mark(x.src)
Base.unmark(x::PcapStreamReader) = unmark(x.src)
Base.ismarked(x::PcapStreamReader) = ismarked(x.src)
Base.reset(x::PcapStreamReader) = reset(x.src)
Base.eof(x::PcapStreamReader) = eof(x.src)

"""
    read(x::PcapStreamReader) -> PcapRecord

Read one record from pcap data. Record is valid until next read().
Throws `EOFError` if no more data available.
"""
function Base.read(x::PcapStreamReader)
    p = pointer(x.record_buffer)
    GC.@preserve x begin
        unsafe_read(x.src, p, sizeof(RecordHeader))
        h = unsafe_load(Ptr{RecordHeader}(p))
        if x.bswapped
            h = bswap(h)
        end
        unsafe_read(x.src, p + sizeof(RecordHeader), h.incl_len)
    end
    t1 = (h.ts_sec + x.header.thiszone) * 1_000_000_000
    t2 = Int64(h.ts_usec) * x.usec_mul
    t = UnixTime(Dates.UTInstant(Nanosecond(t1 + t2)))
    PcapRecord(h, t, x.record_buffer, 0)
end


text:
### DataStreams interface ###

using DataStreams

immutable DataStreamRelation{T<:Data.Source} <: AbstractRelation
    source::T
    relation::Relation
end

function convert_fieldset(source, name=:_)
    sch = Data.schema(source)
    Relation(name, ([Field{sch.types[i]}(symbol(sch.header[i]))
        for i=1:length(sch.header)]...,))
end

Relation(name::Symbol, source::Data.Source) = DataStreamRelation(source, convert_fieldset(source, name))
Relation(source::Data.Source) = DataStreamRelation(source, convert_fieldset(source))
fieldset(x::DataStreamRelation) = x.relation
Base.convert(::Type{Relation}, source::Data.Source) = DataStreamRelation(source)


text:
#
# grad_tools.jl
#
# Provide the RHAPSODIE methods to reconstruct the polarimetric
# parameters from a dataset using:
#
# - VMLM-B (when "mixed" polarimetric parameters are used)
# [Thiébaut, 2002] Thiébaut, E. (2002). Optimization issues in
# blind deconvolution algorithms. In Astronomical Data Analysis II,
# volume 4847, pages 174–183. International Society for Optics and Photonics.
#
# - Forward-Backward with backtracking (when "stokes" parameters are used)
# [Beck and Teboulle, 2009] Beck, A. and Teboulle, M. (2009).
# A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems.
# SIAM J. Imaging Sci., 2(1) :183–202. (TODO: clean implementation)
#
#
# ------------------------------------------------
#
# This file is part of Rhapsodie
#
#
# Copyright (c) 2017-2021 Laurence Denneulin (see LICENCE.md)
#

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

function apply_gradient!(X::PolarimetricMap, A::D, g::Array{T,3}, d::Array{data_table,1}, μ::Array{Tuple{H,T},1}) where {T <: AbstractFloat, H<:AbstractCost, D <:Mapping}

    n1, n2, n3 = size(g)
    @assert (n1,n2) == size(X)
    @assert n3 == 3
    if X.parameter_type == "intensities"
        error("Global reconstruction not implemented on Iu, Ip, and θ.
               Only use 'stokes' or 'mixed' parameters.")
    end

    Ay = cat(A*X.I[:,:], A*X.Q[:,:], A*X.U[:,:], dims=3)
    # Compute data fidelity term and gradient. (As gradient is initially set to
    # zero, we can recycle it between x and y.)
    @assert size(g) == size(Ay)
    vfill!(g, 0)
    local f::Float64 = 0.0;
    for k = 1:length(d)
        f += fg!(Ay, g, d[k])
    end

    # Convert gradient w.r.t. y into gradient w.r.t. x.  Nothing has to be done
    # for the 2nd and 3rd fields (Q and U) or if Ip = 0.
    @assert size(g) == (n1,n2,3)
    @inbounds for i3 in 1:3
        g[:,:,i3].= A'*view(g,:,:,i3)[:,:];
    end

    if X.parameter_type == "mixed"
        @inbounds for i2 in 1:n2
            for i1 in 1:n1
                if X.Ip[i1,i2] > 0
                    g[i1,i2,1] += (X.Q[i1,i2]*g[i1,i2,2] +
                                   X.U[i1,i2]*g[i1,i2,3])/X.Ip[i1,i2]
                end
            end
        end
 	    f+=cost!(μ[1][2] , μ[1][1], X.Iu[:,:], view(g,:,:,1), false);
     elseif X.parameter_type == "stokes"
 	    f+=cost!(μ[1][2] , μ[1][1], X.I[:,:], view(g,:,:,1), false);
     end

	f+=cost!(μ[2][2] , μ[2][1], cat(X.Q[:,:], X.U[:,:], dims=3), view(g,:,:,2:3), false);

	return f
end

function apply_rhapsodie(x0::PolarimetricMap, A::D, d::Array{data_table,1}, par::Array{T,1}; mem=3, maxeval=50, maxiter=50, xtol=(1e-3,1e-8), gtol=(1e-3,1e-8), ftol=(1e-3,1e-8)) where {T <: AbstractFloat, D <:Mapping}
    n1,n2 = size(x0)
    X0 = convert(Array{T,3},x0);
    μ=[(HyperbolicEdgePreserving(10. ^par[3], (0.5,0.5)),10. ^par[1]);
       (HyperbolicEdgePreserving(10. ^par[4],(0.5,0.5,0.)),10. ^par[2])];

    lower_born=vcreate(X0);
    vfill!(view(lower_born,:,:,1),0.0)
    vfill!(view(lower_born,:,:,2:3),-Inf)

    g=vcreate(X0);
    fg!(x,g)=apply_gradient!(PolarimetricMap(x0.parameter_type,x), A, g, d, μ)
    x = vmlmb(fg!, X0, mem=mem, maxeval=maxeval, maxiter=maxiter, lower=lower_born, xtol=xtol,  gtol=gtol, ftol=ftol);
    return PolarimetricMap(x0.parameter_type,x)
end

max_stars_repo_name: LaurenceDenneulin/Rhapsodie


text:
using Submodular
using FactCheck

Tol = 1e-3

facts("Continuous Sets") do

  context("permutation atom") do
    p = Perm(4)
    d1 = randperm(4)
    @fact in(d1, p) --> true
    d2 = collect(2:5)
    @fact in(d2, p) --> false
    w = [.4,.2,.6,.3]
    @fact fenchel(p, w) --> roughly([2, 4, 1, 3], Tol)
  end

  context("associated polyhedra and the Fujishige-minimum-norm-point algorithm") do
    S = SetVariable(4)
    p(z) = -0.5*z^2 + 4*z + 0.5 * z
    perm_func = compose(p, card(S))

    P₁ = SubmodPoly(perm_func)
    @fact in([1, 2, 5, 3], P₁) --> false
    @fact in([1, 2, 4, 3], P₁) --> true
    @fact in([1, 2, 3, 3], P₁) --> true
    @fact in([-1, -2, 0, 0], P₁) --> true
    @fact in([-2, -2, -5, -3], P₁) --> true

    P₂ = BasePoly(perm_func)
    @fact in([1, 2, 5, 3], P₂) --> false
    @fact in([1, 2, 4, 3], P₂) --> true
    @fact in([1, 2, 3, 3], P₂) --> false
    @fact in([-1, -2, 0, 0], P₂) --> false
    @fact in([-2, -2, -5, -3], P₂) --> false

    P₃ = PosPoly(perm_func)
    @fact in([1, 2, 5, 3], P₃) --> false
    @fact in([1, 2, 4, 3], P₃) --> true
    @fact in([1, 2, 3, 3], P₃) --> true
    @fact in([-1, -2, 0, 0], P₃) --> false
    @fact in([-2, -2, -5, -3], P₃) --> false

    P₄ = SymPoly(perm_func)
    @fact in([1, 2, 5, 3], P₄) --> false
    @fact in([1, 2, 4, 3], P₄) --> true
    @fact in([1, 2, 3, 3], P₄) --> true
    @fact in([-1, -2, 0, 0], P₄) --> true
    @fact in([-2, -2, -5, -3], P₄) --> false
  end

end


text:
struct SystemMatrix{T,Z}
    rows::Vector{Z}
    cols::Vector{Z}
    vals::Vector{T}
    function SystemMatrix(rows::Vector{Z},cols::Vector{Z},
        vals::Vector{T}) where {T<:Real,Z<:Integer}

        @assert length(rows) == length(cols)
        @assert length(cols) == length(vals)

        new{T,Z}(rows,cols,vals)
    end
end

function SystemMatrix()
    Z = default_integer_type()
    R = default_float_type()
    rows = Z[]
    cols = Z[]
    vals = R[]
    return SystemMatrix(rows,cols,vals)
end

function assemble!(matrix::SystemMatrix,rows,cols,vals)
    @assert length(rows) == length(cols)
    @assert length(cols) == length(vals)

    append!(matrix.rows,rows)
    append!(matrix.cols,cols)
    append!(matrix.vals,vals)
end

struct SystemRHS{T,Z}
    rows::Vector{Z}
    vals::Vector{T}
    function SystemRHS(rows::Vector{Z},vals::Vector{T}) where {Z<:Integer,T<:Real}
        @assert length(rows) == length(vals)
        new{T,Z}(rows,vals)
    end
end

function SystemRHS()
    Z = default_integer_type()
    R = default_float_type()
    rows = Z[]
    vals = R[]
    return SystemRHS(rows,vals)
end

function assemble!(rhs::SystemRHS,rows,vals)
    @assert length(rows) == length(vals)
    append!(rhs.rows,rows)
    append!(rhs.vals,vals)
end

function element_dof_start(elid,dofsperelement)
    return (elid-1)*dofsperelement+1
end

function element_dof_stop(elid,dofsperelement)
    return elid*dofsperelement
end

function element_dofs(elid,dofsperelement)
    start = element_dof_start(elid,dofsperelement)
    stop = element_dof_stop(elid,dofsperelement)
    return start:stop
end

function operator_dofs(row_dofs,col_dofs)
    lr = length(row_dofs)
    lc = length(col_dofs)
    rows = repeat(row_dofs,outer=lc)
    cols = repeat(col_dofs,inner=lr)
    return rows,cols
end

function assemble!(system_matrix::SystemMatrix,vals,rowelids,
    colelids,dofsperelement)

    rowdofs = vcat([element_dofs(r,dofsperelement) for r in rowelids]...)
    coldofs = vcat([element_dofs(c,dofsperelement) for c in colelids]...)

    oprows,opcols = operator_dofs(rowdofs,coldofs)
    assemble!(system_matrix,oprows,opcols,vals)
end

function assemble!(system_rhs::SystemRHS,vals,rowelid,dofsperelement)
    rowdofs = element_dofs(rowelid,dofsperelement)
    assemble!(system_rhs,rowdofs,vals)
end

function assemble_displacement_face!(system_matrix::SystemMatrix,sbasis,
    facequad,facescale,elid,dofsperelement)

    HH = HHop(sbasis,facequad,facescale)
    assemble!(system_matrix,vec(HH),elid,elid,dofsperelement)
end

function assemble_displacement_face!(system_matrix::SystemMatrix,
    dgmesh::DGMesh,ufs::UniformFunctionSpace,phaseid,cellid,
    faceid,rowelid,dofsperelement)

    sbasis = ufs.sbasis
    facequad = ufs.fquads[phaseid,cellid][faceid]
    facescale = dgmesh.facescale[faceid]
    assemble_displacement_face!(system_matrix,sbasis,facequad,
        facescale,rowelid,dofsperelement)
end

function assemble_traction_face!(system_matrix::SystemMatrix,HL,iLLxLH,HH,
    rowelid,colelids,dofsperelement)

    tractionop = HL*iLLxLH
    assemble!(system_matrix,tractionop,rowelid,colelids,dofsperelement)
    assemble!(system_matrix,-HH,rowelid,rowelid,dofsperelement)
end

function assemble_traction_face!(system_matrix::SystemMatrix,sbasis,facequad,
    facescale,stabilization,HL,iLLxLH,rowelid,colelids,dofsperelement)

    HH = stabilization*HHop(sbasis,facequad,facescale)
    assemble_traction_face!(system_matrix,HL,iLLxLH,HH,rowelid,
        colelids,dofsperelement)
end

function assemble_traction_face!(system_rhs::SystemRHS,sbasis,facequad,
    facescale,rhsval,rowelid,dofsperelement)

    rhs = -facescale*linear_form(rhsval,sbasis,facequad)
    assemble!(system_rhs,rhs,rowelid,dofsperelement)
end

function assemble_traction_face!(system_rhs::SystemRHS,dgmesh::DGMesh,
    ufs::UniformFunctionSpace,rhsval,facetohelid,phaseid,cellid,faceid)

    rowelid = facetohelid[phaseid,cellid][faceid]
    facequad = ufs.fquads[phaseid,cellid][faceid]
    facescale = dgmesh.facescale[faceid]
    dofsperelement = ufs.dofsperelement
    assemble_traction_face!(system_rhs,ufs.sbasis,facequad,facescale,
        rhsval,rowelid,dofsperelement)
end

function assemble_displacement_face!(system_rhs::SystemRHS,dgmesh,
    ufs,rhsval,facetohelid,phaseid,cellid,faceid)

    assemble_traction_face!(system_rhs,dgmesh,ufs,-rhsval,facetohelid,
        phaseid,cellid,faceid)
end

function assemble_mixed_face!(system_rhs::SystemRHS,dgmesh,ufs,dcomp,dval,
    tcomp,tval,facetohelid,phaseid,cellid,faceid)

    assemble_displacement_face!(system_rhs,dgmesh,ufs,dval*dcomp,
        facetohelid,phaseid,cellid,faceid)
    assemble_traction_face!(system_rhs,dgmesh,ufs,tval*tcomp,
        facetohelid,phaseid,cellid,faceid)
end

function assemble_traction_face!(system_matrix::SystemMatrix,
    dgmesh::DGMesh,ufs::UniformFunctionSpace,cellsolvers::CellSolvers,
    phaseid,cellid,faceid,rowelid,colelids,dofsperelement)

    sbasis = ufs.sbasis
    facequad = ufs.fquads[phaseid,cellid][faceid]
    facescale = dgmesh.facescale[faceid]
    cellsolver = cellsolvers[phaseid,cellid]
    stabilization = cellsolvers.stabilization
    iLLxLH = cellsolver.iLLxLH
    facetosolverid = cellsolver.facetosolverid
    HL = cellsolver.fLH[facetosolverid[faceid]]'
    assemble_traction_face!(system_matrix,sbasis,facequad,
        facescale,stabilization,HL,iLLxLH,rowelid,colelids,dofsperelement)
end

function assemble_traction_coherent_interface!(system_matrix::SystemMatrix,
    HL1,iLLxLH1,HL2,iLLxLH2,HH,hid1,colelids1,hid2,colelids2,dofsperelement)

    op1 = HL1*iLLxLH1
    op2 = HL2*iLLxLH2

    assemble!(system_matrix,op1,hid1,colelids1,dofsperelement)
    assemble!(system_matrix,-HH,hid1,hid1,dofsperelement)
    assemble!(system_matrix,op2,hid1,colelids2,dofsperelement)
    assemble!(system_matrix,-HH,hid1,hid2,dofsperelement)

    assemble!(system_matrix,op1,hid2,colelids1,dofsperelement)
    assemble!(system_matrix,-HH,hid2,hid1,dofsperelement)
    assemble!(system_matrix,op2,hid2,colelids2,dofsperelement)
    assemble!(system_matrix,-HH,hid2,hid2,dofsperelement)
end

function assemble_displacement_coherent_interface!(system_matrix::SystemMatrix,
    HH,rowelid,colelid,dofsperelement)

    vals = vec(HH)
    assemble!(system_matrix,vals,rowelid,rowelid,dofsperelement)
    assemble!(system_matrix,-vals,rowelid,colelid,dofsperelement)
    assemble!(system_matrix,-vals,colelid,rowelid,dofsperelement)
    assemble!(system_matrix,vals,colelid,colelid,dofsperelement)
end

function assemble_coherent_interface!(system_matrix::SystemMatrix,
    dgmesh::DGMesh,ufs::UniformFunctionSpace,cellsolvers::CellSolvers,
    cellid,hid1,colelids1,hid2,colelids2,dofsperelement)

    sbasis = ufs.sbasis
    iquad = ufs.iquad
    imap = ufs.imap
    inormals = ufs.inormals[cellid]
    stabilization = cellsolvers.stabilization
    update!(imap,ufs.icoeffs[cellid])
    HH = HHop(sbasis,iquad,imap,inormals,dgmesh.cellmap)

    cs1 = cellsolvers[1,cellid]
    iLLxLH1 = cs1.iLLxLH
    facetosolverid = cs1.facetosolverid
    HL1 = cs1.fLH[facetosolverid[5]]'

    cs2 = cellsolvers[2,cellid]
    iLLxLH2 = cs2.iLLxLH
    facetosolverid = cs2.facetosolverid
    HL2 = cs2.fLH[facetosolverid[5]]'

    assemble_displacement_coherent_interface!(system_matrix,HH,
        hid1,hid2,dofsperelement)
    assemble_traction_coherent_interface!(system_matrix,
        HL1,iLLxLH1,HL2,iLLxLH2,stabilization*HH,hid1,colelids1,
        hid2,colelids2,dofsperelement)
end

function assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    hybrid_element_numbering)

    facetohelid = hybrid_element_numbering.facetohelid
    interfacehelid = hybrid_element_numbering.interfacehelid
    cellids = findall(x->x==0,dgmesh.cellsign)
    dofsperelement = ufs.dofsperelement

    for cellid in cellids
        hid1 = interfacehelid[1,cellid]
        hid2 = interfacehelid[2,cellid]

        colelids1 = filter(x->x!=0,facetohelid[1,cellid])
        push!(colelids1,hid1)
        colelids2 = filter(x->x!=0,facetohelid[2,cellid])
        push!(colelids2,hid2)

        assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
            cellid,hid1,colelids1,hid2,colelids2,dofsperelement)
    end
end

function assemble_displacement_continuity_incoherent_interface!(system_matrix,
    sbasis,iquad,imap,inormals,cellmap,rowelid,colelid,dofsperelement)

    HH = HHop(sbasis,iquad,imap,inormals,inormals,cellmap)
    assemble_displacement_coherent_interface!(system_matrix,HH,
        rowelid,colelid,dofsperelement)
end

function assemble_free_slip_incoherent_interface!(system_matrix,vbasis,
    sbasis,iquad,imap,inormals,Dhalf,stabilization,cellmap,iLLxLH,
    rowelid,colelids,dofsperelement)

    facescale = scale_area(cellmap,inormals)
    itangents = tangents(inormals)
    HLt = hybrid_local_operator_traction_components(sbasis,vbasis,
        iquad,imap,inormals,itangents,Dhalf,stabilization,facescale)
    HHt = stabilization*HHop(sbasis,iquad,imap,inormals,itangents,cellmap)
    assemble_traction_face!(system_matrix,HLt,iLLxLH,HHt,rowelid,colelids,
        dofsperelement)
end

function assemble_traction_continuity_incoherent_interface!(system_matrix,
    vbasis,sbasis,iquad,imap,inormals,D1,D2,stabilization,cellmap,
    iLLxLH1,iLLxLH2,hid1,colelids1,hid2,colelids2,dofsperelement)

    facescale = scale_area(cellmap,inormals)
    HL1 = hybrid_local_operator_traction_components(sbasis,vbasis,iquad,
        imap,-inormals,inormals,D1,stabilization,facescale)
    HL2 = hybrid_local_operator_traction_components(sbasis,vbasis,iquad,
        imap,inormals,inormals,D2,stabilization,facescale)

    HH = stabilization*HHop(sbasis,iquad,imap,inormals,inormals,cellmap)

    assemble_traction_coherent_interface!(system_matrix,HL1,iLLxLH1,HL2,
        iLLxLH2,HH,hid1,colelids1,hid2,colelids2,dofsperelement)
end

function assemble_incoherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    cellid,hid1,colelids1,hid2,colelids2,dofsperelement)

    vbasis = ufs.vbasis
    sbasis = ufs.sbasis
    iquad = ufs.iquad
    imap = ufs.imap
    inormals = ufs.inormals[cellid]
    dofsperelement = ufs.dofsperelement
    stabilization = cellsolvers.stabilization
    D1,D2 = cellsolvers.stiffness
    cellmap = dgmesh.cellmap
    update!(imap,ufs.icoeffs[cellid])

    ls1 = cellsolvers[1,cellid]
    ls2 = cellsolvers[2,cellid]

    iLLxLH1 = ls1.iLLxLH
    iLLxLH2 = ls2.iLLxLH

    assemble_displacement_continuity_incoherent_interface!(system_matrix,
        sbasis,iquad,imap,inormals,cellmap,hid1,hid2,dofsperelement)

    assemble_free_slip_incoherent_interface!(system_matrix,vbasis,sbasis,
        iquad,imap,-inormals,D1,stabilization,cellmap,iLLxLH1,hid1,colelids1,
        dofsperelement)
    assemble_free_slip_incoherent_interface!(system_matrix,vbasis,sbasis,
        iquad,imap,inormals,D2,stabilization,cellmap,iLLxLH2,hid2,colelids2,
        dofsperelement)

    assemble_traction_continuity_incoherent_interface!(system_matrix,vbasis,
        sbasis,iquad,imap,inormals,D1,D2,stabilization,cellmap,
        iLLxLH1,iLLxLH2,hid1,colelids1,hid2,colelids2,dofsperelement)
end

function assemble_mixed_face!(system_matrix::SystemMatrix,vbasis,sbasis,
    facequad,facemap,facenormal,dcomp,tcomp,Dhalf,stabilization,facescale,
    iLLxLH,rowelid,colelids,dofsperelement)

    HHd = HHop(sbasis,facequad,dcomp,facescale)
    assemble!(system_matrix,vec(HHd),rowelid,rowelid,dofsperelement)

    HLt = hybrid_local_operator_traction_components(sbasis,vbasis,facequad,
        facemap,facenormal,tcomp,Dhalf,stabilization,facescale)
    HHt = stabilization*HHop(sbasis,facequad,tcomp,facescale)
    assemble_traction_face!(system_matrix,HLt,iLLxLH,HHt,rowelid,colelids,
        dofsperelement)
end

function assemble_mixed_face!(system_matrix::SystemMatrix,
    dgmesh::DGMesh,ufs::UniformFunctionSpace,cellsolvers::CellSolvers,
    phaseid,cellid,faceid,dcomp,tcomp,rowelid,colelids,dofsperelement)

    vbasis = ufs.vbasis
    sbasis = ufs.sbasis
    facequad = ufs.fquads[phaseid,cellid][faceid]
    facemap = dgmesh.facemaps[faceid]
    facenormal = ufs.fnormals[faceid]
    facescale = dgmesh.facescale[faceid]
    iLLxLH = cellsolvers[phaseid,cellid].iLLxLH
    D = cellsolvers.stiffness[phaseid]
    stabilization = cellsolvers.stabilization
    assemble_mixed_face!(system_matrix,vbasis,sbasis,facequad,facemap,
        facenormal,dcomp,tcomp,D,stabilization,facescale,iLLxLH,
        rowelid,colelids,dofsperelement)
end

function SparseArrays.sparse(system_matrix::SystemMatrix,ndofs)
    return sparse(system_matrix.rows,system_matrix.cols,system_matrix.vals,
        ndofs,ndofs)
end

function rhs(system_rhs::SystemRHS,ndofs)
    return Array(sparsevec(system_rhs.rows,system_rhs.vals,ndofs))
end

function assign_cell_hybrid_element_ids!(facetohelid,visited,
    isactiveface,helid,nfaces)

    for faceid in 1:nfaces
        if !visited[faceid] && isactiveface[faceid]
            facetohelid[faceid] = helid
            visited[faceid] = true
            helid += 1
        end
    end
    return helid
end

function assign_neighbor_cell_hybrid_ids!(facetohelid,visited,phaseid,cellid,
    connectivity,nfaces)

    for faceid = 1:nfaces
        nbrcellid,nbrfaceid = connectivity[cellid][faceid]
        if nbrcellid != 0 && !visited[phaseid,nbrcellid][nbrfaceid]
            facetohelid[phaseid,nbrcellid][nbrfaceid] = facetohelid[phaseid,cellid][faceid]
            visited[phaseid,nbrcellid][nbrfaceid] = true
        end
    end
end

function cell_hybrid_element_ids(cellsign,isactiveface,connectivity,nfaces)
    ncells = length(cellsign)
    facetohelid = [zeros(Int,nfaces) for i = 1:2, j = 1:ncells]
    interfacehelid = zeros(Int,2,ncells)
    visited = [zeros(Bool,nfaces) for i = 1:2, j = 1:ncells]
    helid = 1
    for cellid = 1:ncells
        s = cellsign[cellid]
        if s == 1
            helid = assign_cell_hybrid_element_ids!(facetohelid[1,cellid],
                visited[1,cellid],isactiveface[1,cellid],helid,nfaces)
            assign_neighbor_cell_hybrid_ids!(facetohelid,visited,1,cellid,
                connectivity,nfaces)
        elseif s == -1
            helid = assign_cell_hybrid_element_ids!(facetohelid[2,cellid],
                visited[2,cellid],isactiveface[2,cellid],helid,nfaces)
            assign_neighbor_cell_hybrid_ids!(facetohelid,visited,2,cellid,
                connectivity,nfaces)
        elseif s == 0
            helid = assign_cell_hybrid_element_ids!(facetohelid[1,cellid],
                visited[1,cellid],isactiveface[1,cellid],helid,nfaces)
            assign_neighbor_cell_hybrid_ids!(facetohelid,visited,1,cellid,
                connectivity,nfaces)
            interfacehelid[1,cellid] = helid
            helid += 1
            helid = assign_cell_hybrid_element_ids!(facetohelid[2,cellid],
                visited[2,cellid],isactiveface[2,cellid],helid,nfaces)
            assign_neighbor_cell_hybrid_ids!(facetohelid,visited,2,cellid,
                connectivity,nfaces)
            interfacehelid[2,cellid] = helid
            helid += 1
        else
            error("Expected cellsign ∈ {-1,0,1}, got cellsign = $s")
        end
    end
    return facetohelid,interfacehelid,helid
end

struct HybridElementNumbering
    facetohelid
    interfacehelid
    number_of_hybrid_elements
end

function HybridElementNumbering(dgmesh,ufs)
    dim = dimension(dgmesh)
    nfaces = number_of_faces(dim)
    facetohelid,interfacehelid,helidstop = cell_hybrid_element_ids(
        dgmesh.cellsign,ufs.isactiveface,dgmesh.connectivity,nfaces
    )
    number_of_hybrid_elements = helidstop - 1
    return HybridElementNumbering(facetohelid,interfacehelid,
        number_of_hybrid_elements)
end

function interior_operator(ls)
    return ls.LH'*ls.iLLxLH
end

function assemble_uniform_interior_faces!(system_matrix::SystemMatrix,dgmesh,
    ufs,cellsolvers,facetohelid)

    cellsign = dgmesh.cellsign
    isinteriorcell = dgmesh.isinteriorcell
    dofsperelement = ufs.dofsperelement
    vals1 = vec(interior_operator(cellsolvers[1]))
    vals2 = vec(interior_operator(cellsolvers[2]))
    ncells = length(cellsign)
    cellids = findall(isinteriorcell)

    for cellid in cellids
        s = cellsign[cellid]
        if s == 1
            colelids = facetohelid[1,cellid]
            @assert isnothing(findfirst(x->x==0,colelids))
            for (faceid,rowelid) in enumerate(colelids)
                assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
                    1,cellid,faceid,rowelid,colelids,dofsperelement)
            end
        elseif s == -1
            colelids = facetohelid[2,cellid]
            @assert isnothing(findfirst(x->x==0,colelids))
            for (faceid,rowelid) in enumerate(colelids)
                assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
                    2,cellid,faceid,rowelid,colelids,dofsperelement)
            end
        end
    end
end

function assemble_cut_interior_faces!(system_matrix::SystemMatrix,dgmesh,ufs,
    cellsolvers,facetohelid,interfacehelid,phaseid,cellid)

    dofsperelement = ufs.dofsperelement
    fhelid = facetohelid[phaseid,cellid]
    ihelid = interfacehelid[phaseid,cellid]
    colelids = filter(x->x!=0,fhelid)
    push!(colelids,ihelid)

    for (faceid,helid) in enumerate(fhelid)
        if helid != 0
            assemble_traction_face!(system_matrix,
                dgmesh,ufs,cellsolvers,phaseid,cellid,faceid,
                helid,colelids,dofsperelement)
        end
    end
end

function assemble_cut_interior_faces!(system_matrix::SystemMatrix,dgmesh,ufs,
    cellsolvers,facetohelid,interfacehelid)

    cellsign = dgmesh.cellsign
    isinteriorcell = dgmesh.isinteriorcell
    cellids = findall([x==0 for x in cellsign] .& isinteriorcell)

    for cellid in cellids
        assemble_cut_interior_faces!(system_matrix,dgmesh,ufs,
            cellsolvers,facetohelid,interfacehelid,1,cellid)
        assemble_cut_interior_faces!(system_matrix,dgmesh,ufs,
            cellsolvers,facetohelid,interfacehelid,2,cellid)
    end
end

function assemble_boundary_face!(system_matrix::SystemMatrix,dgmesh,ufs,
    cellsolvers,phaseid,cellid,faceid,midpoint,rowelid,colelids,
    bc_classifier,bc_data)

    dofsperelement = ufs.dofsperelement
    bdrytype = bc_classifier(midpoint)
    if bdrytype == :displacement
        assemble_displacement_face!(system_matrix,dgmesh,ufs,
            phaseid,cellid,faceid,rowelid,dofsperelement)
    elseif bdrytype == :traction
        assemble_traction_face!(system_matrix,dgmesh,ufs,
            cellsolvers,phaseid,cellid,faceid,rowelid,
            colelids,dofsperelement)
    elseif bdrytype == :mixed
        dcomp,tcomp = bc_data[:mixed_components](midpoint)
        assemble_mixed_face!(system_matrix,dgmesh,ufs,
            cellsolvers,phaseid,cellid,faceid,dcomp,tcomp,
            rowelid,colelids,dofsperelement)
    else
        error("Unexpected boundary type encountered")
    end
end

function assemble_uniform_boundary_faces!(system_matrix::SystemMatrix,dgmesh,
    ufs,cellsolvers,facetohelid,phaseid,cellid,bc_classifier,bc_data)

    cellmap = CellMap(dgmesh.domain[cellid])
    midpoints = face_midpoints(cellmap)
    colelids = facetohelid[phaseid,cellid]
    connectivity = dgmesh.connectivity[cellid]
    dofsperelement = ufs.dofsperelement

    for (faceid,xm) in enumerate(midpoints)
        rowelid = colelids[faceid]
        if connectivity[faceid][1] != 0
            assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
                phaseid,cellid,faceid,rowelid,colelids,dofsperelement)
        else
            assemble_boundary_face!(system_matrix,dgmesh,ufs,cellsolvers,
                phaseid,cellid,faceid,xm,rowelid,colelids,bc_classifier,
                bc_data)
        end
    end
end

function assemble_uniform_boundary_faces!(system_matrix::SystemMatrix,dgmesh,
    ufs,cellsolvers,facetohelid,bc_classifier,bc_data)

    cellids = findall(.!dgmesh.isinteriorcell)
    cellsign = dgmesh.cellsign
    for cellid in cellids
        if cellsign[cellid] == 1
            assemble_uniform_boundary_faces!(system_matrix,dgmesh,
            ufs,cellsolvers,facetohelid,1,cellid,bc_classifier,bc_data)
        elseif cellsign[cellid] == -1
            assemble_uniform_boundary_faces!(system_matrix,dgmesh,
            ufs,cellsolvers,facetohelid,2,cellid,bc_classifier,bc_data)
        end
    end
end

function assemble_cut_boundary_faces!(system_matrix::SystemMatrix,dgmesh,
    ufs,cellsolvers,facetohelid,interfacehelid,phaseid,cellid,
    bc_classifier,bc_data)

    dofsperelement = ufs.dofsperelement
    fhelid = facetohelid[phaseid,cellid]
    ihelid = interfacehelid[phaseid,cellid]
    connectivity = dgmesh.connectivity[cellid]

    colelids = filter(x->x!=0,fhelid)
    push!(colelids,ihelid)
    cellmap = CellMap(dgmesh.domain[cellid])
    midpoints = face_midpoints(cellmap)

    for (faceid,rowelid) in enumerate(fhelid)
        if rowelid != 0
            if connectivity[faceid][1] != 0
                assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
                    phaseid,cellid,faceid,rowelid,colelids,dofsperelement)
            else
                assemble_boundary_face!(system_matrix,dgmesh,ufs,cellsolvers,
                    phaseid,cellid,faceid,midpoints[faceid],rowelid,colelids,
                    bc_classifier,bc_data)
            end
        end
    end
end

function assemble_cut_boundary_faces!(system_matrix::SystemMatrix,dgmesh,
    ufs,cellsolvers,facetohelid,interfacehelid,bc_classifier,bc_data)

    cellids = findall(.!dgmesh.isinteriorcell)
    cellsign = dgmesh.cellsign
    for cellid in cellids
        if cellsign[cellid] == 0
            assemble_cut_boundary_faces!(system_matrix,dgmesh,ufs,cellsolvers,
                facetohelid,interfacehelid,1,cellid,bc_classifier,bc_data)
            assemble_cut_boundary_faces!(system_matrix,dgmesh,ufs,cellsolvers,
                facetohelid,interfacehelid,2,cellid,bc_classifier,bc_data)
        end
    end
end

function assemble_boundary_face!(system_rhs::SystemRHS,dgmesh,ufs,
    facetohelid,phaseid,cellid,faceid,midpoint,bc_classifier,bc_data)

    bdrytype = bc_classifier(midpoint)
    if bdrytype == :displacement
        dispval = bc_data[:displacement](midpoint)
        assemble_displacement_face!(system_rhs,dgmesh,ufs,dispval,
            facetohelid,phaseid,cellid,faceid)
    elseif bdrytype == :traction
        tractionval = bc_data[:traction](midpoint)
        assemble_traction_face!(system_rhs,dgmesh,ufs,tractionval,
            facetohelid,phaseid,cellid,faceid)
    elseif bdrytype == :mixed
        dcomp,tcomp = bc_data[:mixed_components](midpoint)
        dval,tval = bc_data[:mixed_values](midpoint)
        assemble_mixed_face!(system_rhs,dgmesh,ufs,dcomp,dval,tcomp,tval,
            facetohelid,phaseid,cellid,faceid)
    else
        error("Unexpected boundary type encountered")
    end
end

function assemble_boundary_faces!(system_rhs::SystemRHS,dgmesh,ufs,
    facetohelid,bc_classifier,bc_data)

    dim = dimension(dgmesh)
    nfaces = number_of_faces(dim)
    connectivity = dgmesh.connectivity
    cellids = findall(.!dgmesh.isinteriorcell)
    for cellid in cellids
        cellmap = CellMap(dgmesh.domain[cellid])
        midpoints = face_midpoints(cellmap)
        for phaseid in 1:2
            for faceid in 1:nfaces
                if facetohelid[phaseid,cellid][faceid] != 0
                    midpoint = midpoints[faceid]
                    assemble_boundary_face!(system_rhs,dgmesh,ufs,
                        facetohelid,phaseid,cellid,faceid,midpoint,
                        bc_classifier,bc_data)
                end
            end
        end
    end
end


text:
if !@isdefined USE_GPU
    const USE_GPU = length(ARGS) > 0 && ARGS[1] == "gpu"
end

using Test
using LinearAlgebra

using ParallelStencil
if !ParallelStencil.is_initialized()
    @static if USE_GPU
        @init_parallel_stencil(CUDA, Float64, 2)
    else
        @init_parallel_stencil(Threads, Float64, 2)
    end
end

include("../scripts-part2/krylov.jl")

@testset "Test Krylov solver" for execution_policy in [parallel, parallel_shmem]
    @info "Runnint with " execution_policy
    n = 66
    h = 1 / (n - 1)
    c = 3.14
    tol = 1e-6
    Nmax = 1000

    b = @ones(n, n)
    b[1, :] .= 0.0
    b[n, :] .= 0.0
    b[:, 1] .= 0.0
    b[:, n] .= 0.0
    xhat = @zeros(n, n)
    res_rms = cg!(xhat, b, h, h, c, tol, Nmax, execution_policy=execution_policy)

    @test res_rms < tol * sqrt(sum(b.^2) / n^2)
end

@ParallelStencil.reset_parallel_stencil()


text:
@testset "1567.maximum-length-of-subarray-with-positive-product.jl" begin
    @test get_max_len([1, -2, -3, 4]) == 4
    @test get_max_len([0, 1, -2, -3, -4]) == 3
    @test get_max_len([-1, -2, -3, 0, 1]) == 2
end


text:
using Documenter, GridapDistributed

makedocs(;
    modules=[GridapDistributed],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/gridap/GridapDistributed.jl/blob/{commit}{path}#L{line}",
    sitename="GridapDistributed.jl",
    authors="Francesc Verdugo <fverdugo@cimne.upc.edu>",
    assets=String[],
)

deploydocs(;
    repo="github.com/gridap/GridapDistributed.jl",
)


text:
function makescatteringnoisetrans(M,ω=1.,basis="Hermite",Nk=M)
#Noise transform: 2-field quadrature to go from basis to k-space
kn,wkn,Tkn = nfieldtrans(Nk,3,ω=1/ω,basis=basis)
Tkn = eigmat(M,kn,ω=1/ω,basis=basis)*diagm((-im).^(0:M-1))

#Noise: position space roots and weights for 3-field quadrature
xn,wxn,Txn = nfieldtrans(M,3,ω=ω,basis=basis)
Txkn = conj(Tkn'*Txn)
return xn,wxn,Txn,kn,wkn,Tkn
end


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

points = DataFrame(index=rand(0:10,30), val=rand(1:10,30))
line = DataFrame(val=rand(1:10,11), index = collect(0:10))
pointLayer = layer(points, x="index", y="val", Geom.point,Theme(default_color="green"))
lineLayer = layer(line, x="index", y="val", Geom.line)
plot(pointLayer, lineLayer, Guide.manual_color_key("Legend", ["Points", "Line"], ["green", "deepskyblue"]))


text:
"""
    Faber(i, j)

change of indexation: from double indexation i,j to single indexation n
"""
function Faber(i::Int64, j::Int64)
    return Int(2<<(i-1) + j)
end

"""
    binlog(convert(UInt32, v))

efficient log_2 of integers
https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog
"""
binlog(v::Integer) = binlog(convert(UInt32, v))
function binlog(v::UInt32)
    r =     (v > 0xFFFF) << 4; v >>= r;
    shift = (v > 0xFF  ) << 3; v >>= shift; r |= shift;
    shift = (v > 0xF   ) << 2; v >>= shift; r |= shift;
    shift = (v > 0x3   ) << 1; v >>= shift; r |= shift;
                                            r |= (v >> 1);
end

"""
    Faber(n)

change of indexation: from single indexation n to double indexation i,j
"""
function Faber(n)
    return (binlog(n)), Int(n - 1<<(binlog(n)))
end

"""
    OrderBasis(i, j)

compare two index (single indexation) and order it from the smaller
"""
function OrderBasis(i, j)
    if i <= j
        return (i, j)
    else
        return (j, i)
    end
end

"""
    OrderBasis(i1, j1, i2, j2)

compare two index (double indexation) and order it from the smaller
"""
function OrderBasis(i1 ,j1 ,i2 ,j2)
    if i1 <= i2
        return (i1, j1, i2, j2)
    else
        return (i2, j2, i1, j1)
    end
end
"""
    ShareDomain(i,j)

check if the basis relative to index i and index j (single indexation) share the same support
return 1 if they do
return 0 if they don't
"""
function ShareDomain(i,j)
    if i == j return 1
    else
        (i, j) = OrderBasis(i, j)
        i1, j1 = Faber(i)
        i2, j2 = Faber(j)
        if j1*2<<(i2 - i1 - 1) <= j2 < (j1 + 1)*2<<(i2 - i1-1)
            return 1
        else
            return 0
        end
    end
end

"""
    SharedSupport(n, L)

Collect all the indexes (single indexation) relative to the basis function sharing
sharing the same support of the function ϕ_n up to level L
"""
function SharedSupport(n,L)
    a = Int64[]
    for i in 1:(2<<(L)) - 1
        if ShareDomain(n,i) == 1 # && n != i
            push!(a, i)
        end
    end
    return a
end

"""
    SharedSupport(i,j, L)

Collect all the indexes (with single indexation) relative to the basis functions
sharing the same support of the function ϕ_{i,j} up to  level L
"""
function SharedSupport(i, j, L)
    n = Faber(i,j)
    a = SharedSupport(n, L)
    return a
end
"""
    NotSharedSupport(i,j, L)

Collect all the indexes (with single indexation) relative to the basis functions
which do the same support of the function ϕ_{i,j} up to  level L
"""
function NotSharedSupport(n, L)
    a = Int64[]
    for i in 1:(2<<(L)) - 1
        if ShareDomain(n,i) == 0
            push!(a, i)
        end
    end
    return a
end

"""
    FS

the strcut
```
struct Fs
    i::Int16    #idenx i (double indexation)
    j::Int16    #index j (double indexation)
    n::Int64    #idenx n (single indexation)
    lb::Float64 #pper bound domain
    ub::Float64 #lower Bound domain [not used?]
    range::Float64 #ub - lb
    nhb::Array{Int64} #neightbourhoods
    nhb::Array{Int64} #not neightbourhoods
    supr::Float64 #suprimum(ϕ)
    δ::Float64
end
```
"""
struct Fs
    i::Int16    #idenx i (double indexation)
    j::Int16    #index j (double indexation)
    n::Int64    #idenx n (single indexation)
    lb::Float64 #Upper bound domain
    ub::Float64 #Lower Bound domain [not used?]
    range::Float64
    nhb::Array{Int64} #neightbourhoods
    notnhb::Array{Int64}
    supr::Float64 #suprimum(ϕ)
    δ::Float64 # range*suprimum(ϕ)
    function Fs(i, j, L, T) #inizialization
        n = Faber(i, j)
        new(i, j, n, T*j/2^i, T*(j+1)/2^i, T/2^i, SharedSupport(n,L), NotSharedSupport(n,L), 2^(-i/2 - 1)*sqrt(T), T/2^i*2^(-i/2 -1)*sqrt(T))
    end
    function Fs(n, L, T) #inizialization
        i, j = Faber(n)
        Fs(i, j, L, T)
    end
end

"""
    generate(L, T)

generate all the set of basis used consistenly
"""
function generate(L, T)
    ϕ = Array{Fs}(undef, 2<<L - 1)
    for  n in 1:2<<L -1
        ϕ[n] = Fs(n, L, T)
    end
    ϕ
end

"""
    Λ(x, T)

function needed for construction of any Faber Schauder function i,j
see   Λ(x, l, k, T). This function should not be visible.
"""
function Λ(x, T)
    return max(0.0, sqrt(T)/2 - abs(x/sqrt(T) - sqrt(T)/2))
end

"""
    Λ(x, l, k, T)

Faber Schauder function `l`, `k` (double indexation) construct as in definition
rescaled by length `T` and evaluated at `x`
"""
function Λ(x, l, k, T)
    x = mod(x, T)   #ok
    2^(-l/2)*Λ((x)*(1<<l) - k*T, T)
end

"""
    IntQuad1(a, b, n1, n2)
integral of quadratic function1 with `a` quadratic term
`b` linear term and`n1` `n2` uppper and lower bound:

∫_n2^n1 a*t^2 + b*t dt
"""
function IntQuad1(a, b, n1, n2)
    a*(n1^3 - n2^3)/3 + b*(n1^2 - n2^2)/2
end

"""
    flip(i,j)

flip index j to the left to ease the integration:
int_0^1 psi_i (t) psi_j(t) dt

ex flip(1,1) = 0 flip(1,0) = 0 (do nothing), flip(3,5)= 2
"""
function flip(i,j)
    bound = ((1/2)*2^i - 1/2)
    if bound < j
        j = 2*bound - j
    end
    return Integer(j)
end
"""
    flip(i,j)

flip index j on the other side to ease the computation
of bound2(n,T)

ex. flip2(1,1) = 0, flip2(1,0) = 1, flip2(2,0) = 3, flip2(2,1) = 2.

"""
function flip2(i,j)
    bound = ((1/2)*2^i - 1/2)
    j = 2*bound - j
    return Integer(j)
end

"""
    int_prod_fs(i,j, T)

compute Φ_{ij} = ∫_0^T ϕ_i ϕ_j dt
ϕ_i being the ith (single indexation) Faber Schauder function
"""
#int_0^1 psi_i (t) psi_j(t) dt with i<j
function int_prod_fs(i,j, T)
    (i, j) = OrderBasis(i, j)
    i1, j1 = Faber(i) #ordering  i < j --> i1 < i2
    i2, j2 = Faber(j)
    i3 = i2 - i1
    j3 = j2 - 2^(i2-i1)*j1
    j3 = flip(i3,j3)
    A = j3/2^i3
    B = (j3 + 1/2)/2^i3
    C = (j3 + 1)/ 2^i3
    return T^2*2^((i2 - 5*i1)/2)*(IntQuad1(1,-j3/2^i3, B, A) - IntQuad1(1,-(j3+1)/2^i3, C, B))
end

"""
    intfs(i, T)

compute Φ_i = ∫_0^T ϕ_{ij} ds
"""
function intfs(n, T)
    i, j = Faber(n)
    return T^(1.5)*2^(-(3i + 4)/2 )
end

"""
    generate_matrix(L, T)
generate vector whose element is ∫_0^T ϕ(i1,j1)*ϕ(i2,j2) dt
for each i1,j1,i2,j2 until truncation level `L`
"""

function generate_matrix(L::Int64, T)
    dim = 2^(L+1) -1
    V = zeros(dim,dim)
    for i in 1:dim
        for j in 1:i
            if i == j
                i1, j1 = Faber(i)
                V[i,i] = T^2* 2.0^(-2*i1 - 2)/3
            else
                if ShareDomain(i,j)==1
                    res = int_prod_fs(i, j, T)
                    V[i,j], V[j,i] = (res, res)
                else
                    V[i,j], V[j,i] = (0 , 0)
                end
            end
        end
    end
    return V
end

"""
    generate_vector(L, T)
generate vector whose element is ∫_0^T ϕ(i,j) dt
for each i,j until truncation level `L`
"""

function generate_vector(L, T)
    [intfs(i, T) for i in 1:(2<<L -1)]
end

"""
    bound2(L, T)
computes ∫_0^T ϕ(n)(t/T) dt
N.B. `n` single inexation
"""
function bound1(n, T)
    i,j = Faber(n)
    A = j / (2^i)
    B = (0.5 + j)/ (2^i)
    C = (1.0 + j)/ (2^i)
    return T^(3/2)*2^(-i/2)*(IntQuad1(2^i,-j, B, A) + IntQuad1(-2^i, j + 1, C, B))
end

"""
    bound2(L, T)
computes ∫_0^T ϕ(n)(1 - t/T) dt
N.B. `n` single inexation
"""
function bound2(n, T)
    i,j = Faber(n)
    j = flip2(i,j)
    A = j / (2^i)
    B = (0.5 + j)/ (2^i)
    C = (1.0 + j)/ (2^i)
    return T^(3/2)*2^(-i/2)*(IntQuad1(2^i,-j, B, A) + IntQuad1(-2^i, j + 1, C, B))
end

"""
    generate_bound1(L, T)
generate vector whose element is ∫_0^T ϕ(i,j)(t/T) dt
for each i,j until truncation level `L`
"""
function generate_bound1(L, T)
    [bound1(n,T) for n in 1:1:2<<L-1]
end

"""
    generate_bound2(L, T)
generate vector whose element is ∫_0^T ϕ(i,j)(1 - t/T) dt
for each i,j until truncation level `L`
"""
function generate_bound2(L, T)
    [bound2(n,T) for n in 1:1:2<<L-1]
end

"""
    MCintegration(k::Fs)

take a random point inside the support of the basis function k
"""
function MCintegration(k::Fs)
    rand()*k.range + k.lb
end

"""
    Sol2E(a::Float64, b::Float64, c::Float64)

Solve a quadratic equation.
TODO add errors when the solution is not real.
solve quadratic equation ax^2 + bx + c = 0
"""
function Sol2E(a::Float64, b::Float64, c::Float64)
    return (-b + sqrt(b^2 - 4a*c))/2a
end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
using Test
using Discretizers

function array_matches(arr::AbstractVector{S}, arr2::AbstractVector{T}) where {S<:Integer, T<:Integer}
    n = length(arr)
    if length(arr2) != n
        return false
    end
    for (a,b) in zip(arr, arr2)
        if a != b
            return false
        end
    end
    return true
end
function array_matches(arr::AbstractVector{S}, arr2::AbstractVector{T}, abs_tolerance::AbstractFloat=eps(Float64)) where {S<:Real, T<:Real}
    n = length(arr)
    if length(arr2) != n
        return false
    end
    for i in 1 : n
        if abs(arr[i] - arr2[i]) > abs_tolerance
            return false
        end
    end
    return true
end

include("test_categorical_discretizer.jl")
include("test_linear_discretizer.jl")
include("test_hybrid_discretizer.jl")

include("test_common.jl")

include("test_disc_uniformwidth.jl")
include("test_disc_uniformcount.jl")
include("test_disc_bayesianblocks.jl")
include("test_disc_MODL.jl")


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

@recipe function f(domain::AbstractDomain, data::AbstractVector) where {N,T}
  @series begin
    PointSet(coordinates(domain)), data
  end
end

@recipe function f(domain::AbstractDomain)
  @series begin
    PointSet(coordinates(domain))
  end
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file does not have any dependencies upon other scripts

"""
    puzzle_n(rank)

Given 2 <= rank <= 15, return a Sudoku integer array sized rank*rank.
"""
function puzzle_n(rank::Integer)
    if (rank < 2) || (rank > 15)
        throw(DomainError("Invalid puzzle rank outside [2,15]"))
    end
    rank_squared = rank*rank
    result = Array{Integer}(undef,(rank_squared,rank_squared))
    for col = 1:rank_squared
        col_shift = Integer(floor((col-1)/rank))
        for row = 1:rank_squared
            result[row,col] = ( (rank*(col-1) + (row-1)) + col_shift) % rank_squared + 1
        end
    end
    return result
end

"""
A maybe-unknown Sudoku value with possible values tracked.

If value > 0, the entry is known and 1<= value <= rank-squared.
If value == 0, the the entry is unknown, and possible values indices into a bit vector for which true is assigned.
"""
struct PuzzleEntry
    possibilities::UInt16 # Bitmask index = 1 iff index is allowed

    """
    Define entry with a value and a BitVector of possibilities.
    """
    function PuzzleEntry(possibilities::BitVector)
        rank = 0
        rank_squared = length(possibilities)
        if rank_squared == 4
            rank = 2
        elseif rank_squared == 9
            rank = 3
        elseif rank_squared == 16
            rank = 4
        else
            throw(DomainError("Possibilities vector length invalid"))
        end
        new(UInt16(possibilities.chunks[1])) # Use local new method
    end
end

"""
    make_entry(rank_squared,value)

Given a squared-rank > 3 and < 17, create a new PuzzleEntry.
A value of 0 is completely unknown.
A value in 1:rank_squared is completely known.
"""
function make_entry(rank_squared::Integer, value::Integer)
    if !(rank_squared in [4,9,16])
        throw(DomainError("Cannot instantiate an entry for rank < 2 or > 4"))
    end
    if value < 0 || value > rank_squared
        throw(DomainError("Cannot instantiate an entry from an value < 0 or > rank*rank"))
    end
    possibilities = BitVector(undef,rank_squared)
    if value == 0
        possibilities .= true
    else
        possibilities .= false
        possibilities[value] = true
    end
    return PuzzleEntry(possibilities)
end

function get_value(entry::PuzzleEntry)
    temp = BitVector(undef,16) # Driven by storage as Int16
    temp.chunks[1] = entry.possibilities # Convert int to BitVector
    possible_count = 0 # Number of possibilities found
    result = 0 # Default to unknown
    for i in 1:length(temp)
        if temp[i]
            possible_count += 1
            result = i
        end
    end
    if possible_count != 1
        result = 0
    end
    return result
end

"""
A square grid of PuzzleEntry instances.

See also [`random_puzzle`](@ref), [`solve_puzzle!`](@ref).
"""
struct SolvablePuzzle
    grid::Array{PuzzleEntry}
    """

    """
    function SolvablePuzzle(rank::Integer)
        solved_puzzle = puzzle_n(rank)
        puzzle = Array{PuzzleEntry}(undef,size(solved_puzzle))
        rank_squared = rank*rank
        for i = 1:length(solved_puzzle)
            v = Integer(solved_puzzle[i])
            puzzle[i] = make_entry(rank_squared,v)
        end
        return new(puzzle)
    end
end

"""
    assign_values!(puzzle,new_values)

Assign a grid of integers into a SolvablePuzzle.
Zero values are set as unknowns.

# Arguments
- puzzle::SolvablePuzzle: grid of Sudoku values to set.
- new_values::Array: input integers of same shape as SolvablePuzzle.grid
"""
function assign_values!(puzzle::SolvablePuzzle,new_values::Array)
    if size(puzzle.grid) != size(new_values)
        throw(DomainError("Mismatched shape of new values for puzzle"))
    end
    rank_squared = size(new_values)[1]
    for row = 1:rank_squared
        for col = 1:rank_squared
            value = new_values[row,col]
            if value < 0 || value > rank_squared
                throw(DomainError("Invalid new value; must be in [0 (unknown), rank-squared]"))
            elseif new_values[row,col] == 0
                # Unknown value
                set_unknown(puzzle,row,col)
            else
                v = Integer(value)
                puzzle.grid[row,col] = make_entry(rank_squared,v)
            end
        end
    end
end

"""
    get_rank(puzzle)

Return a Sudoku rank from the shape of a 2D aray.
"""
function get_rank(puzzle::Array)
    grid_shape = size(puzzle)
    if grid_shape == (4,4) || grid_shape == (16,)
        return 2
    elseif grid_shape == (9,9) || grid_shape == (81,)
        return 3
    elseif grid_shape == (16,16) || grid_shape == (256,)
        return 4
    else
        throw(DimensionMismatch("Puzzles must of shape (rank*rank,rank*rank), 1 < rank < 5"))
    end
end

"""
    get_rank(puzzle)

Return a Sudoku rank from a SolvablePuzzle.
"""
function get_rank(puzzle::SolvablePuzzle)
    return get_rank(puzzle.grid)
end

"""
    set_unknown(puzzle,row,col)

Assign a position in the grid of a SolvablePuzzle to known.
"""
function set_unknown(puzzle::SolvablePuzzle,row::Integer,col::Integer)
    p = BitVector(undef,size(puzzle.grid)[1])
    p .= true
    puzzle.grid[row,col] = PuzzleEntry(p)
    return nothing
end

"""
    as_text(puzzle)

Return the contents of the puzzle as a grid of strings.

Known values are integers greater than 0.
Unknown values are 0.
"""
function as_text(puzzle::SolvablePuzzle)
    result = fill(" ",size(puzzle.grid)) # unknowns!
    for i = 1:length(puzzle.grid)
        value = get_value(puzzle.grid[i])
        if value != 0 # fill in the knowns
            result[i] = string(value)
        end
    end
    return result
end

"""
    as_values(puzzle)

Return the contents of the puzzle as a grid of integers.

Known values are integers greater than 0.
Unknown values are 0.
"""
function as_values(puzzle::SolvablePuzzle)
    result = Array{Integer}(undef, size(puzzle.grid)) # invalid data
    for i = 1:length(puzzle.grid)
        result[i] = get_value(puzzle.grid[i])
    end
    return result
end

"""
    as_possibilities(puzzle)

Return the contents of the puzzle as a grid of integers representing possible values.

Known values are have one value set to true at that values index.
Unknown values are have more than one bit set to true.
"""
function as_possibilities(puzzle::SolvablePuzzle)
    result = Array{UInt16}(undef, size(puzzle.grid)) # invalid data
    for i = 1:length(puzzle.grid)
        result[i] = puzzle.grid[i].possibilities
    end
    return result
end

"""
    uncertainty(puzzle)

Assess the degree of uncertainty in a SolvablePuzzle.

Zero uncertainty is associated with solved puzzles.
Unsolved puzzles have uncertainty greater than zero.
"""
function uncertainty(puzzle::SolvablePuzzle)
    result::Integer = 0
    rank = get_rank(puzzle.grid)
    rank_squared = rank*rank
    error_detected = false
    for i = 1:length(puzzle.grid)
        temp = BitVector(undef,rank_squared)
        temp.chunks[1] = puzzle.grid[i].possibilities # Convert int to BitVector
        possibility_count = sum(temp)
        error_detected = error_detected || (possibility_count == 0)
        if possibility_count > 1
            result += sum(temp)
        end
    end
    if error_detected
        throw(DomainError("Puzzle entry with zero valid possibilities detected."))
    end
    return result
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia