text:
using Base.Dates
using Base.Test

using MarketData

using TimeSeries

@testset "timearray" begin

@testset "field extraction methods work" begin
    @testset "timestamp, values, colnames and meta" begin
        @test typeof(timestamp(cl)) == Array{Date,1}
        @test typeof(values(cl))    == Array{Float64,1}
        @test typeof(colnames(cl))  == Array{String,1}
        @test meta(mdata)           == "Apple"
    end
end

@testset "type constructors allow views" begin
    source_rows = 101:121
    source_cols = 1:size(AAPL.values)[2]
    tstamps = view(AAPL.timestamp, source_rows)
    tvalues = view(AAPL.values, source_rows, source_cols)

    AAPL1 = TimeArray(AAPL.timestamp[source_rows],
                      AAPL.values[source_rows, source_cols],
                      AAPL.colnames, AAPL.meta)

    AAPL2 = TimeArray(tstamps, tvalues, AAPL.colnames, AAPL.meta)

    @testset "match first date" begin
        @test AAPL1[1].timestamp == AAPL2[1].timestamp
    end

    @testset "match first values" begin
        @test AAPL1[1].values == AAPL2[1].values
    end

    @testset "match all values" begin
        @test AAPL1.values == AAPL2.values
    end
end

@testset "type constructors enforce invariants" begin
    mangled_stamp = vcat(cl.timestamp[200:end], cl.timestamp[1:199])
    dupe_stamp    = vcat(cl.timestamp[1:499], cl.timestamp[499])
    dupe_cnames   = rename(AAPL, ["a", "b", "c", "a", "a", "b", "d", "e", "e", "e", "e", "f"])

    @testset "unequal length between values and timestamp fails" begin
        @test_throws(
            DimensionMismatch,
            TimeArray(cl.timestamp, cl.values[2:end], ["Close"]))
    end

    @testset "unequal length between colnames and array width fails" begin
        @test_throws(
            DimensionMismatch,
            TimeArray(cl.timestamp, cl.values, ["Close", "Open"]))
    end

    @testset "duplicate timestamp values fails" begin
        @test_throws(
            ArgumentError,
            TimeArray(dupe_stamp, cl.values, ["Close"]))
    end

    @testset "mangled order of timestamp values fails" begin
        @test_throws(
            ArgumentError,
            TimeArray(mangled_stamp, cl.values, ["Close"]))
    end

    @testset "flipping occurs when needed" begin
        @test TimeArray(flipdim(cl.timestamp, 1), flipdim(cl.values, 1),  ["Close"]).timestamp[1] == Date(2000,1,3)
        @test TimeArray(flipdim(cl.timestamp, 1), flipdim(cl.values, 1),  ["Close"]).values[1]    == 111.94
    end

    @testset "duplicate column names are enumerated by inner constructor" begin
        @test dupe_cnames.colnames[1]  == "a"
        @test dupe_cnames.colnames[2]  == "b"
        @test dupe_cnames.colnames[3]  == "c"
        @test dupe_cnames.colnames[4]  == "a_1"
        @test dupe_cnames.colnames[5]  == "a_2"
        @test dupe_cnames.colnames[6]  == "b_1"
        @test dupe_cnames.colnames[7]  == "d"
        @test dupe_cnames.colnames[8]  == "e"
        @test dupe_cnames.colnames[9]  == "e_1"
        @test dupe_cnames.colnames[10] == "e_2"
        @test dupe_cnames.colnames[11] == "e_3"
        @test dupe_cnames.colnames[12] == "f"
    end
end

@testset "construction without colnames" begin
    no_colnames_one   = TimeArray(cl.timestamp, cl.values)
    no_colnames_multi = TimeArray(AAPL.timestamp, AAPL.values)

    @testset "default colnames to empty String vector" begin
        @test no_colnames_one.colnames   == String[""]
        @test no_colnames_multi.colnames == String["_1", "_2", "_3", "_4", "_5", "_6", "_7", "_8", "_9", "_10", "_11", "_12"]
    end

    @testset "empty colnames forces meta to nothing" begin
        @test no_colnames_one.meta   == nothing
        @test no_colnames_multi.meta == nothing
    end
end

@testset "conversion methods" begin
    @testset "convert works " begin
        @test isa(convert(TimeArray{Float64,1}, (cl.>op)), TimeArray{Float64,1})                == true
        @test isa(convert(TimeArray{Float64,2}, (merge(cl.<op, cl.>op))), TimeArray{Float64,2}) == true
        @test isa(convert(cl.>op), TimeArray{Float64,1})                                        == true
        @test isa(convert(merge(cl.<op, cl.>op)), TimeArray{Float64,2})                         == true
    end
end

@testset "index by integer works with both 1d and 2d time array" begin
    @testset "1d time array" begin
        @test cl[1].timestamp == [Date(2000,1,3)]
        @test cl[1].values    == [111.94]
        @test cl[1].colnames  == ["Close"]
        @test cl[1].meta      == "AAPL"
    end

    @testset "2d time array" begin
        @test ohlc[1].timestamp == [Date(2000,1,3)]
        @test ohlc[1].values    == [104.88 112.5 101.69 111.94]
        @test ohlc[1].colnames  == ["Open", "High", "Low","Close"]
        @test ohlc[1].meta      == "AAPL"
    end
end

@testset "ordered collection methods" begin
    @testset "iterator protocol is valid" begin
        @test !isempty(op)
        @test isempty(op[op .< 0])
        @test start(op)              == 1
        @test next(op, 1)            == ((op.timestamp[1], op.values[1,:]), 2)
        @test done(op, length(op)+1) == true
    end

    @testset "end keyword returns correct index" begin
        @test ohlc[end].timestamp[1] == ohlc.timestamp[end]
    end

    @testset "getindex on single Int and Date" begin
        @test ohlc[1].timestamp              == [Date(2000,1,3)]
        @test ohlc[Date(2000,1,3)].timestamp == [Date(2000,1,3)]
    end

    @testset "getindex on array of Int and Date" begin
        @test ohlc[[1,10]].timestamp                           == [Date(2000,1,3), Date(2000,1,14)]
        @test ohlc[[Date(2000,1,3),Date(2000,1,14)]].timestamp == [Date(2000,1,3), Date(2000,1,14)]
    end

    @testset "getindex on range of Int and Date" begin
        @test ohlc[1:2].timestamp                                  == [Date(2000,1,3), Date(2000,1,4)]
        @test ohlc[1:2:4].timestamp                                == [Date(2000,1,3), Date(2000,1,5)]
        @test ohlc[Int8(1):Int8(2):Int8(4)].timestamp              == [Date(2000,1,3), Date(2000,1,5)]
        @test ohlc[Date(2000,1,3):Day(1):Date(2000,1,4)].timestamp == [Date(2000,1,3), Date(2000,1,4)]
    end

    @testset "getindex on range of DateTime when only Date is in timestamp" begin
        @test_throws(
            MethodError,
            ohlc[DateTime(2000,1,3,0,0,0)])
        @test_throws(
            MethodError,
            ohlc[[DateTime(2000,1,3,0,0,0),DateTime(2000,1,14,0,0,0)]])
        @test_throws(
            MethodError,
            ohlc[DateTime(2000,1,3,0,0,0):Day(1):DateTime(2000,1,4,0,0,0)])
    end

    @testset "getindex on range of Date" begin
        @test length(cl[Date(2000,1,1):Date(2001,12,31)]) == 500
    end

    @testset "getindex on single column name" begin
        @test size(ohlc["Open"].values, 2)                                        == 1
        @test size(ohlc["Open"][Date(2000,1,3):Day(1):Date(2000,1,14)].values, 1) == 10
    end

    @testset "getindex on multiple column name" begin
        @test ohlc["Open", "Close"].values[1]   == 104.88
        @test ohlc["Open", "Close"].values[2]   == 108.25
        @test ohlc["Open", "Close"].values[501] == 111.94
    end

    @testset "getindex on 1d returns 1d object" begin
        @test isa(cl[1], TimeArray{Float64,1})   == true
        @test isa(cl[1:2], TimeArray{Float64,1}) == true
    end

    @testset "getindex on a 1d Boolean TimeArray returns appropriate rows" begin
        @test ohlc[op .> cl][2].values             == ohlc[4].values
        @test ohlc[op[300:end] .> cl][2].timestamp == ohlc[303].timestamp
        # MethodError, Bool must be 1D-TimeArray
        @test_throws MethodError ohlc[merge(op.>cl, op.<cl)]
    end
end

@testset "Base.size" begin
    @test size(ohlc) == (500, 4)
    @test size(ohlc, 1) == 500
    @test size(ohlc, 2) == 4

    @test size(cl) == (500,)
    @test size(cl, 1) == 500
    @test size(cl, 2) == 1
end

@testset "show methods don't throw errors" begin
    show(ohlc)
    show(ohlc[1:4])
    show(ohlc[1:0])
    show(lag(cl[1:2], padding=true))
end

end  # @testset "timearray"


text:
export newtonmodar

using ForwardDiff
using JSOSolverTemplate
using LinearAlgebra
using CUTEst
using NLPModels, LinearOperators, Krylov, SolverTools, SolverBenchmark

function newtonmodar(nlp ;x :: AbstractVector=copy(nlp.meta.x0),
    atol :: Real=√eps(eltype(x)), rtol :: Real=√eps(eltype(x)), max_time = 30.0, max_iter = 1_000_000, η₁ = 1e-2)
    t₀ = time()
    Δt = time() - t₀
    iter = 0

    f(x) = obj(nlp,x)
    ∇f(x) = grad(nlp,x)
    H(x) = Symmetric(hess(nlp, x), :L)
    fx = f(x)
    gx = ∇f(x)

    x⁺ = similar(x)

    ϵ = atol + rtol * norm(gx)

    status =:unknown

    resolvido = norm(gx) < ϵ
    cansado = Δt > max_time || iter > max_iter
    Hx = H(x)
    β = norm(Hx,2)
    v = Vector{Float64}()
    for i in 1:size(hess(nlp,x))[1]
        Δt = time() - t₀
        if Δt >max_time
            status = :max_time
            break
        end
        push!(v,hess(nlp,x)[i,i])
    end
    if minimum(v) > 0
        τ = 0
    else
        τ = β/2
    end
    F = cholesky(Hx + I*τ, check=false)
    while !(cansado || resolvido)
        while (!issuccess(F))
            Δt = time() - t₀
            if Δt >max_time
                status = :max_time
                break
            end
            τ = max(2τ,β/2)
            F = cholesky(Hx + I*τ, check=false)
        end
        d = -(F \ gx)
        slope = dot(gx,d)

        α = 1.0
        x⁺ = x + α * d
        f⁺ = f(x⁺)
        while f⁺ ≥ fx + η₁ * α * slope
            Δt = time() - t₀
            if Δt >max_time
                status = :max_time
                break
            end
            α = α / 2
            x⁺ = x + α * d
            f⁺ = f(x⁺)
            if α < 1e-8
                status = :small_step
                break
            end
        end
        if status != :unknown
            break
        end
        x .= x⁺
        fx = f⁺
        gx = ∇f(x)
        Hx = H(x)

        resolvido = norm(gx) < ϵ
        Δt = time() - t₀
        iter += 1
        cansado = Δt > max_time || iter > max_iter
    end

    if resolvido
        status = :first_order
    elseif cansado
        if Δt >max_time
            status = :max_time
        elseif iter > max_iter
            status = :max_iter
        end
    end

    return GenericExecutionStats(status,nlp,objective=fx,solution=x,dual_feas=norm(gx),iter = iter, elapsed_time = Δt)
end

ext: jl


text:
using Test
using CategoricalDistributions
using CategoricalArrays
using Distributions
using Random
using StableRNGs # for RNGs stable across all julia versions
rng = StableRNGs.StableRNG(123)

import CategoricalDistributions: classes, decoder, int

@testset "utilities" begin
     @test include("utilities.jl")
end

@testset "types" begin
     @test include("types.jl")
end

@testset "methods" begin
     @test include("methods.jl")
end

@testset "arrays.jl" begin
     @test include("arrays.jl")
end

@testset "arithmetic.jl" begin
     @test include("arithmetic.jl")
end


text:
module SimpleTraits
const curmod = module_name(current_module())

using Compat

# This is basically just adding a few convenience functions & macros
# around Holy Traits.

export Trait, istrait, @traitdef, @traitimpl, @traitfn, Not

# All traits are concrete subtypes of this trait.  SUPER is not used
# but present to be compatible with Traits.jl.
## @doc """
## `abstract Trait{SUPER}`

"""
All Traits are subtypes of abstract type Trait.

A concrete Trait will look like:
```julia
immutable Tr1{X,Y} <: Trait end
```
where X and Y are the types involved in the trait.

(SUPER is not used here but in Traits.jl, thus retained for possible
future compatibility.)
"""
abstract Trait{SUPER}

"""
The set of all types not belonging to a trait is encoded by wrapping
it with Not{}, e.g.  Not{Tr1{X,Y}}
"""
abstract Not{T<:Trait} <: Trait

# Helper to strip an even number of Not{}s off: Not{Not{T}}->T
stripNot{T<:Trait}(::Type{T}) = T
stripNot{T<:Trait}(::Type{Not{T}}) = Not{T}
stripNot{T<:Trait}(::Type{Not{Not{T}}}) = stripNot(T)

