text:
using OnlineStats
export groupreduce, groupby, aggregate, aggregate_vec, summarize, ApplyColwise

"""
`reduce(f, t::Table; select::Selection)`

Reduce `t` by applying `f` pair-wise on values or structs
selected by `select`.

`f` can be:

1. A function
2. An OnlineStat
3. A tuple of functions and/or OnlineStats
4. A named tuple of functions and/or OnlineStats
5. A named tuple of (selector => function or OnlineStat) pairs

```jldoctest reduce
julia> t = table([0.1, 0.5, 0.75], [0,1,2], names=[:t, :x])
Table with 3 rows, 2 columns:
t     x
───────
0.1   0
0.5   1
0.75  2
```

When `f` is a function, it reduces the selection as usual:

```jldoctest reduce
julia> reduce(+, t, select=:t)
1.35
```

If `select` is omitted, the rows themselves are passed to reduce as tuples.

```jldoctest reduce
julia> reduce((a, b) -> @NT(t=a.t+b.t, x=a.x+b.x), t)
(t = 1.35, x = 3)
```

If `f` is an OnlineStat object from the [OnlineStats](https://github.com/joshday/OnlineStats.jl) package, the statistic is computed on the selection.

```jldoctest reduce
julia> using OnlineStats

julia> reduce(Mean(), t, select=:t)
Mean: n=3 | value=0.45
```

# Reducing with multiple functions

Often one needs many aggregate values from a table. This is when `f` can be passed as a tuple of functions:

```jldoctest reduce
julia> y = reduce((min, max), t, select=:x)
(min = 0, max = 2)

julia> y.max
2

julia> y.min
0
```

Note that the return value of invoking reduce with a tuple of functions
will be a named tuple which has the function names as the keys. In the example, we reduced using `min` and `max` functions to obtain the minimum and maximum values in column `x`.

If you want to give a different name to the fields in the output, use a named tuple as `f` instead:

```jldoctest reduce
julia> y = reduce(@NT(sum=+, prod=*), t, select=:x)
(sum = 3, prod = 0)
```

You can also compute many OnlineStats by passing tuple or named tuple of OnlineStat objects as the reducer.

```jldoctest reduce
julia> y = reduce((Mean(), Variance()), t, select=:t)
(Mean = Mean: n=3 | value=0.45, Variance = Variance: n=3 | value=0.1075)

julia> y.Mean
Mean: n=3 | value=0.45

julia> y.Variance
Variance: n=3 | value=0.1075
```

# Combining reduction and selection

In the above section where we computed many reduced values at once, we have been using the same selection for all reducers, that specified by `select`. It's possible to select different inputs for different reducers by using a named tuple of `slector => function` pairs:

```jldoctest reduce
julia> reduce(@NT(xsum=:x=>+, negtsum=(:t=>-)=>+), t)
(xsum = 3, negtsum = -1.35)

```

See [`Selection`](@ref) for more on what selectors can be specified. Here since each output can select its own input, `select` keyword is unsually unnecessary. If specified, the slections in the reducer tuple will be done over the result of selecting with the `select` argument.

"""
function reduce(f, t::Dataset; select=valuenames(t))
    fs, input, T = init_inputs(f, rows(t, select), reduced_type, false)
    acc = init_first(fs, input[1])
    _reduce(fs, input, acc, 2)
end

function reduce(f, v0, t::Dataset; select=valuenames(t))
    fs, input, T = init_inputs(f, rows(t, select), reduced_type, false)
    _reduce(fs, input, v0, 1)
end

@deprecate reduce(f, t::Dataset, v0; select=valuenames(t)) reduce(f, v0, t::Dataset; select=select)

function _reduce(fs, input, acc, start)
    @inbounds @simd for i=start:length(input)
        acc = _apply(fs, acc, input[i])
    end
    acc
end

## groupreduce

addname(v, name) = v
addname(v::Tup, name::Type{<:NamedTuple}) = v
addname(v, name::Type{<:NamedTuple}) = name(v)

struct GroupReduce{F, S, T, P, N}
    f::F
    key::S
    data::T
    perm::P
    name::N
    n::Int

    GroupReduce(f::F, key::S, data::T, perm::P; name::N = nothing) where{F, S, T, P, N} =
        new{F, S, T, P, N}(f, key, data, perm, name, length(key))
end

Base.iteratorsize(::Type{<:GroupReduce}) = Base.SizeUnknown()

Base.start(iter::GroupReduce) = 1

function Base.next(iter::GroupReduce, i1)
    f, key, data, perm, n, name = iter.f, iter.key, iter.data, iter.perm, iter.n, iter.name
    val = init_first(f, data[perm[i1]])
    i = i1+1
    while i <= n && roweq(key, perm[i], perm[i1])
        val = _apply(f, val, data[perm[i]])
        i += 1
    end
    (key[perm[i1]] => addname(val, name)), i
end

Base.done(iter::GroupReduce, state) = state > iter.n

"""
`groupreduce(f, t[, by::Selection]; select::Selection)`

Group rows by `by`, and apply `f` to reduce each group. `f` can be a function, OnlineStat or a struct of these as described in [`reduce`](@ref). Recommended: see documentation for [`reduce`](@ref) first. The result of reducing each group is put in a table keyed by unique `by` values, the names of the output columns are the same as the names of the fields of the reduced tuples.

# Examples

```jldoctest groupreduce
julia> t=table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6],
               names=[:x,:y,:z]);

julia> groupreduce(+, t, :x, select=:z)
Table with 2 rows, 2 columns:
x  +
─────
1  6
2  15

julia> groupreduce(+, t, (:x, :y), select=:z)
Table with 4 rows, 3 columns:
x  y  +
────────
1  1  3
1  2  3
2  1  11
2  2  4

julia> groupreduce((+, min, max), t, (:x, :y), select=:z)
Table with 4 rows, 5 columns:
x  y  +   min  max
──────────────────
1  1  3   1    2
1  2  3   3    3
2  1  11  5    6
2  2  4   4    4
```

If `f` is a single function or a tuple of functions, the output columns will be named the same as the functions themselves. To change the name, pass a named tuple:

```jldoctest groupreduce
julia> groupreduce(@NT(zsum=+, zmin=min, zmax=max), t, (:x, :y), select=:z)
Table with 4 rows, 5 columns:
x  y  zsum  zmin  zmax
──────────────────────
1  1  3     1     2
1  2  3     3     3
2  1  11    5     6
2  2  4     4     4
```

Finally, it's possible to select different inputs for different reducers by using a named tuple of `slector => function` pairs:

```jldoctest groupreduce
julia> groupreduce(@NT(xsum=:x=>+, negysum=(:y=>-)=>+), t, :x)
Table with 2 rows, 3 columns:
x  xsum  negysum
────────────────
1  3     -4
2  6     -4

```

"""
function groupreduce(f, t::Dataset, by=pkeynames(t);
                     select = t isa AbstractIndexedTable ? Not(by) : valuenames(t),
                     cache=false)

    if f isa ApplyColwise
        if !(f.functions isa Union{Function, Type})
            error("Only functions are supported in ApplyColwise for groupreduce")
        end
        return groupby(grp->colwise_group_fast(f.functions, grp), t, by; select=select)
    end

    isa(f, Pair) && (f = (f,))

    data = rows(t, select)

    by = lowerselection(t, by)

    if !isa(by, Tuple)
        by=(by,)
    end
    key  = rows(t, by)
    perm = sortpermby(t, by, cache=cache)

    fs, input, T = init_inputs(f, data, reduced_type, false)

    name = isa(t, NextTable) ? namedtuple(nicename(f)) : nothing
    iter = GroupReduce(fs, key, input, perm, name=name)
    convert(collectiontype(t), collect_columns(iter),
            presorted=true, copy=false)
end

colwise_group_fast(f, grp::Union{Columns, Dataset}) = map(c->reduce(f, c), columns(grp))
colwise_group_fast(f, grp::AbstractVector) = reduce(f, grp)

## GroupBy

_apply_with_key(f::Tup, data::Tup, process_data) = _apply(f, map(process_data, data))
_apply_with_key(f::Tup, data, process_data) = _apply_with_key(f, columns(data), process_data)
_apply_with_key(f, data, process_data) = _apply(f, process_data(data))

_apply_with_key(f::Tup, key, data::Tup, process_data) = _apply(f, map(t->key, data), map(process_data, data))
_apply_with_key(f::Tup, key, data, process_data) = _apply_with_key(f, key, columns(data), process_data)
_apply_with_key(f, key, data, process_data) = _apply(f, key, process_data(data))

struct GroupBy
    f
    key
    data
    perm
    usekey::Bool
    name
    n::Int

    GroupBy(f, key, data, perm; usekey = false, name = nothing) =
        new(f, key, data, perm, usekey, name, length(key))
end

Base.iteratorsize(::Type{<:GroupBy}) = Base.SizeUnknown()

Base.start(::GroupBy) = 1

function Base.next(iter::GroupBy, i1)
    f, key, data, perm, usekey, n, name = iter.f, iter.key, iter.data, iter.perm, iter.usekey, iter.n, iter.name
    i = i1+1
    while i <= n && roweq(key, perm[i], perm[i1])
        i += 1
    end
    process_data = t -> view(t, perm[i1:(i-1)])
    val = usekey ? _apply_with_key(f, key[perm[i1]], data, process_data) :
                   _apply_with_key(f, data, process_data)
    (key[perm[i1]] => addname(val, name)), i
end

Base.done(iter, state) = state > iter.n

collectiontype(::Type{<:NDSparse}) = NDSparse
collectiontype(::Type{<:NextTable}) = NextTable
collectiontype(t::Dataset) = collectiontype(typeof(t))

"""
`groupby(f, t[, by::Selection]; select::Selection, flatten)`

Group rows by `by`, and apply `f` to each group. `f` can be a function or a tuple of functions. The result of `f` on each group is put in a table keyed by unique `by` values. `flatten` will flatten the result and can be used when `f` returns a vector instead of a single scalar value.

# Examples

```jldoctest groupby
julia> t=table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6],
               names=[:x,:y,:z]);

julia> groupby(mean, t, :x, select=:z)
Table with 2 rows, 2 columns:
x  mean
───────
1  2.0
2  5.0

julia> groupby(identity, t, (:x, :y), select=:z)
Table with 4 rows, 3 columns:
x  y  identity
──────────────
1  1  [1, 2]
1  2  [3]
2  1  [5, 6]
2  2  [4]

julia> groupby(mean, t, (:x, :y), select=:z)
Table with 4 rows, 3 columns:
x  y  mean
──────────
1  1  1.5
1  2  3.0
2  1  5.5
2  2  4.0
```

multiple aggregates can be computed by passing a tuple of functions:

```jldoctest groupby
julia> groupby((mean, std, var), t, :y, select=:z)
Table with 2 rows, 4 columns:
y  mean  std       var
──────────────────────────
1  3.5   2.38048   5.66667
2  3.5   0.707107  0.5

julia> groupby(@NT(q25=z->quantile(z, 0.25), q50=median,
                   q75=z->quantile(z, 0.75)), t, :y, select=:z)
Table with 2 rows, 4 columns:
y  q25   q50  q75
──────────────────
1  1.75  3.5  5.25
2  3.25  3.5  3.75
```

Finally, it's possible to select different inputs for different functions by using a named tuple of `slector => function` pairs:

```jldoctest groupby
julia> groupby(@NT(xmean=:z=>mean, ystd=(:y=>-)=>std), t, :x)
Table with 2 rows, 3 columns:
x  xmean  ystd
─────────────────
1  2.0    0.57735
2  5.0    0.57735
```

By default, the result of groupby when `f` returns a vector or iterator of values will not be expanded. Pass the `flatten` option as `true` to flatten the grouped column:

```jldoctest
julia> t = table([1,1,2,2], [3,4,5,6], names=[:x,:y])

julia> groupby((:normy => x->Iterators.repeated(mean(x), length(x)),),
                t, :x, select=:y, flatten=true)
Table with 4 rows, 2 columns:
x  normy
────────
1  3.5
1  3.5
2  5.5
2  5.5
```

The keyword option `usekey = true` allows to use information from the indexing column. `f` will need to accept two
arguments, the first being the key (as a `Tuple` or `NamedTuple`) the second the data (as `Columns`).

```jldoctest
julia> t = table([1,1,2,2], [3,4,5,6], names=[:x,:y])

julia> groupby((:x_plus_mean_y => (key, d) -> key.x + mean(d),),
                              t, :x, select=:y, usekey = true)
Table with 2 rows, 2 columns:
x  x_plus_mean_y
────────────────
1  4.5
2  7.5
```

"""
function groupby end

