text:
## W O R K B O O K
# Empty workbook
Workbook() = Workbook(
    Dict{String,Sheet}(), # sheets
    OrderedDict{GlobalRef, Cell}(), # map of positions to cells
    SimpleDiGraph(), # mapping dependencies
    Int[] # topological ordering
)

# getters
# get a cell
get(wb::Workbook, p::GlobalPosition, default::Any = EmptyCell(p)) = Base.get(wb.map, p, default)
get(wb::Workbook, R::GlobalRange) = get.([wb], collect(R))
get!(wb::Workbook, p::GlobalPosition, default::Cell = EmptyCell(p)) = Base.get!(wb.map, p, default)
get!(wb::Workbook, R::GlobalRange) = get!.([wb], collect(R))

# get a sheet
Base.get(wb::Workbook, s::AbstractString) = getindex(wb, s)
function Base.getindex(wb::Workbook, s::AbstractString, default = missing)
    # synonym of wb[s], e.g. wb["Sheet1"] or wb["Sheet1!A1"]
    m = match(sheetregex,s) # when s = e.g. "Sheet1"
    if !isnothing(m)
        if !(s in keys(wb.sheets))
            if ismissing(default)
                wb.sheets[s] = Sheet(s, wb)
            else
                @assert (default.wb === wb) "Sheet must belong to this workbook"
                wb.sheets[s] = default
            end
        end
        return wb.sheets[s]
    else
        return get(wb, GlobalPosition(s))
    end
end

## calculation agents
# get values
function value(wb::Workbook, p::GlobalPosition)
    c = wb.map[p]
    c.value
end
function value(wb::Workbook, r::GlobalRange)
    map(p -> value(wb,p), collect(r))
end

# calculate the value of the cell based on current state of sheet
refresh!(p::AbstractString, wb::Workbook) = refresh!(GlobalPosition(p), wb)
refresh!(p::GlobalPosition, wb::Workbook) = refresh!(get(wb,p), wb)
function refresh!(_::ConstCell,_::Workbook) end # nothing to do
function refresh!(c::DynamicCell,wb::Workbook)
    c.value = c.func(value.([wb],c.ancestors)...)
end
function refresh_latest!(_::ConstCell,_::Workbook) end
function refresh_latest!(c::DynamicCell,wb::Workbook)
    c.value = Base.invokelatest(c.func, value.([wb],c.ancestors)...)
end

function orderedFromSuperset(xx,xxx)
    xx[sortperm(findfirst.(isequal.(xx),[xxx]))]
end
linear_location(c::Cell, wb::Workbook) = linear_location(c.position, wb)
linear_location(p::GlobalPosition, wb::Workbook) = findfirst(isequal(p),wb.map.keys)
function cascade(p::GlobalPosition,wb::Workbook,latest::Bool=false)
    loc = linear_location(p,wb)
    refreshNeeded = orderedFromSuperset(neighbors(wb.graph, loc),wb.order)
    for key in refreshNeeded
        latest ? refresh_latest!(wb.map[wb.map.keys[key]],wb) : refresh!(wb.map[wb.map.keys[key]],wb)
        push!(refreshNeeded, orderedFromSuperset(neighbors(wb.graph, key),wb.order)...) # add cells that are further downstream to the list
    end
end
# get the cell into the overall weave of the thing
function extendgraph(_::ConstCell, _::Workbook) true end # nothing to do
function extendgraph(c::DynamicCell, wb::Workbook)
    # basically: make sure cell's ancestors exist & are recorded

    # create nonexistent ancestors
    for anc in c.ancestors # iterating over GlobalPositions
        get!(wb, anc) # creates empty cells where necessary
    end
    # extend graph
    incrementalNodes = length(wb.map) - nv(wb.graph)
    add_vertices!(wb.graph,incrementalNodes)
    # add edges
    cyclic = false # flag for our loop
    added_ancestors = GlobalPosition[]
    for ancestor in c.ancestors # iterating over GlobalRefs to add edges to graph
        cyclic |= any(add_ancestor(wb,ancestor,c.position))
         # record the edges we have been adding
        if isa(ancestor,GlobalPosition)
            push!(added_ancestors,ancestor)
        elseif isa(ancestor,GlobalRange)
            push!(added_ancestors, collect(ancestor)...)
        else
            error()
        end
        if cyclic # creates cyclic dependencies
            # unset all of the previously added edges
            for ancestor in added_ancestors
                remove_ancestor(wb, ancestor, c.position)
            end
            return false # escalate this problem!
        end
    end

    # update total ordering
    wb.order = topological_sort_by_dfs(wb.graph)
    true # extended w/o problems
end
function add_ancestor(wb::Workbook, ancestor::GlobalPosition, descendant::GlobalPosition)
    loc = linear_location(descendant,wb)
    ancloc = findfirst(isequal(ancestor),wb.map.keys)
    add_edge!(wb.graph, ancloc, loc) # ancestor -> descendant
    if is_cyclic(wb.graph) || has_self_loops(wb.graph)
        return true
    end
    false
end
function remove_ancestor(wb::Workbook, ancestor::GlobalPosition, descendant::GlobalPosition)
    loc = linear_location(descendant,wb)
    ancloc = findfirst(isequal(ancestor),wb.map.keys)
    rem_edge!(wb.graph, ancloc, loc) # ancestor -> descendant
end
add_ancestor(wb::Workbook, ancestor::GlobalRange, descendant::GlobalPosition) = add_ancestor.([wb], collect(ancestor), [descendant])

# write a cell
Base.setindex!(wb::Workbook, f::AbstractString, p) = set!(wb, Cell(p,f))
Base.setindex!(wb::Workbook, f::Number, p) = set!(wb, ConstCell(p,f))
Base.setindex!(wb::Workbook, f, x) = setindex!(wb,string(f),x)

function setcellbasics(wb::Workbook, c::Cell)
    p = c.position
    # add the cell to the wb
    push!(wb.map, p => c)
    # add the cell to the graph
    add_vertices!(wb.graph,length(wb.map) - nv(wb.graph))
    # add the cell to the sheet
    push!(wb[p.sheet].map, p.lref => c)
end
# remove all ancestors
function cleanup(wb::Workbook, p::GlobalPosition)
    loc = linear_location(p, wb)
    for anc in vertices(wb.graph)
        rem_edge!(wb.graph, anc, loc) # many false positives
    end
end
function set!(wb::Workbook, c::Cell)
    p = c.position
    is_new_node = !haskey(wb.map,p)
    if !is_new_node
        oldCell = wb.map[p]
        cleanup(wb,p)
    end
    setcellbasics(wb,c)
    graph_extended = extendgraph(c,wb)
    if graph_extended
        refresh_latest!(c,wb)
        # update downstream
        is_new_node || cascade(c.position,wb,true) # only if an existing cell has been replaced
    else
        if !is_new_node
            set!(wb, oldCell)
        else
            set!(wb, EmptyCell(p))
        end
        error("The new cell was not set because it creates a cyclic reference")
    end
end

## O U T P U T
const minPos = LocalPosition(1,1)
const maxPos = LocalPosition(typemax(UInt),typemax(UInt))

function Base.string(wb::Workbook)
    if isempty(wb.sheets)
        s = "Empty workbook"
    else
        s = "Workbook with "
        s *= string(length(wb.sheets)) * " sheet" * ((length(wb.sheets) > 1) ? "s" : "")
        s *= ": " * join(string.(keys(wb.sheets)),", "," and ")
    end
    s
end
Base.show(io::IO, wb::Workbook) = Base.print(io, string(wb))

hassheet(wb::Workbook, sname::String) = haskey(wb.sheets,sname)


text:
function add(obj::AbstractRealMatrix, arg0::RealMatrix)
    return jcall(obj, "add", RealMatrix, (RealMatrix,), arg0)
end