"""
A trait is defined as full-filled if this function is the identity
function for that trait.  Otherwise it returns the trait wrapped in
`Not`.

Example:
```
trait(IsBits{Int}) # returns IsBits{Int}
trait(IsBits{Array}) # returns Not{IsBits{Array}}
```

Usually this function is defined when using the `@traitimpl` macro.

However, instead of using `@traitimpl` one can define a method for
`trait` to implement a trait, see the README.
"""
trait{T<:Trait}(::Type{T}) = Not{T}
trait{T<:Trait}(::Type{Not{T}}) = trait(T)

## Under the hood, a trait is then implemented for specific types by
## defining:
#   trait(::Type{Tr1{Int,Float64}}) = Tr1{Int,Float64}
# or
#   trait{I<:Integer,F<:FloatingPoint}(::Type{Tr1{I,F}}) = Tr1{I,F}
#
# Note due to invariance, this does probably not the right thing:
#   trait(::Type{Tr1{Integer,FloatingPoint}}) = Tr1{Integer, FloatingPoint}

"""
This function checks whether a trait is fulfilled by a specific
set of types.
```
istrait(Tr1{Int,Float64}) => return true or false
```
"""
istrait(::Any) = error("Argument is not a Trait.")
istrait{T<:Trait}(tr::Type{T}) = trait(tr)==stripNot(tr) ? true : false # Problem, this can run into issue #265
                                                                        # thus is redefine when traits are defined
"""

Used to define a trait.  Traits, like types, are camel cased.  I
suggest to start them with a verb, e.g. `IsImmutable`, to distinguish
them from actual types, which are usually nouns.

Traits need to have one or more (type-)parameters to specify the type
to which the trait is applied.  For instance `IsImmutable{Int}`
signifies that `Int` is part of `IsImmutable` (although whether that
is true needs to be checked with the `istrait` function).  Most traits
will be one-parameter traits, however, several parameters are useful
when there is a "contract" between several types.

Examples:
```julia
@traitdef IsFast{X}
@traitdef IsSlow{X,Y}
```
"""
macro traitdef(tr)
    :(immutable $(esc(tr)) <: Trait end)
end

"""
Used to add a type or type-tuple to a trait.  By default a type does
not belong to a trait.

Example:
```julia
@traitdef IsFast{X}
@traitimpl IsFast{Array{Int,1}}
```
"""
macro traitimpl(tr)
    # makes
    # trait{X1<:Int,X2<:Float64}(::Type{Tr1{X1,X2}}) = Tr1{X1,X2}
    typs = tr.args[2:end]
    trname = esc(tr.args[1])
    curly = Any[]
    paras = Any[]
    for (ty,v) in zip(typs, GenerateTypeVars{:upcase}())
        push!(curly, Expr(:(<:), esc(v), esc(ty)))  #:($v<:$ty)
        push!(paras, esc(v))
    end
    arg = :(::Type{$trname{$(paras...)}})
    fnhead = :($curmod.trait{$(curly...)}($arg))
    isfnhead = :($curmod.istrait{$(curly...)}($arg))
    quote
        $fnhead = $trname{$(paras...)}
        $isfnhead = true # Add the istrait definition as otherwise
                         # method-caching can be an issue.
    end
end

# Defining a function dispatching on the trait (or not)
# @traitfn f{X,Y;  Tr1{X,Y}}(x::X,y::Y) = ...
# @traitfn f{X,Y; !Tr1{X,Y}}(x::X,y::Y) = ... # which is just sugar for:
# @traitfn f{X,Y; Not{Tr1{X,Y}}}(x::X,y::Y) = ...
let dispatch_cache = Set()  # to ensure that the trait-dispatch function is defined only once per pair
    global traitfn
    function traitfn(tfn)
        # Need
        # f{X,Y}(x::X,Y::Y) = f(trait(Tr1{X,Y}), x, y)
        # f(::False, x, y)= ...
        if tfn.head==:macrocall
            hasmac = true
            mac = tfn.args[1]
            tfn = tfn.args[2]
        else
            hasmac = false
        end
        fhead = tfn.args[1]
        fbody = tfn.args[2]
        fname = fhead.args[1].args[1]
        args0 = fhead.args[2:end]  # version before adding gensyms
        args = insertdummy(fhead.args[2:end])
        typs = fhead.args[1].args[3:end]
        trait_opposite = fhead.args[1].args[2].args[1]
        trait = trait_opposite
        if isnegated(trait)
            trait = trait.args[2]
            val = :(::Type{$curmod.Not{$trait}})
        else
            val = :(::Type{$trait})
            trait_opposite = Expr(:call, :!, trait_opposite)  # generate the opposite
        end
        # Get line info for better backtraces
        ln = findline(tfn)
        if isa(ln, Expr)
            pushloc = Expr(:meta, :push_loc, ln.args[2], fname, ln.args[1])
            poploc = Expr(:meta, :pop_loc)
        else
            pushloc = poploc = nothing
        end
        retsym = gensym()
        if hasmac
            fn = :(@dummy $fname{$(typs...)}($val, $(args...)) = ($pushloc; $retsym = $fbody; $poploc; $retsym))
            fn.args[1] = mac # replace @dummy
        else
            fn = :($fname{$(typs...)}($val, $(args...)) = ($pushloc; $retsym = $fbody; $poploc; $retsym))
        end
        ex = fn
        key = (current_module(), fname, typs, args0, trait_opposite)
        if !(key ∈ dispatch_cache)
            ex = quote
                $fname{$(typs...)}($(args...)) = (Base.@_inline_meta(); $fname($curmod.trait($trait), $(strip_tpara(args)...)))
                $ex
            end
            push!(dispatch_cache, key)
        else
            delete!(dispatch_cache, key) # permits function redefinition if that's what we want
        end
        ex
    end
end

"""
Defines a function dispatching on a trait. Examples:

```julia
@traitfn f{X;  Tr1{X}}(x::X,y) = ...
@traitfn f{X; !Tr1{X}}(x::X,y) = ...

@traitfn f{X,Y;  Tr2{X,Y}}(x::X,y::Y) = ...
@traitfn f{X,Y; !Tr2{X,Y}}(x::X,y::Y) = ...
```

Note that the second example is just syntax sugar for `@traitfn f{X,Y; Not{Tr1{X,Y}}}(x::X,y::Y) = ...`.
"""
macro traitfn(tfn)
    esc(traitfn(tfn))
end

######
## Helpers
######

# true if :(!(Tr{x}))
isnegated(t::Expr) = t.head==:call

# [:(x::X)] -> [:x]
# also takes care of :...
strip_tpara(args::Vector) = Any[strip_tpara(a) for a in args]
strip_tpara(a::Symbol) = a
function strip_tpara(a::Expr)
    if a.head==:(::)
        return a.args[1]
    elseif a.head==:...
        return Expr(:..., strip_tpara(a.args[1]))
    else
        error("Cannot parse argument: $a")
    end
end

# insert dummy: ::X -> gensym()::X
# also takes care of :...
insertdummy(args::Vector) = Any[insertdummy(a) for a in args]
insertdummy(a::Symbol) = a
function insertdummy(a::Expr)
    if a.head==:(::) && length(a.args)==1
        return Expr(:(::), gensym(), a.args[1])
    elseif a.head==:...
        return Expr(:..., insertdummy(a.args[1]))
    else
        return a
    end
end

# generates: X1, X2,... or x1, x2.... (just symbols not actual TypeVar)
type GenerateTypeVars{CASE}
end
Base.start(::GenerateTypeVars) = 1
Base.next(::GenerateTypeVars{:upcase}, state) = (Symbol("X$state"), state+1) # X1,..
Base.next(::GenerateTypeVars{:lcase}, state) = (Symbol("x$state"), state+1)  # x1,...
Base.done(::GenerateTypeVars, state) = false

####
# Annotating the source location
####

function findline(ex::Expr)
    ex.head == :line && return ex
    for a in ex.args
        ret = findline(a)
        isa(ret, Expr) && return ret
    end
    nothing
end
findline(arg) = nothing

####
# Extras
####

include("base-traits.jl")

end # module


text:
const COL_NUM_PARTICLES = :num_particles
const COL_NUM_METAINFERENCE = :num_metainference
const COL_PROPOSAL_NAME = :proposal_name
const COL_AIDE_ESTIMATE = :aide_estimate
const COL_AIDE_STDERR = :aide_stderr
const COL_GOLD_STANDARD_NAME = :gold_standard_name
const PRIOR_PROPOSAL_NAME = "prior"
const OPTIMAL_PROPOSAL_NAME = "optimal"
const EXACT_GOLD_STANDARD = "exact"
const APPROXIMATE_GOLD_STANDARD = "approximate"


text:
using HumanReadableIDs
using Test
using Random
using UUIDs

# Monkey patch uuid4
_rand(T) = rand(T)
_rand(::Type{UUID}) = uuid4()

@testset "proquints" begin
    # Literal test cases
    for (x,str) in [ # all letters for first 6 bit positions
                    (0x0000, "babab"),
                    (0x0001, "babad"),
                    (0x0002, "babaf"),
                    (0x0003, "babag"),
                    (0x0004, "babah"),
                    (0x0005, "babaj"),
                    (0x0006, "babak"),
                    (0x0007, "babal"),
                    (0x0008, "babam"),
                    (0x0009, "baban"),
                    (0x000a, "babap"),
                    (0x000b, "babar"),
                    (0x000c, "babas"),
                    (0x000d, "babat"),
                    (0x000e, "babav"),
                    (0x000f, "babaz"),
                    (0x0010, "babib"),
                    (0x0020, "babob"),
                    (0x0030, "babub"),
                    # Other bits
                    (0x03c0, "bazab"),
                    (0x0c00, "bubab"),
                    (0xf000, "zabab"),
                    # Larger sizes
                    (0xf0000000, "zabab-babab"),
                    (0x000f0000, "babaz-babab"),
                   ]
        @test str == encode_proquint(x)
    end
    @testset "round-trip type tests" for T in [UInt8, UInt16, UInt32,
                                               UInt64, UInt128, UUID]
        for i=1:10
            x = _rand(T)
            y = decode_proquint(T, encode_proquint(x))
            @test y isa T
            @test x == y
        end
    end
    @testset "Exhaustive round-trip test for all 16 bit numbers" begin
        for i=0x0000:0xffff
            @test i == decode_proquint(UInt16, encode_proquint(i))
        end
    end
end


text:
using Base.Test
using MIPVerify: ReLU, mip

@testset "relu.jl" begin

@testset "ReLU" begin
    @testset "Initialize without tightening algorithm" begin
        p = ReLU()
        @test isnull(p.tightening_algorithm)
    end

    @testset "Initialize with tightening algorithm" begin
        p = ReLU(mip)
        @test !isnull(p.tightening_algorithm)
        @test get(p.tightening_algorithm) == mip
    end

    @testset "Base.show" begin
        p = ReLU()
        io = IOBuffer()
        Base.show(io, p)
        @test String(take!(io)) == "ReLU()"
    end
end

end


text:
@testset "Linear algebra" begin
  include("LinearAlgebra/Sparse.jl")
end