function groupby(f, t::Dataset, by=pkeynames(t);
            select = t isa AbstractIndexedTable ? Not(by) : valuenames(t),
            flatten=false, usekey = false)

    isa(f, Pair) && (f = (f,))
    data = rows(t, select)
    f = init_func(f, data)
    by = lowerselection(t, by)
    if !(by isa Tuple)
        by = (by,)
    end

    key = by == () ? fill((), length(t)) : rows(t, by)

    fs, input, S = init_inputs(f, data, reduced_type, true)

    if by == ()
        res = usekey ? _apply_with_key(fs, (), input, identity) : _apply_with_key(fs, input, identity)
        res_tup = addname(res, namedtuple(nicename(f)))
        return flatten ? res_tup[end] : res_tup
    end

    perm = sortpermby(t, by)
    # Note: we're not using S here, we'll let _groupby figure it out
    name = isa(t, NextTable) ? namedtuple(nicename(f)) : nothing
    iter = GroupBy(fs, key, input, perm, usekey = usekey, name = name)

    t = convert(collectiontype(t), collect_columns(iter), presorted=true, copy=false)
    t isa NextTable && flatten ?
        IndexedTables.flatten(t, length(columns(t))) : t
end

struct ApplyColwise{T}
    functions::T
    names::Vector{Symbol}
    stack::Bool
    variable::Symbol
end

ApplyColwise(f; stack = false, variable = :variable) = ApplyColwise(f, [Symbol(f)], stack, variable)
ApplyColwise(t::Tuple; stack = false, variable = :variable) = ApplyColwise(t, [map(Symbol,t)...], stack, variable)
ApplyColwise(t::NamedTuple; stack = false, variable = :variable) = ApplyColwise(Tuple(values(t)), keys(t), stack, variable)

init_func(f, t) = f
init_func(ac::ApplyColwise{<:Tuple}, t::AbstractVector) =
    Tuple(Symbol(n) => f for (f, n) in zip(ac.functions, ac.names))
init_func(ac::ApplyColwise{<:Tuple}, t::Columns) =
    ac.stack ? dd -> Columns(colnames(t), ([f(x) for x in columns(dd)] for f in ac.functions)...; names = vcat(ac.variable, ac.names)) :
        Tuple(Symbol(s, :_, n) => s => f for s in colnames(t), (f, n) in zip(ac.functions, ac.names))
init_func(ac::ApplyColwise, t::Columns) =
    ac.stack ? dd -> Columns(colnames(t), [ac.functions(x) for x in columns(dd)]; names = vcat(ac.variable, ac.names)) :
        Tuple(s => s => ac.functions for s in colnames(t))
init_func(ac::ApplyColwise, t::AbstractVector) = ac.functions

"""
`summarize(f, t, by = pkeynames(t); select = excludecols(t, by), stack = false, variable = :variable)`

Apply summary functions column-wise to a table. Return a `NamedTuple` in the non-grouped case
and a table in the grouped case. Use `stack=true` to stack results of the same summary function for different columns.

# Examples

```jldoctest colwise
julia> t = table([1, 2, 3], [1, 1, 1], names = [:x, :y]);

julia> summarize((mean, std), t)
(x_mean = 2.0, y_mean = 1.0, x_std = 1.0, y_std = 0.0)

julia> s = table(["a","a","b","b"], [1,3,5,7], [2,2,2,2], names = [:x, :y, :z], pkey = :x);

julia> summarize(mean, s)
Table with 2 rows, 3 columns:
x    y    z
─────────────
"a"  2.0  2.0
"b"  6.0  2.0

julia> summarize(mean, s, stack = true)
Table with 4 rows, 3 columns:
x    variable  mean
───────────────────
"a"  :y        2.0
"a"  :z        2.0
"b"  :y        6.0
"b"  :z        2.0
```

Use a `NamedTuple` to have different names for the summary functions:

```jldoctest colwise
julia> summarize(@NT(m = mean, s = std), t)
(x_m = 2.0, y_m = 1.0, x_s = 1.0, y_s = 0.0)
```

Use `select` to only summarize some columns:

```jldoctest colwise
julia> summarize(@NT(m = mean, s = std), t, select = :x)
(m = 2.0, s = 1.0)
```

"""
function summarize(f, t, by = pkeynames(t); select = t isa AbstractIndexedTable ? excludecols(t, by) : valuenames(t), stack = false, variable = :variable)
    flatten = stack && !(select isa Union{Int, Symbol})
    s = groupby(ApplyColwise(f, stack = stack, variable = variable), t, by, select = select, flatten = flatten)
    s isa Columns ? table(s, copy = false, presorted = true) : s
end

Base.@deprecate aggregate(f, t;
                          by=pkeynames(t),
                          with=valuenames(t)) groupreduce(f, t, by; select=with)

Base.@deprecate aggregate_vec(
    fs::Function, x;
    names=nothing,
    by=pkeynames(x),
    with=valuenames(x)) groupby(names === nothing ? fs : (names => fs,), x; select=with)

Base.@deprecate aggregate_vec(
    fs::AbstractVector, x;
    names=nothing,
    by=pkeynames(x),
    with=valuenames(x)) groupby(names === nothing ? (fs...) : (map(=>, names, fs)...,), x; select=with)

Base.@deprecate aggregate_vec(t; funs...) groupby(namedtuple(first.(funs)...)(last.(funs)...), t)

"""
`convertdim(x::NDSparse, d::DimName, xlate; agg::Function, vecagg::Function, name)`

Apply function or dictionary `xlate` to each index in the specified dimension.
If the mapping is many-to-one, `agg` or `vecagg` is used to aggregate the results.
If `agg` is passed, it is used as a 2-argument reduction function over the data.
If `vecagg` is passed, it is used as a vector-to-scalar function to aggregate
the data.
`name` optionally specifies a new name for the translated dimension.
"""
function convertdim(x::NDSparse, d::DimName, xlat; agg=nothing, vecagg=nothing, name=nothing, select=valuenames(x))
    ks = setcol(pkeys(x), d, d=>xlat)
    if name !== nothing
        ks = renamecol(ks, d, name)
    end

    if vecagg !== nothing
        y = convert(NDSparse, ks, rows(x, select))
        return groupby(vecagg, y)
    end

    if agg !== nothing
        return convert(NDSparse, ks, rows(x, select), agg=agg)
    end
    convert(NDSparse, ks, rows(x, select))
end

convertdim(x::NDSparse, d::Int, xlat::Dict; agg=nothing, vecagg=nothing, name=nothing, select=valuenames(x)) = convertdim(x, d, i->xlat[i], agg=agg, vecagg=vecagg, name=name, select=select)

convertdim(x::NDSparse, d::Int, xlat, agg) = convertdim(x, d, xlat, agg=agg)

sum(x::NDSparse) = sum(x.data)

"""
`reducedim(f, x::NDSparse, dims)`

Drop `dims` dimension(s) and aggregate with `f`.

```jldoctest
julia> x = ndsparse(@NT(x=[1,1,1,2,2,2],
                        y=[1,2,2,1,2,2],
                        z=[1,1,2,1,1,2]), [1,2,3,4,5,6])
3-d NDSparse with 6 values (Int64):
x  y  z │
────────┼──
1  1  1 │ 1
1  2  1 │ 2
1  2  2 │ 3
2  1  1 │ 4
2  2  1 │ 5
2  2  2 │ 6

julia> reducedim(+, x, 1)
2-d NDSparse with 3 values (Int64):
y  z │
─────┼──
1  1 │ 5
2  1 │ 7
2  2 │ 9

julia> reducedim(+, x, (1,3))
1-d NDSparse with 2 values (Int64):
y │
──┼───
1 │ 5
2 │ 16

```
"""
function reducedim(f, x::NDSparse, dims)
    keep = setdiff([1:ndims(x);], map(d->fieldindex(x.index.columns,d), dims))
    if isempty(keep)
        throw(ArgumentError("to remove all dimensions, use `reduce(f, A)`"))
    end
    groupreduce(f, x, (keep...))
end

reducedim(f, x::NDSparse, dims::Symbol) = reducedim(f, x, [dims])

"""
`reducedim_vec(f::Function, arr::NDSparse, dims)`

Like `reducedim`, except uses a function mapping a vector of values to a scalar instead
of a 2-argument scalar function.
"""
function reducedim_vec(f, x::NDSparse, dims; with=valuenames(x))
    keep = setdiff([1:ndims(x);], map(d->fieldindex(x.index.columns,d), dims))
    if isempty(keep)
        throw(ArgumentError("to remove all dimensions, use `reduce(f, A)`"))
    end
    idxs, d = collect_columns(GroupBy(f, keys(x, (keep...)), rows(x, with), sortpermby(x, (keep...)))).columns
    NDSparse(idxs, d, presorted=true, copy=false)
end

reducedim_vec(f, x::NDSparse, dims::Symbol) = reducedim_vec(f, x, [dims])


text:
__precompile__(false)

module Game

## INCLUDES
include(joinpath(@__DIR__,"libs.jl"))

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

####################
# JULIA
####################
using InteractiveUtils #versioninfo
using Dates
using Distributed
using SharedArrays

####################
# PACKAGES
####################
using DataStructures #SortedDict
using Images
using ImageMagick
using JLD2
using StaticArrays
using ModernGL
using GLFW

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

using CodeManager
using RessourceManager
#using ThreadManager
using TimeManager
using FileManager
using LogManager
using MathManager
using WindowManager
using GraphicsManager
using CameraManager

################################################################################
# PATHS
RessourceManager.setPath(:LOGS,"logs")
RessourceManager.setPath(:ASSETS,"assets")
RessourceManager.setPath(:SCRIPTS,"assets/scripts")
RessourceManager.setPath(:MODELS,"assets/models")
RessourceManager.setPath(:TEXTURES,"assets/textures")
RessourceManager.setPath(:SHADERS,"assets/shaders")

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

const Window = WindowManager.Window

WINDOW_FOCUS = true

function OnFocus(window::Window,focus::Number)
  global WINDOW_FOCUS = focus > 0 ? true : false
end

function createWindow()
  println("Create Window...")
  #GraphicsManager.cleanUp()

  window=WindowManager.Window()

  WindowManager.open(window,[
    :OPENGL_DEBUG_CONTEXT=>true,
    :CONTEXT_VERSION_MAJOR=>4,
    :CONTEXT_VERSION_MINOR=>6,
    #:SAMPLES=>4,
  ])

  WindowManager.addListener(window, :OnWindowFocus, OnFocus)
  WindowManager.addListener(window, :OnMousePos, OnCursorPos)
  WindowManager.addListener(window, :OnKey, OnKey)
  WindowManager.addListener(window, :OnMouseKey, OnMouseKey)
  WindowManager.addListener(window, :OnWindowClose, OnWindowClose)

  GraphicsManager.init()

  GraphicsManager.glDebug(true) # set debugging

  glinfo = GraphicsManager.createcontextinfo()

  println("OpenGL $(stringColor(glinfo[:gl_version];color=:red))")
  println("GLSL $(stringColor(glinfo[:glsl_version];color=:red))")
  println("Vendor $(stringColor(glinfo[:gl_vendor];color=:red))")
  println("Renderer $(stringColor(glinfo[:gl_renderer];color=:red))")
  file=RessourceManager.getPath(:LOGS,"gl_exentsioninfo.txt")
  open(file, "w") do f
    for name in glinfo[:gl_extensions]
      write(f,string(name,"\n"))
    end
  end
  info("$(length(glinfo[:gl_extensions])) exensions saved in $file")
  println("---------------------------------------------------------------------")
  sleep(0)

  window
end

SLEEP=0 #1f0/200
TITLE = "My Game"
STARTTIME = time()
PREVTIME = STARTTIME
FRAMES = 0
MAX_FRAMES = 0
FPS = 0
MAX_FPS = 0
ITERATION = 0
REF_PREVTIME = Ref(0.0)

"""
TODO
"""
function showFrames(this::Window)
  global FRAMES, ITERATION, MAX_FRAMES, REF_PREVTIME, BLOCK_COUNT
  currenttime = time() #GetTimer("FRAME_TIMER")

  ITERATION +=1
  if !OnTime(1.0, REF_PREVTIME; time=currenttime) FRAMES += 1; return end

  #FPS = FRAMES/(currenttime - PREVTIME)
  #PREVTIME = currenttime
  #if MAX_FPS < FPS MAX_FPS = FPS end
  #if FPS > 15 COUNT += 1 end
  #fpms = FPS > 0 ? (1000.0 / FPS) : 0
  #max_fmps = MAX_FPS > 0 ? (1000.0 / MAX_FPS) : 0
  #norm_fps = FPS/MAX_FPS

  if MAX_FRAMES < FRAMES MAX_FRAMES = FRAMES end
  fps = FRAMES
  max_fps = MAX_FRAMES
  fpms = FRAMES > 0 ? (1000.0 / FRAMES) : 0
  max_fmps = MAX_FRAMES > 0 ? (1000.0 / MAX_FRAMES) : 0
  norm_fps = FRAMES / MAX_FRAMES

  WindowManager.title(this,"$(TITLE) - FPS $(round(fps; digits=2))[$(round(max_fps; digits=2))] | FMPS $(round(fpms; digits=2))[$(round(max_fmps; digits=2))] - IT $ITERATION")
  FRAMES = 0