function add_to_entry(obj::AbstractRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "addToEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function copy(obj::AbstractRealMatrix)
    return jcall(obj, "copy", RealMatrix, ())
end

function copy_sub_matrix(obj::AbstractRealMatrix, arg0::Vector{jint}, arg1::Vector{jint}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "copySubMatrix", void, (Vector{jint}, Vector{jint}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function copy_sub_matrix(obj::AbstractRealMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::Vector{Vector{jdouble}})
    return jcall(obj, "copySubMatrix", void, (jint, jint, jint, jint, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3, arg4)
end

function create_matrix(obj::AbstractRealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", RealMatrix, (jint, jint), arg0, arg1)
end

function equals(obj::AbstractRealMatrix, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_column(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "getColumn", Vector{jdouble}, (jint,), arg0)
end

function get_column_dimension(obj::AbstractRealMatrix)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_column_matrix(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "getColumnMatrix", RealMatrix, (jint,), arg0)
end

function get_column_vector(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "getColumnVector", RealVector, (jint,), arg0)
end

function get_data(obj::AbstractRealMatrix)
    return jcall(obj, "getData", Vector{Vector{jdouble}}, ())
end

function get_entry(obj::AbstractRealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", jdouble, (jint, jint), arg0, arg1)
end

function get_frobenius_norm(obj::AbstractRealMatrix)
    return jcall(obj, "getFrobeniusNorm", jdouble, ())
end

function get_norm1(obj::RealMatrix)
    return jcall(obj, "getNorm1", jdouble, ())
end

function get_norm_infty(obj::RealMatrix)
    return jcall(obj, "getNormInfty", jdouble, ())
end

function get_row(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "getRow", Vector{jdouble}, (jint,), arg0)
end

function get_row_dimension(obj::AbstractRealMatrix)
    return jcall(obj, "getRowDimension", jint, ())
end

function get_row_matrix(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "getRowMatrix", RealMatrix, (jint,), arg0)
end

function get_row_vector(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "getRowVector", RealVector, (jint,), arg0)
end

function get_sub_matrix(obj::AbstractRealMatrix, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(obj, "getSubMatrix", RealMatrix, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function get_sub_matrix(obj::AbstractRealMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint)
    return jcall(obj, "getSubMatrix", RealMatrix, (jint, jint, jint, jint), arg0, arg1, arg2, arg3)
end

function get_trace(obj::AbstractRealMatrix)
    return jcall(obj, "getTrace", jdouble, ())
end

function hash_code(obj::AbstractRealMatrix)
    return jcall(obj, "hashCode", jint, ())
end

function is_square(obj::AbstractRealMatrix)
    return jcall(obj, "isSquare", jboolean, ())
end

function is_transposable(obj::RealLinearOperator)
    return jcall(obj, "isTransposable", jboolean, ())
end

function map(obj::RealMatrix, arg0::UnivariateFunction)
    return jcall(obj, "map", RealMatrix, (UnivariateFunction,), arg0)
end

function map_to_self(obj::RealMatrix, arg0::UnivariateFunction)
    return jcall(obj, "mapToSelf", RealMatrix, (UnivariateFunction,), arg0)
end

function multiply(obj::AbstractRealMatrix, arg0::RealMatrix)
    return jcall(obj, "multiply", RealMatrix, (RealMatrix,), arg0)
end

function multiply_entry(obj::AbstractRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "multiplyEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function multiply_transposed(obj::RealMatrix, arg0::RealMatrix)
    return jcall(obj, "multiplyTransposed", RealMatrix, (RealMatrix,), arg0)
end

function operate(obj::AbstractRealMatrix, arg0::RealVector)
    return jcall(obj, "operate", RealVector, (RealVector,), arg0)
end

function operate(obj::AbstractRealMatrix, arg0::Vector{jdouble})
    return jcall(obj, "operate", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function operate_transpose(obj::RealLinearOperator, arg0::RealVector)
    return jcall(obj, "operateTranspose", RealVector, (RealVector,), arg0)
end

function power(obj::AbstractRealMatrix, arg0::jint)
    return jcall(obj, "power", RealMatrix, (jint,), arg0)
end

function pre_multiply(obj::AbstractRealMatrix, arg0::RealMatrix)
    return jcall(obj, "preMultiply", RealMatrix, (RealMatrix,), arg0)
end

function pre_multiply(obj::AbstractRealMatrix, arg0::RealVector)
    return jcall(obj, "preMultiply", RealVector, (RealVector,), arg0)
end

function pre_multiply(obj::AbstractRealMatrix, arg0::Vector{jdouble})
    return jcall(obj, "preMultiply", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function scalar_add(obj::AbstractRealMatrix, arg0::jdouble)
    return jcall(obj, "scalarAdd", RealMatrix, (jdouble,), arg0)
end

function scalar_multiply(obj::AbstractRealMatrix, arg0::jdouble)
    return jcall(obj, "scalarMultiply", RealMatrix, (jdouble,), arg0)
end

function set_column(obj::AbstractRealMatrix, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "setColumn", void, (jint, Vector{jdouble}), arg0, arg1)
end

function set_column_matrix(obj::AbstractRealMatrix, arg0::jint, arg1::RealMatrix)
    return jcall(obj, "setColumnMatrix", void, (jint, RealMatrix), arg0, arg1)
end

function set_column_vector(obj::AbstractRealMatrix, arg0::jint, arg1::RealVector)
    return jcall(obj, "setColumnVector", void, (jint, RealVector), arg0, arg1)
end

function set_entry(obj::AbstractRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "setEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function set_row(obj::AbstractRealMatrix, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "setRow", void, (jint, Vector{jdouble}), arg0, arg1)
end

function set_row_matrix(obj::AbstractRealMatrix, arg0::jint, arg1::RealMatrix)
    return jcall(obj, "setRowMatrix", void, (jint, RealMatrix), arg0, arg1)
end

function set_row_vector(obj::AbstractRealMatrix, arg0::jint, arg1::RealVector)
    return jcall(obj, "setRowVector", void, (jint, RealVector), arg0, arg1)
end

function set_sub_matrix(obj::AbstractRealMatrix, arg0::Vector{Vector{jdouble}}, arg1::jint, arg2::jint)
    return jcall(obj, "setSubMatrix", void, (Vector{Vector{jdouble}}, jint, jint), arg0, arg1, arg2)
end

function subtract(obj::AbstractRealMatrix, arg0::RealMatrix)
    return jcall(obj, "subtract", RealMatrix, (RealMatrix,), arg0)
end

function to_string(obj::AbstractRealMatrix)
    return jcall(obj, "toString", JString, ())
end

function transpose(obj::AbstractRealMatrix)
    return jcall(obj, "transpose", RealMatrix, ())
end

function transpose_multiply(obj::RealMatrix, arg0::RealMatrix)
    return jcall(obj, "transposeMultiply", RealMatrix, (RealMatrix,), arg0)
end

function walk_in_column_order(obj::AbstractRealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_column_order(obj::AbstractRealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_column_order(obj::AbstractRealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_column_order(obj::AbstractRealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::AbstractRealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::AbstractRealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::AbstractRealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::AbstractRealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::AbstractRealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_row_order(obj::AbstractRealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::AbstractRealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_row_order(obj::AbstractRealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

ext: jl


text:
using AutoGrid

func = function (w)
    s = im*w
    tf1 = 0.1*exp(-s)/(s^2 + 1e-4*s +(0.1)^2)
    tf2 = (10^2/(s^2 + s +(10)^2))*(10^2/(s^2 + s +(100)^2))
    return (abs.(tf1), atan.(imag.(tf1),real.(tf1)), abs.(tf2), atan.(imag.(tf2),real.(tf2)))
end
lims = (1e-2, 1e3)
xscale = log10
xscale_inv = exp10
funcConfigs = (FuncConfig(log10, find_peaks=true), WrappingFuncConfig((0.0, 2pi), find_peaks=false), FuncConfig(log10, find_peaks=true), WrappingFuncConfig((0.0, 2pi), find_peaks=false))
agc = AutoGridConfig(func, lims, xscale, xscale_inv, funcConfigs, abstol_dd=1e-2, reltol_dd=1.0, min_eps=1e-4)

values, grid = auto_grid(agc);
@show length(grid)
using Plots
plot(grid, values[1], m=:o, layout=(2,1), subplot=1, scale=:log10, yscale=:log10)
plot!(grid, values[3], m=:o, layout=(2,1), subplot=1, scale=:log10, yscale=:log10)
plot!(grid, values[2]*180/pi, m=:o, xscale=:log10, subplot=2)
plot!(grid, values[4]*180/pi, m=:o, xscale=:log10, subplot=2)

grid2 = 10 .^LinRange(-2,3,length(grid))
plot!(grid2, 0.01.*getindex.(func.(grid2), 1), l=(:green, :dash), subplot=1, m=:o)

values, grid = auto_grid(agc)
using BenchmarkTools
@btime values, grid = auto_grid($agc);

function bench1(grid)
    func.(grid)
end
@btime bench1(grid);

function bench2(grid)
    vec = [func(grid[1]),]
    for i in 2:length(grid)
        push!(vec, func(grid[i]))
    end
    return vec
end
@btime bench2(grid)
# No peaks
# 180.647 μs (2412 allocations: 91.73 KiB) Length: 85
# Both peaks
# 35.296 μs (61 allocations: 20.11 KiB) Length: 196
# 68.038 μs (59 allocations: 29.14 KiB) Length: 376

# Bench 1:
# 1.836 μs (4 allocations: 2.88 KiB)
# Bench 2:
# 8.815 μs (350 allocations: 16.31 KiB)


text:
# using PathDistribution, Gadfly
#
# using LsqFit, Gadfly
# include("../src/misc.jl")
# include("../src/monte_carlo.jl")
# include("../src/fitting.jl")

using PathDistribution

function do_exp(filename, origin, destination, start_node, end_node, link_length, beta_fit)

    N1 = 5000
    N2 = 20000

    samples = monte_carlo_path_sampling(origin, destination, start_node, end_node, link_length, N1, N2)

    option = :unique
    x_data, y_data = estimate_cumulative_count(samples, option)

    writecsv("$(filename)-$origin-$destination.csv", [x_data y_data])

    println("$filename, origin=$origin, destination=$destination")
    println("number of paths estimate: $(y_data[end])")

end

filename = "Buffalo-Data"
data = readcsv("$(filename).csv", header=true) # 405,094
origin = 1
destination = 84
beta_fit = [401103.602284; 61.890097; 9.811411]
start_node = round(Int64, data[1][:,1])
end_node = round(Int64, data[1][:,2])
link_length = data[1][:,3]
@time do_exp(filename, origin, destination, start_node, end_node, link_length, beta_fit)

#
# filename = "Albany-Data"
# data = readcsv("$(filename).csv", header=true) #3,551
# origin = 1
# destination = 12
# beta_fit = [3504.287463; 71.552699; 6.485460]
# start_node = round(Int64, data[1][:,1])
# end_node = round(Int64, data[1][:,2])
# link_length = data[1][:,3]
# @time do_exp(filename, origin, destination, start_node, end_node, link_length, beta_fit)
#
#
#
# filename = "Sioux-Falls"
# data = readcsv("$(filename).csv", header=true) #3,705
# start_node = round(Int64, data[1][:,1])
# end_node = round(Int64, data[1][:,2])
# link_length = data[1][:,3]
#
# origin = 2
# destination = 13
# beta_fit = [4559.889822; 235.971248; 6.032994]
# @time do_exp(filename, origin, destination, start_node, end_node, link_length, beta_fit)
#
# origin = 3
# destination = 19
# beta_fit = [3627.177263; 220.070292; 5.327236]
# @time do_exp(filename, origin, destination, start_node, end_node, link_length, beta_fit)
#
# origin = 10
# destination = 21
# beta_fit = [2777.234877; 239.599660; 4.416381]
# @time do_exp(filename, origin, destination, start_node, end_node, link_length, beta_fit)

#
# start_node_un = round(Int64, data[1][:,1])
# end_node_un = round(Int64, data[1][:,2])
# link_length_un = data[1][:,3]
# start_node = [start_node_un; end_node_un]
# end_node = [end_node_un; start_node_un]
# link_length = [link_length_un; link_length_un]


text:
using RBNF
using MLStyle

# rmlines = @λ begin
#            e :: Expr           -> Expr(e.head, filter(x -> x !== nothing, map(rmlines, e.args))...)
#              :: LineNumberNode -> nothing
#            a                   -> a
# end

# node = quote
# RBNF.@parser TestLang begin
#     @token
#     a := r"\d"
#     @grammar
#     b := [(k = k.?), (seqa :: Vector = Many(a => parse_int))]
# end
# end

# @info :show rmlines(macroexpand(Main, node))

# struct TestLang end
# parse_int(x :: RBNF.Token) = parse(Int, x.str)

# RBNF.@parser TestLang begin
#     @token
#     a := r"\G\d"
#     k := "kkk"
#     @grammar
#     b := [(k = k.?), (seqa :: Vector = Many(a => parse_int))]
# end

# res1, _ = RBNF.runparser(b, RBNF.runlexer(TestLang, "123"))
# println(res1)

# res2, _ = RBNF.runparser(b, RBNF.runlexer(TestLang, "kkk123"))
# println(res2)

# struct TestLang2 end

# RBNF.@parser TestLang2 begin
#     @token
#     num := r"\G\d+"
#     space := r"\G\s+"
#     @grammar
#     manynumbers := @direct_recur begin
#         init = num
#         prefix = [recur, space, num]
#     end
# end

# res3, _ = RBNF.runparser(manynumbers, RBNF.runlexer(TestLang2, "52 123 123 14312 213"))
# println(res3)
include("qasm.jl")
include("reml.jl")


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

using .LRUExample

TestLRU = LRUExample.UnboundedLRU{String, String}()
TestBLRU = LRUExample.BoundedLRU{String, String}(1000)

get_str(i) = String(vcat(map(x->[x>>4; x&0x0F], reinterpret(UInt8, [Int32(i)]))...))

isbounded{L<:LRUExample.LRU}(::Type{L}) = any(map(n->n==:maxsize, fieldnames(L)))
isbounded{L<:LRUExample.LRU}(l::L) = isbounded(L)

nmax = round.(Int, logspace(2, 5, 4))

function lrutest()
    #println("LRU consistency tests")
    for lru in (TestLRU,TestBLRU)
        for n in nmax
            empty!(lru)
            #@printf("  %s, %d items\n", lru, n)
            #print("    Simple eviction: ")
            for i in 1:n
                str = get_str(i)
                lru[str] = str
                @assert lru.q[1].v == str
                if isbounded(lru) && length(lru) >= lru.maxsize
                    tailstr = get_str(i-lru.maxsize+1)
                    @assert lru.q[end].v == tailstr
                end
            end
            #println("pass")

            #print("    Lookup, random access: ")
            for i in 1:n
                str = get_str(rand(1:n))
                if haskey(lru, str) # the bounded LRUs can have cache misses
                    blah = lru[str]
                    @assert lru.q[1].v == blah
                end
            end
            #println("pass")
        end
        empty!(lru)
    end
end

lrutest()

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
#=
Compute 1D NFFT on each column/row of a matrix.
# TODO: Not very elegant, try to mimick the mapreduce approach of e.g. sum
=#

function NFFT.nfft!{T}(p::NFFTPlan{1}, f::AbstractMatrix{T}, fHat::StridedMatrix{T}, ::Type{Val{1}})
	( Nrows = size(f,2) ) == size(fHat,2) || throw(DimensionMismatch())

	for rowidx in 1:Nrows
		frow = slice(f, :, rowidx)
		fhatrow = slice(fHat, :, rowidx)

		NFFT.nfft!(p, frow, fhatrow)
	end

	return fHat
end

function NFFT.nfft!{T}(p::NFFTPlan{1}, f::AbstractMatrix{T}, fHat::StridedMatrix{T}, ::Type{Val{2}})
	( Nrows = size(f,1) ) == size(fHat,2) || throw(DimensionMismatch())

	for rowidx in 1:Nrows
		frow = slice(f, rowidx, :)
		fhatrow = slice(fHat, :, rowidx)

		NFFT.nfft!(p, frow, fhatrow)
	end

	return fHat
end

function NFFT.nfft_adjoint!{T}(p::NFFTPlan{1}, fHat::AbstractMatrix{T}, f::StridedMatrix{T}, ::Type{Val{1}})
	( Nrows = size(f,2) ) == size(fHat,2) || throw(DimensionMismatch())

	for rowidx in 1:Nrows
		frow = slice(f, :, rowidx)
		fhatrow = slice(fHat, :, rowidx)

		NFFT.nfft_adjoint!(p, fhatrow, frow)
	end

	return f
end

function NFFT.nfft_adjoint!{T}(p::NFFTPlan{1}, fHat::AbstractMatrix{T}, f::StridedMatrix{T}, ::Type{Val{2}})
	( Nrows = size(f,1) ) == size(fHat,2) || throw(DimensionMismatch())

	for rowidx in 1:Nrows
		frow = slice(f, rowidx, :)
		fhatrow = slice(fHat, :, rowidx)

		NFFT.nfft_adjoint!(p, fhatrow, frow)
	end

	return f
end


text:
struct Framed{T}
    frame::CartesianFrame3D
    val::T
end

function Base.show(io::IO, x::Framed)
    print(io, x.val)
    print(io, " (in frame: $(x.frame))")
end

in_frame(frame::CartesianFrame3D, val) = Framed(frame, val)
in_frame(::Nothing, val) = val

unwrap(x::Framed) = x.val
unwrap(x) = x

@inline function same_frame(args::Tuple, frame::Union{Missing, CartesianFrame3D}=missing)
    isempty(args) && return frame
    head = args[1]
    tail = Base.tail(args)
    if head isa Framed
        if frame !== missing && head.frame !== frame
            throw(ArgumentError("Frame mismatch.")) # TODO: nicer message
        end
        return same_frame(tail, head.frame)
    else
        return same_frame(tail, frame)
    end
end

function framechecked(f, args...)
    frame = same_frame(args)
    in_frame(frame, f(map(unwrap, args)...))
end

macro framechecked(expr)
    expr.head == :call || return :(error("Expected a function call"))
    f = esc(expr.args[1])
    args = map(esc, expr.args[2 : end])
    return quote
        framechecked($f, $(args...))
    end
end


text:
println("-------------------")
println(" Testing JuLIP.ASE")
println("-------------------")

using JuLIP

# ======================================================================

println("Check that the cubic unit cell for Al is reproduced correctly")
a0 = 2.025
p = [ [0.0;0.0;0.0] [0.0;a0;a0] [a0;0.0;a0] [a0;a0;0.0] ]
at = Atoms("Al", cubic=true)
@test (positions(at) |> mat) == p

# ======================================================================

println("Check neighbourlist without periodicity")
# TODO: implement a test with periodicity!!!
println("   TODO: implement test with periodicity as well?")
println("   ... assemble neighbour list ...")
at = Atoms("Al", cubic=true, pbc = (false,false,false), repeatcell=(3,3,3))
cutoff = 1.7 * a0
nlist = neighbourlist(at, cutoff)
# create a neighbourlist via a naive double-loop
simple = zeros(length(at), length(at))
X = positions(at)
for n = 2:length(at), m = 1:n-1
   if norm(X[m]-X[n]) <= cutoff
      simple[n,m] = simple[m,n] = 1
   end
end
println("   ... check the bond-iterator ... ")
for (i,j,r,R,S) in bonds(nlist)
   @test simple[i,j] == 1
   @test_approx_eq_eps(norm(X[i]-X[j]), r, 1e-12)
   @test_approx_eq_eps(X[j]-X[i], R, 1e-12)
   @test norm(S) == 0
   # switch the flag
   simple[i,j] = -1
end
# check that all pairs have been found
@test maximum(simple) == 0
# revert to original
simple *= -1
println("   ... check the site iterator ... ")
for (i,j,r,R,S) in sites(nlist)
   for n = 1:length(j)
      @test simple[i,j[n]] == 1
      simple[i,j[n]] = -1
   end
   @test maximum(simple[i,:]) == 0
end
@test maximum(simple) == 0

# ======================================================================


text:
abstract type LTISystem <: AbstractSystem end
+(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...)

# Fallback number
+(sys1::LTISystem, sys2::Number) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::Number) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::Number) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::Number) = /(promote(sys1, sys2)...)

+(sys1::Number, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::Number, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::Number, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::Number, sys2::LTISystem) = /(promote(sys1, sys2)...)

# Fallback Matrix
+(sys1::LTISystem, sys2::AbstractMatrix) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::AbstractMatrix) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::AbstractMatrix) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::AbstractMatrix) = /(promote(sys1, sys2)...)

+(sys1::AbstractMatrix, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::AbstractMatrix, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::AbstractMatrix, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::AbstractMatrix, sys2::LTISystem) = /(promote(sys1, sys2)...)

# TODO We should have proper fallbacks for matrices too

feedback(sys1::Union{LTISystem,Number,AbstractMatrix{<:Number}},
         sys2::Union{LTISystem,Number,AbstractMatrix{<:Number}}; kwargs...) = feedback(promote(sys1, sys2)...; kwargs...)

Base.inv(G::LTISystem) = 1/G

"""`issiso(sys)`

Returns `true` if `sys` is SISO, else returns `false`."""
function issiso(sys::LTISystem)
    return ninputs(sys) == 1 && noutputs(sys) == 1
end

"""`iscontinuous(sys)`

Returns `true` if `sys` is continuous, else returns `false`."""
function iscontinuous(sys::LTISystem)
    return sys.Ts == 0
end

"""`isstable(sys)`

Returns `true` if `sys` is stable, else returns `false`."""
function isstable(sys::LTISystem)
    if iscontinuous(sys)
        if any(real.(pole(sys)).>=0)
            return false
        end
    else
        if any(abs.(pole(sys)).>=1)
            return false
        end
    end
    return true
end

function _check_consistent_sampling_time(systems::AbstractVector{LTISystem})
    if !all(s.Ts == Ts for s in systems)
       error("Sampling time mismatch")
   end
end
function _check_consistent_sampling_time(sys1::LTISystem, sys2::LTISystem)
    if sys1.Ts != sys2.Ts
       error("Sampling time mismatch")
   end
end

# Fallback since LTISystem not AbstractArray
Base.size(sys::LTISystem, i::Integer) = size(sys)[i]


text:
include("interface.jl")
include("MGH.jl")


text:
@deprecate einsumexp(args...) loop_einsum(args...)
@deprecate einsumexp!(args...) loop_einsum!(args...)

@deprecate einsumopt einsum


text:
using Benchmark

################################################################################
# im2col is a bottleneck according to profiling. In this benchmark we try to
# compare the performance of C vs. julia in im2col. It turns out that C is a
# little bit faster, plus creating copy of sub-array could be avoided when
# processing im2col for each image in a mini-batch.
#
# | Row | Function    | Average     | Relative | Replications |
# |-----|-------------|-------------|----------|--------------|
# | 1   | "im2col_jl" | 0.000590041 | 1.77187  | 50           |
# | 2   | "im2col_c"  | 0.000333004 | 1.0      | 50           |
#
# Note if we add omp parallel for to the outer-most for loop, the performance
# deteriorate significantly.
# | Row | Function    | Average     | Relative | Replications |
# |-----|-------------|-------------|----------|--------------|
# | 1   | "im2col_jl" | 0.000831314 | 1.0      | 50           |
# | 2   | "im2col_c"  | 0.00514862  | 6.19335  | 50           |
################################################################################

function im2col{T}(img::Array{T}, col::Array{T}, width::Int, height::Int, channels::Int, kernel::NTuple{2,Int}, pad::NTuple{2,Int}, stride::NTuple{2,Int})
  kernel_w, kernel_h = kernel
  pad_w, pad_h = pad
  stride_w, stride_h = stride

  height_col = div(height + 2pad_h - kernel_h, stride_h) + 1
  width_col = div(width + 2pad_w - kernel_w, stride_w) + 1
  channels_col = channels * kernel_h * kernel_w

  for c = 0:channels_col-1
    w_offset = c % kernel_w
    h_offset = div(c, kernel_w) % kernel_h
    c_im = div(c, kernel_h * kernel_w) # channel
    for h = 0:height_col-1
      for w = 0:width_col-1
        h_pad = h*stride_h - pad_h + h_offset
        w_pad = w*stride_w - pad_w + w_offset
        if (h_pad >= 0 && h_pad < height && w_pad >= 0 && w_pad < width)
            @inbounds col[1 + (c*height_col+h) * width_col + w] =
                img[1 + (c_im * height + h_pad) * width + w_pad]
        else
          @inbounds col[1 + (c*height_col+h) * width_col + w] = 0
        end
      end
    end
  end
end

library = dlopen("./libextim2col.so")
func_handle = dlsym(library, :im2col)
function im2col_native(img::Array{Float64}, col::Array{Float64}, width::Int, height::Int, channels::Int, kernel::NTuple{2,Int}, pad::NTuple{2,Int}, stride::NTuple{2,Int})
  kernel_w, kernel_h = kernel
  pad_w, pad_h = pad
  stride_w, stride_h = stride

  ccall(func_handle, Void,
      (Ptr{Float64},Ptr{Float64}, Cint, Cint, Cint,
      Cint, Cint, # kernel
      Cint, Cint, # pad
      Cint, Cint, # stride
      ), img, col, width, height, channels, kernel_w, kernel_h, pad_w, pad_h, stride_w, stride_h)
end

############################################################
# Benchmark data preparation
############################################################
width, height, channels = (28, 28, 50)
kernel = (5,5)
pad = (2,2)
stride = (2,2)
img = rand(width, height, channels)

width_out  = div(width  + 2*pad[1]-kernel[1], stride[1]) + 1
height_out = div(height + 2*pad[2]-kernel[2], stride[2]) + 1
col_buffer = Array(Float64, width_out, height_out, channels*prod(kernel))
col_buffer2 = zeros(size(col_buffer))

im2col_jl() = im2col(img, col_buffer, width, height, channels, kernel, pad, stride)
im2col_c() = im2col_native(img, col_buffer2, width, height, channels, kernel, pad, stride)

im2col_jl()
im2col_c()

@assert all(abs(col_buffer-col_buffer2) .< 1e-10)

df = compare([im2col_jl, im2col_c], 50)
println("$df")


text:
using Pkg
Pkg.activate(".")

using AtomicLevels
using AtomicPotentials

using HTTP
using Gumbo
using AbstractTrees

const BASE_URL = "http://www.tc.uni-koeln.de/cgi-bin/pp.pl"

function ask_cologne(;kwargs...)
    kwargs = Dict(kwargs...)
    kwargs[:format] = get(kwargs, :format, "molpro")
    kwargs[:language] = get(kwargs, :language, "en")

    url = "$(BASE_URL)?"*join(["$(k)=$(w)" for (k,w) in kwargs], ",")
    r = HTTP.request("GET", url)
    if r.status == 200
        parsehtml(String(r.body))
    else
        throw(ArgumentError("Unexpected HTTP status: $(r.status)"))
    end
end

function get_ecps_list(element::String)
    doc = ask_cologne(job="listecps", element=element)

    ecps = String[]
    for elem in PreOrderDFS(doc.root)
        if elem isa HTMLElement{:a}
            a = attrs(elem)
            if "href" in keys(a) && occursin("getecp",a["href"])
                push!(ecps, first(Gumbo.children(elem)).text)
            end
        end
    end
    ecps
end

function core_model(element::String, Q::Integer)
    Z = element_number(Symbol(element))
    gst = ground_state(PointCharge(Z))
    o = 0
    for i in 1:length(gst)
        o += gst.occupancy[i]
        gst.states[i] = Z - o < Q ? :open : :closed
    end
    gst
end

function get_ecp(element::String, ecp::String)
    doc = ask_cologne(job="getecp", element=element, ecp=ecp)
    data = String[]
    for elem in PreOrderDFS(doc.root)
        if elem isa HTMLElement{:pre}
            push!(data, strip(first(Gumbo.children(elem)).text))
        end
    end
    m = match(r"Q=(.*?),", data[1])
    if !isnothing(m)
        Q = parse(Float64, m[1])
        if isinteger(Q)
            c = core_model(element, Int(Q))
            prepend!(data, ["! $(ascii(c))"])
        else
            @warn "Can't generate model configuration for non-integer Q = $Q"
        end
    end
    data = join(data, "\n")
    data
end

function get_ecps(element::String)
    ecps = get_ecps_list(element)
    if isempty(ecps)
        println("No ECPs available for $(element)")
        return
    end
    println(repeat("=", 100))
    println("""ECPnXY
n = number of core electrons
X = S/M: single/multi electron fit
Y = HF/WB/DF: non/quasi/fully relativistic""")
    for ecp in ecps
        println(repeat("=", 100))
        println("==> $(ecp)")
        println(get_ecp(element, ecp))
    end
    println(repeat("=", 100))
end


text:
module PKG

import Pkg

import ..Utils: current, download_cache, cachefile, replace_in_file

const CACHEDICT = Dict{String, Any}()

function get_registry_list(baseurl, update=false)
    registries_file = cachefile("registries.txt")
    if update || !ispath(registries_file)
        @info "Updating mirror information..."
        download_cache(joinpath(baseurl, "registries", "list.txt"), "registries.txt")
    end
    registries = open(registries_file) do io
        readlines(io)
    end
    return registries
end

function activate()
    baseurl = current().url
    registry_list = get_registry_list(baseurl, true)
    default_list = Pkg.Types.DEFAULT_REGISTRIES
    if typeof(default_list) <: Function
        default_list = default_list()
    end
    CACHEDICT["default_registries"] = deepcopy(default_list)
    for registry in default_list
        if registry.name in registry_list
            registry.url = joinpath(baseurl, "registries", "$(registry.name).git")
        end
    end
    for registry in Pkg.Types.collect_registries()
        git_config_file = joinpath(registry.path, ".git", "config")
        if registry.name in registry_list && isfile(git_config_file)
            i = findfirst(x -> x.name == registry.name, default_list)
            if i == nothing
                continue
            end
            replace_in_file(git_config_file, registry.url => default_list[i].url)
            toml_file = joinpath(registry.path, "Registry.toml")
            replace_in_file(toml_file, registry.url => default_list[i].url)
        end
    end

    @eval function Pkg.Operations.get_archive_url_for_version(url::String, ref)
        if (m = match(r"https://github.com/(.*?)/(.*?).jl.git", url)) != nothing
            mirror = current()
            # TODO: include registry info.
            pkgname = m.captures[2]
            return "$(mirror.url)/packages/$pkgname/General/$pkgname-$ref.tar.gz"
        end
        return nothing
    end

end

function deactivate()
    baseurl = current().url
    registry_list = get_registry_list(baseurl, false)
    registries = Pkg.Types.collect_registries()
    original_list = CACHEDICT["default_registries"]
    for registry in registries
        if registry.name in registry_list
            i = findfirst(x -> x.name == registry.name, original_list)
            if i == nothing
                continue
            end
            git_config_file = joinpath(registry.path, ".git", "config")
            replace_in_file(git_config_file, registry.url => original_list[i].url)
            toml_file = joinpath(registry.path, "Registry.toml")
            replace_in_file(toml_file, registry.url => original_list[i].url)
        end
    end
    for (i, x) in enumerate(original_list)
        Pkg.Types.DEFAULT_REGISTRIES[i] = x
    end

    @eval function Pkg.Operations.get_archive_url_for_version(url::String, ref)
        if (m = match(r"https://github.com/(.*?)/(.*?).git", url)) != nothing
            return "https://api.github.com/repos/$(m.captures[1])/$(m.captures[2])/tarball/$(ref)"
        end
        return nothing
    end
end

end


text:
module OpenCL

using LLVM
using LLVM.Interop
using Core: LLVMPtr

export cl
module cl

abstract type CLObject end

Base.hash(x::CLObject) = hash(pointer(x))
Base.isequal(x :: T, y :: T) where {T <: CLObject} = Base.hash(x) == Base.hash(y)
Base.:(==)(x :: T, y :: T) where {T <: CLObject} = Base.hash(x) == Base.hash(y)

# OpenCL Types
include("types.jl")

# The arrays contain a nullbyte that we pop first
function CLString(v :: Array{CL_char})
    pop!(v)
    String(reinterpret(UInt8, v))
end

# OpenCL Constants
include("constants.jl")

# OpenCL low level api
include("api.jl")

# Errors
include("error.jl")

# Macros
include("macros.jl")

# OpenCL Platform
include("platform.jl")

# OpenCL Device
include("device.jl")

# OpenCL Context
include("context.jl")

# OpenCL Queue
include("queue.jl")

# OpenCL Event
include("event.jl")

# OpenCL MemObject
include("memory.jl")

# OpenCL Buffer
include("buffer.jl")

# OpenCL Program
include("program.jl")

# OpenCL Kernel
include("kernel.jl")

# Util functions
include("util.jl")

get_global_id(dimindx::Integer=0) = @builtin_ccall("get_global_id", Csize_t, (Cuint,), dimindx) % Int + 1

# Multidimensional array
include("array.jl")

@deprecate release! finalize
end # cl

function __init__()
    if cl.api.libopencl == ""
        @warn "Could not locate an OpenCL library\nOpenCL API calls will be unavailable"
    end
end

end # module


text:
using ChaosTools
using ChaosTools.DynamicalSystemsBase
using Test
using StatsBase
using Statistics
using Entropies

test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax

println("\nTesting generalized entropy (genentropy) & linear scaling...")
@testset "Generalized Dimensions" begin
    @testset "Henon Map" begin
        ds = Systems.henon()
        ts = trajectory(ds, 200000)
        es = 10 .^ range(-0, stop = -3, length = 7)
        dd = zero(es)
        for q in [0,2,1, 2.56]
            for (i, ee) in enumerate(es)
                dd[i] = genentropy(ts, ee; q)
            end
            linr, dim = linear_region(-log.(es), dd)
            test_value(dim, 1.1, 1.3)
        end
    end
    @testset "Lorenz System" begin
        ds = Systems.lorenz()
        ts = trajectory(ds, 5000)
        es = 10 .^ range(1, stop = -3, length = 11)
        dd = zero(es)
        for q in [0,1,2]
            for (i, ee) in enumerate(es)
                dd[i] = genentropy(ts, ee; q)
            end
            linr, dim = linear_region(-log.(es), dd)
            if q == 0
                test_value(dim, 1.7, 1.9)
            else
                test_value(dim, 1.85, 2.0)
            end
        end
    end
end

println("\nTesting generalized entropy using Molteno's boxing method...")
@testset "Molteno's boxing method" begin
    @testset "Henon Map" begin
        ds = Systems.henon()
        ts = trajectory(ds, 200000)
        boxes, ϵs = molteno_boxing(ts)
        boxes = Probabilities.(boxes)
        for q in [0,2,1, 2.56]
            dd = genentropy.(boxes; q)
            linr, dim = linear_region(-log.(ϵs), dd)
            test_value(dim, 1.1, 1.3)
        end
    end
    @testset "Lorenz System" begin
        ds = Systems.lorenz()
        ts = trajectory(ds, 5000)
        boxes, ϵs = molteno_boxing(ts)
        boxes = Probabilities.(boxes)
        for q in [0,2,1, 2.56]
            dd = genentropy.(boxes; q)
            linr, dim = linear_region(-log.(ϵs), dd)
            test_value(dim, 1.85, 2.2)
        end
    end
end

println("\nTesting dimension calls (all names)...")
@testset "Dimension calls" begin
    ds = Systems.henon()
    ts = trajectory(ds, 50000)
    test_value(generalized_dim(ts; q = 1.32), 1.0, 1.3)
    test_value(generalized_dim(ts; base = 2, q = 1.32), 1.0, 1.3)
end


text:
module Chatterjee

    include("correlation_coefficient.jl")
    export chatterjee_correlation_symmetric, chatterjee_correlation, max_value_chatterjee_correlation, min_value_chatterjee_correlation, p_value_asymptotic_null, p_value_asymptotic_null_non_continuous, chatterjee_std_dev_asymptotic, chatterjee_std_dev_asymptotic_non_continuous

end


text:
using Revise
using Jtac

G = MetaTac

#chain = @chain G Conv(256) Pool() Deconv(512) Conv(512) Conv(512) Dense(82)
#chain = @chain G Dense(1000) Dense(1000) Dense(1000) Dense(1000) Dense(Jtac.policy_length(G)+1)

model = ShallowConv(G, 200)
#model = NeuralModel(G, chain)
amodel = Async(model, max_batchsize = 100)

games = map(1:10000) do _
  game = G()
  for i in 1:rand(1:20)
    if !Jtac.is_over(game)
      Jtac.random_turn!(game)
    end
  end
  game
end

println("CPU performance for 10000 games")
@time model(games);
@time map(model, games);
@time asyncmap(amodel, games, ntasks = Jtac.ntasks(amodel));

gmodel = model |> to_gpu
gamodel = Async(gmodel, max_batchsize = 100)

println()
println("GPU performance for 10000 games")
@time gmodel(games)
@time map(gmodel, games)
@time asyncmap(gamodel, games, ntasks = Jtac.ntasks(gamodel))

max_stars_repo_name: roSievers/jtac


text:
mordat = [ population_deaths(popdata;alter=a, geschlecht=g, jahr=j, deaths=(kw=kw,) )
           # or switch to interpolate deaths for kw:
           # population_deaths(alter=a, geschlecht=g, jahr=j, kw=kw )
           for a in DeutschlandDaten.Altersgruppen_KW,
               g in ["Männlich", "Weiblich"],
               j in 2000:2021,
               kw in 1:52 ] |> DataFrame;
filter!(x-> x.D !==missing && x.D>0, mordat)

@transform! mordat :date = endofweek.(:jahr,:kw)
dataprep!(mordat)

sort!(mordat, [:alter, :geschlecht, :date])
@transform! mordat :PD_12months = gliding_mean(:PD, div(52, 2))
@transform! mordat :PD_quarter = gliding_mean(:PD, div(52, 4))

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
module _Combinatorics

using Brainstorm._Math._Combinatorics
using Base.Test

include("permutation.jl")
include("combination.jl")
include("multipermutation.jl")
include("multicombination.jl")

function test_all()
    test_permutation_all()
    test_combination_all()
    test_multipermutation_all()
    test_multicombination_all()
end

end


text:
# OptimalTransport.jl -- optimal transportation algorithms for Julia
# See prettyprinted documentation at https://juliaoptimaltransport.github.io/OptimalTransport.jl/dev/
#

module OptimalTransport

using Reexport

@reexport using ExactOptimalTransport:
    emd, emd2, ot_cost, ot_plan, wasserstein, squared2wasserstein, discretemeasure

using LinearAlgebra
using IterativeSolvers
using LogExpFunctions: LogExpFunctions
using NNlib: NNlib

export SinkhornGibbs, SinkhornStabilized, SinkhornEpsilonScaling
export SinkhornBarycenterGibbs
export QuadraticOTNewton

export sinkhorn, sinkhorn2
export sinkhorn_stabilized, sinkhorn_stabilized_epsscaling, sinkhorn_barycenter
export sinkhorn_unbalanced, sinkhorn_unbalanced2
export sinkhorn_divergence
export quadreg

include("utils.jl")

include("entropic/sinkhorn.jl")
include("entropic/sinkhorn_divergence.jl")
include("entropic/sinkhorn_gibbs.jl")
include("entropic/symmetric.jl")
include("entropic/sinkhorn_stabilized.jl")
include("entropic/sinkhorn_epsscaling.jl")
include("entropic/sinkhorn_unbalanced.jl")
include("entropic/sinkhorn_barycenter.jl")
include("entropic/sinkhorn_barycenter_gibbs.jl")
include("entropic/sinkhorn_solve.jl")

include("quadratic.jl")
include("quadratic_newton.jl")

include("dual/entropic_dual.jl")

end


text:
function mark_angles(r)
    num_steps = through_rectangles_into_angle(r,(Sud,West))
    for side in (Nord,Ost,Sud,West)
        moves!(r,side)
        putmarker!(r)
    end
    move!(r,(Ost,Nord),num_steps)
end

function through_rectangles_into_angle(r,angle::NTuple{2,HorizonSide})
    num_steps=[]
    while isborder(r,angle[1])==false || isborder(r,angle[2])
        push!(num_steps, moves!(r, angle[2]))
        push!(num_steps, moves!(r, angle[1]))
    end
    return num_steps
end

function moves!(r,sides,num_steps::Vector{Int64})
    for (i,n) in enumerate(reverse!(num_steps))
        moves!(r, sides[mod(i-1, length(sides))+1], n)
    end
end

max_stars_repo_name: Anpilogova/-

max_issues_repo_name: Anpilogova/-

max_forks_repo_name: Anpilogova/-


text:
module LaTeXChanges

export apply_changes

"""
    apply_changes(input_file, output_file = input_file)

Read the `input_file`, apply the changes from the [LaTeX changes](https://ctan.org/pkg/changes) package and write the result to `output_file`.
"""
function apply_changes(input_file, output_file = input_file)
    str_in = join(readlines(input_file), '\n')
    str_clean = _apply_changes(str_in)
    of = open(output_file, "w")
    write(of, str_clean)
    close(of)
end

function _apply_changes(str)
    @info "Applying \\added[]{}"
    str = apply_added(str)
    @info "Applying \\deleted[]{}"
    str = apply_deleted(str)
    @info "Applying \\replaced[]{}{}"
    str = apply_replaced(str)
end

function apply_added(str)
    newstr = ""

    inds = findfirst(r"\\added", str)

    while inds !== nothing
        i_open = opening_index(str, nextind(str, last(inds)))
        i_close = closing_index(str, nextind(str, i_open))

        before = str[firstindex(str):prevind(str, first(inds))]
        between = str[nextind(str, i_open):prevind(str, i_close)]
        after = str[nextind(str, i_close):lastindex(str)]

        newstr *= before * between
        str = after

        inds = findfirst(r"\\added", str)
    end

    return newstr * str
end

function apply_deleted(str)
    newstr = ""

    inds = findfirst(r"\\deleted", str)

    while inds !== nothing
        i_open = opening_index(str, nextind(str, last(inds)))
        i_close = closing_index(str, nextind(str, i_open))

        before = str[firstindex(str):prevind(str, first(inds))]
        after = str[nextind(str, i_close):lastindex(str)]

        newstr *= before
        str = after

        inds = findfirst(r"\\deleted", str)
    end

    return newstr * str
end

function apply_replaced(str)
    newstr = ""

    inds = findfirst(r"\\replaced", str)

    while inds !== nothing
        i_open_add = opening_index(str, nextind(str, last(inds)))
        i_close_add = closing_index(str, nextind(str, i_open_add))

        i_open_del = nextind(str, i_close_add)
        i_close_del = closing_index(str, nextind(str, i_open_del))

        before = str[firstindex(str):prevind(str, first(inds))]
        between = str[nextind(str, i_open_add):prevind(str, i_close_add)]
        after = str[nextind(str, i_close_del):lastindex(str)]

        newstr *= before * between
        str = after

        inds = findfirst(r"\\replaced", str)
    end

    return newstr * str
end

function opening_index(str, i = firstindex(str))
    n = 0
    while i <= lastindex(str)
        if str[i] == '['
            n += 1
        elseif str[i] == ']' && n > 0
            n -= 1
        elseif str[i] == '{' && n == 0
            return i
        end
        i = nextind(str, i)
    end
    error("No opening index found in `str`.")
end

function closing_index(str, i = firstindex(str))
    n = 0
    lastind = lastindex(str)
    while i <= lastind
        if str[i] == '{'
            n += 1
        elseif str[i] == '}'
            if n == 0
                return i
            else
                n -= 1
            end
        end
        i = nextind(str, i)
    end
    error("No closing index found in `str`.")
end

end # module


text:
import GeometryBasics: Point

function get_dims(box::S.Box)
    [box.sizeX, box.sizeY, box.sizeZ]
end

function axis_aligned_bbox_from_point_cloud(point_cloud::Matrix)
    x_bounds = min_max(point_cloud[1, :])
    y_bounds = min_max(point_cloud[2, :])
    z_bounds = min_max(point_cloud[3, :])

    return Pose([mean(x_bounds), mean(y_bounds), mean(z_bounds)], IDENTITY_ORIENTATION),
        S.Box(x_bounds[2] - x_bounds[1],
                            y_bounds[2] - y_bounds[1],
                            z_bounds[2] - z_bounds[1])
end

function get_bbox_corners(bbox::S.Box)
    x,y,z = get_dims(bbox) ./ 2
    nominal_corners = collect([
        -x -y -z;
        -x y -z;
        x y -z;
        x -y -z;
        -x -y  z;
        -x y  z;
        x y  z;
        x -y  z;
    ]')
    nominal_corners
end

function get_bbox_segments_point_list(bbox::S.Box, pose::Pose)
    x,y,z = get_dims(bbox) ./ 2
    nominal_corners = collect([
        # bottom square
        x y -z;
        x -y -z;

        x -y -z;
        -x -y -z;

        -x -y -z;
        -x y -z;

        -x y -z;
        x y -z;

        # top square
        x y z;
        x -y z;

        x -y z;
        -x -y z;

        -x -y z;
        -x y z;

        -x y z;
        x y z;

        # top to bottom
        -x -y z;
        -x -y -z;

        x -y z;
        x -y -z;

        -x y z;
        -x y -z;

        x y z;
        x y -z;
    ]')
    corners = move_points_to_frame_b(nominal_corners, pose)
    Point.(corners[1,:],corners[2,:], corners[3,:])
end

export axis_aligned_bbox_from_point_cloud, get_bbox_corners, get_bbox_segments_point_list


text:
"""
    regression(data::Dataset, labels::Labels; kwargs...)

perform learning of logistic regression classifier on given `data` and correct `labels`

outputs vector `w` which detrmines a hyperplane optimally separating the data samples

- `kwargs` can contain `epsilon` to specify accuracy of optimization method and initial `step` size
"""
function regression(data::Dataset, labels::Labels; kwargs...)
    trainX = data()
    trainy = labels()

    w = gd(trainX, trainy; kwargs...)

    return w
end

"""
    classify(data::Dataset, w::Vector{<:Real})

classify the `data` using separating hyperplane determined by `w`
"""
function classify(data::Dataset, w::Vector{<:Real})
    X = data()
    _, n = size(X)
    return Labels(data, [dot(X[:, i], w) >= 0 ? 1 : -1 for i in 1:n])
end

"""
    computeE(X::Matrix{<:Real}, k::Vector{<:Integer}, w::Vector{<:Real})

returns energy of hyperplane `w` with respect to data `X` and labels `k`; this objective is to be minimized
"""
function computeE(X::Matrix{<:Real}, k::Vector{<:Integer}, w::Vector{<:Real})
    N = length(k)
    return sum([log(1 + exp(-k[i] * dot(X[:, i], w))) for i in 1:N]) ./ N
end

"""
    gradE(X::Matrix{<:Real}, k::Vector{<:Integer}, w::Vector{<:Real})

returns gradient of energy with respect to vector w
"""
function gradE(X::Matrix{<:Real}, k::Vector{<:Integer}, w::Vector{<:Real})
    N = length(k)
    return - sum([(k[i] * X[:, i]) / (1 + exp(k[i] * dot(X[:, i], w))) for i in 1:N]; dims=1)[1] ./ N
end

"""
    l2norm(x)

computes l2-norm of given vector `x`
"""
function l2norm(x)
    return sqrt.(dot(x, x))
end

"""
    gd(X::Matrix{<:Real}, k::Vector{<:Integer}; epsilon=1e-4, step=1.0)

performs gradient descent method to find optimal separating hyperplane `w` on data `X` and labels `k`

keyword arguments can be used to alter the optimization:
- `epsilon` to specify accuracy/terminal condition
- `step` to adjust initial step size
"""
function gd(X::Matrix{<:Real}, k::Vector{<:Integer}; epsilon=1e-4, step=1.0)
    dim, _ = size(X)

    w = zeros(Float64, dim)
    lastw = fill(Inf64, dim)
    E = Inf64

    while l2norm(w - lastw) > epsilon
        grad = gradE(X, k, w)
        neww = w - step * grad
        newE = computeE(X, k, neww)

        if newE < E
            step *= 2
            E = newE
            lastw = w
            w = neww
        else
            step /= 2
        end

    end

    return w
end


text:
using IscaTools
using Test

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


text:
for N in [Float64, Rational{Int}, Float32]
    # overapproximating a set of type T1 with an unsupported type T2 is the
    # identity if T1 = T2
    @test_throws MethodError overapproximate(ZeroSet{N}(2), EmptySet)
    e = EmptySet{N}()
    @test overapproximate(e, EmptySet) == e

    # HPolygon approximation with box directions
    c = N[0, 0]
    b = Ball1(c, N(1))
    p = overapproximate(b, HPolygon)
    for d in [N[1, 0], N[-1, 0]]
        @test σ(d, p)[1] ≈ σ(d, b)[1]
    end
    for d in [N[0, 1], N[0, -1]]
        @test σ(d, p)[2] ≈ σ(d, b)[2]
    end

    # Hyperrectangle approximation
    c = N[0, 0]
    b = Ball1(c, N(1))
    p = overapproximate(b, Hyperrectangle)
    for d in [N[1, 0], N[-1, 0]]
        @test σ(d, p)[1] ≈ σ(d, b)[1]
    end
    for d in [N[0, 1], N[0, -1]]
        @test σ(d, p)[2] ≈ σ(d, b)[2]
    end
    @test p.center ≈ c
    @test p.radius ≈ N[1, 1]

    # Interval approximation
    b = Ball1(N[0], N(1))
    p = overapproximate(b, LazySets.Interval)
    for d in [N[1], N[-1]]
        @test σ(d, p)[1] ≈ σ(d, b)[1]
    end

    # approximation with an axis-aligned hyperrectangle
    Z = Zonotope(N[-1.0, -1.0], N[-1/2 0; -1.2 -1])
    Zoa = overapproximate(Z, Hyperrectangle) # faster o.a.
    Zba = box_approximation(Z) # default o.a. implementation that uses supp function
    @test Zoa.center ≈ Zba.center && Zoa.radius ≈ Zba.radius

    # same but using cartesian product
    h1 = Hyperrectangle(N[1/2],  N[1/2])
    h2 = Hyperrectangle(N[2.5, 4.5],  N[1/2, 1/2])
    H = overapproximate(h1 × h2, Hyperrectangle) # defaults to convert method
    @test low(H) == N[0, 2, 4] && high(H) == N[1, 3, 5]

    # overapproximation of the lazy linear map of a hyperrectangular set
    H = Hyperrectangle(N[0, 0], N[1/2, 1])
    M = Diagonal(N[2, 2])
    OA = overapproximate(M*H, Hyperrectangle)
    @test OA isa Hyperrectangle && OA.center == N[0, 0] && OA.radius == N[1, 2]
end

# tests that do not work with Rational{Int}
for N in [Float64, Float32]
    # useful for benchmarking overapproximate and LinearMap's support vector
    # (see #290)
    function overapproximate_lmap(n)
        B = BallInf(ones(N, n), N(2))
        π = sparse(N[1, 2], N[1, 2], ones(N, 2), 2, n)
        return Approximations.overapproximate(π*B)
    end
    o = overapproximate_lmap(50)
    @test o.center == N[1, 1] && o.radius == N[2, 2]

    # Approximation of a 2D centered unit ball in norm 1
    # All vertices v should be like this:
    # ‖v‖ >= 1 and ‖v‖ <= 1+ε
    # Where ε is the given error bound
    b = Ball1(N[0, 0], N(1))
    ε = N(0.01)
    p = tovrep(overapproximate(b, ε))
    for v in vertices_list(p)
    @test norm(v) >= N(1)
    @test norm(v) <= N(1 + ε)
    end

    # Check that there are no redundant constraints for a ballinf
    b = BallInf(N[0.5, 0.5], N(0.1))
    lcl = overapproximate(b, N(0.001)).constraints
    @test length(lcl) == 4
    @test lcl[1].a == N[1, 0]
    @test lcl[1].b == N(0.6)
    @test lcl[2].a == N[0, 1]
    @test lcl[2].b == N(0.6)
    @test lcl[3].a == N[-1, 0]
    @test lcl[3].b == N(-0.4)
    @test lcl[4].a == N[0, -1]
    @test lcl[4].b == N(-0.4)

    # Check that there are no redundant constraints for a HPolygon (octagon)
    p = HPolygon{N}()
    addconstraint!(p, LinearConstraint(N[1, 0], N(1)))
    addconstraint!(p, LinearConstraint(N[0, 1], N(1)))
    addconstraint!(p, LinearConstraint(N[-1, 0], N(1)))
    addconstraint!(p, LinearConstraint(N[0, -1], N(1)))
    addconstraint!(p, LinearConstraint(N[sqrt(2)/2, sqrt(2)/2], N(1)))
    addconstraint!(p, LinearConstraint(N[-sqrt(2)/2, sqrt(2)/2], N(1)))
    addconstraint!(p, LinearConstraint(N[sqrt(2)/2, -sqrt(2)/2], N(1)))
    addconstraint!(p, LinearConstraint(N[-sqrt(2)/2, -sqrt(2)/2], N(1)))
    lcl = overapproximate(p, N(.001)).constraints
    @test length(lcl) == 8
    @test lcl[1].a ≈ N[1, 0]
    @test lcl[1].b ≈ N(1)
    @test lcl[2].a ≈ N[sqrt(2)/2, sqrt(2)/2]
    @test lcl[2].b ≈ N(1)
    @test lcl[3].a ≈ N[0, 1]
    @test lcl[3].b ≈ N(1)
    @test lcl[4].a ≈ N[-sqrt(2)/2, sqrt(2)/2]
    @test lcl[4].b ≈ N(1)
    @test lcl[5].a ≈ N[-1, 0]
    @test lcl[5].b ≈ N(1)
    @test lcl[6].a ≈ N[-sqrt(2)/2, -sqrt(2)/2]
    @test lcl[6].b ≈ N(1)
    @test lcl[7].a ≈ N[0, -1]
    @test lcl[7].b ≈ N(1)
    @test lcl[8].a ≈ N[sqrt(2)/2, -sqrt(2)/2]
    @test lcl[8].b ≈ N(1)

    # Zonotope approximation
    Z1 = Zonotope(ones(N, 2), [N[1, 0], N[0, 1], N[1, 1]])
    Z2 = Zonotope(-ones(N, 2), [N[0.5, 1], N[-0.1, 0.9], N[1, 4]])
    Y = ConvexHull(Z1, Z2)
    Y_polygon = overapproximate(Y, N(1e-3)) # overapproximate with a polygon
    Y_zonotope = overapproximate(Y, Zonotope) # overapproximate with a zonotope
    @test Y_polygon ⊆ Y_zonotope
    @test !(Y_zonotope ⊆ Y_polygon)
end


text:
################################################################################

"""Abstract type representing the ordering of pixels in a Healpix map.
See also `Ring` and `Nest`.
"""
abstract type Order end

"""The `Ring` type should be used when creating `Map` types in
order to specify that the pixels in the map are sorted in ring
ordering. (See also `Nest`.)
"""
abstract type Ring <: Order end

"""The `Nest` type should be used when creating `Map` types in
order to specify that the pixels in the map are sorted in ring
ordering. (See also `Ring`.)
"""
abstract type Nest <: Order end

"""
    AbstractMap{T, O} <: AbstractArray{T, 1}

An abstract type representing an Healpix map without a specified
ordering. This can be used to implement multiple dispatch when you
don't care about the ordering of a map."""
abstract type AbstractMap{T, O} <: AbstractArray{T, 1} end
abstract type OrderedMap{T, O <: Order} <: AbstractMap{T, O} end

"""
    struct Map{T, O <: Order} <: AbstractMap{T, O}

A Healpix map. The type `T` is used for the value of the pixels in
a map, and it can be anything (even a string!). The type `O` is used
to specify the ordering of the pixels, and it can either be
`Ring` or `Nest`.

A `Map` type contains the following fields:

- `pixels`: array of pixels
- `resolution`: instance of a `Resolution` object

"""
struct Map{T, O} <: OrderedMap{T, O}
    pixels::Array{T}
    resolution::Resolution
end

"""
    Map{T, O <: Order}(nside) -> Map{T, O}

Create an empty map with the specified NSIDE.
"""
function Map{T, O}(nside::Number) where {T, O}
	Map{T,O}(zeros(T, nside2npix(nside)), Resolution(nside))
end

function Map{T, O}(healpixels::Vector{T}) where {T, O}
	Map{T, O}(healpixels, Resolution(npix2nside(length(healpixels))))
end

"""
Create a map with the specified array of pixels, without redundancy of above
syntax.
"""
function Map{O}(healpixels::Vector{T}) where {T, O}
	Map{T, O}(healpixels, Resolution(npix2nside(length(healpixels))))
end

"""
Convenience constructor--default Nest
"""
Map(healpixels::Vector{T}) where T = Map{Nest}(healpixels)

###########################################################################
# Map & MaskedMap initialization functions

struct MaskedMap{T, O} <: OrderedMap{T, O}
	pixels::Vector{T}
	inds::Vector{Int64}
	resolution::Resolution
end

# MaskedMap to Map
function Map(mm::MaskedMap{T, O}) where {T, O}
	healpixels = fill(NaN, mm.resolution.numOfPixels)
	healpixels[mm.inds] = mm.pixels
	Map{O}(healpixels)
end

# Map to MaskedMap
function MaskedMap(m::Map{T, O}, θlims, ϕlims) where {T,O}
	inds = Vector{Int64}()
	for i in eachindex(m.pixels)
		θ, ϕ  = get_latlong_deg(m, i)
		if (θlims[1] <= θ <= θlims[2]) & (ϕlims[1] <= ϕ <= ϕlims[2])
			push!(inds, i)
		end
	end
	MaskedMap{T,O}(m.pixels[inds], inds, m.resolution)
end

# Pixel vector to MaskedMap
function MaskedMap{O}(healpixels::Vector{T}, θlims, ϕlims) where {T, O}
	MaskedMap(Map{O}(healpixels), θlims, ϕlims)
end

# Pixel vector to defaults to ringOrdering
function MaskedMap(healpixels::Vector{T}, θlims, ϕlims) where T <: Number
	MaskedMap{Nested}(healpixels, θlims, ϕlims)
end

###########################################################################
# reordering and coordinate querying

Base.eachindex(mm::MaskedMap) = mm.inds
function reorder_inds(m::AbstractMap{T,O}) where {T,O}
    reorder_func = (O ≡ Ring ? ring2nest : nest2ring)
    sortperm([reorder_func(m.resolution, i) for i in eachindex(m)])
end

function reorder(m::Map{T,O}) where {T,O}
    Map{T, (O ≡ Ring ? Nest : Ring)}(m.pixels[reorder_inds(m)])
end
function reorder(mm::MaskedMap{T,O}) where {T,O}
    inds = reorder_inds(mm)
    MaskedMap{T, (O ≡ Ring ? Nest : Ring)}(mm.pixels[inds], inds, mm.resolution)
end

function get_latlong_deg(m::OrderedMap, ind)
	colat, lon = pix2ang(m, ind)
	rad2deg.((colat2lat(colat), lon))
end

function getcoords(m::OrderedMap)
	coords = zeros((2, length(m)))
	for (i, healpixel_ind) in enumerate(eachindex(m))
		coords[:,i] .= get_latlong_deg(m, healpixel_ind)
	end
	coords
end

function resize(m::Map{T,O}, nside) where {T, O}
    order =
	Map{O}(healpy.ud_grade(m.pixels, nside, order_in=(O ≡ Ring ? "RING" : "NEST")))
end
function resize(mm::MaskedMap{T,O}, nside) where {T, O}
	θlims, ϕlims = extrema(getcoords(mm), dims=2)
	MaskedMap(resize(Map(mm), nside), θlims, ϕlims)
end

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

import Base: +, -, *, /

+(a::AbstractMap{T, O}, b::AbstractMap{T, O}) where {T <: Number, O} = AbstractMap{T, O}(a.pixels .+ b.pixels)
-(a::AbstractMap{T, O}, b::AbstractMap{T, O}) where {T <: Number, O} = AbstractMap{T, O}(a.pixels .- b.pixels)
*(a::AbstractMap{T, O}, b::AbstractMap{T, O}) where {T <: Number, O} = AbstractMap{T, O}(a.pixels .* b.pixels)
/(a::AbstractMap{T, O}, b::AbstractMap{T, O}) where {T <: Number, O} = AbstractMap{T, O}(a.pixels ./ b.pixels)

+(a::AbstractMap{T, O}, b::Number) where {T <: Number, O} = AbstractMap{T, O}(a.pixels .+ b)
-(a::AbstractMap{T, O}, b::Number) where {T <: Number, O} = a + (-b)
*(a::AbstractMap{T, O}, b::Number) where {T <: Number, O} = AbstractMap{T, O}(a.pixels .* b)
/(a::AbstractMap{T, O}, b::Number) where {T <: Number, O} = AbstractMap{T, O}(a.pixels ./ b)

+(a::Number, b::AbstractMap{T, O}) where {T <: Number, O} = b + a
-(a::Number, b::AbstractMap{T, O}) where {T <: Number, O} = b + (-a)
*(a::Number, b::AbstractMap{T, O}) where {T <: Number, O} = b * a
/(a::Number, b::AbstractMap{T, O}) where {T <: Number, O} = AbstractMap{T, O}(a ./ b.pixels)

################################################################################
# Iterator interface

Base.size(m::AbstractMap{T, O}) where {T, O} = (length(m.pixels),)

Base.IndexStyle(::Type{<:AbstractMap{T, O}}) where {T, O} = IndexLinear()

function getindex(m::AbstractMap{T, O}, i::Integer) where {T, O}
    1 ≤ i ≤ m.resolution.numOfPixels || throw(BoundsError(m, i))
    m.pixels[i]
end

function setindex!(m::AbstractMap{T, O}, val, i::Integer) where {T, O}
    1 ≤ i ≤ m.resolution.numOfPixels || throw(BoundsError(m, i))
    m.pixels[i] = val
end

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

"""
    conformables{T, S, O1 <: Order, O2 <: Order}(map1::Map{T, O1},
                                                 map2::Map{S, O2}) -> Bool

Determine if two Healpix maps are "conformables", i.e., if their
shape and ordering are the same.
"""
conformables(map1::OrderedMap{T, Ring}, map2::OrderedMap{S, Ring}) where {T, S} =
    map1.resolution.nside == map2.resolution.nside

conformables(map1::OrderedMap{T, Nest}, map2::OrderedMap{S, Nest}) where {T, S} =
    map1.resolution.nside == map2.resolution.nside

conformables(map1::OrderedMap{T, O1},
             map2::OrderedMap{S, O2}) where {T, S, O1 <: Order, O2 <: Order} = false


text:
module Periodic_Cov
using Statistics, LinearAlgebra
using GaussianRandomFields
import GaussianRandomFields: shortname, apply, apply_symmetric!, apply_non_symmetric!

abstract type Distance end

shortname(::Distance) = "custom distance"
function apply(::Distance, x::NTuple{d, Real}, y::NTuple{d, Real}) where d end

struct CustomDistanceCovarianceStructure{T} <: CovarianceStructure{T}
    distance::Distance
    cov::IsotropicCovarianceStructure{T}
end

apply(s::CustomDistanceCovarianceStructure, ρ::Real) = apply(s.cov, ρ)

shortname(s::CustomDistanceCovarianceStructure) = shortname(s.cov) * " on " * shortname(s.distance)

Statistics.std(c::IsotropicCovarianceStructure) = c.σ
Statistics.std(s::CustomDistanceCovarianceStructure) = Statistics.std(s.cov)
Statistics.std(f::AbstractCovarianceFunction) = Statistics.std(f.cov)

function apply_symmetric!(
    C::Matrix,
    s::CustomDistanceCovarianceStructure,
    x::NTuple{d, AbstractVector},
    y::NTuple{d, AbstractVector}
) where d
    xiterator = enumerate(Iterators.product(x...))
    for (j, idy) in enumerate(Iterators.product(y...))
        for (i, idx) in Iterators.take(xiterator, j)
            @inbounds C[i, j] = apply(s.cov, apply(s.distance, idx, idy))
        end
    end
    Symmetric(C, :U)
end

function apply_non_symmetric!(
    C::Matrix,
    s::CustomDistanceCovarianceStructure,
    x::NTuple{d, AbstractVector},
    y::NTuple{d, AbstractVector}
) where d
    xiterator = enumerate(Iterators.product(x...))
    for (j, idy) in enumerate(Iterators.product(y...))
        for (i, idx) in xiterator
            @inbounds C[i, j] = apply(s.cov, apply(s.distance, idx, idy))
        end
    end
    C
end

struct SphericalDistance <: Distance end

shortname(::SphericalDistance) = "spherical distance"

function apply(::SphericalDistance, (λ₁, φ₁)::Tuple{Float64, Float64}, (λ₂, φ₂)::Tuple{Float64, Float64})
    Δλ = λ₂ - λ₁  # longitude difference
    Δφ = φ₂ - φ₁  # latitude difference
    # haversine formula
    a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2
    # distance
    2 * ( asin( min(√a, one(a)) ))
end

struct EuclideanDistance <: Distance end
shortname(::EuclideanDistance) = "Euclidean distance"
apply(::EuclideanDistance, x::NTuple{d, Real}, y::NTuple{d, Real}) where d = norm(x .- y)
end

max_stars_repo_name: Team-RADDISH/ParticleDA-UseCases

max_issues_repo_name: Team-RADDISH/ParticleDA-UseCases

max_forks_repo_name: Team-RADDISH/ParticleDA-UseCases


text:
import GenDirectionalStats: uniform_rot3
import Rotations

struct UniformPose <: Gen.Distribution{Pose} end
const uniformPose = UniformPose()
(::UniformPose)(xmin, xmax, ymin, ymax, zmin, zmax) = Gen.random(
    uniformPose, xmin, xmax, ymin, ymax, zmin, zmax)

function Gen.random(::UniformPose,
                    xmin::Real, xmax::Real,
                    ymin::Real, ymax::Real,
                    zmin::Real, zmax::Real,
                   )::Pose
  Pose([Gen.uniform(xmin, xmax),
        Gen.uniform(ymin, ymax),
        Gen.uniform(zmin, zmax)],
        uniform_rot3())
end
function Gen.logpdf(::UniformPose,
                    val::Pose,
                    xmin::Real, xmax::Real,
                    ymin::Real, ymax::Real,
                    zmin::Real, zmax::Real)
  +(Gen.logpdf(Gen.uniform, val.pos[1], xmin, xmax),
    Gen.logpdf(Gen.uniform, val.pos[2], ymin, ymax),
    Gen.logpdf(Gen.uniform, val.pos[3], zmin, zmax),
    Gen.logpdf(uniform_rot3, Rotations.UnitQuaternion(val.orientation)))
end
Gen.has_output_grad(::UniformPose) = false
Gen.has_argument_grads(::UniformPose) = fill(false, 6) |> Tuple
function Gen.logpdf_grad(::UniformPose,
                         xmin::Real, xmax::Real,
                         ymin::Real, ymax::Real,
                         zmin::Real, zmax::Real)
  fill(nothing, 7) |> Tuple
end

export uniformPose


text:
export AbstractLinearOperator,
  LinearOperator,
  LinearOperatorException,
  hermitian,
  ishermitian,
  symmetric,
  issymmetric,
  nprod,
  ntprod,
  nctprod,
  reset!,
  shape

mutable struct LinearOperatorException <: Exception
  msg::AbstractString
end

# when indexing, Colon() is treated separately
const LinearOperatorIndexType{I} = Union{UnitRange{I}, StepRange{I, I}, AbstractVector{I}} where {I<:Integer}

# import methods we overload
import Base.eltype, Base.isreal, Base.size, Base.show
import LinearAlgebra.issymmetric, LinearAlgebra.ishermitian

abstract type AbstractLinearOperator{T} end
OperatorOrMatrix = Union{AbstractLinearOperator, AbstractMatrix}

eltype(A::AbstractLinearOperator{T}) where {T} = T
isreal(A::AbstractLinearOperator{T}) where {T} = T <: Real

"""
Base type to represent a linear operator.
The usual arithmetic operations may be applied to operators
to combine or otherwise alter them. They can be combined with
other operators, with matrices and with scalars. Operators may
be transposed and conjugate-transposed using the usual Julia syntax.
"""
mutable struct LinearOperator{T, S, I<:Integer, F, Ft, Fct} <: AbstractLinearOperator{T}
  nrow::I
  ncol::I
  symmetric::Bool
  hermitian::Bool
  prod!::F
  tprod!::Ft
  ctprod!::Fct
  Mv::S # storage vector for prod!
  Mtu::S # storage vector for tprod!
  Maw::S # storage vector for ctprod!
  nprod::I
  ntprod::I
  nctprod::I
end

LinearOperator{T}(nrow::I, ncol::I, symmetric::Bool, hermitian::Bool,
                  prod!::F, tprod!::Ft, ctprod!::Fct,
                  Mv::S, Mtu::S, Maw::S,
                  nprod::I, ntprod::I, nctprod::I
                  ) where {T,S,I<:Integer,F,Ft,Fct} = LinearOperator{T,S,I,F,Ft,Fct}(nrow, ncol, symmetric, hermitian,
                                                                                prod!, tprod!, ctprod!,
                                                                                Mv, Mtu, Maw,
                                                                                nprod, ntprod, nctprod)

LinearOperator{T}(
  nrow::I,
  ncol::I,
  symmetric::Bool,
  hermitian::Bool,
  prod!,
  tprod!,
  ctprod!,
  Mv::S,
  Mtu::S,
  Maw::S
) where {T,S,I<:Integer} = LinearOperator{T}(nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, Mv, Mtu, Maw, 0, 0, 0)

nprod(op::AbstractLinearOperator) = op.nprod
ntprod(op::AbstractLinearOperator) = op.ntprod
nctprod(op::AbstractLinearOperator) = op.nctprod

increase_nprod(op::AbstractLinearOperator) = (op.nprod += 1)
increase_ntprod(op::AbstractLinearOperator) = (op.ntprod += 1)
increase_nctprod(op::AbstractLinearOperator) = (op.nctprod += 1)

"""
  reset!(op)

Reset the product counters of a linear operator.
"""
function reset!(op::AbstractLinearOperator)
  op.nprod = 0
  op.ntprod = 0
  op.nctprod = 0
  return op
end

"""
    m, n = size(op)

Return the size of a linear operator as a tuple.
"""
size(op::AbstractLinearOperator) = (op.nrow, op.ncol)

"""
    m = size(op, d)

Return the size of a linear operator along dimension `d`.
"""
function size(op::AbstractLinearOperator, d::Integer)
  nrow, ncol = size(op)
  if d == 1
    return nrow
  end
  if d == 2
    return ncol
  end
  throw(LinearOperatorException("Linear operators only have 2 dimensions for now"))
end

"""
    m, n = shape(op)

An alias for size.
"""
shape(op::AbstractLinearOperator) = size(op)

"""
    hermitian(op)
    ishermitian(op)

Determine whether the operator is Hermitian.
"""
hermitian(op::AbstractLinearOperator) = op.hermitian
ishermitian(op::AbstractLinearOperator) = op.hermitian

"""
    symmetric(op)
    issymmetric(op)

Determine whether the operator is symmetric.
"""
symmetric(op::AbstractLinearOperator) = op.symmetric
issymmetric(op::AbstractLinearOperator) = op.symmetric

"""
    show(io, op)

Display basic information about a linear operator.
"""
function show(io::IO, op::AbstractLinearOperator)
  s = "Linear operator\n"
  nrow, ncol = size(op)
  s *= @sprintf("  nrow: %s\n", nrow)
  s *= @sprintf("  ncol: %d\n", ncol)
  s *= @sprintf("  eltype: %s\n", eltype(op))
  s *= @sprintf("  symmetric: %s\n", symmetric(op))
  s *= @sprintf("  hermitian: %s\n", hermitian(op))
  s *= @sprintf("  nprod:   %d\n", nprod(op))
  s *= @sprintf("  ntprod:  %d\n", ntprod(op))
  s *= @sprintf("  nctprod: %d\n", nctprod(op))
  s *= "\n"
  print(io, s)
end

"""
    A = Matrix(op)

Materialize an operator as a dense array using `op.ncol` products.
"""
function Base.Matrix(op::AbstractLinearOperator{T}) where {T}
  (m, n) = size(op)
  A = Array{T}(undef, m, n)
  ei = zeros(T, n)
  for i = 1:n
    ei[i] = one(T)
    A[:, i] = op * ei
    ei[i] = zero(T)
  end
  return A
end


text:
function reinsert_gcmarker_pass!(fn::LLVM.Function)
    reinsert_gcmarker!(fn)
    return true
end


text:
using PIC: runtests

runtests()

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
__precompile__()

module POMDPBounds

using POMDPs

export
    # Abstract type
    Bound,

    # Bound functions
    upper_bound,
    lower_bound,
    initial_upper_bound,
    initial_lower_bound

abstract Bound

"""
    upper_bound(bound::Bound, prob::Union{POMDP,MDP}, b)

Provide an upper bound on the value function for belief or state b.
"""
function upper_bound end

"""
    lower_bound{B}(bound::Bound, prob::Union{POMDP,MDP}, b::B)

Provide a lower bound on the value function for belief or state b.
"""
function lower_bound end

initial_upper_bound(bound::Bound, prob::Union{POMDP,MDP}) = nothing
initial_lower_bound(bound::Bound, prob::Union{POMDP,MDP}) = nothing

end # module


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666],
     y=[0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512],
     color=[1.17353, 3.78735, 6.1502, 1.78478, 4.76983, 3.82285, 1.17567, 5.83429, 1.19491, 3.15933],
     Scale.color_log2(minvalue=1, maxvalue=8))


text:
module MongocUtils

using Mongoc,Dates,InteractiveUtils
import Base.setindex!
import Mongoc.BSON

export @BSON,@BSON_setindex
export as_struct

const BSON_PRIMITIVE=Union{Mongoc.BSONObjectId,Number,AbstractString,DateTime,Mongoc.BSON,Type,Dict,Mongoc.BSONCode,Date,Vector{UInt8},Nothing,Enum}
const BSON_VALUE_PRIMITIVE=Union{Mongoc.BSONObjectId,Number,AbstractString,DateTime,Mongoc.BSONCode,Date,Vector{UInt8},Nothing}

include("BSON.jl")
include("as_struct.jl")

end


text:
# Unit tests for api/md_simulator.jl

@testset "api/md_simulator.jl" begin
    struct MockSystem <: AbstractSystem{3} end
    struct MockSimulator <: MolecularDynamicsSimulator end
    struct MockPotential <: InteratomicPotentials.AbstractPotential end

    @test_throws UnimplementedError{MockSimulator} simulate(MockSystem(), MockSimulator(), MockPotential())
end


text:
const ContactImplicitMPC = Main
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
model = s.model
env = s.env
nq = model.dim.q
nu = model.dim.u
nc = model.dim.c

function run_policy(s::Simulation; H_sim::Int = 4000, verbose = verbose,
		offset::AbstractVector = zeros(s.model.dim.q))
	model = s.model
	env = s.env
	ref_traj = deepcopy(get_trajectory(model, env,
	    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_in_place.jld2"),
	    load_type = :joint_traj))
	# time
	H = ref_traj.H
	h = ref_traj.h
	N_sample = 5
	H_mpc = 10
	h_sim = h / N_sample

	# barrier parameter
	κ_mpc = 1.0e-4

	obj = TrackingVelocityObjective(model, env, H_mpc,
		q = [[Diagonal(1.0e-2 * [0.1,3,1,3])   for t = 1:H_mpc-2]; [Diagonal(1.0e-1 * [0.1,3,1,3])   for t = 1:2]],
		v = [[Diagonal(1.0e+1 * [0.1,3,1,3])   for t = 1:2]; [Diagonal(1.0e-3 * [0.1,3,1,3])   for t = 1:H_mpc-2]],
	    u = [Diagonal(1.0e-0 * [3e-3, 1e0]) for t = 1:H_mpc],
	    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
	    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

	p = linearized_mpc_policy(ref_traj, s, obj,
	    H_mpc = H_mpc,
	    N_sample = N_sample,
	    κ_mpc = κ_mpc,
		mode = :configuration,
		ip_type = :mehrotra,
	    n_opts = NewtonOptions(
			r_tol = 3e-4,
			max_iter = 5,
			max_time = ref_traj.h, # HARD REAL TIME
			),
	    mpc_opts = LinearizedMPCOptions(
	        ),
		ip_opts = InteriorPointOptions(
			max_iter = 100,
			verbose = false,
			r_tol = 1.0e-4,
			κ_tol = 1.0e-4,
			diff_sol = true,
			# κ_reg = 1e-3,
			# γ_reg = 1e-1,
			solver = :empty_solver,
			),
	    )

	q1_ref = copy(ref_traj.q[2]) + offset
	q0_ref = copy(ref_traj.q[1]) + offset
	q1_sim = SVector{model.dim.q}(q1_ref)
	q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
	@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

	sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
	    p = p,
	    ip_opts = ContactImplicitMPC.InteriorPointOptions(
	        r_tol = 1.0e-8,
	        κ_init = 1.0e-6,
	        κ_tol = 2.0e-6,
	        diff_sol = true),
	    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

	status = ContactImplicitMPC.simulate!(sim, verbose = false)
	return deepcopy(sim.traj)
end

function collect_runs(s::Simulation; n::Int = 1, H_sim::Int = 4000, verbose::Bool = false)
	trajs = []
	Random.seed!(100)
	offset_max = [+0.50, 0.30, +0.30, +0.00]
	offset_min = [-0.50, 0.00, -0.30, -0.30]
	offset_Δ = offset_max .- offset_min
	for i = 1:n
		verbose && println("sample = $i/$n")
		offset = offset_min .+ offset_Δ .* rand(s.model.dim.q)
		traj = run_policy(s, H_sim = H_sim, offset = offset, verbose = verbose)
		push!(trajs, traj)
	end
	return trajs
end

function visualize_runs!(vis::Visualizer, model::ContactModel, trajs::AbstractVector;
		sample=max(1, Int(floor(trajs[1].H / 100))), h=trajs[1].h*sample,  α=1.0,
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=model_name(model))

	for (i, traj) in enumerate(trajs)
		name_i = Symbol(string(name) * string(i))
		anim = visualize_robot!(vis, model, traj, sample = sample, h = h, α = α, anim = anim, name = name_i)
	end
	return anim
end

H_sim = 5000
trajs = collect_runs(s, n = 100, H_sim = H_sim, verbose = true)
anim = visualize_runs!(vis, s.model, trajs, α=0.1, sample=5)
rep_ref_traj = repeat_ref_traj(ref_traj, Int(ceil(H_sim/N_sample/H)))
anim = visualize_robot!(vis, model, rep_ref_traj, name=:ref, anim=anim, sample=1, α=1.0)
plot_surface!(vis, s.env, xlims=[-2,2], ylims=[-0.2,0.2])

plt = plot(layout=(3,1), legend=false)
plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[1,1], hcat(Vector.(sim.traj.q)...)', color=:blue, linewidth=1.0)
plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
    color=:red, linewidth=3.0)
plot!(plt[2,1], hcat(Vector.([u[1:nu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)

rep_ref_traj = repeat_ref_traj(ref_traj, Int(ceil(H_sim/N_sample/H)))
plot_surface!(vis, env, n=200)
plot_lines!(vis, model, sim.traj.q[1:10:end])
anim = visualize_robot!(vis, model, sim.traj, sample=5)
anim = visualize_robot!(vis, model, rep_ref_traj, name=:ref, anim=anim, sample=1, α=0.3)

plot(hcat([u[[1,2]] for u in Vector.(sim.traj.u)[1:end]]...)')

filename = "quadruped_slow_drop"
MeshCat.convert_frames_to_video(
    "/home/simon/Downloads/$filename.tar",
    "/home/simon/Documents/$filename.mp4", overwrite=true)

convert_video_to_gif(
    "/home/simon/Documents/$filename.mp4",
    "/home/simon/Documents/$filename.gif", overwrite=true)

# const ContactImplicitMPC = Main


text:
"""
`GAPA(α=1.0, β=0.0; kwargs...)``

GAP Adaptive
Same as GAP but with adaptive `α1,α2` set to optimal `αopt=2/(1+sinθ)`
according to the estimate of the Friedrischs angle `θ`.
`β` is averaging factor between `αopt` and `2`: `α1=α2= (1-β)*αopt + β*2.0`.
"""
mutable struct GAPA  <: FOSAlgorithm
    α::Float64
    β::Float64
    direct::Bool
    options
end
GAPA(α=1.0, β=0.0; direct=false, kwargs...) = GAPA(α, β, direct, kwargs)

mutable struct GAPAData{T1,T2} <: FOSSolverData
    α12::Float64
    tmp1::Array{Float64,1}
    tmp2::Array{Float64,1}
    constinit::Base.RefValue{Bool} # If constant term has been calculated (linesearch feature)
    S1::T1
    S2::T2
end
#GAPAData{T1,T2}(α12, tmp1, tmp2, S1::T1, S2::T2) = GAPAData{T1,T2}(α12, tmp1, tmp2, S1, S2)

function init_algorithm!(alg::GAPA, model::AbstractFOSModel)
    S1, S2, n, status_generator = get_sets_and_status(alg, model)
    data = GAPAData(2.0, Array{Float64,1}(undef, n), Array{Float64,1}(undef, n),
            Ref(false), S1, S2)
    return data, status_generator
end

""" normedScalar(x1,x2,y1,y2)
Calculate |<x1-x2,y1-y2>|/(||x1-x2||*||y1-y2||)"""
function normedScalar(x1,x2,y1,y2)
    sum, norm1, norm2 = 0.0, 0.0, 0.0
    #@inbounds @simd
    for i in 1:length(x1)
        d1 = x1[i] - x2[i]
        d2 = y1[i] - y2[i]
        sum   += d1*d2
        norm1 += d1^2
        norm2 += d2^2
    end
    return abs(sum)/sqrt(norm1*norm2)
end

function S1constant!(y, alg::GAPA, data::GAPAData, mem, i)
    α12 = data.α12
    if !data.constinit.x
        y .= 0.0
        #TODO high accuracy?
        prox!(mem, S1, y)
        data.constinit.x = true
    end
    y .= α12.*mem# .+ (1-α1).*0.0
    return
end

function S1!(y, alg::GAPA, data::GAPAData, x, i, status::AbstractStatus, longstep=nothing)
    α12, S1 =  data.α12, data.S1
    prox!(y, S1, x)
    addprojeq(longstep, y, x)
    y .= α12.*y .+ (1-α12).*x
end

function S2!(y, alg::GAPA, data::GAPAData, x, i, status::AbstractStatus, longstep=nothing)
    α12, S2 =  data.α12, data.S2
    prox!(y, S2, x)
    checkstatus(status, y)
    addprojineq(longstep, y, x)
    y .= α12.*y .+ (1-α12).*x
end

function Base.step(alg::GAPA, data::GAPAData, x, i, status::AbstractStatus, longstep=nothing)
    α,β = alg.α,alg.β
    α12,tmp1,tmp2,S1,S2 = data.α12,data.tmp1,data.tmp2,data.S1,data.S2
    # Relaxed projection onto S1
    S1!(tmp1, alg, data, x,    i, status, longstep)
    # prox!(tmp1, S1, x)
    # addprojeq(longstep, tmp1, x)
    # tmp1 .= α12.*tmp1 .+ (1-α12).*x
    # Relaxed projection onto S2
    S2!(tmp2, alg, data, tmp1, i, status, longstep)
    # prox!(tmp2, S2, tmp1)
    # checkstatus(status, tmp2)
    # #println("α12: $α12")
    # addprojineq(longstep, tmp2, tmp1)
    # tmp2 .= α12.*tmp2 .+ (1-α12).*tmp1
    #Calculate θ_F estimate, normedScalar(x1,x2,y1,y2) = |<x1-x2,y1-y2>|/(||x1-x2||*||y1-y2||)
    scl = clamp(normedScalar(tmp2,tmp1,tmp1,x), 0.0, 1.0)
    s = sqrt(1-scl^2) #Efficient way to sin(acos(scl))
    # Update α1, α2
    αopt = 2/(1+s) # α1,α2 = 2/(1+sin(Θ_F))
    data.α12 = (1-β)*αopt + β*2.0
    #Set output
    x .= α.*tmp2 .+ (1-α).*x
    return
end

function getsol(::GAPA, data::GAPAData, x)
    tmp1,tmp2,S1,S2 = data.tmp1,data.tmp2,data.S1,data.S2
    prox!(tmp1, S1, x)
    prox!(tmp2, S2, tmp1)
    return tmp2
end

support_longstep(::GAPA) = true
projections_per_step(::GAPA) = (1,1)

support_linesearch(::GAPA) = Val{:Fast}


text:
function threedBilinear( kernel, ϕ, test, w_g )
#= """
  threedBilinear.m - routine to compute \int{ kernel*ϕ*test }
                     (same as twodBilinear)

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions.jl

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    M = threedBilinear(kernel, ϕ, test, w_g)
  ```

  Arguments:
  - `kernel`: Kernel function in the integral evaluated at the quadrature points
              (dim: nQuadrature, nDOF)
  - `ϕ`: Matrix of element test functions evaluated at the quadrature points
              (dim: nQuadrature, nDOF)
  - `test`: Matrix of test functions evaluated at the quadrature points
              (dim: nQuadrature, nDOF)
  - `wg`: Row vector of quadrature weights

  Output argument:
  - `M`: Numerical integral of kernel*ϕ*test
""" =#
#   [n_quadrature,n_row] = size(test)
#   [n_g1   ,n_col] = size(phi )
#
#   M = zeros(n_row,n_col)
#   for i=1:n_row
#     for j=1:n_col
#        M(i,j) = ( w_g'    .* test(:,i)' ) * ( kernel .* phi(:,j) )
#     end
#   end

  wk = wg.*kernel

  M = test'*(LinearAlgebra.Diagonal(wk)*ϕ)

  return M
end