text:
## We will start with vectors
v = [1; 2; 3; 4]
##
v = [1 2 3 4]
##
v = [1,2,3,4]
## If we want floats
v = [1.0,2,3,4]
## If you want to convert to Float64
w = map(Float64, [1,2,3,4]) # this will build an array, then convert to Float64
##
v[1]
##
v[3]
##
v[end]
##
v[end-1] # this is one before 3nd...
## This gives an error
v[0]
##
w = [1 2 3 4]
v = [1,2,3,4]
w+v
## Aside: This is an Atom thing
## Command enter shows the result of the current line
1+1
## If we run alt-enter then it runs the entire block
a = 5
a+1
##
v+w' # this should work
##
vec(v+w') # this gives a 1d vector
## This is as expected
v.+w'
## This is interesting ...
A = v.+w
## This gives a 4x4 array, not a 1x4 or 4x1 array of element-wise products.
A[3,4] == v[3] + w[4] # w becomes columns, v gets mapped to rows.
##
function vecsum(x::Vector)
    s = zero(eltype(x)) # this is a variable that is zero in the same
    # element type of as x.
    s = 0.0 # this only works for floats, not BigFloats, etc...
    for i=1:length(x)
        s += x[i]
    end
    return s
end
vecsum(v)
##
vecsum(vec(w))
## A word of warning about vec.
w[3]
##

z = vec(w)
z[3]
##
z[3] = 10
##
w[3]
## If you want to de-couple the memory ... just make a copy!
z = copy(vec(w))
z[3] = 50
w[3]

## Now on to Matrices.
A = [1 2 3; 4 5 6; 7 8 0]
##
b = [0; 1; 2] # this is a 3x1 vector, but really a 1d array

## We can solve linear systems of equations!
x = A\b # this computes x such that A*x = b
## It's always good to check solutions
A*x - b
## Other operations with matrices
using LinearAlgebra # adds lots of linear algebra routines to the Julia Session.
##
rank(A)
##
eigen(A)
##
##
svd(A)
##
qr(A)
##
F = qr(A)
F.Q # this is the matrix you get by doing Gram-Schmit

##
Q = F.Q
Q'*Q
##
Q'*Q - I
##
norm(Q'*Q - I) # this computes the Frobenius norm

##
F = lu(A) # this computes an LU factorization
##
F.L
##
F.U

##
# You can use F to solve systems with A
F\b - A\b
## A comment on Norms (this is for Matlab people)
norm(A) # Julia uses the element-wise matrix norms by default.
# you get this in Matlab with "norm(A,'fro')"
## This is another way to compute the norm.
sqrt(sum(A.^2))
## Matlab's default matrix norm is what Julia calls "opnorm"
# These are operator norms. A slightly different concept.
opnorm(A) #

## Let's chat about Comments!
# Comments start with a "#"
#= You can do multiline Comments
these keep going
until they end =#
1+1

##
s = """ These are not comments,
these are strings """
typeof(s)
## Let's chat about printing!
print(x)
##
println(x)
println("Hi")
##
using Printf
##
@printf("%.18e", x[1])
##
println(x[1])
## Ranges
z = 1:5
## This creates a range! It doesn't actually build a vector, but it acts like it does.
for i=z # this means i will range from 1 to 5
    println(i)
end

## Ranges have types
z = 1.0:5 # this is going to be a FLoat64 Ranges
for i=z # this means i will range from 1 to 5
    println(i)
end
##
for i=1.0:5
    println(i)
end
##
x
##
w = x .+ (1:3) # this will add elements to x

##
r = 1:3
x.+r

##
r = 3:-1:1
r[1]

##
r[end]
##
r
##


text:
module SparseFEM

# assemble the finite-difference laplacian
function fdlaplacian(N)
    # create a 1D laplacian and a sparse identity
    fdl1 = spdiagm((ones(N-1),-2*ones(N),ones(N-1)), [-1,0,1])
    # laplace operator on the full grid
    return kron(speye(N), fdl1) + kron(fdl1, speye(N))
end

# get the list of boundary dof-indices
function get_free(N)
    L = zeros(Int, N, N)
    L[2:N-1, 2:N-1] = 1
    return find(L)
end

# timing of assembly, slice and solve
function perf_sparse_fem(N)
    Ifree = get_free(N)
    # assembly
    A = fdlaplacian(N)
    # boundary condition
    B = A[Ifree, Ifree]
    # solver
    return lufact(B)
end

end # module


text:
using Documenter
using LatexSVG

makedocs(
    modules = [LatexSVG],
    sitename = "LatexSVG.jl",
    format = Documenter.HTML(
        assets = ["assets/svgtheme.css"]
    ),
    doctest = false,
    pages = Any[
        "Introduction" => "index.md",
        "Guide" => [
            "man/installation.md",
            "man/usage.md",
            "man/htmldisplay.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/EHBaozi/LatexSVG.jl.git"
)


text:
using SoftSquishyMatter # load simulation package
using Random # useful if you want to shuffle arrays
cd(dirname(@__FILE__)) # set directory to this file's folder

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize empty simulation with a description for future reference
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation = Simulation()
simulation.descriptor = "Poly-dispersed colloids with Lennard-Jones interactions"

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize hot and cold particles
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
R_small = 1e-6 # radius of particle
γ_trans_small = 6 * pi * 8.9e-4 * R_small # friction coefficient of particle assuming Stokes flow
D_trans_small = 1.38e-23 * 300 / γ_trans_small # diffusivity of particle

R_large = 5e-6 # radius of large particles
γ_trans_large = 6 * pi * 8.9e-4 * R_large
D_trans_large = 1.38e-23 * 300 / γ_trans_large

# initialize small and large particles on triangular lattice
# use the lattice dimensions of the small particles for the simulation dimensions (L_x, L_y)
positions_small, simulation.L_x, simulation.L_y = triangular_lattice(s = 3e-6, M_x = 30, M_y = 20)
positions_large, _, _ = triangular_lattice(s = 20e-6, M_x = 4, M_y = 3)

# remove small particles that overlap with the large ones
remove_overlaps!(positions_small; fixed_positions = positions_large, minimum_distance = 2^(1.0 / 6.0) * (R_small + R_large), period_x = simulation.L_x, period_y = simulation.L_y)

# create the particles with properties
for (x, y) in positions_small
    push!(simulation.particles, Particle(ptype = :small, x = x, y = y, R = R_small, γ_trans = γ_trans_small, D_trans = D_trans_small))
end
for (x, y) in positions_large
    push!(simulation.particles, Particle(ptype = :large, x = x, y = y, R = R_large, γ_trans = γ_trans_large, D_trans = D_trans_large))
end
# particle groups for small and large
# these will be useful for the interactions between small-small, small-large, and large-large
pgroup_small = group_by_type(simulation.particles; ptype = :small)
pgroup_large = group_by_type(simulation.particles; ptype = :large)
pgroup_all = simulation.particles

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize cell lists and Lennard-Jones interactions
# We here use different cell lists for the different particles.  This speeds up
# the simulation compared to if we used one cell list (see comment below).
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cell_list_ss = CellList(particles = pgroup_small, L_x = simulation.L_x, L_y = simulation.L_y, cutoff = 2^(1.0 / 6.0) * 2 * R_small)
cell_list_sl = CellList(particles = pgroup_small, L_x = simulation.L_x, L_y = simulation.L_y, cutoff = 2^(1.0 / 6.0) * (R_small + R_large))
cell_list_ll = CellList(particles = pgroup_large, L_x = simulation.L_x, L_y = simulation.L_y, cutoff = 2^(1.0 / 6.0) * 2 * R_large)
lj_ss = LennardJones(particles = pgroup_small, cell_list = cell_list_ss, ϵ = 1.38e-23 * 300, multithreaded = true)
lj_sl = LennardJones(particles = pgroup_large, cell_list = cell_list_sl, ϵ = 1.38e-23 * 300, multithreaded = false, use_newton_3rd = true)
lj_ll = LennardJones(particles = pgroup_large, cell_list = cell_list_ll, ϵ = 1.38e-23 * 300, multithreaded = false)
append!(simulation.cell_lists, [cell_list_ss, cell_list_sl, cell_list_ll])
append!(simulation.interactions, [lj_ss, lj_sl, lj_ll])
# this works too, but is slower since small particles can use a smaller cutoff
#=cell_list = CellList(particles = pgroup_all, L_x = simulation.L_x, L_y = simulation.L_y, cutoff = 2^(1.0 / 6.0) * 2 * R_large)
lj = LennardJones(particles = pgroup_all, cell_list = cell_list, ϵ = 1.38e-23 * 300, multithreaded = true)
push!(simulation.cell_lists, cell_list)
push!(simulation.interactions, lj)=#

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Initialize Brownian integrator for overdamped dynamics
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation.dt = 0.001
# we apply the same integrator to all the particles
brownian = Brownian(particles = pgroup_all, dt = simulation.dt, multithreaded = true)
push!(simulation.integrators, brownian)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set the number of steps and which particles to save periodically
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation.num_steps = trunc(Int64, 10 / simulation.dt)
simulation.save_interval = trunc(Int64, 1 / simulation.dt)
simulation.particles_to_save = pgroup_all

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Run simulation and save data
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
run_simulation!(simulation, save_as = "out/Example_PolyDispersed_data.out")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# OPTIONAL: Load simulation data and animate the simulation.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
simulation = load_simulation(file = "out/Example_PolyDispersed_data.out")
visualize!(simulation; save_as = (:gif, "frames/Example_PolyDispersed.gif"), fps = 1, particle_colors = Dict(:small => "black", :large => "blue"))


text:
using BetaDispersion, Distances, StatsBase,BenchmarkTools, RCall
using Test
R"library(vegan)"

   x = rand(1000,100)
   y =rand(1:20,1000)
   j = @benchmark dispersion($x,$y,BrayCurtis)
   r = @benchmark R"betadisper(vegdist($x),as.factor($y))"
   @test mean(j.times) < mean(r.times)

   Dj = pairwise(BrayCurtis(),x,dims = 1)
   Dr =R"vegdist($x)"
   dispj = dispersion(Dj,y)
   dispr = R"betadisper($Dr,as.factor($y))"
   j = @benchmark permutest(dispj,9999)
   r = @benchmark R"permutest($dispr,pairwise = TRUE, permutations = 9999)"
   @test mean(j.times) < mean(r.times)


text:
using Fetch
using Documenter

DocMeta.setdocmeta!(Fetch, :DocTestSetup, :(using Fetch); recursive=true)

makedocs(;
    modules=[Fetch],
    authors="JingYu Ning <foldfelis@gmail.com> and contributors",
    repo="https://github.com/foldfelis/Fetch.jl/blob/{commit}{path}#{line}",
    sitename="Fetch.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://foldfelis.github.io/Fetch.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/foldfelis/Fetch.jl",
    devbranch="main"
)


text:
using JuMP, SCIP, DataFrames, CSV

struct Technician
    name::String
    cap::Int64
    depot::String
end

struct Job
    name::String
    priority::Int64
    duration::Int64
    coveredBy::Array{String,1}
end

struct Customer
    name::String
    loc::String
    job::Job
    tStart::Int64
    tEnd::Int64
    tDue::Int64
end

function loadproblem()

    ####################
    # Technicians data #
    ####################

    # Read technicians data
    technicians_file = CSV.read("data/technicians.txt",DataFrame)
    #println(technicians_file)

    technicians = [Technician(row["Name"],
                              row["Capacity"],
                              row["Depot"]) for row in eachrow(technicians_file)]

    #@show technicians
    #return technicians,0

    #################
    # Jobs data #
    #################

    # Read services data
    jobs_file = CSV.read("data/jobs.txt",DataFrame)
    #println(jobs_file)

    # Read technician-services data
    technicians_jobs_file = CSV.read("data/technicians_jobs.txt",DataFrame)
    #println(technicians_jobs_file)

    cobertura_servicos = Dict([row["Job Type"] => [] for row in eachrow(jobs_file)])
    for row in eachrow(technicians_jobs_file)
        job = row["Job"]
        technician = row["Technician"]
        tempo = cobertura_servicos[job]
        push!(tempo, technician)
        cobertura_servicos[job] = tempo
    end

    #@show cobertura_servicos
    #return technician_services_file

    jobs = [Job(row["Job Type"],
                row["Priority"],
                row["Duration (min)"],
                cobertura_servicos[ row["Job Type"] ])
            for row in eachrow(jobs_file)]

    #@show jobs
    #return jobs, 0

    ##################
    # Locations data #
    ##################
    # Read locations data
    locations_file = CSV.read("data/locations.txt",DataFrame)
    #println(locations_file)

    # Build useful data structures
    L = String.(unique(locations_file[:,"Location Start"]))
    locations_End = unique(locations_file[:,"Location End"])
    for l in locations_End
        if !(l in L)
            push!(L,String(l))
        end
    end

    #return locations,0

    dist = Dict([(l, l) => 0 for l in L])
    for row in eachrow(locations_file)
        saida = row["Location Start"]
        chegada = row["Location End"]
        tempo = row["Tempo Viagem"]
        dist[(saida,chegada)] = tempo
        dist[(chegada,saida)] = tempo
    end

    #@show dist
    #return dist,0

    ##################
    # Customers data #
    ##################
    # Read customers data
    customers_file = CSV.read("data/customers.txt",DataFrame)
    #println(customers_file)

    customers = Customer[]
    for row in eachrow(customers_file)
        for job in jobs
            if row["Job Type"] == job.name
                thisCustomer = Customer(row["Name"],
                                        row["Location"],
                                        job,
                                        row["Time window Start"],
                                        row["Time window end"],
                                        row["Due time"])
                push!(customers,thisCustomer)
            end
        end
    end

    #@show customers
    #return customers, 0

    K = [k.name for k in technicians]
    C = [j.name for j in customers]
    J = [j.loc for j in customers]
    D = String.(unique(technicians_file.Depot))
    cap = Dict([k.name => k.cap for k in technicians])
    loc = Dict([j.name => j.loc for j in customers])
    depot = Dict([k.name => k.depot for k in technicians])
    canCover = Dict([j.name => j.job.coveredBy for j in customers])
    dur = Dict([j.name => j.job.duration for j in customers])
    tStart = Dict([j.name => j.tStart for j in customers])
    tEnd = Dict([j.name => j.tEnd for j in customers])
    tDue = Dict([j.name => j.tDue for j in customers])
    priority = Dict([j.name => j.job.priority for j in customers])

    #=
    @show K
    @show C
    @show J
    @show L
    @show D
    @show cap
    @show loc
    @show depot
    @show canCover
    @show dur
    @show tStart
    @show tEnd
    @show tDue
    @show priority
    =#

    return technicians,customers,dist,
           K,C,J,L,D,
           cap,loc,depot,canCover,dur,tStart,tEnd,tDue,priority
end

function experimento()

    technicians,customers,dist,K,C,J,L,D,
    cap,loc,depot,canCover,dur,tStart,tEnd,tDue,priority = loadproblem()

    model = Model(SCIP.Optimizer)

    ######################
    # Decision variables #
    ######################

    # Customer-technician assignment
    @variable(model,x[C,K],Bin)

    # Technician assignment
    @variable(model,u[K],Bin)

    # Edge-route assignment to technician
    @variable(model,y[L,L,K],Bin)

    # Technician cannot leave or return to a depot that is not its base
    for k in technicians
        for d in D
            if k.depot != d
                for i in L
                    JuMP.set_upper_bound(y[i,d,k.name],0)
                    JuMP.set_upper_bound(y[d,i,k.name],0)
                end
            end
        end
    end

    # Start time of service
    @variable(model, 0 <= t[L] <= 600,Int)

    # Lateness of service
    @variable(model,z[C] >= 0,Int)

    # Artificial variables to correct time window upper and lower limits
    @variable(model,xa[C] >= 0,Int)
    @variable(model,xb[C] >= 0,Int)

    # Unfilled jobs
    @variable(model,g[C],Bin)

    ###############
    # Constraints #
    ###############

    # A technician must be assigned to a job, or a gap is declared (1)
    for j in C
        constraint = @constraint(model,sum(x[j,k] for k in canCover[j]) + g[j] == 1)
        JuMP.set_name(constraint, "assignToJob[$j]")
    end

    # At most one technician can be assigned to a job (2)
    for j in C
        constraint = @constraint(model, sum(x[j,k] for k in K) <= 1)
        JuMP.set_name(constraint, "assignOne[$j]")
    end

    print(model)

end

experimento()
println("ok")
#technicians,customers,dist,K,C,J,L,D,
#cap,loc,depot,canCover,dur,tStart,tEnd,tDue,priority = loadproblem()

max_stars_repo_name: LuizHNLorena/workforceILP

max_issues_repo_name: LuizHNLorena/workforceILP

max_forks_repo_name: LuizHNLorena/workforceILP


text:
function test_holdout_1()

    rng = StableRNG(0)
    N = 10
    T = 200
    λvals = [1, 2]
    ρvals = [1, 2]
    weights = [1]
    xDistsTest = [Uniform(100, 100.000000001)]
    νDistsTest = [Uniform(30, 30.000000001)]
    xDistsTrain = [Uniform(10, 500)]
    νDistsTrain = νDistsTest
    noiseDist = Normal(0, 0.01)
    signalModels = [(x, ν) -> exp(-ν / x)]
    kernelgenerator = Λ -> GaussianKernel(Λ)
    (λ, ρ, Ψ) = PERK.holdout(rng, N, T, λvals, ρvals, weights, xDistsTest,
                             νDistsTest, xDistsTrain, νDistsTrain, noiseDist,
                             signalModels, kernelgenerator, showprogress = false)

    return λ == 1 && ρ == 1

end

function test_holdout_2()

    rng = StableRNG(0)
    N = 10
    T = 200
    λvals = [1, 2]
    ρvals = [1, 2]
    weights = [1]
    xDistsTest = [Uniform(100, 100.000000001)]
    xDistsTrain = [Uniform(10, 500)]
    noiseDist = Normal(0, 0.01)
    signalModels = [x -> exp(-30 / x)]
    kernelgenerator = Λ -> GaussianKernel(Λ)
    (λ, ρ, Ψ) = PERK.holdout(rng, N, T, λvals, ρvals, weights, xDistsTest,
                             xDistsTrain, noiseDist, signalModels,
                             kernelgenerator, showprogress = false)

    return λ == 1 && ρ == 1

end

function test_holdout_3()

    rng = StableRNG(0)
    N = 10
    T = 200
    λvals = [1, 2]
    ρvals = [1, 2]
    xDistsTest = Uniform(100, 100.000000001)
    νDistsTest = Uniform(30, 30.000000001)
    xDistsTrain = Uniform(10, 500)
    νDistsTrain = νDistsTest
    noiseDist = Normal(0, 0.01)
    signalModels = [(x, ν) -> exp(-ν / x)]
    kernelgenerator = Λ -> GaussianKernel(Λ)
    (λ, ρ, Ψ) = PERK.holdout(rng, N, T, λvals, ρvals, xDistsTest,
                             νDistsTest, xDistsTrain, νDistsTrain, noiseDist,
                             signalModels, kernelgenerator, showprogress = false)

    return λ == 1 && ρ == 1

end

function test_holdout_4()

    rng = StableRNG(0)
    N = 10
    T = 200
    λvals = [1, 2]
    ρvals = [1, 2]
    xDistsTest = Uniform(100, 100.000000001)
    xDistsTrain = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = x -> [exp(-30 / x), x]
    kernelgenerator = Λ -> GaussianKernel(Λ)
    (λ, ρ, Ψ) = PERK.holdout(rng, N, T, λvals, ρvals, xDistsTest,
                             xDistsTrain, noiseDist, signalModels,
                             kernelgenerator, showprogress = false)

    return λ == 1 && ρ == 1

end

@testset "Holdout" begin

    @test test_holdout_1()
    @test test_holdout_2()
    @test test_holdout_3()
    @test test_holdout_4()

end


text:
# NOTE: Based on https://rdrr.io/cran/sdols/ by David B. Dahl.

"""
Create pairwise allocation matrix
"""
function pairwiseallocmat(Z, W, i)
  J, K = size(Z)
  A = zeros(J, J)
  for r in 1:J
    for c in 1:J
      for k in 1:K
        if Z[r, k] == 1 && Z[c, k] == 1
          A[r, c] += W[i, k]
          A[c, r] += W[i, k]
        end
      end
    end
  end

  return A
end

function estimate_ZWi_index(Zs, Ws, i)
  As = [pairwiseallocmat(z, w, i) for (z, w) in zip(Zs, Ws)]

  Amean = mean(As)
  mse = [mean((A - Amean) .^ 2) for A in As]

  return argmin(mse)
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
include("signalobject.jl")
include("peakobject.jl")
include("config/ecgconfig.jl")
include("utils/utils.jl")

using Wavelets
using PyCall
using Statistics

using .signalobject
using .peakobject
import .ecgconfig as cfg
using .utils
pywt = pyimport("pywt")

function preprocess(signals::Vector{T}, loss_threshold = 6000,
    noise_threshold = 100) where T<:Real

    edge = sec_to_sample(2, cfg.FS, true)
    question_segments = utils.mask_segment(signals, loss_threshold, edge)
    question_segments = utils.connect_segment(question_segments)
    utils.replace_zeros!(signals, question_segments)

    signals::Vector{T} = signals - utils.ma(signals, msec_to_sample(600, cfg.FS, true))

    question_segments2 = utils.mask_segment(signals, noise_threshold, edge)
    question_segments2 = utils.connect_segment(question_segments2)
    utils.replace_zeros!(signals, question_segments2)

    question_segments = utils.union_array_of_section([question_segments, question_segments2])

    signals = utils.butter_lowpass(signals, 25, cfg.FS, 2)
    return signalobject.SignalObject(signals), question_segments

end

function test(signals_object::SignalObject)
    print("OK")
end

function detect_peak(signals_object::SignalObject{T},
    detect_algo::String = "above_mean_fix", enhance_algo::String = "abs_dwt_enhancer",
    return_enhanced::Bool = false) where T

    enhance_signals::Vector{T} = enhancer(signals_object.signals, enhance_algo)

    index_array = peak_detector(enhance_signals, detect_algo)

    if return_enhanced
        return index_array, enhance_signals
    end

    return index_array
end

function detect_r_peak(signals_object::SignalObject{<:Real}, index_array::Vector{<:Integer})

    peakarray = PeakArray([Event(0)])
    peakarray.fs = cfg.FS

    r_spacing = msec_to_sample(cfg.DETECTOR_RPEAK_SPACING_MSEC, peakarray.fs)
    r_search_range = msec_to_sample(50, peakarray.fs)

    n_signals = length(signals_object.signal_mean_absolute)

    for index in index_array
        if index - peakarray[end].idx > r_spacing
            start_idx = max(index - r_search_range, 0)
            end_idx = min(index + r_search_range, n_signals)

            _, max_idx = findmax((@view signals_object.signal_mean_absolute[start_idx:end_idx] ))
            push!(peakarray, QRS(max_idx+ start_idx -1))
        end
    end

    return peakarray
end

# function detect_peak_above_mean_fix2(signals::Vector{T}) where T <:Real
#     n_signals = length(signals)

#     smooth::Vector{T} = utils.ma(signals, msec_to_sample(cfg.DETECTOR_AMF_SMOOTH_MSEC, cfg.FS))

#     stride = 2

#     index_array = Vector{Int}(undef, 0)

#     width_threshold = msec_to_sample( cfg.DETECTOR_AMF_WIDTH_THRESHOLD_MSEC,cfg.FS, false)/stride

#     smooth = abs.(smooth) .+  cfg.DETECTOR_AMF_THRESHOLD

#     starting::Bool = false
#     start::Int = 0

#     for i in 1:stride:n_signals

#         if signals[i] > smooth[i]
#             if !starting
#                 start = i
#                 starting = true
#             end
#             continue
#         elseif i-start > width_threshold && starting
#             window = @view signals[start:i]
#             max_val, max_idx = findmax(window)

#             if max_val - smooth[(start + max_idx -1) ]<0.2*max_val
#                 starting = false
#                 continue
#             end
#             push!(index_array, start + max_idx -1)
#         end
#         starting = false

#     end

#     return index_array
# end

# function detect_peak_above_mean_fix(signals::Vector{T}) where T <:Real
#     n_signals = length(signals)

#     smooth::Vector{T} = utils.ma(signals, msec_to_sample(cfg.DETECTOR_AMF_SMOOTH_MSEC, cfg.FS))

#     stride = 2

#     window = Vector{T}(undef, 0)
#     index_array = Vector{Integer}(undef, 0)

#     width_threshold = msec_to_sample( cfg.DETECTOR_AMF_WIDTH_THRESHOLD_MSEC,cfg.FS, false)/stride

#     smooth = abs.(smooth) .+  cfg.DETECTOR_AMF_THRESHOLD

#     for i in 1:stride:n_signals

#         if signals[i] > smooth[i]
#             push!(window, signals[i])
#             continue
#         elseif length(window) > width_threshold
#             max_val, max_idx = findmax(window)

#             if max_val - smooth[i - length(window) + max_idx]<0.2*max_val
#                 window = Vector{T}(undef, 0)
#                 continue
#             end
#             push!(index_array, i - length(window) + max_idx)
#         end
#         window = Vector{T}(undef, 0)

#     end

#     return index_array
# end

function detect_peak_above_mean_fix(sig::SignalObject{T}) where T <:Real
    signals = sig.signal_enhance
    n_signals = length(signals)

    smooth::Vector{T} = utils.ma(signals, msec_to_sample(cfg.DETECTOR_AMF_SMOOTH_MSEC, cfg.FS))

    stride = 2

    index_array = Vector{Int}(undef, 0)

    width_threshold = msec_to_sample( cfg.DETECTOR_AMF_WIDTH_THRESHOLD_MSEC,cfg.FS, false)/stride

    smooth = abs.(smooth) .+  cfg.DETECTOR_AMF_THRESHOLD

    starting::Bool = false
    start::Int = 0

    for i in 1:stride:n_signals

        if signals[i] > smooth[i]
            if !starting
                start = i
                starting = true
            end
            continue
        elseif i-start > width_threshold && starting
            window = @view signals[start:i]
            max_val, max_idx = findmax(window)
            max_idx *= stride

            if max_val - smooth[(start + max_idx -1) ]<0.2*max_val
                starting = false
                continue
            end
            push!(index_array, start + max_idx -1)
        end
        starting = false

    end

    return index_array
end

function detect_t_second_derivative!(signals::SignalObject{T}, peakarray::PeakArray{PeakT}) where {T <:Real, PeakT <: Integer}

    start_from_r_peak = msec_to_sample(cfg.DETECTOR_T_FROM_R_MSEC ,peakarray.fs)
    n_signal = length(signals)
    interval_ratio = cfg.DETECTOR_T_INTERVAL_RATIO
    signal_mean_absolute::Vector{T} = signals.signal_mean_absolute
    signal_second_derivative::Vector{T} = signals.signal_second_derivative

    for peak in peakarray.array
        if (!(peak isa QRS) || isnothing(peak.next_interval) || !peak.normal) continue end
        end_from_r_peak = convert(PeakT, floor(peak.next_interval*interval_ratio))

        if end_from_r_peak <= start_from_r_peak continue end

        start_idx = min(peak.idx+start_from_r_peak, n_signal)
        end_idx = min(peak.idx+end_from_r_peak, n_signal)

        landmark = @view signal_mean_absolute[start_idx:end_idx]

        if length(landmark) > 2
            t_peak_idx = argmax(landmark) + peak.idx + start_from_r_peak - 2
            peak.t = t_peak_idx
        else
            continue
        end
        t_peak = peak.t
        end_idx = min(convert( PeakT ,floor(0.3*peak.next_interval))+t_peak, end_from_r_peak+peak.idx)

        landmark = @view signal_second_derivative[t_peak: end_idx]
        if length(landmark) >2
            peak.t_end = argmax(landmark) + t_peak -1
        end
    end
end

function detect_qrs_onset_offset_second_derivative!(signals::SignalObject{T}, peakarray::PeakArray{PeakT}) where {T <:Real, PeakT <: Integer}
    search_near_r_long =  utils.msec_to_sample(100, peakarray.fs)
    search_near_r_short =  utils.msec_to_sample(20, peakarray.fs)

    second_dif_signals::Vector{T} =  signals.signal_second_derivative
    n_signals = length(second_dif_signals)

    for peak in peakarray.array
        if (!(peak isa QRS) || isnothing(peak.interval) || isnothing(peak.next_interval)) continue end

        start_idx = max(peak.idx-search_near_r_long, 0)
        end_idx = max(peak.idx-search_near_r_short, 0)

        landmask = @view second_dif_signals[start_idx: end_idx]
        if length(landmask) != 0
            onset_idx = argmax(landmask) + start_idx -1
            peak.onset = onset_idx
        end

        start_idx = min(peak.idx + search_near_r_short, n_signals )
        end_idx = min(peak.idx + search_near_r_long, n_signals)

        landmask = @view second_dif_signals[start_idx: end_idx]
        if length(landmask) != 0
            offset_idx = argmax(landmask) + start_idx-1
            peak.offset = offset_idx
        end

    end
end

# function detect_peak_above_mean_fix_mod(sig::SignalObject{T}) where T <:Real
#     signals = sig.signal_enhance
#     n_signals = length(signals)

#     smooth::Vector{T} = utils.ma(signals, msec_to_sample(cfg.DETECTOR_AMF_SMOOTH_MSEC, cfg.FS))

#     stride = 2

#     window = Vector{T}(undef, 0)
#     index_array = Vector{Integer}(undef, 0)

#     width_threshold = msec_to_sample( cfg.DETECTOR_AMF_WIDTH_THRESHOLD_MSEC,cfg.FS, false)/stride

#     smooth = abs.(smooth) .+  cfg.DETECTOR_AMF_THRESHOLD

#     for i in 1:stride:n_signals

#         if signals[i] > smooth[i]
#             push!(window, signals[i])
#             continue
#         elseif length(window) > width_threshold
#             max_val, max_idx = findmax(window)

#             if max_val - smooth[i - length(window) + max_idx]<0.2*max_val
#                 window = Vector{T}(undef, 0)
#                 continue
#             end
#             push!(index_array, i - length(window) + max_idx)
#         end
#         window = Vector{T}(undef, 0)

#     end

#     return index_array
# end

function abs_swt_enhancer!(signals::SignalObject{<:Real})

    level = 3
    check = 2^level
    left = length(signals.signals)%check
    if left !=0
        signals_en = utils.pad_end(signals.signals, check - left, 0)
    else
        signals_en = signals.signals
    end

    signals_en = abs.(pywt.swt(signals_en, "db3", 3)[1][2])
    signals_en = utils.butter_bandpass(signals_en, 0.01, 8.0, cfg.FS).*20
    signals_en = signals_en - utils.ma(signals_en, msec_to_sample(600, cfg.FS))

    signals.signal_enhance = signals_en

end

function abs_dwt_enhancer!(signals::SignalObject{<:Real})
    level = 3
    check = 2^level
    left = length(signals.signals)%check
    if left !=0
        signals_en = utils.pad_end(signals.signals, check - left, 0)
    else
        signals_en = signals.signals
    end

    signals_en = abs.( dwt(signals_en, wavelet(WT.db3), level))
    signals_en = utils.butter_bandpass(signals_en, 0.01, 8.0, cfg.FS).*20
    signals_en = signals_en - utils.ma(signals_en, msec_to_sample(600, cfg.FS))
    signals.signal_enhance = signals_en
end

function tachy_brady!(peakarray::peakobject.PeakArray,segment_index::Vector{<:Integer} ,algo::String = "time")
    if algo == "time"
        tachy_brady_time!(peakarray, segment_index)
    #TODO
    # elseif algo == "peak"
    #     tachy_brady_pause_peak(peakarray)
    end
end

function tachy_brady_time!(peakarray::peakobject.PeakArray, segment_index::Vector{<:Integer})

    start_idx = 1
    threshold_tachy = cfg.HR_TACHY
    threshold_brady = cfg.HR_BRADY
    fs = peakarray.fs

    for end_idx in segment_index
        peakarray_landmark = peakarray[start_idx: end_idx]
        interval, _ = peakobject.get_interval(peakarray_landmark,false, false, true)
        if length(interval) < 3 continue end

        hr = sample_to_hr(mean(interval), fs)

        for peak in peakarray_landmark.array
            if peak isa QRS peakobject.assigned_hr!(peak, hr, threshold_tachy, threshold_brady) end
        end
        start_idx = end_idx

    end
end

function tachy_brady_peak!(peakarray::peakobject.PeakArray)
end

function ectopic!(peakarray::peakobject.PeakArray, segment_index::Vector{<:Integer}, algo::String = "time")
    if algo == "time"
        ectopic_time!(peakarray, segment_index )
    else
        ectopic_time!(peakarray, segment_index )
    end

end

function ectopic_time!(peakarray::peakobject.PeakArray,segment_index::Vector{<:Integer} )
    start_idx = 1
    ectopic_ratio = cfg.ECTOPIC_RATIO
    qrs_width_sample = msec_to_sample(cfg.ECTOPIC_WIDTH_MSEC, peakarray.fs)

    for end_idx in segment_index

        peakarray_landmark = peakarray[start_idx: end_idx]
        interval_array, _ = peakobject.get_interval(peakarray_landmark,false, false, false)
        if length(interval_array) == 0 continue end
        m_rr = median(interval_array)

        for peak in peakarray_landmark.array
            if !(peak isa QRS) continue end
            if isnothing(peak.interval) continue end

            if peak.interval < m_rr * ectopic_ratio
                pvc_pac!(peak, qrs_width_sample)
            end
        end
        start_idx = end_idx
    end
end

function pvc_pac!(peak::peakobject.QRS, qrs_width_sample)

    if isnothing(peak.onset) || isnothing(peak.offset)
        peakobject.add_diagnosis!(peak, peakobject.PAC)
        return
    end

    qrs_width = peak.offset - peak.onset
    if qrs_width > qrs_width_sample peakobject.add_diagnosis!(peak,peakobject.PVC)
    else peakobject.add_diagnosis!(peak, peakobject.PAC) end
end


text:
# This is a Julia implementation of Empirical Bayes model selection and estimation
# in sparse high-dimensional Gaussian linear regression models as described in the
# paper:
#    Martin, Mess, and Walker, "Empirical Bayes posterior concentration in sparse
#    high-dimensional linear models," arXiv:1406:7718
#

# We need a lasso solution to start from so use GLMNet.
using GLMNet

# Get the lasso solution from a cross-validation run. Since it is only for starting
# the MH chain in empirical bayes below we only run a 5-fold CV.
function lasso_cv(x, y)
  cv = glmnetcv(x, y)
  cv.path.betas[:, indmin(cv.meanloss)]
end

# Prior for complexity/size of a model.
function dcomplex(x, n, p, a, b, calclog = true)
  o = -x * (log(b) + a * log(p)) + log(float(x .<= n))
  return (calclog ? o : exp(o))
end

# Empirical Bayes model selection and estimation Run Metropolis-Hastings chain to sample
# from marginal posterior.
#
#   S, s, i, beta = empirical_bayes_linreg(x, y; options...)
#
# Inputs:
#  x = predictor variables, n*p, with n samples of p variables
#  y = response variables, n*1, with n samples of response
#
# Outputs:
#
# Options:
#  sigma = error variance, if nothing then it is estimated from data
#  alpha = likelihood fraction (0 < alpha < 1)
#  gamma = conditional prior parameter (0 <= gamma, should be close to 0)
#  numSamples = number of Monte Carlo samples
#  burnIn = fraction of additional burn in samples taken first but not saved
#  sampleCoefs = true iff samples of the coefficients should be returned
#  logF = function to calc log of prior for the model size
#
function empirical_bayes_linreg(X, y;
  sigma2 = nothing,
  alpha = 0.999,
  gamma = 0.001,
  numSamples = 10000,
  burnIn = 0.25,
  logF = dcomplex
)

  # Check inputs and prep.
  n, p = size(X)
  @assert size(y, 1) == n
  @assert size(y, 2) == 1

  numBurnin = iround(burnIn * numSamples)

  # We will start the chain from the cross-validated Lasso solution so get that.
  lasso_coefs = lasso_cv(X, y)

  if sigma2 == nothing
    z = y - X * lasso_coefs
    sigma2 = sum(z .^ 2) / max(1, n - sum(S))
  end

  twosigma2 = 2 * sigma2
  v = gamma + alpha / sigma2
  logvhalf = log(v) / 2

  # We will use an array of active coefficients to represent a model, S (selected).
  S = find(lasso_coefs .!= 0.0)

  # Get the log priors so we need not recalc them all the time
  logpriors = map((len) -> logF(len) - lchoose(p, len), 0:n)

  # Save the models and which model was used at each step
  models = Any[]        # Each model is unique here and saved as (model, length, betahat, sse, U)
  modelhash = Dict{Vector{Int64}, Int64}() # Map from a model to its index in models array
  sampledmodels = zeros(Int64, numSamples) # Index of model sampled in each step of the chain

  addmodel(model) = begin
    if haskey(modelhash, model)
      entry = models[modelhash[model]]
    else
      # Do linear regression on predictor variables selected by the model
      betas, residuals, sse = linreg_for_subset_of_columns(X, y, model)

      # Create entry for models array
      i = 1 + length(models) # will be at next free index
      modelhash[model] = i
      s = length(model)
      logposterior = logpriors[s] - alpha * sse / twosigma2 - s * logvhalf
      entry = (i, s, model, betas, sse, logposterior)

      # Save entry in models array
      push!(models, entry)
    end
    return entry
  end

  # Add our first model, the one from the Lasso.
  i, s, S, coefs, sse, logposterior = addmodel(S)

  # Now take the Monte Carlo samples
  for j in 1:(numBurnin + numSamples)

    # Get a new model proposal based on random sampling from current model
    Snew = random_proposal_model(S, n, p)

    # Add the model after calculating its implications/variables.
    inew, snew, Snew, coefs, sse, logposteriornew = addmodel(Snew)

    println("model ", j, ": ", snew, ", ", sse, ", ", logposteriornew, "\n")

    # Metropolis-Hastings accept step: Randomly accept proposal dependent on
    # its posterior compared to the posterior of the current model.
    if rand() <= exp(logposteriornew - logposterior)
      i = inew
      S = Snew
      logposterior = logposteriornew
    end

    # Save sample if we are past the burnin
    if j > numBurnin
      sampledmodels[j - numBurnin] = i
    end

  end

  return models, sampledmodels

end

bigfloat_to_float(bf, digits = 7) = int(bf*(10^digits))/(10^digits)
lchoose(n, m) = bigfloat_to_float(log(binomial(BigInt(n), BigInt(m))))

function linreg_for_subset_of_columns(X, y, model)
  xselected = X[:, model]
  betas = linreg(xselected, y)
  residuals = y .- (betas[1] + xselected * betas[2:end])
  sse = sum( residuals.^2 )
  return betas, residuals, sse
end

# Helper functions for proposal function
delete_at(a, i) = [a[1:(i-1)], a[(i+1):end]]
sample(ary) = ary[rand(1:length(ary))]
add_new_from_range(a, maxval, minval = 1) = push!(copy(a), sample(collect(setdiff(Set(minval:maxval), Set(a)))))

# Proposal function that randomly samples a new model close to given one.
function random_proposal_model(model, maxlen, numvars)
  s = length(model)
  if s == maxlen # If we are at the max length we can only take away variables
    return delete_at(model, rand(1:s))
  elseif s == 1 # If we are at length 1 we can only add
    return sort(add_new_from_range(model, numvars))
  else
    if rand() <= 0.5
      return delete_at(model, rand(1:s))
    else
      return sort(add_new_from_range(model, numvars))
    end
  end
end

# To sample coefficients given a model with a restricted set of
# predictor variables we must first calculate the U as defined
# in the paper.
function calc_U(xselected, model)
  xtxselected = xselected' * xselected
  V = pinv(xtxselected)
  U = chol(V)
  return U
end

sample_coefficients(s, U, sqrtv, betas) = betas + sqrtv * U * randn(s)


text:
import Base: eof, close, position
export BedReader, close, value, eof, advance!, eachposition

type BedReader <: SequenceReader
    bedStream
    done::Bool
    position::Int64
    currRegionStart::Int64
    currRegionEnd::Int64
    currRegionValue::Float64
    contigs::ReferenceContigs
end

function BedReader(bedFileName::ASCIIString, contigs)
    f = open(bedFileName)
    r = BedReader(f, false, 0, 0, 0, 0, contigs)
    advance!(r)
    r
end
function BedReader(f::IO, contigs)
    r = BedReader(f, false, 0, 0, 0, 0, contigs)
    advance!(r)
    r
end
close(reader::BedReader) = close(reader.bedStream)
value(reader::BedReader) = r.currRegionValue
position(reader::BedReader) = reader.position
eof(reader::BedReader) = reader.position == -1

function advance!(r::BedReader)

    while !r.done
        if peek(r.bedStream) == -1 # eof does not work on .gz files (gzip issue?)
            r.done = true
            r.position = -1
            return
        end

        r.position += 1

        # read until we have the right BED line
        while r.position > r.currRegionEnd && !eof(r.bedStream)
            line = readline(r.bedStream)
            parts = split(line, '\t', limit=6)
            r.currRegionValue = 1.0
            if length(parts) >= 5
                r.currRegionValue = parse(Float64, parts[5])
            end
            offset = contig_offset(r.contigs, parts[1])
            r.currRegionStart = offset + parse(Float64, parts[2]) + 1
            r.currRegionEnd = offset + parse(Float64, parts[3])
        end

        # skip blank areas
        if r.position < r.currRegionStart
            r.position = r.currRegionStart
        end

        # If we found a region then
        if r.position <= r.currRegionEnd
            return
        end
    end
end

# here we want to update the reader
eachposition(r::BedReader) = BedReaderIterator(r)
immutable BedReaderIterator
	reader::BedReader
end
Base.start(it::BedReaderIterator) = it.reader.position
Base.done(it::BedReaderIterator, position) = position == -1
function Base.next(it::BedReaderIterator, position)
	pos = it.reader.position
	advance!(it.reader)
	pos,it.reader.position
end


text:
function opisreduced(op::Operation)
  for rdep ∈ reduceddependencies(op)
    rdep ∉ loopdependencies(op) && return true
  end
  false
end
function storeinstr_preprend(op::Operation, vloopsym::Symbol)
    # defaultstoreop = :vstore!
    # defaultstoreop = :vnoaliasstore!
    isvectorized(op) && return Symbol("")
    vloopsym ∉ reduceddependencies(op) && return Symbol("")
    # vectorized is not a loopdep, but is a reduced dep
    opp::Operation = first(parents(op))
    # while vectorized ∉ loopdependencies(opp)
    while ((!isvectorized(opp)) || opisreduced(opp))
        oppold = opp
        for oppp ∈ parents(opp)
            if vloopsym ∈ reduceddependencies(oppp)
                @assert opp !== oppp "More than one parent is a reduction over the vectorized variable."
                opp = oppp
            end
        end
        @assert opp !== oppold "Failed to find any parents "
    end
    reduction_to_scalar(reduction_instruction_class(instruction(opp)))
end

function reduce_expr_u₂(toreduct::Symbol, instr::Instruction, u₂::Int, suffix::Symbol)
    t = Expr(:tuple)
    for u ∈ 0:u₂-1
        push!(t.args, Symbol(toreduct, u, suffix))
    end
    Expr(:call, lv(:reduce_tup), reduce_to_onevecunroll(instr), t)
end
function reduce_expr!(q::Expr, toreduct::Symbol, instr::Instruction, u₁::Int, u₂::Int, isu₁unrolled::Bool, isu₂unrolled::Bool)
    if isu₂unrolled# u₂ != -1
        _toreduct = Symbol(toreduct, 0)
        push!(q.args, Expr(:(=), _toreduct, reduce_expr_u₂(toreduct, instr, u₂, Symbol(""))))
    else#if u₂ == -1
        _toreduct = Symbol(toreduct, '_', u₁)
    # else
        # _toreduct = Symbol(toreduct, 0)
    end
    # @show toreduct, _toreduct, u₁, u₂, isu₁unrolled, isu₂unrolled
    if (u₁ == 1) | (~isu₁unrolled)
        push!(q.args, Expr(:(=), Symbol(toreduct, "##onevec##"), _toreduct))
    else
        push!(q.args, Expr(:(=), Symbol(toreduct, "##onevec##"), Expr(:call, lv(reduction_to_single_vector(instr)), _toreduct)))
        # push!(q.args, :(@show $_toreduct))
        # push!(q.args, Expr(:(=), Symbol(toreduct, "##onevec##"), :(@show $(Expr(:call, lv(reduction_to_single_vector(instr)), _toreduct)))))
    end
    nothing
end

function lower_store_collection!(
    q::Expr, ls::LoopSet, op::Operation, ua::UnrollArgs, mask::Bool, inds_calc_by_ptr_offset::Vector{Bool}
)
    omop = offsetloadcollection(ls)
    batchid, bopind = omop.batchedcollectionmap[identifier(op)]
    collectionid, copind = omop.opidcollectionmap[identifier(op)]
    opidmap = offsetloadcollection(ls).opids[collectionid]
    idsformap = omop.batchedcollections[batchid]

    @unpack u₁, u₁loop, u₁loopsym, u₂loopsym, vloopsym, vloop, u₂max, suffix = ua
    ops = operations(ls)
    # __u₂max = ls.unrollspecification.u₂
    nouter = length(idsformap)

    t = Expr(:tuple)
    for (i,(opid,_)) ∈ enumerate(idsformap)
        opp = first(parents(ops[opidmap[opid]]))

        isu₁, isu₂ = isunrolled_sym(opp, u₁loopsym, u₂loopsym, vloopsym, ls)#, __u₂max)
        u = Core.ifelse(isu₁, u₁, 1)
        mvar = Symbol(variable_name(opp, ifelse(isu₂, suffix, -1)), '_', u)
        # mvar = Symbol(variable_name(_op, suffix), '_', u)
        # mvar = Symbol(variable_name(_op, ifelse(isu₂, suffix, -1)), '_', u)
        # @show mvar, isu₂, isu₂unrolled(opp)
        push!(t.args, mvar)
    end

    offset_dummy_loop = Loop(first(getindices(op)), MaybeKnown(1), MaybeKnown(1024), MaybeKnown(1), Symbol(""), Symbol(""))
    unrollcurl₂ = unrolled_curly(op, nouter, offset_dummy_loop, vloop, mask, 1)
    inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, false, 0, ls)
    falseexpr = Expr(:call, lv(:False));
    aliasexpr = falseexpr;
    # trueexpr = Expr(:call, lv(:True));
    rs = staticexpr(reg_size(ls));
    manualunrollu₁ = if isu₁unrolled(op) && u₁ > 1 # both unrolled
        if isknown(step(u₁loop)) && sum(Base.Fix2(===,u₁loopsym), getindicesonly(op)) == 1
            if first(getindices(op)) === vloopsym
                interleaveval = -nouter
            else
                interleaveval = 0
            end
            unrollcurl₁ = unrolled_curly(op, u₁, ua.u₁loop, vloop, mask, interleaveval)
            inds = Expr(:call, unrollcurl₁, inds)
            false
        else
            true
        end
    else
        false
    end
    uinds = Expr(:call, unrollcurl₂, inds)
    vp = vptr(op)
    sptrsym = sptr!(q, op)
    storeexpr = Expr(:call, lv(:_vstore!), sptrsym, Expr(:call, lv(:VecUnroll), t), uinds)
    # not using `add_memory_mask!(storeexpr, op, ua, mask, ls)` because we checked `isconditionalmemop` earlier in `lower_load_collection!`
    u₁vectorized = u₁loopsym === vloopsym
    if mask# && isvectorized(op))
        if !(manualunrollu₁ & u₁vectorized)
            push!(storeexpr.args, MASKSYMBOL)
        end
    end
    push!(storeexpr.args, falseexpr, aliasexpr, falseexpr, rs)
    if manualunrollu₁
        masklast = mask & u₁vectorized
        gf = GlobalRef(Core,:getfield)
        tv = Vector{Symbol}(undef, length(t.args))
        for i ∈ eachindex(tv)
            s = gensym!(ls, "##tmp##collection##store##")
            tv[i] = s
            push!(q.args, Expr(:(=), s, Expr(:call, gf, t.args[i], 1)))
        end
        # @show u₁, t
        for u ∈ 0:u₁-1
            lastiter = (u+1) == u₁
            storeexpr_tmp = if lastiter
                storeexpr
                (((u+1) == u₁) & masklast) && insert!(storeexpr.args, length(storeexpr.args)-3, MASKSYMBOL) # 3 for falseexpr, aliasexpr, falseexpr
                storeexpr
            else
                copy(storeexpr)
            end
            vut = Expr(:tuple)
            for i ∈ eachindex(tv)
                push!(vut.args, Expr(:call, gf, tv[i], u+1, false))
            end
            storeexpr_tmp.args[3] = Expr(:call, lv(:VecUnroll), vut)
            if u ≠ 0
                storeexpr_tmp.args[4] = Expr(:call, unrollcurl₂, mem_offset_u(op, ua, inds_calc_by_ptr_offset, false, u, ls))
            end
            push!(q.args, storeexpr_tmp)
        end
    else
        push!(q.args, storeexpr)
    end
    nothing
end
gf(s::Symbol, n::Int) = Expr(:call, GlobalRef(Core,:getfield), s, n, false)
function lower_store!(
  q::Expr, ls::LoopSet, op::Operation, ua::UnrollArgs, mask::Bool,
  reductfunc::Symbol = storeinstr_preprend(op, ua.vloop.itersymbol), inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)
)
  @unpack u₁, u₁loopsym, u₂loopsym, vloopsym, vloop, u₂max, suffix = ua
  omop = offsetloadcollection(ls)
  batchid, opind = omop.batchedcollectionmap[identifier(op)]
  if ((batchid ≠ 0) && isvectorized(op)) && (!rejectinterleave(op))
    (opind == 1) && lower_store_collection!(q, ls, op, ua, mask, inds_calc_by_ptr_offset)
    return
  end
  falseexpr = Expr(:call, lv(:False));
  aliasexpr = falseexpr;
  # trueexpr = Expr(:call, lv(:True));
  rs = staticexpr(reg_size(ls));
  opp = first(parents(op))
  if ((opp.instruction.instr === reductfunc) || (opp.instruction.instr === :identity))
    parents_opp = parents(opp)
    opppstate = Base.iterate(parents_opp)
    if opppstate ≢ nothing
      oppp, state = opppstate
      if (Base.iterate(parents_opp, state) === nothing) && isu₂unrolled(op) == isu₂unrolled(oppp)
        opp = oppp
      end
    end
  end
  # __u₂max = ls.unrollspecification.u₂
  isu₁, isu₂ = isunrolled_sym(opp, u₁loopsym, u₂loopsym, vloopsym, ls)#, __u₂max)
  # @show isu₁, isu₂, u₁loopsym, u₂loopsym
  # @show isu₁, isu₂, opp, u₁loopsym, u₂loopsym, vloopsym
  u = isu₁ ? u₁ : 1
  mvar = Symbol(variable_name(opp, ifelse(isu₂, suffix, -1)), '_', u)
  if all(op.ref.loopedindex)
    inds = unrolledindex(op, ua, mask, inds_calc_by_ptr_offset, ls)
    storeexpr = if reductfunc === Symbol("")
      Expr(:call, lv(:_vstore!), sptr(op), mvar, inds)
    else
      Expr(:call, lv(:_vstore!), lv(reductfunc), sptr(op), mvar, inds)
    end
    add_memory_mask!(storeexpr, op, ua, mask, ls)
    push!(storeexpr.args, falseexpr, aliasexpr, falseexpr, rs)
    push!(q.args, storeexpr)
  else
    parents_op = parents(op)
    data_u₁ = isu₁ & (u₁ > 1)

    indices_u₁ = data_u₁
    if !data_u₁ & (length(parents_op) > 1)
      indices_u₁ = first(isunrolled_sym(op, u₁loopsym, u₂loopsym, vloopsym, ls))
    end
    if indices_u₁
      mvard = Symbol(mvar, "##data##")
      # isu₁ &&
      data_u₁ && push!(q.args, Expr(:(=), mvard, Expr(:call, lv(:data), mvar)))
      sptrsym = sptr!(q, op)
      for u ∈ 1:u₁
        inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, true, u-1, ls)
        # @show isu₁unrolled(opp), opp
        storeexpr = if data_u₁
          if reductfunc === Symbol("")
            Expr(:call, lv(:_vstore!), sptrsym, gf(mvard,u), inds)
          else
            Expr(:call, lv(:_vstore!), lv(reductfunc), sptrsym, gf(mvard,u), inds)
          end
        elseif reductfunc === Symbol("")
          Expr(:call, lv(:_vstore!), sptrsym, mvar, inds)
        else
          Expr(:call, lv(:_vstore!), lv(reductfunc), sptrsym, mvar, inds)
        end
        domask = mask && (isvectorized(op) & ((u == u₁) | (vloopsym !== u₁loopsym)))
        add_memory_mask!(storeexpr, op, ua, domask, ls)# & ((u == u₁) | isvectorized(op)))
        push!(storeexpr.args, falseexpr, aliasexpr, falseexpr, rs)
        push!(q.args, storeexpr)
      end
    else
      inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, true, 0, ls)
      storeexpr = if reductfunc === Symbol("")
        Expr(:call, lv(:_vstore!), sptr(op), mvar, inds)
      else
        Expr(:call, lv(:_vstore!), lv(reductfunc), sptr(op), mvar, inds)
      end
      add_memory_mask!(storeexpr, op, ua, mask, ls)
      push!(storeexpr.args, falseexpr, aliasexpr, falseexpr, rs)
      push!(q.args, storeexpr)
    end
  end
  nothing