end

function OnWindowClose(window::Window)
  GraphicsManager.freeMemory()
end

function cleanUp()
  info("cleanUp")
  WindowManager.cleanUp()
  CodeManager.safe_clean!(@__MODULE__, :SCRIPT)
end

const EMPTY_FUNCTION = (args...) -> nothing

#script_call(name::Symbol, args...) = safe_eval(:(if isdefined(SCRIPT, Symbol($(string(name)))) SCRIPT.$name($(args...)); end))
function script_function(name::Symbol)
  result = safe_eval(:(isdefined(SCRIPT, Symbol($(string(name)))) ? SCRIPT.$name : nothing))
  if !result[1] || result[2] == nothing result = EMPTY_FUNCTION
  else result = result[2]
  end
  result
end

#script_function(name::Symbol) = eval(Meta.parse("SCRIPT.$name"))
const script_call = Base.invokelatest

script_OnUpdate = EMPTY_FUNCTION
script_OnRender = EMPTY_FUNCTION
OnUpdate = () -> script_call(script_OnUpdate)
OnRender = () -> script_call(script_OnRender)

function script_init(window::Window)
  println("Include Script...")
  include_module(RessourceManager.getPath(:SCRIPTS,"scenes/mcWorld.jl"))
  args = Dict{Symbol,Any}()
  args[:WINDOW] = window.ref

  println("Init Script...")
  global script_OnUpdate = script_function(:OnUpdate)
  global script_OnRender = script_function(:OnRender)

  script_call(script_function(:main), args)
  script_call(script_function(:OnInit))

  #SCRIPT.main(args)
  #SCRIPT.OnInit()
  #global script_OnUpdate = SCRIPT.OnUpdate
  #global script_OnRender = SCRIPT.OnRender
end

function OnUpdatedKeys(window::Window)
  keyValue, keyPressed = getKey()
  if keyPressed
    if keyValue == 294 #F5 #88=x
      reload(window)
      return false
    end
  end
end

function reload(window::Window)
  println("Reload Script...")
  #script_call(script_function(:OnDestroy))
  #GraphicsManager.cleanUp() #remove all
  #CodeManager.safe_clean!(@__MODULE__, :SCRIPT)
  #safe_invoke(:(SCRIPT.free!))

  #waitTime = programTime() + 3
  #print("Wait...")
  #while programTime() <= waitTime
  #  GLFW.PollEvents()
  #  #glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT)
  #  GLFW.SwapBuffers(window)
  #  if OnTime(1.0, REF_PREVTIME) print(".") end
  #end #wait
  #println()

  script_init(window)
  script_call(script_function(:OnReload))
end

function run2()
  println("---------------------------------------------------------------------")
  println("Start Program @ ", Dates.date())
  InteractiveUtils.versioninfo()

  window = createWindow()
  script_init(window)

  while WindowManager.isOpen(window)
    CameraManager.resetKeys()
    WindowManager.OnUpdateEvents() # Poll for and process events
    UpdateTimers()
    showFrames(window)
    OnUpdatedKeys(window)
  end
end

function run()
  println("---------------------------------------------------------------------")
  println("Start Program @ ", Dates.time())
  InteractiveUtils.versioninfo()

  window = createWindow()
  script_init(window)

  WindowManager.loop(window, (this)->begin
    UpdateTimers()
    showFrames(this)
    OnUpdatedKeys(this)
    OnUpdate()
    OnRender()
    # Swap front and back buffers
    #WindowManager.swap(window)
    #if SLEEP>0 Libc.systemsleep(SLEEP) end
    #sleep(SLEEP)
    CameraManager.resetKeys()
  end)

  # while WindowManager.isOpen(window)
  #   CameraManager.resetKeys()
  #   WindowManager.OnUpdateEvents() # Poll for and process events
  #   UpdateTimers()
  #   showFrames(window)
  #   if !WINDOW_FOCUS
  #     sleep(0.1)
  #   else
  #     OnUpdatedKeys(window)
  #     OnUpdate()
  #     OnRender()
  #     # Swap front and back buffers
  #     WindowManager.swap(window)
  #     #if SLEEP>0 Libc.systemsleep(SLEEP) end
  #     #sleep(SLEEP)
  #   end
  # end

  #GLFW.DestroyWindow(window) #might bug sometimes
  cleanUp()
end

end #Game

ext: jl

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
include("../imports.jl")

@testset ExtendedTestSet "`StopOnNanLoss`" begin
    learner = testlearner(coeff = NaN, Recorder(), StopOnNaNLoss())
    # Epoch will be cancelled
    @test_throws CancelFittingException fit!(learner, 1)
    @test learner.cbstate.history[TrainingPhase()].epochs == 0
end


text:
StatsBase.@weights PriorityWeights

const DirectlyWeighted = Union{Bonferroni,BenjaminiHochberg}
#TODO: assumes weights sum to length(pvals)
function adjust(pvals, ws::Union{PriorityWeights,UnitWeights}, method::DirectlyWeighted)
    weighted_pvals = copy(pvals)
    @inbounds for i = 1:length(pvals)
        weighted_pvals[i] = pvals[i] == 0.0 ? 0.0 : min(1.0, pvals[i] / ws[i])
    end
    adjust(weighted_pvals, method)
end


text:
# Day Count (adapted from Ito.jl and InterestRates.jl)
using Base.Dates

abstract DayCount

type Actual360 <:DayCount ; end
type Actual365 <: DayCount ; end

abstract Thirty360 <:DayCount

type BondThirty360 <: Thirty360; end
type EuroBondThirty360 <: Thirty360; end
type ItalianThirty360 <: Thirty360; end

typealias USAThirty360 BondThirty360
typealias EuroThirty360 EuroBondThirty360

abstract ActualActual <: DayCount

type ISMAActualActual <: ActualActual; end
type ISDAActualActual <: ActualActual; end
type AFBActualActual <: ActualActual; end

typealias ActualActualBond ISMAActualActual

type SimpleDayCount <: DayCount end

# Day Counting
# default day count method
function day_count(c::BondThirty360, d_start::Date, d_end::Date)
  dd1 = day(d_start)
  dd2 = day(d_end)

  mm1 = month(d_start)
  mm2 = month(d_end)

  yy1 = year(d_start)
  yy2 = year(d_end)

  if dd2 == 31 && dd1 < 30
    dd2 = 1
    mm2 += 1
  end

  return 360.0 * (yy2 - yy1) + 30.0 * (mm2 - mm1 - 1) + max(0, 30 - dd1) + min(30, dd2)
end

function day_count(c::EuroBondThirty360, d_start::Date, d_end::Date)
  dd1 = day(d_start)
  dd2 = day(d_end)

  mm1 = month(d_start)
  mm2 = month(d_end)

  yy1 = year(d_start)
  yy2 = year(d_end)

  return 360.0 * (yy2 - yy1) + 30.0 * (mm2 - mm1 - 1) + max(0, 30 - dd1) + min(30, dd2)
end

day_count(c::DayCount, d_start::Date, d_end::Date) = Int(d_end - d_start)

# days per year
days_per_year(::Union{Actual360, Thirty360}) = 360.0
days_per_year(::Actual365) = 365.0

# year fractions
# default
year_fraction(c::SimpleDayCount, d_start::Date, d_end::Date) = year_fraction(c, d_start, d_end, Date(), Date())

year_fraction(c::DayCount, d_start::Date, d_end::Date) = day_count(c, d_start, d_end) / days_per_year(c)

# add'l methods
# year_fraction(c::Union{Actual360, Thirty360, Actual365}, d_start::Date, d_end::Date) = year_fraction(c, d_start, d_end, Date(), Date())
year_fraction(c::Union{Actual360, Thirty360, Actual365}, d_start::Date, d_end::Date, ::Date, ::Date) = year_fraction(c, d_start, d_end)

function year_fraction(::SimpleDayCount, d_start::Date, d_end::Date, ::Date, ::Date)
  dm_start = Dates.Day(d_start)
  dm_end = Dates.Day(d_end)

  if dm_start == dm_end || (dm_start > dm_end && Dates.lastdayofmonth(d_end) == d_end) || (dm_start < dm_end && Dates.lastdayofmonth(d_start) == d_start)
    return Int(Dates.Year(d_end) - Dates.Year(d_start)) +  (Int(Dates.Month(d_end)) - Int(Dates.Month(d_start))) / 12.0
  else
    return year_fraction(BondThirty360(), d_start, d_end)
  end
end

function year_fraction(dc::ISDAActualActual, d1::Date, d2::Date, ::Date = Date(), ::Date = Date())
  if d1 == d2
    return 0.0
  end

  if d1 > d2
    return -year_fraction(dc, d2, d1, Date(), Date())
  end

  y1 = year(d1)
  y2 = year(d2)

  dib1 = daysinyear(d1)
  dib2 = daysinyear(d2)
  # println(y1)
  # println(y2)
  # println(dib1)
  # println(dib2)

  sum = y2 - y1 - 1

  sum += day_count(dc, d1, Date(y1+1, 1, 1)) / dib1
  # println(d2)
  # println(day_count(dc, Date(y2, 1, 1), d2))
  sum += day_count(dc, Date(y2, 1, 1), d2) / dib2
  return sum
end

function year_fraction(dc::ISMAActualActual, d1::Date, d2::Date, d3::Date = Date(), d4::Date = Date())
  if d1 == d2
    return 0.0
  end

  if d1 > d2
    return -year_fraction(dc, d2, d1, d3, d4)
  end

  ref_period_start = d3 != Date() ? d3 : d1
  ref_period_end = d4 != Date() ? d4 : d2

  months = floor(Int, 0.5 + 12 * Int(ref_period_end - ref_period_start) / 365)

  if months == 0
    ref_period_start = d1
    ref_period_end = d1 + Year(1)
    months = 12
  end

  period = months / 12.0

  if d2 <= ref_period_end
    if d1 >= ref_period_start
      return period * day_count(dc, d1, d2) / day_count(dc, ref_period_start, ref_period_end)
    else
      previous_ref = ref_period_start - Month(months)
      if d2 > ref_period_start
        return year_fraction(dc, d1, ref_period_start, previous_ref, ref_period_end) + year_fraction(dc, d1, d2, previous_ref, ref_period_start)
      else
        return year_fraction(dc, d1, d2, previous_ref, ref_period_start)
      end
    end
  else
    sum = year_fraction(dc, d1, ref_period_end, ref_period_start, ref_period_end)
    i = 0
    new_ref_start = new_ref_end = Date()
    while true
      new_ref_start = ref_period_end + Month(i * months)
      new_ref_end = ref_period_end + Month((i + 1) * months)
      if d2 < new_ref_end
        break
      else
        sum += period
        i += 1
      end
    end
    sum += year_fraction(dc, new_ref_start, d2, new_ref_start, new_ref_end)
    return sum
  end
end


text:
using Test
using SparseIR

@testset "matsubara.jl" begin
    @testset "single pole with stat = $stat, Λ = $Λ" for stat in (fermion, boson),
        Λ in (10, 42, 10_000)

        ε = (@isdefined Float64x2) ? nothing : 1e-15
        wmax = 1.0
        pole = 0.1 * wmax
        β = Λ / wmax
        basis = FiniteTempBasis(stat, β, wmax, ε; sve_result=sve_logistic[Λ])

        stat_shift = (stat == fermion) ? 1 : 0
        spr = SparsePoleRepresentation(basis, [pole])
        gl = to_IR(spr, [1.0])

        func_G(n) = 1 / (im * (2n + stat_shift) * π / β - pole)

        # Compute G(iwn) using unl
        matsu_test = Int[-1, 0, 1, 1e2, 1e4, 1e6, 1e8, 1e10, 1e12]
        prj_w = transpose(basis.uhat(2matsu_test .+ stat_shift))
        Giwn_t = prj_w * gl

        # Compute G(iwn) from analytic expression
        Giwn_ref = func_G.(matsu_test)

        magnitude = maximum(abs, Giwn_ref)
        diff = abs.(Giwn_t - Giwn_ref)
        tol = max(10 * last(basis.s) / first(basis.s), 1e-10)

        # Absolute error
        @test maximum(diff ./ magnitude) < tol

        # Relative error
        @test maximum(abs, diff ./ Giwn_ref) < tol
    end
end


text:
# tests for 1d example
using Test