end

function lower_tiled_store!(
    blockq::Expr, opsv1::Vector{Operation}, opsv2::Vector{Operation}, ls::LoopSet, unrollsyms::UnrollSymbols, u₁::Int, u₂::Int, mask::Bool
)
    ua = UnrollArgs(ls, u₁, unrollsyms, u₂, 0)
    for opsv ∈ (opsv1, opsv2)
        for op ∈ opsv
            lower_tiled_store!(blockq, op, ls, ua, u₁, u₂, mask)
        end
    end
end

function donot_tile_store(ls::LoopSet, op::Operation, vloop::Loop, reductfunc::Symbol, u₂::Int)
    ((!((reductfunc === Symbol("")) && all(op.ref.loopedindex))) || (u₂ ≤ 1) || isconditionalmemop(op)) && return true
    rejectcurly(op) && return true
    omop = offsetloadcollection(ls)
    batchid, opind = omop.batchedcollectionmap[identifier(op)]
    return ((batchid ≠ 0) && isvectorized(op)) && (!rejectinterleave(op))
end

# VectorizationBase implements optimizations for certain grouped stores
# thus we group stores together here to allow for these possibilities.
# (In particular, it tries to replace scatters with shuffles when there are groups
#   of stores offset from one another.)
function lower_tiled_store!(blockq::Expr, op::Operation, ls::LoopSet, ua::UnrollArgs, u₁::Int, u₂::Int, mask::Bool)
    @unpack u₁loopsym, u₂loopsym, vloopsym, u₁loop, u₂loop, vloop = ua
    reductfunc = storeinstr_preprend(op, vloopsym)
    inds_calc_by_ptr_offset = indices_calculated_by_pointer_offsets(ls, op.ref)

    if donot_tile_store(ls, op, vloop, reductfunc, u₂)
        # If we have a reductfunc, we're using a reducing store instead of a contiuguous or shuffle store anyway
        # so no benefit to being able to handle that case here, vs just calling the default `lower_store!` method
        @unpack u₁, u₂max = ua
        for t ∈ 0:u₂-1
            unrollargs = UnrollArgs(u₁loop, u₂loop, vloop, u₁, u₂max, t)
            lower_store!(blockq, ls, op, unrollargs, mask, reductfunc, inds_calc_by_ptr_offset)
        end
        return
    end
    opp = first(parents(op))
    if (opp.instruction.instr === reductfunc) && isone(length(parents(opp)))
        throw("Operation $opp's instruction is $reductfunc, shouldn't be able to reach here.")
        # opp = only(parents(opp))
    end
    isu₁, isu₂ = isunrolled_sym(opp, u₁loopsym, u₂loopsym, vloopsym, ls)#, u₂)
    # It's reasonable for this to be `!isu₁`
    u = Core.ifelse(isu₁, u₁, 1)
    tup = Expr(:tuple)
    for t ∈ 0:u₂-1
        push!(tup.args, Symbol(variable_name(opp, ifelse(isu₂, t, -1)), '_', u))
    end
    vut = Expr(:call, lv(:VecUnroll), tup) # `VecUnroll` of `VecUnroll`s
    inds = mem_offset_u(op, ua, inds_calc_by_ptr_offset, false, 0, ls)
    unrollcurl₂ = unrolled_curly(op, u₂, u₂loop, vloop, mask)
    falseexpr = Expr(:call, lv(:False));
    aliasexpr = falseexpr;
    # trueexpr = Expr(:call, lv(:True));
    rs = staticexpr(reg_size(ls));
    if isu₁ && u₁ > 1 # both unrolled
        unrollcurl₁ = unrolled_curly(op, u₁, u₁loop, vloop, mask)
        inds = Expr(:call, unrollcurl₁, inds)
    end
    uinds = Expr(:call, unrollcurl₂, inds)
    storeexpr = Expr(:call, lv(:_vstore!), sptr(op), vut, uinds)
    if mask && isvectorized(op)
        # add_memory_mask!(storeexpr, op, ua, mask, ls)
        # we checked for `isconditionalmemop` earlier, so we skip this check
        # and just directly take the branch in `add_memory_mask!`
        push!(storeexpr.args, MASKSYMBOL)
    end
    push!(storeexpr.args, falseexpr, aliasexpr, falseexpr, rs)
    push!(blockq.args, storeexpr)
    nothing
end


text:
"""
    analytic_histogram(d::Distributions.UnivariateDistribution, G::LinRange; normalize=true)

Compute the probability mass of `d` in equally-sized intervals with midpoints given by `G`.

If `normalize`, then the probabilities are re-scaled such that their sum is 1.
"""
function analytic_histogram(d::Distributions.UnivariateDistribution,
                            G::LinRange;
                            normalize=true)
    cdf(x) = Distributions.cdf(d, x)

    mps = get_midpoints(G)
    hist = zeros(length(G))
    for i in eachindex(hist)
        hist[i] = cdf(mps[i + 1]) - cdf(mps[i])
    end

    if normalize
        s = sum(hist)
        hist .= hist ./ s
    end

    return hist
end

# aux for analytic_histogram
function get_midpoints(G::LinRange)
    half_step = step(G) / 2
    return LinRange(first(G) - half_step,
                    last(G) + half_step,
                    length(G) + 1)
end


text:
# Find optimal bin width for a histogram using the method described in:
#  Shimazaki H. and Shinomoto S., A method for selecting the bin size of a
#  time histogram Neural Computation (2007) Vol. 19(6), 1503-1527.

# We need the biased variance and not the unbiased which is implemented in Julia.
function biased_var(x)
  n = length(x)
  var(x) * (n - 1) / n