push!(LOAD_PATH, "/home/felix/Work/Research/Code")
cd("/home/felix/Work/Research/Code/BSE_k_ISDF/experiments")

using Revise # remove after debugging
using BSE_k_ISDF, LinearAlgebra, FFTW

l = 1.5
V_sp = r -> 20 * cos(4π / l * r) + 0.2 * sin(2π / l * r)

V_func = (r_1, r_2) -> 1 /  sqrt((r_1 - r_2)^2 + 0.01)
W_func = (r_1, r_2, l, L) -> 0.0625 * (3.0 + sin(2π / l * r_1)) * (3.0 + cos(4π / l * r_2)) *
    exp(-abs(BSE_k_ISDF.supercell_difference(r_1, r_2, L))^2 / (2 * (4 * l)^2)) * V_func(BSE_k_ISDF.supercell_difference(r_1, r_2, L), 0.)

# fixed parameters
N_r = 128

N_core = 0
N_v = 4
N_c = 10
N_k = 256

N_μ_cc = 22
N_μ_vv = 16
N_μ_vc = 21

sp_prob = BSE_k_ISDF.SPProblem1D(V_sp, l, N_r, N_k);
prob = BSE_k_ISDF.BSEProblem1D(sp_prob, N_core, N_v, N_c, V_func, W_func);
isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc);

ir = 4100
jr = 34
r = prob.prob.r_super[ir]
rp = prob.prob.r_unit[jr]
W_func(r, rp, l, N_k * l)

1 / (N_k * l) * sum(exp(im * (q + G) * r) * prob.w_hat[iG, iGp, iq] * exp(-im * (q + Gp) * rp)
    for (iq, q) in enumerate(2 * pi / l .* (prob.prob.k_bz .- 0.5)),
    (iG, G) in enumerate(2 * pi / l .* (0:127)), (iGp, Gp) in enumerate(2 * pi / l .* (0:127)))