end

# Assumes the data is sorted.
function count_events_per_bin(delta, sorteddata)
  minv = sorteddata[1]
  numbins = convert(Int64, floor((sorteddata[end] - minv) / delta)) + 1
  counts = zeros(Int64, numbins)

  current_bin = 1
  next_bin_start = minv + delta

  for(data in sorteddata)
    if data >= next_bin_start
      next_bin_start += delta
      current_bin += 1
    end
    counts[current_bin] += 1
  end

  counts
end

# The function to be optimized is a function taking the bin-width, delta, as
# the input, but it also needs the data to calc C(delta).
function c_objective(delta, data)
  bincounts = count_events_per_bin(delta, data)
  k = mean(bincounts)
  v = biased_var(bincounts)
  (2*k - v) / delta
end

using BlackBoxOptim

function find_optimal_bin_width(data; minwidth = nothing, maxwidth = nothing)
  sorted = sort(data)
  minvalue = sorted[1]
  maxvalue = sorted[end]
  range = maxvalue - minvalue

  min_width = (minwidth == nothing) ? range/200 : minwidth
  max_width = (maxwidth == nothing) ? range/2   : maxwidth

  best, fitness = bboptimize((delta) -> c_objective(delta[1], sorted), (min_width, max_width);
    dimensions = 1, iterations = 1e4)

  best[1], count_events_per_bin(best[1], sorted)
end


text:
"""
    AbstractDiagnosticConsumer <: Any
Supretype for DiagnosticConsumers.
"""
abstract type AbstractDiagnosticConsumer end

"""
    struct DiagnosticConsumer <: AbstractDiagnosticConsumer
"""
struct DiagnosticConsumer <: AbstractDiagnosticConsumer
    ptr::CXDiagnosticConsumer
end

Base.unsafe_convert(::Type{CXDiagnosticConsumer}, x::DiagnosticConsumer) = x.ptr
Base.cconvert(::Type{CXDiagnosticConsumer}, x::DiagnosticConsumer) = x

"""
    struct IgnoringDiagConsumer <: AbstractDiagnosticConsumer
Hold a pointer to a `clang::IgnoringDiagConsumer` object.
"""
struct IgnoringDiagConsumer <: AbstractDiagnosticConsumer
    ptr::CXDiagnosticConsumer
end

Base.unsafe_convert(::Type{CXDiagnosticConsumer}, x::IgnoringDiagConsumer) = x.ptr
Base.cconvert(::Type{CXDiagnosticConsumer}, x::IgnoringDiagConsumer) = x

"""
    struct DiagnosticsEngine <: Any
Hold a pointer to a `clang::DiagnosticsEngine` object.
"""
struct DiagnosticsEngine
    ptr::CXDiagnosticsEngine
end

Base.unsafe_convert(::Type{CXDiagnosticsEngine}, x::DiagnosticsEngine) = x.ptr
Base.cconvert(::Type{CXDiagnosticsEngine}, x::DiagnosticsEngine) = x


text:
using Plots
using JLD

plotly()

N=100
for f in ARGS
    ext = last(split(f, '.'))
    if ext == "jld"
        rewards = load(f, "rdict")["despot"]
    elseif ext == "txt"
        rewards = readdlm(f)
    else
        warn("Unrecognized file type: $f")
    end
    histogram!(rewards[1:N], title=f)
    println(f)
    println("rewards < -19: $(find(rewards.<-19.0))")
end

gui()


text:
module Indexer

using DocStringExtensions: METHODLIST
using PkgUtility: read_nc, size_nc

# export public functions
export lat_ind, lon_ind, read_LUT

# functions to calculate the lat and lon index
"""
    lat_ind(lat::Number; res::Number=1)

Round the latitude and return the index in a matrix, given
- `lat` Latitude
- `res` Resolution in latitude

---
# Examples
```julia
ilat = lat_ind(0.3);
ilat = lat_ind(0.3; res=0.5);
```
"""
function lat_ind(lat::Number; res::Number=1)
    @assert -90 <= lat <= 90;

    return Int(fld(lat + 90, res)) + 1
end

"""
    lon_ind(lon::Number; res::Number=1)

Round the longitude and return the index in a matrix, given
- `lon` Longitude
- `res` Resolution in longitude

---
# Examples
```julia
ilon = lon_ind(90.3);
ilon = lon_ind(90.3; res=0.5);
```
"""
function lon_ind(lon::Number; res::Number=1)
    @assert -180 <= lon <= 180;

    return Int(fld(lon + 180, res)) + 1
end

# read data and std from given NetCDF file
"""
This function reads look-up-table (LUT) for whole dataset or selected pieces. Supported methods are

$(METHODLIST)

"""
function read_LUT end

"""
    read_LUT(fn::String, FT::DataType = Float32)

Read the entire look-up-table from collection, given
- `fn` Path to the target file
- `FT` Float number type, default is Float32

---
# Examples
```julia
read_LUT(query_collection(vcmax_collection()));
read_LUT(query_collection(vcmax_collection()), Float64);
```
"""
read_LUT(fn::String, FT::DataType = Float32) = (
    @assert isfile(fn);

    return read_nc(FT, fn, "data"), read_nc(FT, fn, "std")
);

"""
    read_LUT(fn::String, cyc::Int, FT::DataType = Float32)

Read the entire look-up-table from collection, given
- `fn` Path to the target file
- `cyc` Cycle number, such as 8 for data in Augest in 1 `1M` resolution dataset
- `FT` Float number type, default is Float32

---
# Examples
```julia
read_LUT(query_collection(gpp_collection()), 8);
read_LUT(query_collection(gpp_collection()), 8, Float64);
```
"""
read_LUT(fn::String, cyc::Int, FT::DataType = Float32) = (
    @assert isfile(fn);
    @assert size_nc(fn, "data")[1] == 3;

    return read_nc(FT, fn, "data", cyc), read_nc(FT, fn, "std", cyc)
);

"""
    read_LUT(fn::String, lat::Number, lon::Number, res::Number, FT::DataType = Float32; interpolation::Bool = false)

Read the selected part of a look-up-table from collection, given
- `fn` Path to the target file
- `lat` Latitude in `°`
- `lon` Longitude in `°`
- `res` Spatial resolution in `°`
- `FT` Float number type, default is Float32
- `interpolation` If true, interpolate the dataset

---
# Examples
```julia
read_LUT(query_collection(vcmax_collection()), 30, 116, 0.5);
read_LUT(query_collection(vcmax_collection()), 30, 116, 0.5, Float64);
read_LUT(query_collection(vcmax_collection()), 30, 116, 0.5, Float64; interpolation=true);
```
"""
read_LUT(fn::String, lat::Number, lon::Number, res::Number, FT::DataType = Float32; interpolation::Bool = false) = (
    @assert isfile(fn);

    # if not at interpolation mode
    _ilat = lat_ind(lat; res=res);
    _ilon = lon_ind(lon; res=res);
    _raw_dat = read_nc(FT, fn, "data", _ilon, _ilat);
    _raw_std = read_nc(FT, fn, "std" , _ilon, _ilat);

    if !interpolation
        return _raw_dat, _raw_std
    end;

    # if at interpolation mode
    _nlat = Int(180 / res);
    _nlon = Int(360 / res);

    # locate the south and north lines of the target pixel
    _ilat_s = Int(fld(lat + 90 - res/2, res) + 1);
    _ilat_n = _ilat_s + 1;
    _dlat_s = lat - (_ilat_s - 1 + 1/2) * res + 90;
    _dlat_n = (_ilat_n - 1 + 1/2) * res - 90 - lat;

    _ilat_s = max(_ilat_s, 1);
    _ilat_n = min(_ilat_n, _nlat);

    # locate the west and east lines of the target pixel
    _ilon_w = Int(fld(lon + 180 - res/2, res) + 1);
    _ilon_e = _ilon_w + 1;
    _dlon_w = lon - (_ilon_w - 1 + 1/2) * res + 180;
    _dlon_e = (_ilon_e - 1 + 1/2) * res - 180 - lon;

    if _ilon_w < 1 _ilon_w = _nlon end;
    if _ilon_e > _nlon _ilon_e = 1 end;

    # interpolate the value
    _val_s = _dlon_w ./ res .* read_nc(FT, fn, "data", _ilon_e, _ilat_s) .+ _dlon_e ./ res .* read_nc(FT, fn, "data", _ilon_w, _ilat_s);
    _val_n = _dlon_w ./ res .* read_nc(FT, fn, "data", _ilon_e, _ilat_n) .+ _dlon_e ./ res .* read_nc(FT, fn, "data", _ilon_w, _ilat_n);
    _val_i = _dlat_s ./ res .* _val_n .+ _dlat_n ./ res .* _val_s;

    # use non-interpolated value if the interpolated one is NaN
    if typeof(_val_i) <: Number
        if isnan(_val_i)
            _val_i = _raw_dat;
        end;
    else
        for _i in eachindex(_val_i)
            if isnan(_val_i[_i])
                _val_i[_i] = _raw_dat[_i];
            end;
        end;
    end;

    # interpolate the standard deviation
    _std_s = _dlon_w ./ res .* read_nc(FT, fn, "std", _ilon_e, _ilat_s) .+ _dlon_e ./ res .* read_nc(FT, fn, "std", _ilon_w, _ilat_s);
    _std_n = _dlon_w ./ res .* read_nc(FT, fn, "std", _ilon_e, _ilat_n) .+ _dlon_e ./ res .* read_nc(FT, fn, "std", _ilon_w, _ilat_n);
    _std_i = _dlat_s ./ res .* _std_n .+ _dlat_n ./ res .* _std_s;

    # use non-interpolated value if the interpolated one is NaN
    if typeof(_std_i) <: Number
        if isnan(_std_i)
            _std_i = _raw_std;
        end;
    else
        for _i in eachindex(_std_i)
            if isnan(_std_i[_i])
                _std_i[_i] = _raw_std[_i];
            end;
        end;
    end;

    return _val_i, _std_i
);