w_realspace = W_func.(prob.prob.r_super, prob.prob.r_unit', l, N_k * l)

D = BSE_k_ISDF.setup_D(prob);
V = BSE_k_ISDF.setup_V(prob, isdf);
W = BSE_k_ISDF.setup_W(prob, isdf);
H = BSE_k_ISDF.setup_H(prob, isdf);

@time V_dense = Matrix(V);
@time V_exact = BSE_k_ISDF.assemble_exact_V(prob)
# @time V_realspace = BSE_k_ISDF.assemble_exact_V_realspace(prob, V_func)
norm(V_dense - V_exact) / norm(V_exact)

@time W_dense = Matrix(W);
@time W_exact = BSE_k_ISDF.assemble_exact_W(prob)
# @time W_realspace = BSE_k_ISDF.assemble_exact_W_realspace(prob, W_func)
norm(W_dense - W_exact) / norm(W_exact)

direction = 1
N_iter = 200
σ = 0.1
g = ω -> 1 / π * σ / (ω^2 + σ^2)
E_min, E_max = 0.0, 30.0
Erange = E_min:0.001:E_max

optical_absorption_lanc = BSE_k_ISDF.lanczos_optical_absorption(prob, isdf, direction, N_iter, g, Erange)

plot!(Erange, optical_absorption_lanc)


text:
# (weighted) sum of L2 norm and L1 norm
# for Group Lasso, this can be used together with src/calculus/slicedSeparableSum
export NormL1plusL2

"""
**L1 norm plus L2 norm**

    NormL1plusL2(λ_1=1, λ_2=1)

With two nonegative scalars λ_1 and λ_2, returns the function
```math
f(x) = λ_1 ∑_{i=1}^{n} |x_i| + λ_2 \\sqrt{x_1^2 + … + x_n^2}.
```
With nonnegative array λ_1 and nonnegative scalar λ_2, returns the function
```math
f(x) =  ∑_{i=1}^{n} {λ_1}_i |x_i| + λ_2 \\sqrt{x_1^2 + … + x_n^2}.
```
"""
struct NormL1plusL2{L1<:NormL1, L2 <: NormL2} <: ProximableFunction
    l1::L1
    l2::L2
end

is_separable(f::NormL1plusL2) = false
is_convex(f::NormL1plusL2) = true
is_positively_homogeneous(f::NormL1plusL2) = true

function NormL1plusL2(lambda1::L=1, lambda2::M=1) where {L <: Union{Real, AbstractArray}, M <: Real}
    NormL1plusL2(NormL1(lambda1), NormL2(lambda2))
end

function (f::NormL1plusL2)(x::AbstractArray{T}) where T <: RealOrComplex
    return f.l1(x) + f.l2(x)
end

function prox!(y::AbstractArray{T}, f::NormL1plusL2, x::AbstractArray{T}, gamma::Real=1) where T <: RealOrComplex
    prox!(y, f.l1, x, gamma)
    prox!(y, f.l2, y, gamma)
    return f(y)
end

fun_name(f::NormL1plusL2) = "L2-norm + L1-norm"
fun_dom(f::NormL1plusL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::NormL1plusL2{L1,L2}) where {L1<:NormL1{<:Real}, L2} = "x ↦ λ_1 ||x||_1 + λ_2||x||_2"
fun_expr(f::NormL1plusL2{L1,L2}) where {L1<:NormL1{<:AbstractArray}, L2} = "x ↦ sum( (λ_1)_i |x_i| ) + λ_2||x||_2"
fun_params(f::NormL1plusL2{L1,L2}) where {L1<:NormL1{<:Real}, L2}  = "λ_1 = $(f.l1.lambda), λ_2 = $(f.l2.lambda)"
fun_params(f::NormL1plusL2{L1,L2}) where {L1<:NormL1{<:AbstractArray}, L2}  = "λ_1 = $(typeof(f.lambda)) of size $(size(f.lambda)), λ_2 = $(f.l2.lambda)"

function prox_naive(f::NormL1plusL2, x::AbstractArray{T}, gamma::Real=1) where T <: RealOrComplex
    y1, v1 = prox_naive(f.l1, x, gamma)
    y2, v2 = prox_naive(f.l2, y1, gamma)
    return y2, f(y2)
end


text:
include("base_utils.jl")
include("comparison_utils.jl")
include("properties_utils.jl")
include("behaviour_utils.jl")


text:
module LowLevelParticleFilters

export KalmanFilter, UnscentedKalmanFilter, ExtendedKalmanFilter, ParticleFilter, AuxiliaryParticleFilter, AdvancedParticleFilter, SigmaFilter, PFstate, index, state, covariance, num_particles, effective_particles, weights, expweights, particles, particletype, smooth, sample_measurement, simulate, loglik, log_likelihood_fun, forward_trajectory, mean_trajectory, mode_trajectory, weigthed_mean, weigthed_cov, update!, predict!, correct!, reset!, metropolis, shouldresample, TupleProduct

export densityplot, debugplot, commandplot, trajectorydensity, dimensiondensity

using StatsBase, Parameters, Lazy, Random, LinearAlgebra, Printf, LoopVectorization
import PDMats # To extend some methods on static arrays
import PDMats: PDMat
using StaticArrays
using Distributions
using RecipesBase
using ForwardDiff

abstract type ResamplingStrategy end
struct ResampleSystematic <: ResamplingStrategy end

abstract type AbstractFilter end

include("PFtypes.jl")
include("kalman.jl")
include("filtering.jl")
include("resample.jl")
include("utils.jl")
include("smoothing.jl")
include("plotting.jl")
include("ekf.jl")

index(f::AbstractFilter) = f.t[]

using Requires
function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots.jl")
    @require MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" include("mcm.jl")
end

end # module


text:
function convert(::Type{NodalData}, S::T) where T <: GphysData
  (T == NodalData) && (return deepcopy(S))
  @assert minimum(S.fs) == maximum(S.fs)
  if T != SeisData
    S = convert(SeisData, S)
  end
  sync!(S, s="last", t="first")

  L = [length(i) for i in S.x]
  L0 = minimum(L)
  @assert L0 == maximum(L)

  # clear everything from buffer TDMS
  reset_tdms()

  ts = S.t[1][1,2]
  data = Array{Float32, 2}(undef, L0, S.n)
  for i in 1:S.n
    data[:, i] .= S.x[i]
  end
  TD = NodalData(data, TDMS.hdr, 1:S.n, ts)
  for f in convertible_fields
    setfield!(TD, f, deepcopy(getfield(S, f)))
  end
  return TD
end
NodalData(S::T) where T <: GphysData = convert(NodalData, S)

function convert(::Type{SeisData}, TD::NodalData)
  S = SeisData(getfield(TD, :n))
  nx = size(TD.data, 1)

  # convertible_fields plus :t are directly copied
  for f in convertible_fields
    setfield!(S, f, deepcopy(getfield(TD, f)))
  end
  setfield!(S, :t, deepcopy(getfield(TD, :t)))

  # :x is set by copying from :data, to prevent GC problems if TD is cleared
  for i in 1:S.n
    S.x[i] = copy(TD.data[:, i])
  end
  return S
end

function convert(::Type{SeisChannel}, D::NodalChannel)
  C = SeisChannel()
  for f in datafields
    setfield!(C, f, deepcopy(getfield(D, f)))
  end
  return C
end

function convert(::Type{NodalChannel}, C::T) where T<:GphysChannel
  (T == NodalChannel) && (return deepcopy(C))
  if T != SeisChannel
    C = convert(SeisChannel, C)
  end

  D = NodalChannel()
  for f in datafields
    setfield!(D, f, deepcopy(getfield(C, f)))
  end
  return D
end
NodalChannel(C::T) where T <: GphysChannel = convert(NodalChannel, C)

push!(TD::NodalData, C::SeisChannel) = push!(TD, convert(NodalChannel, C))


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

function generate_symmetric_sparse_matrix(m)
    A = SparseMatrixCSC{Cdouble, Cint}(sprand(m, m, 0.1))
    A = A*A'
    return A
end

function generate_symmetric_nonzerodiagonal_sparse_matrix(m)
    A = generate_symmetric_sparse_matrix(m)

    # Make digonal nonzero
    for i = 1:m
        if abs(A[i, i]) < 0.5
            A[i, i] = 1.0
        end
    end

    return A
end

function check_symmetry(A)
    println("**** Checking symmetry")
    m = size(A, 1)
    for i = 1:m
        for j = 1:m
            if A[i, j] != A[j, i]
                println("Matrix is asymmetric!")
                println("A[", i, ",", j, "] != A[", j, ",", i, "]")
                break
            end
        end
    end
    println("Done checking")
    flush(STDOUT::IO)
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
using DataFrames

export FeatureExtractor

"""
    FeatureExtractor

Abstract type for features that may be reported about the solution to a PC-TAPF
    or sequential task assignment problem
"""
abstract type FeatureExtractor{T} end

export feature_type
feature_type(extractor::FeatureExtractor{T}) where {T} = T

export
    RunTime,
    MemAllocs,
    GCTime,
    ByteCount,
    SolutionCost,
    PrimaryCost,
    OptimalityGap,
    OptimalFlag,
    FeasibleFlag,
    NumConflicts,
    FinalTimeFeature,
    IterationCount,
    LowLevelIterationCount

export
    TimeOutStatus,
    IterationMaxOutStatus

export
    RobotPaths,
    RobotSeparation

struct RunTime          <: FeatureExtractor{Float64} end
struct MemAllocs        <: FeatureExtractor{Float64} end
struct GCTime           <: FeatureExtractor{Float64} end
struct ByteCount        <: FeatureExtractor{Float64} end
struct SolutionCost     <: FeatureExtractor{Vector{Float64}} end
struct PrimaryCost      <: FeatureExtractor{Float64} end
struct OptimalityGap    <: FeatureExtractor{Float64} end
struct OptimalFlag      <: FeatureExtractor{Bool} end
struct FeasibleFlag     <: FeatureExtractor{Bool} end
struct NumConflicts     <: FeatureExtractor{Int} end
struct FinalTimeFeature <: FeatureExtractor{Int} end
struct IterationCount   <: FeatureExtractor{Int} end
struct LowLevelIterationCount <: FeatureExtractor{Int} end
struct TimeOutStatus    <: FeatureExtractor{Bool} end
struct IterationMaxOutStatus <: FeatureExtractor{Bool} end
struct RobotPaths       <: FeatureExtractor{Vector{Int}} end
struct RobotSeparation  <: FeatureExtractor{Vector{Int}} end

export
    extract_feature,
    load_feature

"""
    extract_feature(solution,extractor,solution,solve_time)

Function must be implemented to extract features from solutions
Args:
- solver (pc_tapf solver)
- extractor <: FeatureExtractor
- solution <: SearchEnv
- timer_results::NamedTuple{t,bytes,gctime,memallocs}
"""
function extract_feature end

extract_feature(solver,feats::Tuple,args...) = map(feat->extract_feature(solver,feat,args...),feats)

extract_feature(solver,::RunTime,       mapf,solution,timer_results) = timer_results.t
extract_feature(solver,::MemAllocs,     mapf,solution,timer_results) = timer_results.memallocs
extract_feature(solver,::GCTime,        mapf,solution,timer_results) = timer_results.gctime
extract_feature(solver,::ByteCount,     mapf,solution,timer_results) = timer_results.bytes
extract_feature(solver,::SolutionCost,  mapf,solution,timer_results) = get_cost(solution)
extract_feature(solver,::PrimaryCost,  args...) = extract_feature(solver,SolutionCost(),args...)[1]
extract_feature(solver,::OptimalityGap, mapf,solution,timer_results) = optimality_gap(solver)
extract_feature(solver,::OptimalFlag,   mapf,solution,timer_results) = optimality_gap(solver) <= 0
extract_feature(solver,::FeasibleFlag,  mapf,solution,timer_results) = best_cost(solver) < typemax(typeof(best_cost(solver)))
extract_feature(solver,::NumConflicts,  mapf,solution,timer_results) = count_conflicts(detect_conflicts(solution))
extract_feature(solver,::FinalTimeFeature,     mapf,solution,timer_results) = maximum(map(length,get_paths(solution)))
extract_feature(solver,::IterationCount, mapf,solution,timer_results) = iterations(solver)
extract_feature(solver,::LowLevelIterationCount, args...) = -1
extract_feature(solver::BiLevelPlanner,::LowLevelIterationCount, args...) = max_iterations(low_level(solver))
extract_feature(solver,::TimeOutStatus, mapf,solution,timer_results) = time_out_status(solver)
extract_feature(solver,::IterationMaxOutStatus, mapf,solution,timer_results) = iteration_max_out_status(solver)
extract_feature(solver,::RobotPaths,    mapf,solution,timer_results) = convert_to_vertex_lists(solution)
function extract_feature(solver,::RobotSeparation, mapf,solution,timer_results)
    separation_distances = Vector{Float64}()
    for (t,nodes) in enumerate(zip(map(p->p.path_nodes,get_paths(solution))...))
        d_min = Inf
        for (j,n1) in enumerate(nodes)
            for (k,n2) in enumerate(nodes)
                if k <= j
                    continue
                end
                d_min = min(get_distance(mapf,get_sp(n1),get_sp(n2)),d_min)
            end
        end
        push!(separation_distances,d_min)
    end
    return separation_distances
end

function load_feature(feat::FeatureExtractor, results)
    feature_key = string(typeof(feat))
    @assert(haskey(results,feature_key))
    return results[feature_key]
end

export compile_results

function compile_results(solver,feats,args...)
    results = Dict{String,Any}(string(typeof(feat))=>extract_feature(solver,feat,args...) for feat in feats)
    for (k,v) in results
        if isa(v,Tuple)
            results[k] = [v...] # because TOML doesn't handle tuples
        end
    end
    return results
end

export profile_solver!

"""
    profile_solver!(solver,mapf)

Profile `solver` on `mapf`:
Returns
- solution, timer_results
"""
function profile_solver!(solver,mapf)
    reset_solver!(solver)
    (solution, cost), t, bytes, gctime, memallocs = @timed solve!(solver,mapf)
    timer_results = (
        t=t,
        bytes=bytes,
        gctime=gctime,
        memallocs=memallocs
        )
    return solution, timer_results
end

export load_problem

"""
    load_problem(loader,filename)

Generic function that allows a problem instance to be retrived/constructed from
`loader` (which may cache information) and `filename::String` (which points to a
problem file).
"""
function load_problem end

export get_results_path, write_results

function get_results_path(loader,solver_config,problem_file)
    problem_name = splitext(splitdir(problem_file)[end])[1]
    outpath = joinpath(solver_config.results_path)
    mkpath(outpath)
    outfile = joinpath(outpath,string(problem_name,".results"))
end
"""
    write_results(loader,solver_config,prob,problem_file,results)

Write results `results` of profiling `solver_config.solver` on problem `prob` of
name `problem_name`. Can be overloaded to dispatch on problem type when applicable.
"""
function write_results(loader,solver_config,prob,problem_file,results)
    # problem_name = splitext(splitdir(problem_file)[end])[1]
    # outpath = joinpath(solver_config.results_path)
    # outfile = joinpath(outpath,string(problem_name,".results"))
    # mkpath(outpath)
    write_results(get_results_path(loader,solver_config,problem_file),results)
end
function write_results(outfile::String,results::Dict)
    open(outfile,"w") do io
        TOML.print(io,results)
    end
end

export run_profiling

"""
    run_profiling(config,loader)

Profile the performance of one or more solvers on a set of problems defined
by files in a directory.
Args:
* `config`: A named tuple or struct with at least the following fields:
    - `problem_dir::String` : path to problem files.
    - `solver_configs` : A vector of objects, each of which has the fields:
        - `solver` : a solver complying with the CRCBS interface.
        - `results_path::String` : path at which to store results.
    - feats : A vector of `Feature` objects, which defines the specific
    features to be extracted and saved.
* `Loader`: Any kind of cache on which `load_problem(loader,problem_file)`
    can be called.
"""
function run_profiling(config,loader)
    for solver_config in config.solver_configs
        mkpath(solver_config.results_path)
    end
    for problem_file in readdir(config.problem_dir;join=true)
        mapf = load_problem(loader,problem_file)
        for solver_config in config.solver_configs
            solver = solver_config.solver
            solution, timer_results = profile_solver!(solver,mapf)
            results_dict = compile_results(
                solver,
                config.feats,
                mapf,
                solution,
                timer_results
                )
            results_dict["problem_file"] = problem_file
            write_results(loader,solver_config,mapf,problem_file,results_dict)
        end
    end
end

################################################################################
############################### Results Analysis ###############################
################################################################################

export init_dataframe

# extract_key_types(feats::Vector{FeatureExtractor}) = [Symbol(typeof(f))=>feature_type(f) for f in feats if feature_type(f) <: Real]
# extract_key_types(dict::Dict) = [Symbol(string(k))=>typeof(v) for (k,v) in dict if isa(v,Union{Real,String})]
extract_key_types(feats::Vector{FeatureExtractor}) = [Symbol(typeof(f))=>feature_type(f) for f in feats]
extract_key_types(dict::Dict) = [Symbol(string(k))=>typeof(v) for (k,v) in dict]

export get_problem_name
get_problem_name(loader,prob_file) = split(prob_file,"/")[end]

export get_config_path
get_config_path(loader,problem_file) = joinpath(problem_file,"config.toml")

export load_results
load_dict(f::String) = Dict(Symbol(k)=>v for (k,v) in TOML.parsefile(f))
load_results(results_file::String) = load_dict(results_file)
function load_results(loader,results_file)
    # results = Dict(Symbol(k)=>v for (k,v) in TOML.parsefile(results_file))
    results = load_results(results_file)
    prob_file = results[:problem_file]
    results[:problem_name] = get_problem_name(loader,prob_file)
    results
end
function load_config(loader,prob_file)
    config_file = get_config_path(loader,prob_file)
    # config = Dict(Symbol(k)=>v for (k,v) in TOML.parsefile(config_file))
    config = load_dict(config_file)
    config[:problem_name] = get_problem_name(loader,prob_file)
    config
end

"""
    init_dataframe(feats::Tuple)

Instantiate an empty dataframe based on the names and types of feature
extractors in `feats`.
"""
function init_dataframe(feats,extras=[:problem_name=>String])
    df = DataFrame()
    for (k,VType) in extract_key_types(feats)
        df[!,k] = VType[]
    end
    for (k,VType) in extras
        df[!,k] = VType[]
    end
    df
end

export construct_results_dataframe
"""
    construct_results_dataframe(loader,solver_config,config_template)

Compile results at `solver_config.results_path` into a `DataFrame` based on the
features stored in `solver_config.feats`
"""
function construct_results_dataframe(loader,solver_config,config_template;
        feats = solver_config.feats,
    )
    results_df = init_dataframe(feats)
    for results_file in readdir(solver_config.results_path;join=true)
        results = load_results(loader,results_file)
        push!(results_df,results;cols=:intersect)
    end
    results_df
end

export construct_config_dataframe
"""
    construct_results_dataframe(loader,solver_config,config_template)

Compile results at `solver_config.results_path` into a `DataFrame` based on the
features stored in `solver_config.feats`
"""
function construct_config_dataframe(loader,problem_path,config_template)
    df = init_dataframe(config_template)
    for problem_file in readdir(problem_path;join=true)
        config = load_config(loader,problem_file)
        push!(df,config;cols=:intersect)
    end
    df
end

# export innerjoin_dataframes
# function innerjoin_dataframes(results_df,config_df)
#     innerjoin(results_df,config_df,on=:problem_name)
# end


text:
using Test, LibCEED, LinearAlgebra, StaticArrays

@testset "LibCEED Development Tests" begin
    @test ceedversion() isa VersionNumber
    @test isrelease() == false

    @testset "CeedVector" begin
        n = 10
        c = Ceed()

        v1 = rand(n)
        v2 = rand(n)
        v3 = rand(n)

        cv1 = CeedVector(c, v1)
        cv2 = CeedVector(c, v2)
        cv3 = CeedVector(c, v3)

        alpha = rand()

        scale!(cv1, alpha)
        v1 .*= alpha
        @test @witharray_read(a = cv1, a == v1)

        pointwisemult!(cv1, cv2, cv3)
        v1 .= v2.*v3
        @test @witharray_read(a = cv1, a == v1)

        axpy!(alpha, cv2, cv1)
        axpy!(alpha, v2, v1)
        @test @witharray_read(a = cv1, a ≈ v1)
    end
end

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


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

include("rich.jl")

# Utils

function withtag(f, io::IO, tag, attrs...)
    print(io, "<$tag")
    for (attr, value) in attrs
        print(io, " ")
        htmlesc(io, attr)
        print(io, "=\"")
        htmlesc(io, value)
        print(io, "\"")
    end
    f === nothing && return print(io, " />")

    print(io, ">")
    f()
    print(io, "</$tag>")
end

tag(io::IO, tag, attrs...) = withtag(nothing, io, tag, attrs...)

const _htmlescape_chars = Dict('<'=>"<",   '>'=>">",
                               '"'=>""", '&'=>"&",
                               # ' '=>" ",
                               )
for ch in "'`!\$\%()=+{}[]"
    _htmlescape_chars[ch] = "&#$(Int(ch));"
end

function htmlesc(io::IO, s::AbstractString)
    # s1 = replace(s, r"&(?!(\w+|\#\d+);)", "&")
    for ch in s
        print(io, get(_htmlescape_chars, ch, ch))
    end
end
function htmlesc(io::IO, s::Symbol)
    htmlesc(io, string(s))
end
function htmlesc(io::IO, xs::Union{AbstractString,Symbol}...)
    for s in xs
        htmlesc(io, s)
    end
end
function htmlesc(s::Union{AbstractString,Symbol})
    sprint(htmlesc, s)
end

# Block elements

function html(io::IO, content::Vector)
    for md in content
        html(io, md)
        println(io)
    end
end

html(io::IO, md::MD) = html(io, md.content)

function html{l}(io::IO, header::Header{l})
    withtag(io, "h$l") do
        htmlinline(io, header.text)
    end
end

function html(io::IO, code::Code)
    withtag(io, :pre) do
        maybe_lang = !isempty(code.language) ? Any[:class=>"language-$(code.language)"] : []
        withtag(io, :code, maybe_lang...) do
            htmlesc(io, code.code)
            # TODO should print newline if this is longer than one line ?
        end
    end
end

function html(io::IO, md::Paragraph)
    withtag(io, :p) do
        htmlinline(io, md.content)
    end
end

function html(io::IO, md::BlockQuote)
    withtag(io, :blockquote) do
        println(io)
        html(io, md.content)
    end
end

function html(io::IO, f::Footnote)
    withtag(io, :div, :class => "footnote", :id => "footnote-$(f.id)") do
        withtag(io, :p, :class => "footnote-title") do
            print(io, f.id)
        end
        html(io, f.text)
    end
end

function html(io::IO, md::Admonition)
    withtag(io, :div, :class => "admonition $(md.category)") do
        withtag(io, :p, :class => "admonition-title") do
            print(io, md.title)
        end
        html(io, md.content)
    end
end

function html(io::IO, md::List)
    maybe_attr = md.ordered > 1 ? Any[:start => string(md.ordered)] : []
    withtag(io, isordered(md) ? :ol : :ul, maybe_attr...) do
        for item in md.items
            println(io)
            withtag(io, :li) do
                html(io, item)
            end
        end
        println(io)
    end
end

function html(io::IO, md::HorizontalRule)
    tag(io, :hr)
end

html(io::IO, x) = tohtml(io, x)

# Inline elements

function htmlinline(io::IO, content::Vector)
    for x in content
        htmlinline(io, x)
    end
end

function htmlinline(io::IO, code::Code)
    withtag(io, :code) do
        htmlesc(io, code.code)
    end
end

function htmlinline(io::IO, md::Union{Symbol,AbstractString})
    htmlesc(io, md)
end

function htmlinline(io::IO, md::Bold)
    withtag(io, :strong) do
        htmlinline(io, md.text)
    end
end

function htmlinline(io::IO, md::Italic)
    withtag(io, :em) do
        htmlinline(io, md.text)
    end
end

function htmlinline(io::IO, md::Image)
    tag(io, :img, :src=>md.url, :alt=>md.alt)
end

function htmlinline(io::IO, f::Footnote)
    withtag(io, :a, :href => "#footnote-$(f.id)", :class => "footnote") do
        print(io, "[", f.id, "]")
    end
end

function htmlinline(io::IO, link::Link)
    withtag(io, :a, :href=>link.url) do
        htmlinline(io, link.text)
    end
end

function htmlinline(io::IO, br::LineBreak)
    tag(io, :br)
end

htmlinline(io::IO, x) = tohtml(io, x)

# API

export html

html(md) = sprint(html, md)

function show(io::IO, ::MIME"text/html", md::MD)
    withtag(io, :div, :class=>"markdown") do
        html(io, md)
    end
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
function OnBoardAntennaRangeModifier(arg0::Vector3D)
    return OnBoardAntennaRangeModifier((Vector3D,), arg0)
end

function get_parameters_drivers(obj::OnBoardAntennaRangeModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::OnBoardAntennaRangeModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
module QJuliaUtils

export gridsize_from_cmdline

gridsize_from_cmdline = Array{Cint, 1}(undef, 4);
rank_order            = Int32(0)
device_id 	      = Int32(0)
rank_offset	      = Int32(0)

function lex_rank_from_coords_t_(comms, fdata)::Cint
	rank = gridsize_from_cmdline[1]

	for i in 1:4
          rank = gridsize_from_cmdline[i] * rank + gridsize_from_cmdline[i]
          println("Communication direction: ", i , ", rank : ", rank , ", grid point ", gridsize_from_cmdline[i] )
	end
	return rank
end #lex_rank_from_coords_t_
lex_rank_from_coords_t_c = @cfunction(lex_rank_from_coords_t_, Cint, (Ptr{Cint}, Ptr{Cvoid}));

function lex_rank_from_coords_x_(comms, fdata)::Cint
	rank = gridsize_from_cmdline[4]
	for i in 4:-1:1
          rank = gridsize_from_cmdline[i] * rank + gridsize_from_cmdline[i]
          println("Communication direction: ", i , ", rank : ", rank , ", grid point ", gridsize_from_cmdline[i] )
	end
	return rank
end #lex_rank_from_coords_x_
lex_rank_from_coords_x_c = @cfunction(lex_rank_from_coords_x_, Cint, (Ptr{Cint}, Ptr{Cvoid}));

function get_rank_order(rorder::String)
  if rorder == "col"
    println("Setting column major rank order t->z->y->x, t is fastest")
    rank_order = Int32(0)
  elseif rorder == "row"
    println("Setting row major rank order x->y->z->t, x is fastest")
    rank_order = Int32(1)
  else
    throw("func 'get_rank_order':: Rank order is not supported.")
  end
end #get_rank_order

@inline function gen_random_spinor!(spin::Vector{Complex{T}}) where T <: AbstractFloat
  for i in 1:length(spin);  spin[i] = Complex{T}(rand(), rand()); end
end #gen_random_spinor!

@inline function gen_unit_spinor!(spin::Vector{Complex{T}}) where T <: AbstractFloat
  field_volume = sqrt(2.0*Float64(length(spin)))
  for i in 1:length(spin);  spin[i] = Complex{T}( 1.0/field_volume , 1.0/field_volume ); end
end #gen_random_spinor!

@inline function gen_const_spinor!(spin::Vector{Complex{T}}, rea = 0.0, img = 0.0) where T <: AbstractFloat
  for i in 1:length(spin);  spin[i] = Complex{T}(rea, img); end
end #gen_random_spinor!

@inline function print_spinor(spin::Vector{Complex{T}}, len) where T <: AbstractFloat
  for i in 1:len; println("Real = ", real(spin[i]), " Imaginary = ", imag(spin[i])); end
end #print_spinor

end #QJuliaUtils

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
"""
```
Kalman{S<:AbstractFloat}
```
### Fields:

- `loglh`: vector of conditional log-likelihoods log p(y_t | y_{1:t-1}), t = 1:T
- `s_T`: state vector in the last period for which data is provided
- `P_T`: variance-covariance matrix for `s_T`
- `s_pred`: `Ns` x `Nt` matrix of s_{t|t-1}, t = 1:T
- `P_pred`: `Ns` x `Ns` x `Nt` array of P_{t|t-1}, t = 1:T
- `s_filt`: `Ns` x `Nt` matrix of s_{t|t}, t = 1:T
- `P_filt`: `Ns` x `Ns` x `Nt` array of P_{t|t}, t = 1:T
- `s_0`: starting-period state vector. If there are presample periods in the
  data, then `s_0` is the state vector at the end of the presample/beginning of
  the main sample
- `P_0`: variance-covariance matrix for `s_0`
- `total_loglh`: log p(y_{1:t})
"""
struct Kalman{S<:AbstractFloat}
    loglh::Vector{S}            # log p(y_t | y_{1:t-1}), t = 1:T
    s_pred::AbstractArray{S}    # s_{t|t-1}, t = 1:T
    P_pred::Array{S, 3}         # P_{t|t-1}, t = 1:T
    s_filt::AbstractArray{S}    # s_{t|t}, t = 1:T
    P_filt::Array{S, 3}         # P_{t|t}, t = 1:T
    s_0::Vector{S}              # s_0
    P_0::AbstractArray{S}       # P_0
    s_T::Vector{S}              # s_{T|T}
    P_T::AbstractArray{S}       # P_{T|T}
    total_loglh::S              # log p(y_{1:t})
end

function Kalman(loglh::Vector{S},
                s_pred::AbstractArray{S}, P_pred::Array{S, 3},
                s_filt::AbstractArray{S}, P_filt::Array{S, 3},
                s_0::Vector{S}, P_0::AbstractArray{S},
                s_T::Vector{S}, P_T::AbstractArray{S}) where S<:AbstractFloat

    return Kalman{S}(loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T, sum(loglh))
end

function Base.getindex(K::Kalman, d::Symbol)
    if d in (:loglh, :s_pred, :P_pred, :s_filt, :P_filt, :s_0, :P_0, :s_T, :P_T, :total_loglh)
        return getfield(K, d)
    else
        throw(KeyError(d))
    end
end

function Base.getindex(kal::Kalman, inds::Union{Int, UnitRange{Int}})
    t0 = first(inds)
    t1 = last(inds)

    return DSGE.Kalman(kal[:loglh][inds],        # loglh
                       kal[:s_pred][:,    inds], # s_pred
                       kal[:P_pred][:, :, inds], # P_pred
                       kal[:s_filt][:,    inds], # filt
                       kal[:P_filt][:, :, inds], # P_filt
                       kal[:s_filt][:,    t0],   # s_0
                       kal[:P_filt][:, :, t0],   # P_0
                       kal[:s_filt][:,    t1],   # s_T
                       kal[:P_filt][:, :, t1],   # P_T
                       sum(kal[:loglh][inds]))   # total_loglh
end

function Base.cat(m::AbstractDSGEModel, k1::Kalman{S},
                  k2::Kalman{S}; allout::Bool = true) where S<:AbstractFloat

    loglh  = cat(k1[:loglh], k2[:loglh], dims = 1)
    s_pred = cat(k1[:s_pred], k2[:s_pred], dims = 2)
    P_pred = cat(k1[:P_pred], k2[:P_pred], dims = 3)
    s_filt = cat(k1[:s_filt], k2[:s_filt], dims = 2)
    P_filt = cat(k1[:P_filt], k2[:P_filt], dims = 3)
    s_0    = k1[:s_0]
    P_0    = k1[:P_0]
    s_T    = k2[:s_T]
    P_T    = k2[:P_T]
    total_loglh = k1[:total_loglh] + k2[:total_loglh]

    return Kalman(loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T, total_loglh)
end

"""
```
zlb_regime_indices(m, data, start_date = date_presample_start(m))
```

Returns a Vector{UnitRange{Int64}} of index ranges for the pre- and post-ZLB
regimes. The optional argument `start_date` indicates the first quarter of
`data`.
"""
function zlb_regime_indices(m::AbstractDSGEModel{S}, data::AbstractArray,
                            start_date::Dates.Date=date_presample_start(m)) where S<:AbstractFloat
    T = size(data, 2)
    if n_mon_anticipated_shocks(m) > 0 && !isempty(data)
        if start_date < date_presample_start(m)
            error("Start date $start_date must be >= date_presample_start(m)")

        elseif 0 < subtract_quarters(date_zlb_start(m), start_date) < T
            n_nozlb_periods = subtract_quarters(date_zlb_start(m), start_date)
            regime_inds::Vector{UnitRange{Int64}} = [1:n_nozlb_periods, (n_nozlb_periods+1):T]
        else
            regime_inds = UnitRange{Int64}[1:T]
        end
    else
        regime_inds = UnitRange{Int64}[1:T]
    end
    return regime_inds
end

"""
```
function zlb_plus_regime_indices(m::AbstractDSGEModel{S}, data::AbstractArray,
                                 start_date::Dates.Date=date_presample_start(m)) where S<:AbstractFloat
```
returns (1) a Vector{UnitRange{Int64}} of index ranges for regime switches with the pre- and post-ZLB
regimes spliced into the regime switches, (2) the index in the vector indicating where the ZLB starts,
and (3) a Boolean for whether the ZLB coincides with that regime or occurs in the middle of the regime.

The optional argument `start_date` indicates the first quarter of
`data`. Use the `Setting` with key `:regime_dates` to set the start dates of different regimes (excluding
the ZLB regime), and use the `Setting` key `:date_zlb_start` to set the start of the post-ZLB regime.
"""
function zlb_plus_regime_indices(m::AbstractDSGEModel{S}, data::AbstractArray,
                                 start_date::Dates.Date=date_presample_start(m)) where S<:AbstractFloat

    T = size(data, 2)
    if !isempty(data)
        # Calculate the number of periods since start date for each regime
        if haskey(m.settings, :n_hist_regimes)
            n_regime_periods = Vector{Int}(undef, get_setting(m, :n_hist_regimes) + get_setting(m, :n_cond_regimes))
            for k in 1:(get_setting(m, :n_hist_regimes)  + get_setting(m, :n_cond_regimes))
                n_regime_periods[k] = subtract_quarters(get_setting(m, :regime_dates)[k], start_date)
            end
        else
            n_regime_periods = Vector{Int}(undef, length(get_setting(m, :regime_dates)))
            for (k, v) in get_setting(m, :regime_dates)
                n_regime_periods[k] = subtract_quarters(v, start_date)
            end
        end

        if start_date < date_presample_start(m)
            error("Start date $start_date must be >= date_presample_start(m)")
        elseif 0 < subtract_quarters(date_zlb_start(m), start_date) < T
            regime_inds = Vector{UnitRange{Int64}}(undef, length(n_regime_periods) + 1)

            n_nozlb_periods = subtract_quarters(date_zlb_start(m), start_date) # number of periods since start date for start of ZLB

            # Get index of next regime after ZLB starts.
            # Note that it cannot be 1 b/c the first regime starts at the start date
            i_splice_zlb = findfirst(n_nozlb_periods .< n_regime_periods)

            if isnothing(i_splice_zlb)
                check_zlb_coincide, i_zlb_start = if n_nozlb_periods == n_regime_periods[end]
                    n_regime_periods[end], length(n_regime_periods)
                else
                    -1, 0
                end
            else
                i_zlb_start = i_splice_zlb - 1
                check_zlb_coincide = n_regime_periods[i_zlb_start]
            end
            if check_zlb_coincide == n_nozlb_periods
                # In this case, the post-ZLB coincides with a specific regime,
                # so we return the regime indices w/out splicing the ZLB in.
                return regime_indices(m, start_date, iterate_quarters(start_date, T - 1)), i_zlb_start, false
            else
                # Populate vector of regime indices
                if isnothing(i_splice_zlb) # post-ZLB is the last regime
                    for reg in 1:(length(n_regime_periods) - 1)
                        regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[end - 1] = (n_regime_periods[end] + 1):n_nozlb_periods
                    regime_inds[end]     = (n_nozlb_periods + 1):T

                    i_zlb_start = length(regime_inds)
                elseif i_splice_zlb > 2 # at least one full regime before ZLB starts
                    regime_inds[1] = 1:n_regime_periods[2]
                    for reg in 2:i_splice_zlb - 2 # if i_splice_zlb == 3, then this loop does not run
                        regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[i_splice_zlb - 1] = (n_regime_periods[i_splice_zlb - 1] + 1):n_nozlb_periods
                    regime_inds[i_splice_zlb]     = (n_nozlb_periods + 1):(n_regime_periods[i_splice_zlb])

                    i_zlb_start = i_splice_zlb

                    # Index reg + 1 b/c we have spliced pre- and post- ZLB regime in
                    for reg in i_splice_zlb:(length(n_regime_periods) - 1)
                        regime_inds[reg + 1] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[end] = (n_regime_periods[end] + 1):T
                else # post- ZLB regimes start in the very first regime specified by user
                    regime_inds[1] = 1:n_nozlb_periods
                    regime_inds[2] = (n_nozlb_periods + 1):n_regime_periods[2]

                    i_zlb_start = 2

                    # Index reg + 1 b/c spliced pre- and post-ZLB regime in
                    for reg in 2:(length(n_regime_periods) - 1)
                        regime_inds[reg + 1] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
                    end
                    regime_inds[end] = (n_regime_periods[end] + 1):T
                end
                splice_zlb_regime = true
            end
        else
            # This is the case that date_zlb_start <= start_date so the first regime is the post-ZLB regime (no pre-ZLB)
            # Since there is no ZLB regime to splice in, we return the regime_inds corresponding to the user-specified regimes
            regime_inds = Vector{UnitRange{Int64}}(undef, length(n_regime_periods))
            for reg in 1:(length(n_regime_periods) - 1)
                regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
            end
            regime_inds[end] = (n_regime_periods[end] + 1):T

            i_zlb_start = 1 # start immediately in post-ZLB

            splice_zlb_regime = date_zlb_start(m) == start_date ? false : true
        end
    else # Empty, so we ignore regime switching
        regime_inds = UnitRange{Int64}[1:T]
        i_zlb_start = 0
        splice_zlb_regime = false
    end
    # Remove regimes that only are for the forecast period (i.e. bigger than number of rows in dataframe)
    filter!(x -> last(x) <= T, regime_inds)
    return regime_inds, i_zlb_start, splice_zlb_regime
end

"""
```
function regime_indices(m::AbstractDSGEModel{S}, start_date::Dates.Date = date_presample_start(m),
                        end_date::Dates.Date = prev_quarter(date_forecast_start(m))) where S<:AbstractFloat
```
returns a `Vector{UnitRange{Int64}}` of index ranges for regime switches
The optional argument `start_date` indicates the first quarter of
`data`. The optimal argument `end_date` indicates the last quarter of `data`.
Use the `Setting` with key `:regime_dates` to set the start dates of different regimes.

This function differs from `zlb_plus_regime_indices` because it does not splice the
pre- and post-ZLB regimes into the output of index ranges for regime switches. When it is not necessary
to explicitly zero out anticipated monetary policy shocks to establish the pre-ZLB regime, there is no need
to distinguish between the pre-ZLB and post-ZLB regimes. For example, with shock decompositions,
the user needs to specify the history of filtered shocks, which should have zeroes already for anticipated shocks
in the pre-ZLB regime. Thus, there is no need to explicitly account for the pre-ZLB regime.
"""
function regime_indices(m::AbstractDSGEModel{S}, start_date::Dates.Date = date_presample_start(m),
                        end_date::Dates.Date = prev_quarter(date_forecast_start(m))) where S<:AbstractFloat
    T = subtract_quarters(end_date, start_date) + 1

    n_regime_periods = Vector{Int}(undef, length(get_setting(m, :regime_dates)))
    last_regime = findfirst(values(get_setting(m, :regime_dates)) .> end_date) # Last regime based on end date
    if isnothing(last_regime)
        last_regime = length(get_setting(m, :regime_dates))
    end
    for (k, v) in get_setting(m, :regime_dates)
        n_regime_periods[k] = subtract_quarters(v, start_date)
    end
    regime_inds = Vector{UnitRange{Int64}}(undef, min(length(n_regime_periods), last_regime))

    for reg in 1:(min(length(n_regime_periods), last_regime) - 1)
        regime_inds[reg] = (n_regime_periods[reg] + 1):n_regime_periods[reg + 1]
    end
    regime_inds[end] = (n_regime_periods[last_regime] + 1):T

    return regime_inds
end

"""
```
zlb_regime_matrices(m, system, start_date = date_presample_start(m))
```
Returns `TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs`, an 8-tuple of
`Vector{AbstractArray{S}}`s and `Vector{Vector{S}}`s of system matrices for the pre-
and post-ZLB regimes. Of these, only `QQ` changes from pre- to post-ZLB: the
entries corresponding to anticipated shock variances are zeroed out pre-ZLB.
"""
function zlb_regime_matrices(m::AbstractDSGEModel{S}, system::System{S},
                             start_date::Dates.Date=date_presample_start(m)) where S<:AbstractFloat
    if n_mon_anticipated_shocks(m) > 0
        if start_date < date_presample_start(m)
            error("Start date $start_date must be >= date_presample_start(m)")

        # TODO: This technically doesn't handle the case where the end_date of the sample
        # is before the start of the ZLB
        elseif date_presample_start(m) <= start_date <= date_zlb_start(m)
            n_regimes = 2

            shock_inds = inds_shocks_no_ant(m)
            QQ_ZLB = system[:QQ]
            QQ_preZLB = zeros(size(QQ_ZLB))
            QQ_preZLB[shock_inds, shock_inds] = QQ_ZLB[shock_inds, shock_inds]
            QQs = Matrix{S}[QQ_preZLB, QQ_ZLB]

        elseif date_zlb_start(m) < start_date
            n_regimes = 1
            QQs = Matrix{S}[system[:QQ]]
        end
    else
        n_regimes = 1
        QQs = Matrix{S}[system[:QQ]]
    end

    TTTs = fill(system[:TTT], n_regimes)
    RRRs = fill(system[:RRR], n_regimes)
    CCCs = fill(system[:CCC], n_regimes)
    ZZs  = fill(system[:ZZ], n_regimes)
    DDs  = fill(system[:DD], n_regimes)
    EEs  = fill(system[:EE], n_regimes)

    return TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs
end

"""
```
zlb_plus_regime_matrices(m, system, start_date = date_presample_start(m);
    n_regimes = 0, ind_zlb_start = 0, splice_zlb_regime = true)
```
Returns `TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs`, an 8-tuple of
`Vector{AbstractArray{S}}`s and `Vector{Vector{S}}`s of system matrices for the regimes specified
by `system`, along with the pre- and post-ZLB regimes spliced into their correct location.
See `?zlb_regime_matrices` to see how the pre- and post-ZLB matrices are handled.

The keyword `n_regimes` indicates how many regimes to return in the case that fewer
"""
function zlb_plus_regime_matrices(m::AbstractDSGEModel{S}, system::RegimeSwitchingSystem{S},
                                  n_reg::Int, start_date::Dates.Date = date_presample_start(m);
                                  ind_zlb_start::Int = 0,
                                  splice_zlb_regime::Bool = true) where S<:AbstractFloat

    if n_reg <= 0 || ind_zlb_start <= 0
        throw(DomainError())
    end

    if start_date < date_presample_start(m)
        error("Start date $start_date must be >= date_presample_start(m)")

        # TODO: This technically doesn't handle the case where the end_date of the sample
        # is before the start of the ZLB
    elseif date_presample_start(m) <= start_date <= date_zlb_start(m)
        TTTs = Vector{Matrix{S}}(undef, n_reg)
        RRRs = Vector{Matrix{S}}(undef, n_reg)
        CCCs = Vector{Vector{S}}(undef, n_reg)
        ZZs  = Vector{Matrix{S}}(undef, n_reg)
        DDs  = Vector{Vector{S}}(undef, n_reg)
        QQs  = Vector{Matrix{S}}(undef, n_reg)
        EEs  = Vector{Matrix{S}}(undef, n_reg)

        shock_inds = inds_shocks_no_ant(m)

        if splice_zlb_regime
            if ind_zlb_start == 2
                if n_mon_anticipated_shocks(m) == 0
                    QQ_preZLB = copy(system[1, :QQ]) # anticipated shocks are zero already
                else
                    QQ_preZLB                         = zeros(size(system[1, :QQ]))
                    QQ_preZLB[shock_inds, shock_inds] = system[1, :QQ][shock_inds, shock_inds]
                end
                QQ_postZLB = copy(system[1, :QQ])

                TTTs[1:2] = [system[1, :TTT], system[1, :TTT]]::Vector{Matrix{S}}
                RRRs[1:2] = [system[1, :RRR], system[1, :RRR]]::Vector{Matrix{S}}
                CCCs[1:2] = [system[1, :CCC], system[1, :CCC]]::Vector{Vector{S}}
                ZZs[1:2]  = [system[1, :ZZ],  system[1, :ZZ]]::Vector{Matrix{S}}
                DDs[1:2]  = [system[1, :DD],  system[1, :DD]]::Vector{Vector{S}}
                EEs[1:2]  = [system[1, :EE],  system[1, :EE]]::Vector{Matrix{S}}
                QQs[1]    = QQ_preZLB
                QQs[2]    = QQ_postZLB

                if n_reg >= 3
                    input = 2:(n_reg - 1) # minus 1 b/c spliced in a ZLB regime

                    TTTs[3:end] = [system[i, :TTT] for i in input]::Vector{Matrix{S}}
                    RRRs[3:end] = [system[i, :RRR] for i in input]::Vector{Matrix{S}}
                    CCCs[3:end] = [system[i, :CCC] for i in input]::Vector{Vector{S}}
                    ZZs[3:end]  = [system[i, :ZZ] for i in input]::Vector{Matrix{S}}
                    DDs[3:end]  = [system[i, :DD] for i in input]::Vector{Vector{S}}
                    QQs[3:end]  = [system[i, :QQ] for i in input]::Vector{Matrix{S}}
                    EEs[3:end]  = [system[i, :EE] for i in input]::Vector{Matrix{S}}
                end
            elseif ind_zlb_start >= 3
                # Populate matrices for regimes before the ZLB split
                pre_zlb_input = 1:(ind_zlb_start - 2) # minus 2 b/c ind_zlb_start - 1 will be the pre-ZLB regime
                splice_reg    = ind_zlb_start - 1

                TTTs[pre_zlb_input] = [system[i, :TTT] for i in pre_zlb_input]::Vector{Matrix{S}}
                RRRs[pre_zlb_input] = [system[i, :RRR] for i in pre_zlb_input]::Vector{Matrix{S}}
                CCCs[pre_zlb_input] = [system[i, :CCC] for i in pre_zlb_input]::Vector{Vector{S}}
                ZZs[pre_zlb_input]  = [system[i, :ZZ] for i in pre_zlb_input]::Vector{Matrix{S}}
                DDs[pre_zlb_input]  = [system[i, :DD] for i in pre_zlb_input]::Vector{Vector{S}}
                EEs[pre_zlb_input]  = [system[i, :EE] for i in pre_zlb_input]::Vector{Matrix{S}}

                if n_mon_anticipated_shocks(m) == 0
                    QQs[pre_zlb_input] = [system[i, :QQ] for i in pre_zlb_input]::Vector{Matrix{S}}
                else
                    for i in pre_zlb_input
                        QQ_preZLB                         = zeros(size(system[i, :QQ]))
                        QQ_preZLB[shock_inds, shock_inds] = system[i, :QQ][shock_inds, shock_inds]
                        QQs[i]                            = QQ_preZLB
                    end
                end

                # Regime in which the ZLB occurs
                if n_mon_anticipated_shocks(m) == 0
                    QQ_preZLB = copy(system[splice_reg, :QQ]) # anticipated shocks are zero already
                else
                    QQ_preZLB                         = zeros(size(system[splice_reg, :QQ]))
                    QQ_preZLB[shock_inds, shock_inds] = system[splice_reg, :QQ][shock_inds, shock_inds]
                end
                QQ_postZLB = copy(system[splice_reg, :QQ])
                TTTs[splice_reg:ind_zlb_start] = [system[splice_reg, :TTT], system[splice_reg, :TTT]]::Vector{Matrix{S}}
                RRRs[splice_reg:ind_zlb_start] = [system[splice_reg, :RRR], system[splice_reg, :RRR]]::Vector{Matrix{S}}
                CCCs[splice_reg:ind_zlb_start] = [system[splice_reg, :CCC], system[splice_reg, :CCC]]::Vector{Vector{S}}
                ZZs[splice_reg:ind_zlb_start]  = [system[splice_reg, :ZZ],  system[splice_reg, :ZZ]]::Vector{Matrix{S}}
                DDs[splice_reg:ind_zlb_start]  = [system[splice_reg, :DD],  system[splice_reg, :DD]]::Vector{Vector{S}}
                EEs[splice_reg:ind_zlb_start]  = [system[splice_reg, :EE],  system[splice_reg, :EE]]::Vector{Matrix{S}}
                QQs[splice_reg]     = QQ_preZLB
                QQs[ind_zlb_start] = QQ_postZLB

                # Handle regimes after the ZLB split
                if n_reg >= ind_zlb_start + 1
                    post_zlb_input = ind_zlb_start:(n_reg - 1) # minus 1 to both b/c spliced in the ZLB regime

                    TTTs[4:end] = [system[i, :TTT] for i in post_zlb_input]::Vector{Matrix{S}}
                    RRRs[4:end] = [system[i, :RRR] for i in post_zlb_input]::Vector{Matrix{S}}
                    CCCs[4:end] = [system[i, :CCC] for i in post_zlb_input]::Vector{Vector{S}}
                    ZZs[4:end]  = [system[i, :ZZ] for i in post_zlb_input]::Vector{Matrix{S}}
                    DDs[4:end]  = [system[i, :DD] for i in post_zlb_input]::Vector{Vector{S}}
                    QQs[4:end]  = [system[i, :QQ] for i in post_zlb_input]::Vector{Matrix{S}}
                    EEs[4:end]  = [system[i, :EE] for i in post_zlb_input]::Vector{Matrix{S}}
                end
             end
        else
            if ind_zlb_start > 1
                pre_zlb_input = 1:(ind_zlb_start - 1) # minus 1 b/c post-ZLB coincides with a regime

                TTTs[pre_zlb_input] = [system[i, :TTT] for i in pre_zlb_input]::Vector{Matrix{S}}
                RRRs[pre_zlb_input] = [system[i, :RRR] for i in pre_zlb_input]::Vector{Matrix{S}}
                CCCs[pre_zlb_input] = [system[i, :CCC] for i in pre_zlb_input]::Vector{Vector{S}}
                ZZs[pre_zlb_input]  = [system[i, :ZZ] for i in pre_zlb_input]::Vector{Matrix{S}}
                DDs[pre_zlb_input]  = [system[i, :DD] for i in pre_zlb_input]::Vector{Vector{S}}
                EEs[pre_zlb_input]  = [system[i, :EE] for i in pre_zlb_input]::Vector{Matrix{S}}

                if n_mon_anticipated_shocks(m) == 0
                    QQs[pre_zlb_input] = Matrix{S}[system[i, :QQ] for i in pre_zlb_input]::Vector{Matrix{S}}
                else
                    for i in pre_zlb_input
                        QQ_preZLB                         = zeros(size(system[i, :QQ]))
                        QQ_preZLB[shock_inds, shock_inds] = system[i, :QQ][shock_inds, shock_inds]
                        QQs[i]                            = QQ_preZLB
                    end
                end
            end

            # Handle regimes after the ZLB split
            post_zlb_input = ind_zlb_start:n_reg # no minus 1 b/c didn't splice in the regime

            TTTs[post_zlb_input] = [system[i, :TTT] for i in post_zlb_input]::Vector{Matrix{S}}
            RRRs[post_zlb_input] = [system[i, :RRR] for i in post_zlb_input]::Vector{Matrix{S}}
            CCCs[post_zlb_input] = [system[i, :CCC] for i in post_zlb_input]::Vector{Vector{S}}
            ZZs[post_zlb_input]  = [system[i, :ZZ] for i in post_zlb_input]::Vector{Matrix{S}}
            DDs[post_zlb_input]  = [system[i, :DD] for i in post_zlb_input]::Vector{Vector{S}}
            QQs[post_zlb_input]  = [system[i, :QQ] for i in post_zlb_input]::Vector{Matrix{S}}
            EEs[post_zlb_input]  = [system[i, :EE] for i in post_zlb_input]::Vector{Matrix{S}}
         end
    elseif date_zlb_start(m) < start_date
        # We always use post-ZLB matrices in this case,
        # hence we just return the matrices of the `RegimeSwitchingSystem`
        input = 1:n_reg
        TTTs = [system[i, :TTT] for i in input]::Vector{Matrix{S}}
        RRRs = [system[i, :RRR] for i in input]::Vector{Matrix{S}}
        CCCs = [system[i, :CCC] for i in input]::Vector{Vector{S}}
        ZZs  = [system[i, :ZZ] for i in input]::Vector{Matrix{S}}
        DDs  = [system[i, :DD] for i in input]::Vector{Vector{S}}
        QQs  = [system[i, :QQ] for i in input]::Vector{Matrix{S}}
        EEs  = [system[i, :EE] for i in input]::Vector{Matrix{S}}
    end

    return TTTs, RRRs, CCCs, QQs, ZZs, DDs, EEs
end


text:
__precompile__()

module Rectangle

export  Line,
            isHorizontal, isVertical, length, parallelogram_area,
            ratio, intersects, merge_axis_aligned, sx, sy, ex, ey,
            intersect_axis_aligned, start, endof, horiz_desc, vert_asc,
            xplot, yplot,
        Rect,
            x, y, lx, ly, rx, ry, w, h, area, perimeter, hlines,
            vlines, lines, olines,
            union, intersects, inside, cg,
            to_plot_shape,
            projectX, projectY,
            visibleX, visibleY,
            has_x_overlap, has_y_overlap,
            avg_min_dist, min_dist,
            vline_xsection, hline_xsection,
        OrderedRectMapX, OrderedRectMapY,
            create_ordered_map, get_intersect_data, insert_rect!, delete_rect!,
        pcTol,
            AbstractBinaryTree, AbstractBST, RBTree, BinarySearchTree, Iterator,
            IntervalTree, Interval, intersects, depth, XYTree, XYData,
            get_values

include("utils.jl")
include("Line.jl")
include("bst.jl")
include("interval.jl")
include("xytree.jl")
include("Rect.jl")

end # module


text:
# Synchronization (B.6)

export sync_threads, sync_warp

"""
    sync_threads()

Waits until all threads in the thread block have reached this point and all global and
shared memory accesses made by these threads prior to `sync_threads()` are visible to all
threads in the block.
"""
@inline sync_threads() = ccall("llvm.nvvm.barrier0", llvmcall, Cvoid, ())

"""
    sync_warp(mask::Integer=0xffffffff)

Waits threads in the warp, selected by means of the bitmask `mask`, have reached this point
and all global and shared memory accesses made by these threads prior to `sync_warp()` are
visible to those threads in the warp. The default value for `mask` selects all threads in
the warp.
"""
sync_warp

if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
    @inline function sync_warp(mask::Integer=0xffffffff)
        @asmcall("bar.warp.sync \$0;", "r", true,
                 Cvoid, Tuple{UInt32}, convert(UInt32, mask))
    end
else
    @inline sync_warp(mask::Integer=0xffffffff) = nothing
end


text:
"""
The aggregate value of an image over a mask can be computed by the aggregate function. The function multiplies the image and the mask (elementwise) and then performs the summation.
```
rand_image = rand(10, 10)
rand_mask = rand(0:1, 10, 10)
aggregate(rand_image, rand_mask)
```
"""
function aggregate(image, mask)
    masked_image = copy(image).*mask
    return sum(masked_image)
end

"""
To find the time series of aggregate values of a datacube over a mask, using the ```aggregate_timeseries``` function

```
rand_datacube = rand(10, 10, 10)
rand_mask = rand(0:1, 10, 10)
aggregate_timeseries(rand_datacube, rand_mask)
```
"""
function aggregate_timeseries(datacube::Array{T, 3}, mask) where T <: Real
    tmp = apply_mask(copy(datacube), mask)
    lights = []
    for i in 1:size(datacube)[3]
        push!(lights, sum(tmp[:,:,i]))
    end
    return Float32.(lights)
end

"""
To find the time series of aggregate values of a datacube over a mask, using the ```aggregate_timeseries``` function

```
rand_datacube = rand(10, 10, 10)
rand_mask = rand(0:1, 10, 10)
aggregate_timeseries(rand_datacube, rand_mask)
```
"""
function aggregate_timeseries(datacube::VectorOfArray{Any,3}, mask)
    tmp = apply_mask(copy(datacube), mask)
    lights = []
    for i in 1:size(datacube)[3]
        push!(lights, sum(tmp[i]))
    end
    return Float32.(lights)
end

"""
Computes the time series of aggregate value of datacube for each row of a shapefile
# Example:
```
rand_datacube = rand(10, 10, 10)
shapefile_df  = load_shapefile("assets/mumbai_map/mumbai_districts.shp")
aggregate_dataframe(MUMBAI_COORDINATE_SYSTEM, rand_datacube, shapefile_df, "District")
```
"""
function aggregate_dataframe(geometry::CoordinateSystem, datacube, shapefile_df, attribute)
    df = DataFrame()
    datacube = sparse_cube(datacube)
    GC.gc()
    @showprogress for i in 1:length(shapefile_df[:, 1])
        shapefile_row = shapefile_df[i, :]
        geom_polygon = polygon_mask(geometry, shapefile_row)
        df[!, shapefile_df[!, attribute][i]] = aggregate_timeseries(datacube, geom_polygon)
        i = i + 1
    end
    return df
end

"""
Computes the time series of aggregate value per area of datacube for each row of a shapefile
# Example:
```
rand_datacube = rand(10, 10, 10)
shapefile_df  = load_shapefile("assets/mumbai_map/mumbai_districts.shp")
aggregate_per_area_dataframe(MUMBAI_COORDINATE_SYSTEM, rand_datacube, shapefile_df, "District")
```
"""
function aggregate_per_area_dataframe(geometry::CoordinateSystem, datacube, shapefile_df, attribute, res = 15)
    df = DataFrame()
    datacube = sparse_cube(datacube)
    GC.gc()
    @showprogress for i in 1:length(shapefile_df[:, 1])
        shapefile_row = shapefile_df[i, :]
        geom_polygon = polygon_mask(geometry, shapefile_row)
        df[!, shapefile_df[!, attribute][i]] = aggregate_timeseries(datacube, geom_polygon) ./mask_area(geometry, geom_polygon, res)
        i = i + 1
    end
    return df
end


text:
using BH_ED
using Test

@testset "BH_ED.jl" begin
    r1 = BH_solver(N   = 6
                , Lp   = 3
                , L    = 3*6
                , μ    = 1.0
                , J    = 1.0
                , V    = 1.0
                , ε    = 1.0
    )
    t_cf_r = [2.1127883913424594,
    1.0465557381373054,
    1.3234461149460939,
    1.270276914830182,
    1.3234461149460939,
    1.0465557381373054,
    2.1127883913424594]

    @test abs(r1[1][1] + 14.60579110615412) <= 0.1
    @test abs(r1[3] - 6.813133331993588) <= 0.1
    @test abs( sum(r1[4] .- t_cf_r) ) <= 0.1
end


text:
using Gtk.ShortNames
using GtkUtilities, GtkUtilities.Graphics

c = Canvas()
win = Window(c, "RubberBandCanvas")
draw(c) do widget
    ctx = getgc(widget)
    h = height(widget)
    w = width(widget)
    # Paint red rectangle
    rectangle(ctx, 0, 0, w, h/2)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)
    # Paint white rectangle
    rectangle(ctx, 0, h/2, w, 3h/4)
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)
    # Paint blue rectangle
    rectangle(ctx, 0, 3h/4, w, h/4)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)
end

# Set up a rubberband-on-click. Note that `panzoom.jl` shows an
# alternative way to set this up, one that doesn't mess with
# any pre-existing `button1press` function that the user
# may have defined.
c.mouse.button1press = (widget, event) -> begin
    if event.event_type == Gtk.GdkEventType.BUTTON_PRESS
        rubberband_start(widget, event.x, event.y, (widget, bb) -> @show bb)
    end
end

showall(c)


text:
x = CertainValue(2.0)

test_constraints = [
    NoConstraint(),
    TruncateLowerQuantile(0.2),
    TruncateUpperQuantile(0.2),
    TruncateQuantiles(0.2, 0.8),
    TruncateMaximum(50),
    TruncateMinimum(-50),
    TruncateRange(-50, 50),
    TruncateStd(1)
]

T = eltype(x)

@test rand(x) isa T
@test rand(x, 10) isa Vector

@test resample(x) isa T
@test resample(x, 10) isa Vector
@test all(resample(x, 10) .== x.value)

for constraint in test_constraints
    @test resample(x, constraint) isa T
    @test resample(x, constraint, 10) isa Vector
    @test all(resample(x, constraint, 10) .== x.value)
end