"""
    read_LUT(fn::String, lat::Number, lon::Number, FT::DataType = Float32; interpolation::Bool = false)

Read the selected part of a look-up-table from collection, given
- `fn` Path to the target file
- `lat` Latitude in `°`
- `lon` Longitude in `°`
- `FT` Float number type, default is Float32
- `interpolation` If true, interpolate the dataset

---
# Examples
```julia
read_LUT(query_collection(vcmax_collection()), 30, 116);
read_LUT(query_collection(vcmax_collection()), 30, 116, Float64);
read_LUT(query_collection(vcmax_collection()), 30, 116, Float64; interpolation=true);
```
"""
read_LUT(fn::String, lat::Number, lon::Number, FT::DataType = Float32; interpolation::Bool = false) = (
    @assert isfile(fn);

    (_,_sizes) = size_nc(fn, "lat");
    _res = 180 / _sizes[1];

    return read_LUT(fn, lat, lon, _res, FT; interpolation=interpolation)
);

"""
    read_LUT(fn::String, lat::Number, lon::Number, cyc::Int, res::Number, FT::DataType = Float32; interpolation::Bool = false)

Read the selected part of a look-up-table from collection, given
- `fn` Path to the target file
- `lat` Latitude in `°`
- `lon` Longitude in `°`
- `cyc` Cycle number, such as 8 for data in Augest in 1 `1M` resolution dataset
- `res` Spatial resolution in `°`
- `FT` Float number type, default is Float32
- `interpolation` If true, interpolate the dataset

---
# Examples
```julia
read_LUT(query_collection(gpp_collection()), 30, 116, 8, 0.5);
read_LUT(query_collection(gpp_collection()), 30, 116, 8, 0.5, Float64);
read_LUT(query_collection(gpp_collection()), 30, 116, 8, 0.5, Float64; interpolation=true);
```
"""
read_LUT(fn::String, lat::Number, lon::Number, cyc::Int, res::Number, FT::DataType = Float32; interpolation::Bool = false) = (
    @assert isfile(fn);

    # if not at interpolation mode
    _ilat = lat_ind(lat; res=res);
    _ilon = lon_ind(lon; res=res);
    _raw_dat = read_nc(FT, fn, "data", _ilon, _ilat, cyc);
    _raw_std = read_nc(FT, fn, "std" , _ilon, _ilat, cyc);

    if !interpolation
        return _raw_dat, _raw_std
    end;

    # if at interpolation mode
    _nlat = Int(180 / res);
    _nlon = Int(360 / res);

    # locate the south and north lines of the target pixel
    _ilat_s = Int(fld(lat + 90 - res/2, res) + 1);
    _ilat_n = _ilat_s + 1;
    _dlat_s = lat - (_ilat_s - 1 + 1/2) * res + 90;
    _dlat_n = (_ilat_n - 1 + 1/2) * res - 90 - lat;

    _ilat_s = max(_ilat_s, 1);
    _ilat_n = min(_ilat_n, _nlat);

    # locate the west and east lines of the target pixel
    _ilon_w = Int(fld(lon + 180 - res/2, res) + 1);
    _ilon_e = _ilon_w + 1;
    _dlon_w = lon - (_ilon_w - 1 + 1/2) * res + 180;
    _dlon_e = (_ilon_e - 1 + 1/2) * res - 180 - lon;

    if _ilon_w < 1 _ilon_w = _nlon end;
    if _ilon_e > _nlon _ilon_e = 1 end;

    # interpolate the value
    _val_s = _dlon_w ./ res .* read_nc(FT, fn, "data", _ilon_e, _ilat_s, cyc) .+ _dlon_e ./ res .* read_nc(FT, fn, "data", _ilon_w, _ilat_s, cyc);
    _val_n = _dlon_w ./ res .* read_nc(FT, fn, "data", _ilon_e, _ilat_n, cyc) .+ _dlon_e ./ res .* read_nc(FT, fn, "data", _ilon_w, _ilat_n, cyc);
    _val_i = _dlat_s ./ res .* _val_n .+ _dlat_n ./ res .* _val_s;

    # use non-interpolated value if the interpolated one is NaN
    if typeof(_val_i) <: Number
        if isnan(_val_i)
            _val_i = _raw_dat;
        end;
    else
        for _i in eachindex(_val_i)
            if isnan(_val_i[_i])
                _val_i[_i] = _raw_dat[_i];
            end;
        end;
    end;

    # interpolate the standard deviation
    _std_s = _dlon_w ./ res .* read_nc(FT, fn, "std", _ilon_e, _ilat_s, cyc) .+ _dlon_e ./ res .* read_nc(FT, fn, "std", _ilon_w, _ilat_s, cyc);
    _std_n = _dlon_w ./ res .* read_nc(FT, fn, "std", _ilon_e, _ilat_n, cyc) .+ _dlon_e ./ res .* read_nc(FT, fn, "std", _ilon_w, _ilat_n, cyc);
    _std_i = _dlat_s ./ res .* _std_n .+ _dlat_n ./ res .* _std_s;

    # use non-interpolated value if the interpolated one is NaN
    if typeof(_std_i) <: Number
        if isnan(_std_i)
            _std_i = _raw_std;
        end;
    else
        for _i in eachindex(_std_i)
            if isnan(_std_i[_i])
                _std_i[_i] = _raw_std[_i];
            end;
        end;
    end;

    return _val_i, _std_i
);

"""
    read_LUT(fn::String, lat::Number, lon::Number, cyc::Int, FT::DataType = Float32; interpolation::Bool = false)

Read the selected part of a look-up-table from collection, given
- `fn` Path to the target file
- `lat` Latitude in `°`
- `lon` Longitude in `°`
- `cyc` Cycle number, such as 8 for data in Augest in 1 `1M` resolution dataset
- `FT` Float number type, default is Float32
- `interpolation` If true, interpolate the dataset

---
# Examples
```julia
read_LUT(query_collection(gpp_collection()), 30, 116, 8);
read_LUT(query_collection(gpp_collection()), 30, 116, 8, Float64);
read_LUT(query_collection(gpp_collection()), 30, 116, 8, Float64; interpolation=true);
```
"""
read_LUT(fn::String, lat::Number, lon::Number, cyc::Int, FT::DataType = Float32; interpolation::Bool = false) = (
    @assert isfile(fn);

    (_,_sizes) = size_nc(fn, "lat");
    _res = 180 / _sizes[1];

    return read_LUT(fn, lat, lon, cyc, _res, FT; interpolation=interpolation)
);

end


text:
module VortexHelperFloorBooster
using ..Ahorn, Maple

@mapdef Entity "VortexHelper/FloorBooster" FloorBooster(
			x::Integer, y::Integer,
			width::Integer=8, height::Integer=8,
			left::Bool = false,
			speed::Integer = 110,
			iceMode::Bool = false,
			noRefillOnIce::Bool = true,
			notAttached::Bool = false)

const placements = Ahorn.PlacementDict(
    "Floor Booster (Right) (Vortex Helper)" => Ahorn.EntityPlacement(
        FloorBooster,
        "rectangle",
        Dict{String, Any}(
            "left" => false
        )
    ),
    "Floor Booster (Left) (Vortex Helper)" => Ahorn.EntityPlacement(
        FloorBooster,
        "rectangle",
        Dict{String, Any}(
            "left" => true
        )
    )
)

Ahorn.minimumSize(entity::FloorBooster) = 8, 8
Ahorn.resizable(entity::FloorBooster) = true, false

function Ahorn.selection(entity::FloorBooster)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FloorBooster, room::Maple.Room)
    left = get(entity.data, "left", false)
	sprite = 	get(entity.data, "iceMode", false) ? "ice" : "fire"

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))
    tileWidth = div(width, 8)

    if left
        for i in 2:tileWidth - 1
            Ahorn.drawImage(ctx, "objects/VortexHelper/floorBooster/" * sprite * "Mid00", (i - 1) * 8, 0)
        end

        Ahorn.drawImage(ctx, "objects/VortexHelper/floorBooster/" * sprite * "Left00", 0, 0)
        Ahorn.drawImage(ctx, "objects/VortexHelper/floorBooster/" * sprite * "Right00", (tileWidth - 1) * 8, 0)

    else
        Ahorn.Cairo.save(ctx)
        Ahorn.scale(ctx, -1, 1)

        for i in 2:tileWidth - 1
            Ahorn.drawImage(ctx, "objects/VortexHelper/floorBooster/" * sprite * "Mid00", i * -8, 0)
        end

        Ahorn.drawImage(ctx, "objects/VortexHelper/floorBooster/" * sprite * "Right00", -8, 0)
        Ahorn.drawImage(ctx, "objects/VortexHelper/floorBooster/" * sprite * "Left00", tileWidth * -8, 0)

        Ahorn.restore(ctx)
    end
end

end

max_stars_repo_name: catapillie/VortexHelper

max_issues_repo_name: catapillie/VortexHelper

max_forks_repo_name: catapillie/VortexHelper


text:
model = Dynamics.pendulum[1]

# Objective
n,m = model.n, model.m
x0 = [0; 0.]
xf = [pi; 0] # (ie, swing up)
Q = 1e-3*Diagonal(I,n)
Qf = 100.0*Diagonal(I,n)
R = 1e-2*Diagonal(I,m)
tf = 5.
obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)

N = 41
U0 = ones(m,N)

# Test original problem
solver = Solver(model, obj_uncon, N=N)
res, stats = solve(solver, U0)
@test norm(res.X[N]-xf) < 1e-3

# Convert LQRObjective to Generic Objective
ell(x,u) = 0.5*x'Q*x + 0.5*u'R*u - xf'Q*x
ell(x) = 0.5*x'Qf*x - xf'Qf*x

costfun = GenericCost(ell,ell,n,m)
obj_generic = UnconstrainedObjective(costfun,tf,x0,xf)
solver_g = Solver(model, obj_generic, N=N)
res_g, stats_g = solve(solver_g, U0)
@test norm(res_g.X - res.X) == 0

# Try a non-linear cost function
ell2(x,u) = cos(x[1]) + u'R*u + x[2]^2*1e-3
ell2(x) = cos(x[1]) + x[2]^2*1e-1

costfun2 = GenericCost(ell2,ell2,n,m)
obj2 = UnconstrainedObjective(costfun2,tf,x0,xf)
solver_g2 = Solver(model, obj2, N=N)
res_g2, stats_g2 = solve(solver_g2, U0*0)
@test norm(res_g2.X[N] - xf) < 1e-6

# Constrained (non-zero initialization)
x_min = [-2*pi,-Inf]
x_max = [2*pi,Inf]
obj2_con = ConstrainedObjective(obj2,x_min=x_min,x_max=x_max)
solver_g2 = Solver(model, obj2_con, N=N)
res_g2_con, stats_g2_con = solve(solver_g2, U0)
@test norm(res_g2_con.X[N] - xf) < 1e-3


text:
facts("operations") do
    const mock_db = DB(joinpath(dirname(@__FILE__), "mock_db"))

    context("getindex") do
        @fact mock_db['n', "'hood"].word --> "'hood"
        @fact mock_db["'hood", 'n'].word --> "'hood"
    end

    context("synsets") do
        ss = synsets(mock_db, mock_db['n', "'hood"])

        @fact length(ss) --> 1
        @fact ss[1].gloss --> "(slang) a neighborhood"
    end

    context("expanded_hypernyms") do
        ss = synsets(mock_db, mock_db['n', "'hood"])[1]
        @fact length(expanded_hypernyms(mock_db, ss)) --> 8
    end

    context("antonyms") do
        ss = synsets(mock_db, mock_db['n', "'hood"])[1]
        @fact length(antonyms(mock_db, ss)) --> 0
    end
end