text:
using Test
using Metal

@testset "MTL" begin

@testset "devices" begin

devs = devices()
@test length(devs) > 0

dev = first(devs)
@test dev == devs[1]

if length(devs) > 1
    @test dev != devs[2]
end

compact_str = sprint(io->show(io, dev))
full_str = sprint(io->show(io, MIME"text/plain"(), dev))

@test dev.name isa String
@test dev.isLowPower isa Bool
@test dev.isRemovable isa Bool
@test dev.hasUnifiedMemory isa Bool
@test dev.registryID isa Integer
@test dev.maxTransferRate isa Integer

@test dev.recommendedMaxWorkingSetSize isa Integer
@test dev.maxThreadgroupMemoryLength isa Integer
@test dev.maxThreadsPerThreadgroup isa MTL.MtSize
@test dev.maxBufferLength isa Integer

@test dev.currentAllocatedSize isa Integer

end

@testset "compile options" begin

opts = MtlCompileOptions()

compact_str = sprint(io->show(io, opts))
full_str = sprint(io->show(io, MIME"text/plain"(), opts))

@test opts.fastMathEnabled isa Bool
val = !opts.fastMathEnabled
opts.fastMathEnabled = val
@test opts.fastMathEnabled == val

@test opts.languageVersion isa VersionNumber
opts.languageVersion = v"1.0"
@test opts.languageVersion == v"1.0"

end

@testset "libraries" begin

dev = first(devices())
opts = MtlCompileOptions()

let lib = MtlLibrary(dev, "", opts)
    @test lib.device == dev
    @test lib.label == nothing
    @test isempty(lib.functionNames)
end

metal_code = read(joinpath(@__DIR__, "dummy.metal"), String)
let lib = MtlLibrary(dev, metal_code, opts)
    @test lib.device == dev
    @test lib.label == nothing
    fns = lib.functionNames
    @test length(fns) == 2
    @test "kernel_1" in fns
    @test "kernel_2" in fns
end

binary_path = joinpath(@__DIR__, "dummy.metallib")
let lib = MtlLibraryFromFile(dev, binary_path)
    @test lib.device == dev
    @test lib.label == nothing
    fns = lib.functionNames
    @test length(fns) == 2
    @test "kernel_1" in fns
    @test "kernel_2" in fns
end

binary_code = read(binary_path)
let lib = MtlLibraryFromData(dev, binary_code)
    @test lib.device == dev
    @test lib.label == nothing
    fns = lib.functionNames
    @test length(fns) == 2
    @test "kernel_1" in fns
    @test "kernel_2" in fns

    compact_str = sprint(io->show(io, lib))
    full_str = sprint(io->show(io, MIME"text/plain"(), lib))
end

end

@testset "functions" begin

dev = first(devices())
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "dummy.metallib"))
fun = MtlFunction(lib, "kernel_1")

compact_str = sprint(io->show(io, fun))
full_str = sprint(io->show(io, MIME"text/plain"(), fun))

@test fun.device == dev
@test fun.label == nothing
@test fun.name == "kernel_1"
@test fun.functionType == MTL.MtFunctionTypeKernel

end

@testset "events" begin

dev = first(devices())

let ev = MtlEvent(dev)
    @test ev.device == dev
    @test ev.label == nothing
end

let ev = MtlSharedEvent(dev)
    @test ev.device == dev
    @test ev.label == nothing
    @test ev.signaledValue == 0
end

end

@testset "fences" begin

dev = first(devices())

let fen = MtlFence(dev)
    @test fen.device == dev
end

end

@testset "heap" begin

dev = first(devices())

let desc = MtlHeapDescriptor()
    @test desc.type == MTL.MtHeapTypeAutomatic
    desc.type = MTL.MtHeapTypePlacement
    @test desc.type == MTL.MtHeapTypePlacement

    @test desc.size == 0
    desc.size = 1024
    @test desc.size == 1024

    @test desc.storageMode == MTL.MtStorageModePrivate
    desc.storageMode = MTL.MtStorageModeShared
    @test desc.storageMode == MTL.MtStorageModeShared

    @test desc.cpuCacheMode == MTL.MtCPUCacheModeDefaultCache
    desc.cpuCacheMode = MTL.MtCPUCacheModeWriteCombined
    @test desc.cpuCacheMode == MTL.MtCPUCacheModeWriteCombined

    @test desc.hazardTrackingMode == MTL.MtHazardTrackingModeDefault
    desc.hazardTrackingMode = MTL.MtHazardTrackingModeUntracked
    @test desc.hazardTrackingMode == MTL.MtHazardTrackingModeUntracked

    @test desc.resourceOptions == MTL.MtResourceStorageModeShared |
                                  MTL.MtResourceCPUCacheModeWriteCombined |
                                  MTL.MtResourceHazardTrackingModeUntracked
    desc.resourceOptions = MTL.MtResourceStorageModePrivate |
                           MTL.MtResourceCPUCacheModeDefaultCache |
                           MTL.MtResourceHazardTrackingModeDefault
    @test desc.resourceOptions == MTL.MtResourceStorageModePrivate |
                                  MTL.MtResourceCPUCacheModeDefaultCache |
                                  MTL.MtResourceHazardTrackingModeDefault

    # setting resource options should be reflected in individual fields
    @test desc.storageMode == MTL.MtStorageModePrivate
    @test desc.cpuCacheMode == MTL.MtCPUCacheModeDefaultCache
    @test desc.hazardTrackingMode == MTL.MtHazardTrackingModeDefault
end

desc = MtlHeapDescriptor()
desc.size = 0x4000 # TODO: use heapBufferSizeAndAlign
let heap = MtlHeap(dev, desc)
    @test heap.label == nothing

    @test heap.type == desc.type

    @test heap.size == desc.size

    # NOTE: these checks are fragile, as the heap options seems to depend on the requested size

    #=@test=# heap.storageMode == desc.storageMode

    #=@test=# heap.cpuCacheMode == desc.cpuCacheMode

    #=@test=# heap.hazardTrackingMode == desc.hazardTrackingMode

    #=@test=# heap.resourceOptions == desc.resourceOptions
end

end

@testset "buffers" begin

dev = first(devices())

buf = MtlBuffer{Int}(dev, 1)

@test buf.length == 8

# MtlResource properties
@test buf.device == dev
@test buf.label == nothing

@test sizeof(buf) == 8

free(buf)

end

@testset "command queue" begin

dev = first(devices())

cmdq = MtlCommandQueue(dev)

@test cmdq.device == dev
@test cmdq.label == nothing

end

@testset "command buffer" begin

dev = first(devices())
cmdq = MtlCommandQueue(dev)

cmdbuf = MtlCommandBuffer(cmdq)

@test cmdbuf.device == dev
@test cmdbuf.commandQueue == cmdq
@test cmdbuf.label == nothing
@test cmdbuf.error === nothing
@test cmdbuf.status == MTL.MtCommandBufferStatusNotEnqueued
@test cmdbuf.kernelStartTime == 0
@test cmdbuf.kernelEndTime == 0
@test cmdbuf.gpuStartTime == 0
@test cmdbuf.gpuEndTime == 0

let ev = MtlSharedEvent(dev)
    @test ev.signaledValue == 0
    encode_signal!(cmdbuf, ev, 42)
    encode_wait!(cmdbuf, ev, 21)
    commit!(cmdbuf)
    wait_completed(cmdbuf)
    @test ev.signaledValue == 42
end

cmdbuf = MtlCommandBuffer(cmdq)
@test cmdbuf.status == MTL.MtCommandBufferStatusNotEnqueued
enqueue!(cmdbuf)
@test cmdbuf.status == MTL.MtCommandBufferStatusEnqueued
commit!(cmdbuf)
@test cmdbuf.status == MTL.MtCommandBufferStatusCommitted
# Completion happens too quickly to test for committed status to be checked
wait_completed(cmdbuf) == MTL.MtCommandBufferStatusCompleted
@test cmdbuf.status == MTL.MtCommandBufferStatusCompleted

# CommandBufferDescriptor tests
desc = MTL.mtNewCommandBufferDescriptor()
@test MTL.mtCommandBufferDescriptorRetainedReferences(desc) == true
MTL.mtCommandBufferDescriptorRetainedReferencesSet(desc,false)
@test MTL.mtCommandBufferDescriptorRetainedReferences(desc) == false

@test MTL.mtCommandBufferDescriptorErrorOptions(desc) == MTL.MtCommandBufferErrorOptionNone
MTL.mtCommandBufferDescriptorErrorOptionsSet(desc,MTL.MtCommandBufferErrorOptionEncoderExecutionStatus)
@test MTL.mtCommandBufferDescriptorErrorOptions(desc) == MTL.MtCommandBufferErrorOptionEncoderExecutionStatus

cmq = MtlCommandQueue(device())
cmdbuf = MtlCommandBuffer(cmq; retainReferences=false, errorOption=MTL.MtCommandBufferErrorOptionEncoderExecutionStatus)
@test cmdbuf.retainedReferences == false
@test cmdbuf.errorOptions == MTL.MtCommandBufferErrorOptionEncoderExecutionStatus

end

@testset "compute pipeline" begin

dev = first(devices())
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "dummy.metallib"))
fun = MtlFunction(lib, "kernel_1")

pipeline = MtlComputePipelineState(dev, fun)

@test pipeline.device == dev
@test pipeline.label === nothing

@test pipeline.maxTotalThreadsPerThreadgroup isa Integer
@test pipeline.threadExecutionWidth isa Integer
@test pipeline.staticThreadgroupMemoryLength == 0

end

@testset "argument encoder" begin

dev = first(devices())
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "vadd.metallib"))
fun = MtlFunction(lib, "vadd")

encoder = MtlArgumentEncoder(fun, 1)

@test encoder.encodedLength == 0
@test encoder.alignment == 1

# TODO: actually encode arguments

end

# TODO: continue adding tests

end

@testset "arrays" begin

mtl_arr = MtlArray{Int}(undef, 1)
arr = Array(mtl_arr)

@test sizeof(arr) == 8
@test length(arr) == 1
@test eltype(arr) == Int

end

@testset "kernels" begin

function tester(A)
    idx = thread_position_in_grid_1d()
    A[idx] = Int(5)
    return nothing
end

bufferSize = 8
bufferA = MtlArray{Int,1}(undef, tuple(bufferSize), storage=Shared)
vecA = unsafe_wrap(Vector{Int}, bufferA.buffer, tuple(bufferSize))

Metal.@sync @metal threads=(bufferSize) tester(bufferA.buffer)
@test all(vecA .== Int(5))

@testset "launch params" begin
    vecA .= 0
    Metal.@sync @metal threads=(2) tester(bufferA.buffer)
    @test all(vecA == Int.([5, 5, 0, 0, 0, 0, 0, 0]))
    vecA .= 0

    Metal.@sync @metal grid=(3) threads=(2) tester(bufferA.buffer)
    @test all(vecA == Int.([5, 5, 5, 5, 5, 5, 0, 0]))
    vecA .= 0

    @test_throws InexactError @metal threads=(-2) tester(bufferA.buffer)
    @test_throws InexactError @metal grid=(-2) tester(bufferA.buffer)
    @test_throws ArgumentError @metal threads=(1025) tester(bufferA.buffer)
    @test_throws ArgumentError @metal threads=(1000,2) tester(bufferA.buffer)
end

@testset "argument buffers" begin
    Metal.@sync @metal threads=(bufferSize) tester(bufferA)
    @test all(vecA .== Int(5))
    vecA .= 0

    function no_intrinsic(A)
        A[1] += Int(5)
        return nothing
    end
    Metal.@sync @metal no_intrinsic(bufferA)
    @test all(vecA == Int.([5, 0, 0, 0, 0, 0, 0, 0]))

    function types_tester(A::MtlDeviceVector{T}) where T
        idx = thread_position_in_grid_1d()
        A[idx] = T(5)
        return nothing
    end
    types = [Float32, Float16, Int64, Int32, Int16, Int8]
    for typ in types
        bufferA = MtlArray{typ,1}(undef, tuple(bufferSize), storage=Shared)
        vecA = unsafe_wrap(Vector{typ}, bufferA.buffer, tuple(bufferSize))
        Metal.@sync @metal threads=(bufferSize) types_tester(bufferA)
        @test all(vecA .== typ(5))
    end
end

@testset "math intrinsics" begin
    a = ones(Float32,1)
    a .* Float32(3.14)
    bufferA = MtlArray(a)
    vecA = unsafe_wrap(Vector{Float32}, bufferA.buffer, 1)

    function intr_test(buf)
        idx = thread_position_in_grid_1d()
        buf[idx] = cos(buf[idx])
        return nothing
    end
    @metal intr_test(bufferA.buffer)
    synchronize()
    @test vecA ≈ cos.(a)

    function intr_test2(buf)
        idx = thread_position_in_grid_1d()
        buf[idx] = Metal.rsqrt(buf[idx])
        return nothing
    end
    @metal intr_test2(bufferA.buffer)
end

@testset "sync barriers" begin
    function sync_test_kernel(buf)
        idx = thread_position_in_grid_1d()
        buf[idx] += UInt8(1)
        return nothing
    end
    buf = MtlArray{UInt8,1}(undef, tuple(1024); storage=Shared)
    vec = unsafe_wrap(Vector{UInt8}, buf.buffer, (1024))
    @metal threads=1024 sync_test_kernel(buf)
    synchronize()
    @test all(vec .== UInt8(1))

    function barrier_test_kernel(buf)
        idx = thread_position_in_grid_1d()
        if thread_position_in_threadgroup_1d() != UInt32(1)
            for i in range(1,threads_per_threadgroup_1d())
                buf[idx] += UInt32(i)
            end
            buf[idx] = 1
        end

        threadgroup_barrier()

        if thread_position_in_threadgroup_1d() == UInt32(1)
            for i in range(1,threads_per_threadgroup_1d())
                buf[idx] += buf[idx+i-1]
            end
        end
        return nothing
    end

    buf = MtlArray{Int,1}(undef, tuple(1024); storage=Shared)
    vec = unsafe_wrap(Vector{Int}, buf.buffer, (1024))
    Metal.@sync @metal threads=1024 barrier_test_kernel(buf)
    @test vec[1] == 992

    # TODO: simdgroup barrier test
end

@testset "values as references" begin
    function kernel(a, val)
        @inbounds a[thread_index_in_threadgroup()] = val
        return
    end

    a = MtlArray{Float32}(undef, 1)
    @metal kernel(a, 42f0)
    @test Array(a) == [42f0]
end

end # End kernels testset

# Examples
# TODO: Do this in a way more similar to the other backends
@testset "examples" begin
    include("../examples/unified_memory.jl")
    include("../examples/vadd.jl")
end


text:
struct DensityDependentQSpec{T} <: SpatQSimSpec
    realization::Int
    densdep_mult::T
    prep_file::String

    function DensityDependentQSpec(realization::Int,
                             densdep_mult::T,
                             prep_file::String = "prep.h5") where T
        new{T}(realization, densdep_mult, prep_file)
    end
end

sim_value(spec::DensityDependentQSpec) = spec.densdep_mult

function comm_catchability(spec::DensityDependentQSpec)
    DensityDependentCatchability(base_catchability(), spec.densdep_mult)
end

ext: jl


text:
export
    SyncBP,
    AsyncBP,
    CoherentBP

"""
    SyncBP(range_size = 10)

A window filter that uses a synchronous window with ThreePassBP with the given range size.
"""
SyncBP(range_size = 10) = WindowFilter(SyncWindow(), ThreePassBP(range_size))

"""
    AsyncBP(range_size = 10, T = Float64)

A window filter that uses an asynchronous window with ThreePassBP with the given range size.
`T` represents the time type and must be the same as used in creation of the runtime.
"""
AsyncBP(range_size = 10, T = Float64) = WindowFilter(AsyncWindow{T}(), ThreePassBP(range_size))

"""
    CoherentBP(range_size = 10, T = Float64)

A window filter that uses a coherent window with ThreePassBP with the given range size.
`T` represents the time type and must be the same as used in creation of the runtime.
"""
CoherentBP(range_size = 10, T = Float64) = WindowFilter(CoherentWindow{T}(), ThreePassBP(range_size))


text:
module ExcelUtilities

using CSV
using Tables
import InteractiveUtils

export xlclip

"""
    xlclip()

Copy Excel-copied data from the clipboard into a Julia vector or matrix (depending on shape). Single column or row are converted to Julia `Vector`s. Two dimensional data is copied into a matrix (2D array).

The data is parsed using [CSV.jl](https://github.com/JuliaData/CSV.jl); if you have questions about how it's being read check the documentation for that package.

Common issues are that the data you are trying to copy has numbers that are formatted in a way that is not natively recognized as numbers in Julia. Either make sure to make your numbers formatted in a standard format in Excel, e.g. `-1235.32` and not `(1,235.32)`.
"""
function xlclip()
    xlclip_reader(InteractiveUtils.clipboard())
end

# use a barrier function to allow precompilation/testing of expensive part without
# messing with user's clipboard
function xlclip_reader(str)
    res = CSV.File(IOBuffer(str),header=false,delim="\t") |> Tables.matrix
    m, n = size(res)
    if m == 1 || n == 1
        return vec(res)
    else
        return res
    end
end

"""
    xlclip(data)

Copy Julia array to the clipboard in an Excel-friendly format.

Vectors will be copied as Excel columns; to copy a vector `v` to a row for Excel, you can transpose it with a single quote: `xlclip(v')`
"""
function xlclip(data)
    InteractiveUtils.clipboard(xlclip_writer(data)) # drop the trailing newline
end

function xlclip_writer(data)
    #Tables.table wants an array, not vector so convert if vector
    if typeof(data) <: Vector
        data = reshape(data, :,1)
    end
    rows = collect(CSV.RowWriter(Tables.table(data),delim = "\t"))
    result = join(rows[2:end]) # don't include the column names
   return result[1:end-1] # drop the trailing newline
end

include("precompile.jl")
_precompile_()
end

max_issues_repo_name: alecloudenback/ExcelUtilities

max_forks_repo_name: alecloudenback/ExcelUtilities


text:
using IntervalArithmetic, IntervalContractors

include("atan2.jl")
include("c-xsc.jl")
include("fi_lib.jl")
include("ieee1788-constructors.jl")
include("ieee1788-exceptions.jl")
include("libieeep1788_bool.jl")
include("libieeep1788_cancel.jl")
include("libieeep1788_class.jl")
include("libieeep1788_elem.jl")
include("libieeep1788_num.jl")
include("libieeep1788_overlap.jl")
include("libieeep1788_rec_bool.jl")
include("libieeep1788_reduction.jl")
include("libieeep1788_set.jl")
include("mpfi.jl")
include("abs_rev.jl")
include("libieeep1788_rev.jl")
include("libieeep1788_mul_rev.jl")
include("pow_rev.jl")


text:
"""
id: an identifier that uniquely identifies the component in a lens
meta_data: any additional meta data for this component
surfaces
"""
struct OpticalComponent{T<:Real}
    id::Symbol
    type::DataType
    meta_data::Any
    surfaces::AbstractVector{OpticalSurface{T}}
    function OpticalComponent(
        id,
        type,
        meta_data,
        surfaces::AbstractVector{OpticalSurface{T}},
    ) where {T}
        if !allunique(filter(s -> !isnothing(s.id), surfaces))
            error("Surface IDs must be nothing or unique")
        end
        new{T}(id, type, meta_data, surfaces)
    end
end

function track_length(oc::OpticalComponent{T}) where {T}
    sum(map(s -> s.t, oc.surfaces))
end

function track_length(os::AbstractVector{OpticalSurface{T}}) where {T}
    sum(map(s -> s.t, os))
end

function track_length(os::AbstractVector{OpticalSurface{T}}, id::Symbol) where {T}
    l = zero(T)
    for s in os
        l += s.t
        if s.id == id
            return l
        end
    end
    nothing
end

function track_length(c::OpticalComponent{T}, id::Symbol) where {T}
    track_length(c.surfaces, id)
end


text:
########## Rolling window ##########

"""
Apply a function on a rolling window with truncation at boundaries.

# Args
- func: function to be applied
- X: input array. The rolling window is applied on the last dimension from small to large index (time arrow). For example, on a vector it runs through the column direction (i.e. vertically), while on a matrix it runs through the row direction (i.e. horizontally, and each column corresponds to an observation).
- w: number of samples on the rolling window
- d: downsampling factor on the rolling window
- p: step of the rolling window

# Kwargs
- mode: :causal or :anticausal
- boundary: truncation of boundary, :hard or :soft
"""
function rolling_apply(func::Function, X::AbstractArray, w::Integer, d::Integer=1, p::Integer=1; mode::Symbol=:causal, boundary::Symbol=:hard)
    return if boundary == :hard
        rolling_apply_hard(func, X, w, d, p; mode=mode)
    elseif boundary == :soft
        rolling_apply_soft(func, X, w, d, p; mode=mode)
    else
        error("Unknown boundary condition: $(boundary)")
    end
end

"""
# Notes
The result is empty if `d*(w-1)+1 > L`
"""
function rolling_apply_hard(func::Function, X::AbstractArray, w::Integer, d::Integer=1, p::Integer=1; mode::Symbol=:causal)
    @assert w>=1 && d>=1 && p>=1

    L = size(X)[end]  # number of samples in X
    res = []
    if mode==:causal
        gs = t -> reverse(t:-d:(t-d*(w-1)))  # grid of samples, ending by t
        for t=L:-p:(1+d*(w-1))
            # V = ndims(X) == 1 ? reverse(view(X,t:-d:t-d*w+1)) : reverse(view(X,:,t:-d:t-d*w+1),dims=2)  # old version, WRONG
            # correct version:
            # V = ndims(X) == 1 ? view(X,gs(t)) : view(X,:,gs(t))
            # pushfirst!(res, (t, func(V)))
            pushfirst!(res, (t, func(selectdim(X, ndims(X), gs(t)))))
        end
    else  # anti-causal
        gs = t -> t:d:(t+d*(w-1))  # grid of samples, starting by t
        for t=1:p:(L-d*(w-1))
            # V = ndims(X) == 1 ? view(X,t:d:t+d*w-1) : view(X,:,t:d:t+d*w-1)  # old version, WRONG
            # correct version
            # V = ndims(X) == 1 ? view(X,gs(t)) : view(X,:,gs(t))
            # push!(res, (t, func(V)))
            push!(res, (t, func(selectdim(X, ndims(X), gs(t)))))
        end
    end
    return res
end

function rolling_apply_soft(func::Function, X::AbstractArray, w::Integer, d::Integer=1, p::Integer=1; mode::Symbol=:causal)
    @assert w>=1 && d>=1 && p>=1

    L = size(X)[end]  # number of samples in X
    res = []
    if mode==:causal
        gs = t -> reverse(t:-d:max(1,t-d*(w-1)))  # grid of samples, ending by t
        for t=L:-p:1
            pushfirst!(res, (t, func(selectdim(X, ndims(X), gs(t)))))
        end
    else  # anti-causal
        gs = t -> t:d:min(L,t+d*(w-1))  # grid of samples, starting by t
        for t=1:p:L
            push!(res, (t, func(selectdim(X, ndims(X), gs(t)))))
        end
    end
    return res
end

"""
Rolling vectorization.

Vectorize the contents of a rolling window and make horizontal concatenation.

# Args
- X: input array. The window rolls in the last dimension from small to large index.
- w: size of rolling window, or number of consecutive samples to be concatenated together
- d: downsampling factor
- p: period
- mode: :causal (proceed from right to left) or :anticausal (proceed from left to right)

# Returns
A named tuple `(index, value)` which are
- index of `X` where the rolling window starts (mode=:anticausal) or ends (mode=:causal, default)
- output matrix that row dimension equals to `w` if parameters are valid, otherwise it is an empty matrix.

# Notes
When `d==1`, two special cases are
- `w==1`: this function is equivalent to `reshape(X, :, size(X)[end])`
- `w==size(X)[end]`: this function is equivalent to `reshape(X, :, 1)`
"""
function rolling_vectorize(X::AbstractArray{T}, w::Integer, d::Integer=1, p::Integer=1; kwargs...) where {T}
    # res = rolling_apply_hard(x->vec(hcat(x...)), X, w, 1, p; kwargs...)
    # res = rolling_apply_hard(x->vec(x), X, w, 1, p; kwargs...)
    res = rolling_apply_hard(x->vec(x), X, w, d, p; kwargs...)

    tidx::Vector{Integer} = [x[1] for x in res]
    # In case of empty array `hcat([]...)` gives `0-element Array{Any,1}`
    # Assure the return type in this case is a numerical, not `any`.
    val::Matrix{T} = length(tidx)==0 ? Array{T}(undef,0,0) : hcat([x[2] for x in res]...)
    return (index=tidx, value=val)
end

"""
Rolling mean in row direction.
"""
function rolling_mean(X::AbstractVecOrMat{<:Number}, w::Integer, d::Integer=1, p::Integer=1; kwargs...)
    res = rolling_apply(x->mean(x, dims=1), X, w, d, p; kwargs...)
    Xm = hcat([x[2] for x in res]...)[:]
    # tidx = [x[1] for x in res]
    return Xm
end

"""
Rolling median in row direction.
"""
function rolling_median(X::AbstractVecOrMat{<:Number}, w::Integer, d::Integer=1, p::Integer=1; kwargs...)
    res = rolling_apply(x->median(x, dims=1), X, w, d, p; kwargs...)
    Xm = hcat([x[2] for x in res]...)[:]
    # tidx = [x[1] for x in res]
    return Xm
end

######## Below is not maintained ########

"""
Make rolling orders based on fGn-MLE.

# Args
- X: raw price
- d: time lag for finite difference operator of the return
- w: number of samples on the rolling window
- d: downsampling factor on the rolling window
- p: period of the rolling window
- s: size of sub window
- l: length of decorrelation
- mode: :causal or :anticausal
- boundary: truncation, :hard or :soft

"""
function rolling_order_fGn_MLE(R::AbstractVector{<:Real}, d::Integer; wsize::Integer, dfct::Integer, pov::Integer, ssize::Integer, dlen::Integer, nspl::Integer, sdfc::Integer, nprd::Integer, dmode::Symbol=:logscale, kwargs...)
    @assert 0 < ssize <= wsize

    func = x -> fGn_MLE_estim_predict(x, d, ssize, dlen, nspl, sdfc, nprd; dmode=dmode)

    # kwargs_rolling = Dict(k=>v for (k,v) in kwargs if k in [:causal, :boundary])
    return rolling_apply(func, R, wsize, dfct, pov; kwargs...)
end

# function rolling_order_fGn_MLE(P::AbstractVector{<:Real}, d::Integer; wsize::Integer, pov::Integer, dmode::Symbol=:regular, kwargs...)
#     R = price2return(P, d, mode=:causal, method=:ori)  # return
#     # println(size(R))
#     func = x -> fGn_MLE_estim_predict(x, 1, wsize, 1, 0, 1, 1; dmode=dmode)

#     # kwargs_rolling = Dict(k=>v for (k,v) in kwargs if k in [:causal, :boundary])
#     return rolling_apply(func, R, wsize, d, pov; kwargs...)
# end

# abstract type AbstractDataStream <: TimeSeries.AbstractTimeSeries end

# function fetch_data(X::AbstractDataStream, )
# end

# abstract type AbstractRollingFunction end

# struct RollingEstimPredictor <: AbstractRollingFunction
#     prepare::Function  # data preparation
#     estim::Function  # estimator
#     predict::Function  # predictor
# end

"""
    rolling_estim(estim::Function, X0::AbstractVecOrMat{T}, (w,s,d)::Tuple{Integer,Integer,Integer}, p::Integer, trans::Function=(x->vec(x)); mode::Symbol=:causal) where {T<:Real}

Rolling window estimator for 1d or multivariate time series.

# Args
- estim: estimator which accepts either 1d or 2d array
- X0: input, 1d or 2d array with each column being one observation
- (w,s,d): size of rolling window; size of sub-window, length of decorrelation
- p: step of the rolling window
- trans: function of transformation which returns a vector of fixed size or a scalar,  optional
- mode: :causal or :anticausal

# Returns
- array of estimations on the rolling window

# Notes
- The estimator is applied on a rolling window every `p` steps. The rolling window is divided into `n` (possibly overlapping) sub-windows of size `w` at the pace `d`, such that the size of the rolling window equals to `(n-1)*d+w`. For q-variates time series, data on a sub-window is a matrix of dimension `q`-by-`w` which is further transformed by the function `trans` into another vector. The transformed vectors of `n` sub-windows are concatenated into a matrix which is finally passed to the estimator `estim`. As example, for `trans = x -> vec(x)` the data on a rolling window is put into a new matrix of dimension `w*q`-by-`n`, and its columns are the column-concatenantions of data on the sub-window. Moreover, different columns of this new matrix are assumed as i.i.d. observations.
- Boundary is treated in a soft way.
"""
function rolling_estim(estim::Function, X0::AbstractVecOrMat{T}, (w,s,d)::Tuple{Integer,Integer,Integer}, p::Integer, trans::Function=(x->vec(x)); mode::Symbol=:causal, nan::Symbol=:ignore) where {T<:Number}
    X = ndims(X0)>1 ? X0 : reshape(X0, 1, :)  # vec to matrix, create a reference not a copy
    L = size(X,2)
    # println(L); println(w); println(s)
    # @assert w >= s && L >= s
    res = []

    if mode == :causal  # causal
        for t = L:-p:1
            # xv = view(X, :, t:-1:max(1, t-w+1))
            xv = view(X, :, max(1, t-w+1):t)
            xs = if nan == :ignore
                idx = findall(.!any(isnan.(xv), dims=1)[:])  # ignore columns containing nan values
                if length(idx) > 0
                    rolling_apply_hard(trans, view(xv,:,idx), s, d; mode=:causal)
                else
                    []
                end
            else
                rolling_apply_hard(trans, xv, s, d; mode=:causal)
            end

            if length(xs) > 0
                # println(size(xs))
                pushfirst!(res, (t,estim(squeezedims(xs, dims=[1,2]))))
            end
        end
    else  # anticausal
        for t = 1:p:L
            xv = view(X, :, t:min(L, t+w-1))
            xs = if nan == :ignore
                idx = findall(.!any(isnan.(xv), dims=1)[:])  # ignore columns containing nan values
                if length(idx) > 0
                    rolling_apply_hard(trans, view(xv,:,idx), s, d; mode=:anticausal)
                else
                    []
                end
            else
                rolling_apply_hard(trans, xv, s, d; mode=:anticausal)
            end

            if length(xs) > 0
                push!(res, (t,estim(squeezedims(xs, dims=[1,2]))))
            end
        end
    end
    return res
end

function rolling_estim(func::Function, X0::AbstractVecOrMat{T}, w::Integer, p::Integer, trans::Function=(x->x); kwargs...) where {T<:Number}
    return rolling_estim(func, X0, (w,w,1), p, trans; kwargs...)
end

function rolling_regress_predict(regressor::Function, predictor::Function, X0::AbstractVecOrMat{T}, Y0::AbstractVecOrMat{T}, X1::AbstractVecOrMat{T}, Y1::AbstractVecOrMat{T}, (w,s,d)::Tuple{Integer,Integer,Integer}, p::Integer, trans::Function=(x->vec(x)); mode::Symbol=:causal) where {T<:Number}
    X = ndims(X0)>1 ? X0 : reshape(X0, 1, :)  # vec to matrix, create a reference not a copy
    Y = ndims(Y0)>1 ? Y0 : reshape(Y0, 1, :)
    @assert size(X,2) == size(Y,2)

    @assert size(X0) == size(X1)
    @assert size(Y0) == size(Y1)
    Xp = ndims(X1)>1 ? X1 : reshape(X1, 1, :)  # vec to matrix, create a reference not a copy
    Yp = ndims(Y1)>1 ? Y1 : reshape(Y1, 1, :)

    (any(isnan.(X)) || any(isnan.(Y))) && throw(ValueError("Inputs cannot contain nan values!"))
    (any(isnan.(Xp)) || any(isnan.(Yp))) && throw(ValueError("Inputs cannot contain nan values!"))
    L = size(X,2)  # total time
    @assert L >= w >= s

    res = []

    if mode == :causal  # causal
        tf, xf, yf = [], [], []  # future time and data used for prediction

        for t = L:-p:w
            # printfmtln("Processing time {}...\r", t)

            xv = rolling_apply_hard(trans, view(X, :, t-w+1:t), s, d; mode=:causal)
            yv0 = view(Y, :, t-w+1:t)[:, end:-d:1]  # time-reversed
            yv = yv0[:,1:size(xv,2)][:,end:-1:1]  # reverse back

            reg = regressor(yv, xv)  # regression

            xf = trans(Xp[:,t-s+1:t])
            yf = Yp[:,t]
            prd = predictor(reg, xf, yf)  # prediction using future data

            pushfirst!(res, (time=t, regression=reg, prediction=prd))
        end
    else  # anticausal: TODO
        for t = 1:p:L
            xv = rolling_apply_hard(trans, view(X, :, t:min(L, t+w-1)), s, d; mode=:anticausal)
            yv = rolling_apply_hard(trans, view(Y, :, t:min(L, t+w-1)), 1, d; mode=:anticausal)

            reg = regressor(yv, xv)
            prd = isempty(res_reg) ? [] : predictor(res_reg[end], xv, yv[:,end])
            # TODO
        end
    end

    return res
end

function rolling_estim_predict(estim::Function, predict::Function, X0::AbstractVecOrMat{T}, (w,s,d)::Tuple{Integer,Integer,Integer}, p::Integer, trans::Function=(x->vec(x)); mode::Symbol=:causal, nan::Symbol=:ignore) where {T<:Number}
    X = ndims(X0)>1 ? X0 : reshape(X0, 1, :)  # vec to matrix, create a reference not a copy
    L = size(X,2)
    # println(L); println(w); println(s)
    # @assert w >= s && L >= s
    res = []

    if mode == :causal  # causal
        for t = L:-p:1
            # xv = view(X, :, t:-1:max(1, t-w+1))
            xv = view(X, :, max(1, t-w+1):t)
            xs = if nan == :ignore
                idx = findall(.!any(isnan.(xv), dims=1)[:])  # ignore columns containing nan values
                if length(idx) > 0
                    rolling_apply_hard(trans, view(xv,:,idx), s, d; mode=:causal)
                else
                    []
                end
            else
                rolling_apply_hard(trans, xv, s, d; mode=:causal)
            end
            xs  = hcat([x[2] for x in xs]...)  # added 2019-01
            if length(xs) > 0
                # println(size(xs))
                Θ = estim(squeezedims(xs, dims=[1,2]))  # estimations
                Π = predict(Θ, xv)
                pushfirst!(res, (t,Π))
            end
        end
    else  # anticausal
        for t = 1:p:L
            xv = view(X, :, t:min(L, t+w-1))
            xs = if nan == :ignore
                idx = findall(.!any(isnan.(xv), dims=1)[:])  # ignore columns containing nan values
                if length(idx) > 0
                    rolling_apply_hard(trans, view(xv,:,idx), s, d; mode=:anticausal)
                else
                    []
                end
            else
                rolling_apply_hard(trans, xv, s, d; mode=:anticausal)
            end

            if length(xs) > 0
                Θ = estim(squeezedims(xs, dims=[1,2]))  # estimations
                Π = predict(Θ, xv)
                push!(res, (t,Π))
            end
        end
    end
    return res
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
function test_tntp()
  data_dir = download_tntp()

  # Test
  @testset "Loading Network Data" begin
    for d in readdir(data_dir)
      if isdir(joinpath(data_dir, d)) && d != "_scripts"
        @testset "$d" begin
          try
            read_ta_network(d)
            td = load_ta_network(d)
            @info "Network '$d' is OK."
            @test td.network_name == d
          catch e
            @show e
            @warn "Directory '$d' does not contain a compatible dataset."
          end
        end
      end
    end
  end
end


text:
# This file was generated, do not modify it. # hide
x = [I; randn(n - k, k)]

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# included in module Wut

struct ScalarEncoder <: AbstractEncoder

    w::Int64             # width of the contiguous 1 bits, must be odd
    minval::Float64      # minimum value of input signal
    maxval::Float64      # maximum value of input signal (strictly less when periodic = true)
    periodic::Bool       # input values wrap around and maxval input is mapped onto minval
    n::Int64             # output bit width (must be greater than w)
    radius::Float64      # input value difference > radius do not overlap bits
    resolution::Float64  # input value difference >= resolution do not give identical bits
    name::AbstractString # optional descriptive string
    verbosity::Int64
    clipInput::Bool      # Non-periodic inputs are clipped to minval / maxval
    forced::Bool         # skip some safety checks
    padding::Int64
    range::Float64

    function ScalarEncoder(;w=0,minval=0.0,maxval=0.0,periodic=false,
            n=0,radius=0.0,resolution=0.0,name="",
            verbosity=0,clipInput=false,forced=false,padding=0,range=0.0)

        res = rad = rag = 0.0

        pad = padding > 0 ? padding : ( periodic ? 0 : convert(Int64,round((w-1)/2)) )

        if n != 0 && radius == 0.0 && resolution == 0.0
           res = periodic ? (maxval-minval)/n : (maxval-minval)/(n-w)
           rad = w*res
           rag = periodic ? maxval-minval : maxval-minval+res

        elseif n == 0 && radius != 0.0 && resolution == 0.0
           rad = radius
           res = radius / w

        elseif n == 0 && radius == 0.0 && resolution != 0.0
           res = resolution
           rad = res * w

        else

        end

        if n == 0
            rag = periodic ? maxval-minval : maxval-minval+res
            nfloat = (w * (rag / rad)) + (2 * pad)
            n = Int32(ceil(nfloat))
        end

        #pad = periodic ? 0 : convert(Int64,round((w-1)/2))
        nm = (name != "") ? name : @sprintf("[%s:%s]" , minval, maxval)

        new(w,minval,maxval,periodic,n,rad,res,nm,verbosity,clipInput,forced,pad,rag)
    end

end
export ScalarEncoder

function getWidth(e::ScalarEncoder)
    e.n
end
export getWidth

function getDescription(e::ScalarEncoder)
    e.name
end
export getDescription

function getName(e::ScalarEncoder)
    e.name
end
export getName

function getBucketIndices(e::ScalarEncoder)

end
export getBucketIndices

function encodeIntoArray(e::ScalarEncoder,n::Number,b::BitPat;learn=true, offset=0, length=0)
    #fill!(b.b,false)
    for i = 1 : e.n
       b.b[i+offset]=false
    end
    t1 = e.periodic ? n % e.maxval : n
    t2 = t1 > e.maxval ? e.maxval : t1
    t = t2 < e.minval ? e.minval : t2

    c = convert(Int64,trunc((((t - e.minval) + e.resolution/2) / e.resolution ) + e.padding))
    d = convert(Int64,(e.w-1)/2)

    for i = 1+c-d : 1+c+d
      if i <= e.n && i > 0
        b.b[i+offset]=true
      end
    end

    b
end
export encodeIntoArray

function decode(e::ScalarEncoder, b::BitPat; parentFieldName="")

    parentName = (parentFieldName == "") ? e.name : @sprintf("%s.%s", parentFieldName, e.name)
    x=Dict(parentName => [])
    r=("concrete",x,[parentName])

    if any(b.b)

      # ------------------------------------------------------------------------
      # First, assume the input pool is not sampled 100%, and fill in the
      #  "holes" in the encoded representation (which are likely to be present
      #  if this is a coincidence that was learned by the SP).

      # Search for portions of the output that have "holes"
      mz = (e.w - 1)/2
      for i = 0 : mz - 1
         sS = new(BitArray(i+3))
         fill!(sS,false)
         sS[1]=true
         sS[i+3]=true
         subLen=i+3

#        # Does this search string appear in the output?
        if e.periodic
          for j = 0 : E.n -1
#            outputIndices = numpy.arange(j, j + subLen)
#            outputIndices %= self.n
#            if numpy.array_equal(searchStr, tmpOutput[outputIndices]):
#              tmpOutput[outputIndices] = 1
           end
        else
          for j = 1 :(E.n - subLen + 1)
#            if numpy.array_equal(searchStr, tmpOutput[j:j + subLen]):
#              tmpOutput[j:j + subLen] = 1
          end
        end
      end
#      # ------------------------------------------------------------------------
#      # Find each run of 1's.
#      nz = tmpOutput.nonzero()[0]
#      runs = []     # will be tuples of (startIdx, runLength)
#      run = [nz[0], 1]
#      i = 1
#      while (i < len(nz)):
#        if nz[i] == run[0] + run[1]:
#          run[1] += 1
#        else:
#          runs.append(run)
#          run = [nz[i], 1]
#        i += 1
#      runs.append(run)

#      # If we have a periodic encoder, merge the first and last run if they
#      #  both go all the way to the edges
#      if self.periodic and len(runs) > 1:
#        if runs[0][0] == 0 and runs[-1][0] + runs[-1][1] == self.n:
#          runs[-1][1] += runs[0][1]
#          runs = runs[1:]

#      # ------------------------------------------------------------------------
#      # Now, for each group of 1's, determine the "left" and "right" edges, where
#      #  the "left" edge is inset by halfwidth and the "right" edge is inset by
#      #  halfwidth.
#      # For a group of width w or less, the "left" and "right" edge are both at
#      #   the center position of the group.
#      ranges = []
#      for run in runs:
#        (start, runLen) = run
#        if runLen <= self.w:
#          left = right = start + runLen / 2
#        else:
#          left = start + self.halfwidth
#          right = start + runLen - 1 - self.halfwidth

#        # Convert to input space.
#        if not self.periodic:
#          inMin = (left - self.padding) * self.resolution + self.minval
#          inMax = (right - self.padding) * self.resolution + self.minval
#        else:
#          inMin = (left - self.padding) * self.range / self.nInternal + self.minval
#          inMax = (right - self.padding) * self.range / self.nInternal + self.minval
#        # Handle wrap-around if periodic
#        if self.periodic:
#          if inMin >= self.maxval:
#            inMin -= self.range
#            inMax -= self.range

#        # Clip low end
#        if inMin < self.minval:
#          inMin = self.minval
#        if inMax < self.minval:
#          inMax = self.minval

#        # If we have a periodic encoder, and the max is past the edge, break into
#        #  2 separate ranges
#        if self.periodic and inMax >= self.maxval:
#          ranges.append([inMin, self.maxval])
#          ranges.append([self.minval, inMax - self.range])
#        else:
#          if inMax > self.maxval:
#            inMax = self.maxval
#          if inMin > self.maxval:
#            inMin = self.maxval
#          ranges.append([inMin, inMax])

#      desc = self._generateRangeDescription(ranges)

    end

    r

end
export decode

function getBucketValues(e::ScalarEncoder)

end
export getBucketValues

function getBucketInfo(e::ScalarEncoder)

end
export getBucketInfo

function topDownCompute(e::ScalarEncoder)

end
export topDownCompute

function closenessScores(e::ScalarEncoder)

end
export closenessScores

# end # module


text:
include("main.jl")

REWRITE_DATA = true

BUSES, LINES, GENERATORS, rPTDF = load_case_data()
ERROR_VARIANCES_TRUE, ERROR_DIST_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

# Test dependency on variance estimation confidence
casename = "risk_level_variation"
case_path = "results/$(casename)"

if REWRITE_DATA && isdir(case_path)
    rm(case_path, recursive=true)
end

if !(isdir(case_path)) || REWRITE_DATA
    # Only run if there are no saved results
    mkdir(case_path)

    testvalues = [0.15, 0.12, 0.1, 0.08, 0.05]

    sample_size = 750
    # Create global set of true observations to make tests comparable
    OBSERVATIONS = rand(ERROR_DIST_TRUE, sample_size)

    cctest_dict = DataFrame(eta=Any[], v_violation=Any[], f_violation=Any[], median_costs=Any[], mean_costs=Any[], min_costs=Any[], max_costs=Any[], quantile10_cost=Any[], quantile90_cost=Any[], median_delta=Any[], mean_delta=Any[])
    ddtest_dict = DataFrame(eta=Any[], v_violation=Any[], f_violation=Any[], median_costs=Any[], mean_costs=Any[], min_costs=Any[], max_costs=Any[], quantile10_cost=Any[], quantile90_cost=Any[], median_delta=Any[], mean_delta=Any[])

    for v in testvalues
        # Set Confidence Intervals
        η_v = v # 1 - Confidence for voltage limit cc
        η_g = v # 1 - Confidence for generation limit cc
        # SND Quantiles for global use
        snd = Normal(0,1)
        z_v = quantile(snd, 1-η_v)
        z_g = quantile(snd, 1-η_g)
        z_g2 = quantile(snd, 1-(η_g/2))

        info("Runnig for η = $v")
        ξ = 0.005
        ERROR_VARIANCES_TRUE, ERROR_DIST_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

        objective, bus_results, line_results, solvetimecc = Solve_DD_DR_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_SAMPLE)
        cctest = test_results(bus_results, line_results, objective; samples = sample_size, message="Standard CC", print_on=false)
        push!(cctest_dict, vcat(v, cctest))

        ddobjective, ddbus_results, ddline_results, solvetimedr = Solve_DD_DR_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_WC)
        ddtest = test_results(ddbus_results, ddline_results, ddobjective; samples = sample_size, message="DD CC", print_on=false)

        push!(ddtest_dict, vcat(v, ddtest))
    end

    writetable("$(case_path)/cctest_dict.csv", cctest_dict)
    writetable("$(case_path)/ddtest_dict.csv", ddtest_dict)

else
    # If there is saved data load it
    info("Data already exists at $case_path")
    cctest_dict = readtable("$(case_path)/cctest_dict.csv")
    ddtest_dict = readtable("$(case_path)/ddtest_dict.csv")
end


text:
#####
##### Abstract tree/trajectory interface
#####

####
#### Directions
####

"Maximum number of iterations [`next_direction`](@ref) supports."
const MAX_DIRECTIONS_DEPTH = 32

"""
Internal type implementing random directions. Draw a new value with `rand`, see
[`next_direction`](@ref).

Serves two purposes: a fixed value of `Directions` is useful for unit testing, and drawing a
single bit flag collection economizes on the RNG cost.
"""
struct Directions
    flags::UInt32
end

Base.rand(rng::AbstractRNG, ::Type{Directions}) = Directions(rand(rng, UInt32))

"""
$(SIGNATURES)

Return the next direction flag and the new state of directions. Results are undefined for
more than [`MAX_DIRECTIONS_DEPTH`](@ref) updates.
"""
function next_direction(directions::Directions)
    @unpack flags = directions
    Bool(flags & 0x01), Directions(flags >>> 1)
end

####
#### Trajectory interface
####

"""
    $(FUNCTIONNAME)(trajectory, z, is_forward)

Move along the trajectory in the specified direction. Return the new position.
"""
function move end

"""
    $(FUNCTIONNAME)(trajectory, τ)

Test if the turn statistics `τ` indicate that the corresponding tree is turning.

Will only be called on nontrivial trees (at least two nodes).
"""
function is_turning end

"""
    $(FUNCTIONNAME)(trajectory, τ₁, τ₂)

Combine turn statistics on trajectory. Implementation can assume that the trees that
correspond to the turn statistics have the same ordering.

When
```julia
τ = combine_turn_statistics(trajectory, τ₁, τ₂)
is_turning(trajectory, τ)
```
the combined turn statistic `τ` is guaranteed not to escape the caller, so it can eg change
type.
"""
function combine_turn_statistics end

"""
    $(FUNCTIONNAME)(trajectory, v₁, v₂)

Combine visited node statistics for adjacent trees trajectory. Implementation should be
invariant to the ordering of `v₁` and `v₂` (ie the operation is commutative).
"""
function combine_visited_statistics end

"""
    $(FUNCTIONNAME)(trajectory, is_doubling::Bool, ω₁, ω₂, ω)

Calculate the log probability if selecting the subtree corresponding to `ω₂`. Being the log
of a probability, it is always `≤ 0`, but implementations are allowed to return and accept
values `> 0` and treat them as `0`.

When `is_doubling`, the tree corresponding to `ω₂` was obtained from a doubling step (this
can be relevant eg for biased progressive sampling).

The value `ω = logaddexp(ω₁, ω₂)` is provided for avoiding redundant calculations.

See [`biased_progressive_logprob2`](@ref) for an implementation.
"""
function calculate_logprob2 end

"""
    $(FUNCTIONNAME)(rng, trajectory, ζ₁, ζ₂, logprob2::Real, is_forward::Bool)

Combine two proposals `ζ₁, ζ₂` on `trajectory`, with log probability `logprob2` for
selecting `ζ₂`.

 `ζ₁` is before `ζ₂` iff `is_forward`.
"""
function combine_proposals end

"""
    ζωτ_or_nothing, v = $(FUNCTIONNAME)(trajectory, z, is_initial)

Information for a tree made of a single node. When `is_initial == true`, this is the first
node.

The first value is either

1. `nothing` for a divergent node,

2. a tuple containing the proposal `ζ`, the log weight (probability) of the node `ω`, the
turn statistics `τ` (never tested with `is_turning` for leaf nodes).

The second value is the visited node information.
"""
function leaf end

####
#### utilities
####

"""
$(SIGNATURES)

Combine turn statistics with the given direction. When `is_forward`, `τ₁` is before `τ₂`,
otherwise after.

Internal helper function.
"""
@inline function combine_turn_statistics_in_direction(trajectory, τ₁, τ₂, is_forward::Bool)
    if is_forward
        combine_turn_statistics(trajectory, τ₁, τ₂)
    else
        combine_turn_statistics(trajectory, τ₂, τ₁)
    end
end

function combine_proposals_and_logweights(rng, trajectory, ζ₁, ζ₂, ω₁::Real, ω₂::Real,
                                          is_forward::Bool, is_doubling::Bool)
    ω = logaddexp(ω₁, ω₂)
    logprob2 = calculate_logprob2(trajectory, is_doubling, ω₁, ω₂, ω)
    ζ = combine_proposals(rng, trajectory, ζ₁, ζ₂, logprob2, is_forward)
    ζ, ω
end

"""
$(SIGNATURES)

Given (relative) log probabilities `ω₁` and `ω₂`, return the log probabiliy of
drawing a sample from the second (`logprob2`).

When `bias`, biases towards the second argument, introducing anti-correlations.
"""
function biased_progressive_logprob2(bias::Bool, ω₁::Real, ω₂::Real, ω = logaddexp(ω₁, ω₂))
    ω₂ - (bias ? ω₁ : ω)
end

####
#### abstract trajectory interface
####

"""
$(SIGNATURES)

Information about an invalid (sub)tree, using positions relative to the starting node.

1. When `left < right`, this tree was *turning*.

2. When `left == right`, this is a *divergent* node.

3. `left == 1 && right == 0` is used as a sentinel value for reaching maximum depth without
encountering any invalid trees (see [`REACHED_MAX_DEPTH`](@ref). All other `left > right`
values are disallowed.
"""
struct InvalidTree
    left::Int
    right::Int
end

InvalidTree(i::Integer) = InvalidTree(i, i)

is_divergent(invalid_tree::InvalidTree) = invalid_tree.left == invalid_tree.right

function Base.show(io::IO, invalid_tree::InvalidTree)
    msg = if is_divergent(invalid_tree)
        "divergence at position $(invalid_tree.left)"
    elseif invalid_tree == REACHED_MAX_DEPTH
        "reached maximum depth without divergence or turning"
    else
        @unpack left, right = invalid_tree
        "turning at positions $(left):$(right)"
    end
    print(io, msg)
end

"Sentinel value for reaching maximum depth."
const REACHED_MAX_DEPTH = InvalidTree(1, 0)

"""
    result, v = adjacent_tree(rng, trajectory, z, i, depth, is_forward)

Traverse the tree of given `depth` adjacent to point `z` in `trajectory`.

`is_forward` specifies the direction, `rng` is used for random numbers in
[`combine_proposals`](@ref). `i` is an integer position relative to the initial node (`0`).

The *first value* is either

1. an `InvalidTree`, indicating the first divergent node or turning subtree that was
encounteted and invalidated this tree.

2. a tuple of `(ζ, ω, τ, z′, i′), with

    - `ζ`: the proposal from the tree.

    - `ω`: the log weight of the subtree that corresponds to the proposal

    - `τ`: turn statistics

    - `z′`: the last node of the tree

    - `i′`: the position of the last node relative to the initial node.

The *second value* is always the visited node statistic.
"""
function adjacent_tree(rng, trajectory, z, i, depth, is_forward)
    i′ = i + (is_forward ? 1 : -1)
    if depth == 0
        z′ = move(trajectory, z, is_forward)
        ζωτ, v = leaf(trajectory, z′, false)
        if ζωτ ≡ nothing
            InvalidTree(i′), v
        else
            (ζωτ..., z′, i′), v
        end
    else
        # “left” tree
        t₋, v₋ = adjacent_tree(rng, trajectory, z, i, depth - 1, is_forward)
        t₋ isa InvalidTree && return t₋, v₋
        ζ₋, ω₋, τ₋, z₋, i₋ = t₋

        # “right” tree — visited information from left is kept even if invalid
        t₊, v₊ = adjacent_tree(rng, trajectory, z₋, i₋, depth - 1, is_forward)
        v = combine_visited_statistics(trajectory, v₋, v₊)
        t₊ isa InvalidTree && return t₊, v
        ζ₊, ω₊, τ₊, z₊, i₊ = t₊

        # turning invalidates
        τ = combine_turn_statistics_in_direction(trajectory, τ₋, τ₊, is_forward)
        is_turning(trajectory, τ) && return InvalidTree(i′, i₊), v

        # valid subtree, combine proposals
        ζ, ω = combine_proposals_and_logweights(rng, trajectory, ζ₋, ζ₊, ω₋, ω₊,
                                                is_forward, false)
        (ζ, ω, τ, z₊, i₊), v
    end
end

"""
$(SIGNATURES)

Sample a `trajectory` starting at `z`, up to `max_depth`. `directions` determines the tree
expansion directions.

Return the following values

- `ζ`: proposal from the tree

- `v`: visited node statistics

- `termination`: an `InvalidTree` (this includes the last doubling step turning, which is
  technically a valid tree) or `REACHED_MAX_DEPTH` when all subtrees were valid and no
  turning happens.

- `depth`: the depth of the tree that was sampled from. Doubling steps that lead to an
  invalid adjacent tree do not contribute to `depth`.
"""
function sample_trajectory(rng, trajectory, z, max_depth::Integer, directions::Directions)
    @argcheck max_depth ≤ MAX_DIRECTIONS_DEPTH
    (ζ, ω, τ), v = leaf(trajectory, z, true)
    z₋ = z₊ = z
    depth = 0
    termination = REACHED_MAX_DEPTH
    i₋ = i₊ = 0
    while depth < max_depth
        is_forward, directions = next_direction(directions)
        t′, v′ = adjacent_tree(rng, trajectory, is_forward ? z₊ : z₋, is_forward ? i₊ : i₋,
                               depth, is_forward)
        v = combine_visited_statistics(trajectory, v, v′)

        # invalid adjacent tree: stop
        t′ isa InvalidTree && (termination = t′; break)

        # extract information from adjacent tree
        ζ′, ω′, τ′, z′, i′ = t′

        # update edges and combine proposals
        if is_forward
            z₊, i₊ = z′, i′
        else
            z₋, i₋ = z′, i′
        end

        # tree has doubled successfully
        ζ, ω = combine_proposals_and_logweights(rng, trajectory, ζ, ζ′, ω, ω′,
                                                is_forward, true)
        depth += 1

        # when the combined tree is turning, stop
        τ = combine_turn_statistics_in_direction(trajectory, τ, τ′, is_forward)
        is_turning(trajectory, τ) && (termination = InvalidTree(i₋, i₊); break)
    end
    ζ, v, termination, depth
end


text:
# logistic regression on mnist data

using MultiStochGrad

using  LinearAlgebra, Test
using Random, Distributions
using Statistics

include("../src/applis/LogisticRegression.jl")
include("../src/mnist.jl")

const IMAGE_FNAME_STR  = "/home.1/jpboth/Data/MNIST/train-images-idx3-ubyte"
const LABEL_FNAME_STR  = "/home.1/jpboth/Data/MNIST/train-labels-idx1-ubyte"

function initMnitObservations()
    # check path and load
    if !isfile(IMAGE_FNAME_STR) || !isfile(LABEL_FNAME_STR)
        @warn("\n\n mnist_logistic_regression : bad path to mnist data file, cannot do tests, download and install files")
        return nothing
    end
    mnistdata = mnist.MnistData(IMAGE_FNAME_STR, LABEL_FNAME_STR)
    # transform data into logistic regression problem
    # labels are numbered from 0 to 9.
    datas = Vector{Vector{Float64}}()
    values = Vector{Float64}()
    images = mnistdata.images
    dims = size(images)
    nbimages = dims[3]
    nbpixels = dims[1] * dims[2]
    for i in 1:nbimages
        v = Vector{Float64}()
        push!(v, 1.)
        # reshape and convert to Float64
        v = cat(v, convert(Array{Float64}, reshape(images[:,:,i], nbpixels)), dims = 1)
        v = v / 256.
        push!(datas, v)
        push!(values, Float64(mnistdata.labels[i]))
    end
    observations = Observations(datas,values)
    Some(observations)
end

function mnist_logistic_regression_scsg()
    @info "\n\n\n begining mnist_logistic_regression_scsg ..."
    mnist = initMnitObservations()
    if isnothing(mnist)
        return true
    end
    #
    observations = something(mnist)
    nbclass = 10
    # struct LogisticRegression takes care of interceptions terms and constraint on first class
    # define TermFunction , TermGradient and Evaluator, dims is one less than number of classes for identifiability constraints
    dims = Dims{2}((length(observations.datas[1]) , nbclass-1))
    nbterms = length(observations.datas)
    #
    logistic_term_value = make_logistic_term_value(observations)
    term_function =  TermFunction{typeof(logistic_term_value)}(logistic_term_value, nbterms, dims)
    #
    logistic_term_gradient = make_logistic_term_gradient(observations)
    term_gradient = TermGradient{typeof(logistic_term_gradient)}(logistic_term_gradient ,nbterms, dims)
    #
    evaluator = Evaluator{typeof(logistic_term_value),typeof(logistic_term_gradient)}(term_function, term_gradient)
    # define parameters for scsg
    scsg_pb = SCSG(0.1, 0.006, 1 , 0.02)
    @debug "scsg_pb" scsg_pb
    # solve.
    nb_iter = 70
    initial_position= fill(0.5, dims)
    initial_error = compute_value(evaluator, initial_position)
    @info "initial error " initial_error
    @time position, value = minimize(scsg_pb, evaluator, nb_iter, initial_position)
    @info "value reached : " value
    @info "\n\n mnist_logistic_regression_scsg done \n"
    value < 0.3 ? true : false
end

function mnist_logistic_regression_svrg()
    @info "\n\n\n begining mnist_logistic_regression_svrg ..."
    mnist = initMnitObservations()
    if isnothing(mnist)
        return true
    end
    #
    observations = something(mnist)
    nbclass = 10
    # struct LogisticRegression takes care of interceptions terms and constraint on first class
    # define TermFunction , TermGradient and Evaluator, dims is one less than number of classes for identifiability constraints
    dims = Dims{2}((length(observations.datas[1]) , nbclass-1))
    nbterms = length(observations.datas)
    #
    logistic_term_value = make_logistic_term_value(observations)
    term_function =  TermFunction{typeof(logistic_term_value)}(logistic_term_value, nbterms, dims)
    #
    logistic_term_gradient = make_logistic_term_gradient(observations)
    term_gradient = TermGradient{typeof(logistic_term_gradient)}(logistic_term_gradient ,nbterms, dims)
    #
    evaluator = Evaluator{typeof(logistic_term_value),typeof(logistic_term_gradient)}(term_function, term_gradient)
    # define parameters for svrg  1000 minibatch , step 0.02
    svrg_pb = SVRG(2000, 0.05)
    # solve.
    nb_iter = 60
    initial_position= fill(0.5, dims)
    initial_error = compute_value(evaluator, initial_position)
    @info "initial error " initial_error
    @time position, value = minimize(svrg_pb, evaluator, nb_iter, initial_position)
    @info "value reached " value
    @info "\n\n mnist_logistic_regression_svrg done \n"
    value < 0.3 ? true : false
end


text:
using OrdinaryDiffEq, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

prob = prob_ode_linear

sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=false)
sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[1/2])

@test sort!(symdiff(sol.t,sol2.t)) == [0.0,1/2,1.0]

sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[0.0,1/2,1.0])

@test symdiff(sol.t,sol2.t) == [1/2]

sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=1/2)

@test sort!(symdiff(sol.t,sol2.t)) == [1/2]

sol3=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[1/2],tstops=[1/2])

@test sol3.t == [1/2]

sol3=solve(prob,DP5(),dt=1//2^(2),saveat=[0.0,1/2,1.0],tstops=[1/2])

@test sol3.t == [0.0,1/2,1.0]

sol3=solve(prob,DP5(),dt=1//2^(2),saveat=1/10,tstops=[1/2])

@test sol3.t == collect(0.0:0.1:1.00)

sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false)
sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true)
sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true)
sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

prob = prob_ode_2Dlinear

sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=true)
sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=true,saveat=[0.0,1/2,1.0])

@test symdiff(sol.t,sol2.t) == [1/2]

sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false)
sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false,saveat=[0.0,.125,.6,.61,.8,1.0])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true)
sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])

@test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=false)
sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),saveat=[.125,.6,.61,.8])

@test sort!(symdiff(sol.t,sol2.t)) == [0.0,.125,.6,.61,.8,1.0]

sol=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[0,.125,.6,.61,.8])

@test !(sol.t[2] ≈ 0)

# Test Iterators

sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=0:1//100:1)

@test sol2.t ≈ collect(0:1//100:1)

sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=range(0,stop=1,length=100))

@test sol2.t ≈ range(0,stop=1,length=100)

f = (du,u,p,t) -> prob.f(du,u,p,t)
prob2 = ODEProblem(f,vec(prob.u0),prob.tspan,1.01)

sol2=solve(prob2,DP5(),dt=1//2^(2),saveat=.1,save_idxs=1:2:5)

for u in sol2.u
  @test length(u) == 3
end

sol2=solve(prob2,DP5(),dt=1//2^(2),saveat=.1,save_idxs=1:2:5,save_everystep=true)

sol=solve(prob2,DP5(),dt=1//2^(2),save_start=false)

@test sol.t[1] == 1//2^(2)

# Test save_on switch
sol = solve(prob, DP5(), save_on=false, save_start=false, save_end=false)
@test isempty(sol.t) && isempty(sol.u)
sol = solve(prob, DP5(), saveat=0.2, save_on=false, save_start=false, save_end=false)
@test isempty(sol.t) && isempty(sol.u)


text:
using LocalBinaryPatterns
using ImageCore
using ImageCore.OffsetArrays
using IntegralArrays
using Interpolations
using Statistics
using Test
using Aqua, Documenter

@testset "LocalBinaryPatterns.jl" begin

@testset "meta-quality" begin
    Aqua.test_ambiguities(LocalBinaryPatterns)
    Aqua.test_all(LocalBinaryPatterns; ambiguities=false)
    if VERSION <= v"1.8.0-DEV.1070"
        # TODO(johnnychen94): remove this version check when `@_inline_meta` deprecation is handled in upstream ecosystems
        doctest(LocalBinaryPatterns; manual = false)
    end
end

include("bitrotate_encoding.jl")
include("uniform_encoding.jl")
include("modes.jl")
include("local_binary_pattern.jl")

end


text:
"""
    branching_order!(mtg; ascend = true)

Compute the topological branching order of the nodes in an mtg.

# Arguments

- `mtg`: the mtg, *e.g.* output from `read_mtg()`
- `ascend`: If `true`, the order is computed from the base (acropetal), if `false`,
it is computed from the tip (basipetal).

# Notes

The order of a node is computed from the maximum order of their children when using the
basipetal computation.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
branching_order!(mtg)
DataFrame(mtg, :branching_order)
# 7×2 DataFrame
#  Row │ tree                        branching_order
#      │ String                      Int64
# ─────┼───────────────────────────────────────────────
#    1 │ / 1: \$                                      1
#    2 │ └─ / 2: Individual                          1
#    3 │    └─ / 3: Axis                             1
#    4 │       └─ / 4: Internode                     1
#    5 │          ├─ + 5: Leaf                       2
#    6 │          └─ < 6: Internode                  1
#    7 │             └─ + 7: Leaf                    2

branching_order!(mtg, ascend = false)
DataFrame(mtg, :branching_order)
# 7×2 DataFrame
#  Row │ tree                        branching_order
#      │ String                      Int64
# ─────┼───────────────────────────────────────────────
#    1 │ / 1: \$                                      2
#    2 │ └─ / 2: Individual                          2
#    3 │    └─ / 3: Axis                             2
#    4 │       └─ / 4: Internode                     2
#    5 │          ├─ + 5: Leaf                       1
#    6 │          └─ < 6: Internode                  2
#    7 │             └─ + 7: Leaf                    1
```
"""
function branching_order!(mtg; ascend = true)

    if ascend
        @mutate_mtg!(mtg, branching_order = branching_order_ascend!(node))
    else
        branching_order_descend!(mtg)
    end
end

function branching_order_ascend!(node)
        if isroot(node)
        return 1
    else
        if node.MTG.link == "+"
            return parent(node)[:branching_order] + 1
        else
            return parent(node)[:branching_order]
        end
    end
end

function branching_order_descend!(node)

    if isleaf(node)
        val = 1
    else
        val_child = Int[]
        for chnode in ordered_children(node)
            push!(val_child, branching_order_descend!(chnode))
        end
        val = maximum(val_child)

        if node.MTG.link == "+"
            val += 1
        end
    end

    node[:branching_order] = val

    return val
end


text:
# example 3: read all at once with readfasta
#            default return type (String)

module FastaIOExample3
using FastaIO

function read_fasta_file(filename::String)
    out = readfasta(filename)

    n = 0
    for (desc, seq) in out
        n += 1
        println("num=$n desc=$desc seq=$seq")
    end
    # here of course n == fr.num_parsed == length(out)
    println("read $n entries")
end

read_fasta_file(joinpath(@__DIR__, "example.fasta.gz"))

end # module


text:
function galilean_search!(m,e,t,llidx,samples,s_limit,cache)
    new_m,cache=galilean_trajectory_sample!(m,e,t,cache)
    m.i==1 && ((m,new_m,cache)=start_trim!(m,new_m,cache))

    while new_m.i==m.i && samples<=s_limit && t.τ > e.GMC_τ_death#failure to find any new step, retune tau and try again until sample limit reached or trajectory's tuner reaches taudeath
        new_m,cache=galilean_trajectory_sample!(m,e,t,cache)
        samples+=1
        m.i==1 && ((m,new_m,cache)=start_trim!(m,new_m,cache))
    end

    if new_m.i==m.i+1
        e.sample_posterior && push!(e.posterior_samples, e.models[llidx])#if sampling posterior, push the model record to the ensemble's posterior samples vector

        rectype=typeof(e.models[1])
        deleteat!(e.models, llidx) #remove least likely model record

        new_model_record = rectype(new_m.trajectory,new_m.i, new_m.pos, string(e.path,'/',new_m.trajectory,'.',new_m.i), new_m.log_Li)
        push!(e.models, new_model_record) #insert new model record
        serialize(new_model_record.path, new_m)

        return true, cache
    else
        return false, cache
    end
end

function start_trim!(m,new_m,cache)
    if new_m.i!=m.i && isapprox(new_m.pos,m.pos)
        params=[getfield(m,pn) for pn in propertynames(m)]
        params[6]=new_m.v
        m=typeof(m)(params...)
        params=[getfield(cache,pn) for pn in propertynames(m)]
        params[2]=2;
        new_m=typeof(m)(params...)
        cache=nothing
    end
    return m,new_m,cache
end

function insert_reflection!(e, cache)
    t_idx=findfirst(t->t==cache.trajectory,[m.trajectory for m in e.models])

    e.sample_posterior && push!(e.posterior_samples, e.models[t_idx])#if sampling posterior, push the model record to the ensemble's posterior samples vector
    rectype=typeof(e.models[1])
    deleteat!(e.models, t_idx) #remove least likely model record

    new_model_record = rectype(cache.trajectory,cache.i, cache.pos, string(e.path,'/',cache.trajectory,'.',cache.i), cache.log_Li)
    push!(e.models, new_model_record) #insert new model record
    serialize(new_model_record.path, cache)
end

function diffusion_search!(e, tuners, llidx, sample_limit, τ_limit)
    sample=0; model_selected=false; new_m=0; src_traj=0
    while !(sample > sample_limit) && !model_selected
        sample +=1
        mrec=rand(e.models)
        model=deserialize(mrec.path)
        randdir=rand(MvNormal(length(model.θ),1.))
        τ=tuners[model.trajectory].τ
        new_pos=model.pos+randdir*τ
        box_bound!(new_pos,e.box)
        new_m=e.model_initλ(e.t_counter, 1, to_prior.(new_pos,e.priors), new_pos, randdir, e.obs, e.constants...)

        while τ*.7>τ_limit && new_m.log_Li<e.contour
            τ*=.7
            randdir=rand(MvNormal(length(model.θ),1.))
            new_pos=model.pos+randdir*τ
            box_bound!(new_pos,e.box)
            new_m=e.model_initλ(e.t_counter, 1, to_prior.(new_pos,e.priors), new_pos, randdir, e.obs, e.constants...)
        end

        new_m.log_Li >= e.contour && (model_selected=true; src_traj=model.trajectory)
    end

    if new_m.log_Li >= e.contour
        e.sample_posterior && push!(e.posterior_samples, e.models[llidx])#if sampling posterior, push the model record to the ensemble's posterior samples vector
        rectype=typeof(e.models[1])
        deleteat!(e.models, llidx) #remove least likely model record

        new_model_record = rectype(new_m.trajectory,new_m.i, new_m.pos, string(e.path,'/',new_m.trajectory,'.',new_m.i), new_m.log_Li)
        push!(e.models, new_model_record) #insert new model record
        serialize(new_model_record.path, new_m)
        tuners[e.t_counter]=τ_PID(e)
        tuners[e.t_counter].τ=tuners[src_traj].τ
        e.t_counter +=1
        return true
    else
        return false
    end
end

function ellipsoid_resample!(e,tdict,llidx,samples,s_limit)
    posmat=zeros(length(e.models[1].pos),length(e.models))
    for (n,rec) in enumerate(e.models)
        posmat[:,n]=rec.pos
    end
    e_ellip=bounding_ellipsoid(posmat)

    new_pos=sample_ellipsoid(e_ellip)
    box_bound!(new_pos,e.box)
    new_m=e.model_initλ(e.t_counter, 1, to_prior.(new_pos,e.priors), new_pos, [0.], e.obs, e.box, e.constants..., v_init=true)
    samples+=1

    while new_m.log_Li < e.contour && samples < s_limit
        new_pos=sample_ellipsoid(e_ellip)
        box_bound!(new_pos,e.box)
        new_m=e.model_initλ(e.t_counter, 1, to_prior.(new_pos,e.priors), new_pos, [0.], e.obs, e.box, e.constants..., v_init=true)
        samples+=1
    end

    if new_m.log_Li > e.contour
        e.sample_posterior && push!(e.posterior_samples, e.models[llidx])#if sampling posterior, push the model record to the ensemble's posterior samples vector
        rectype=typeof(e.models[1])
        deleteat!(e.models, llidx) #remove least likely model record

        new_model_record = rectype(new_m.trajectory,new_m.i, new_m.pos, string(e.path,'/',new_m.trajectory,'.',new_m.i), new_m.log_Li)
        push!(e.models, new_model_record) #insert new model record
        serialize(new_model_record.path, new_m)
        tdict[e.t_counter]=τ_PID(e)
        e.t_counter +=1
        return true
    else
        return false
    end
end


text:
#=
- glider http://jsxgraph.org/docs/symbols/Glider.html
- group http://jsxgraph.org/docs/symbols/Group.html
- grid http://jsxgraph.org/docs/symbols/Grid.html
- integral http://jsxgraph.org/docs/symbols/Integral.html
=#


text:
##
# This file is a part of Dagger.jl. License is MIT
#
# Based upon https://github.com/google/benchmark, which is licensed under Apache v2:
# https://github.com/google/benchmark/blob/master/LICENSE
#
# In compliance with the Apache v2 license, here are the original copyright notices:
# Copyright 2015 Google Inc. All rights reserved.
##

struct TimeSpec
    tv_sec  :: UInt64 # time_t
    tv_nsec :: UInt64
end

maketime(ts) = ts.tv_sec * UInt(1e9) + ts.tv_nsec

# From bits/times.h on a Linux system
# Check if those are the same on BSD
if Sys.islinux()
    const CLOCK_MONOTONIC          = Cint(1)
    const CLOCK_PROCESS_CPUTIME_ID = Cint(2)
    const CLOCK_THREAD_CPUTIME_ID  = Cint(3)
elseif Sys.KERNEL == :FreeBSD # atleast on FreeBSD 11.1
    const CLOCK_MONOTONIC          = Cint(4)
    const CLOCK_PROCESS_CPUTIME_ID = Cint(14)
elseif Compat.Sys.isapple() # Version 10.12 required
    const CLOCK_MONOTONIC          = Cint(6)
    const CLOCK_PROCESS_CPUTIME_ID = Cint(12)
else
    error("""
          BenchmarkTools doesn't currently support your operating system.
          Please file an issue, your kernel is $(Sys.KERNEL)
          """)
end

@inline function clock_gettime(cid)
    ts = Ref{TimeSpec}()
    ccall(:clock_gettime, Cint, (Cint, Ref{TimeSpec}), cid, ts)
    return ts[]
end

@inline function realtime()
    maketime(clock_gettime(CLOCK_MONOTONIC))
end

@inline function cputime()
    maketime(clock_gettime(CLOCK_PROCESS_CPUTIME_ID))
end

@inline function cputhreadtime()
    maketime(clock_gettime(CLOCK_THREAD_CPUTIME_ID))
end

struct Measurement
    realtime::TimeSpec
    cputime::TimeSpec
    function Measurement()
        rtime = clock_gettime(CLOCK_MONOTONIC)
        ctime = clock_gettime(CLOCK_PROCESS_CPUTIME_ID)
        return new(rtime, ctime)
    end
end

struct MeasurementDelta
    realtime::Float64
    cpuratio::Float64
    function MeasurementDelta(t1::Measurement, t0::Measurement)
        rt0 = maketime(t0.realtime)
        ct0 = maketime(t0.cputime)
        rt1 = maketime(t1.realtime)
        ct1 = maketime(t1.cputime)
        realtime = rt1 - rt0
        cputime = ct1 - ct0
        return new(realtime, cputime/realtime)
    end
end


text:
module Gym

using Flux
using Flux.Tracker

import Base.show

export step!, reset!, Ctx, render

#Spaces
include("Spaces/Space.jl")
using .Space
export sample,
    Box, Discrete, TupleSpace, DictSpace, MultiBinary, MultiDiscrete
#Environments
export CartPoleEnv, PendulumEnv, Continuous_MountainCarEnv

include("vis/utils.jl")
include("CartPole.jl")
include("Pendulum.jl")
include("Continuous-MountainCar.jl")

end #module


text:
#create a random transition matrix of n size with rational probabilities with denominator k
function RandomTransitionMatrix(n::Int64,k::Int64)
  RandM=rand(DiscreteUniform(1,n),(n,k))
  RandK=sortcols(RandM[1,1:k])
  for i in 2:n
    RandK=vcat(RandK,sortcols(RandM[i,1:k]))
  end

  Rnd=Array(Float64,n,n)
  for i in 1:n
    for j in 1:n
      Rnd[i,j]=count(y->y==j,[x::Int64 for x in sortrows(RandK)[i,1:k]])
    end
  end
  return Rnd/k
end

function RandomGraphMatrix(n::Int64)
  Mat=eye(Int64,n)
  Mat[1,1:n]=ones(Int64,n)
  Mat[1:n,n]=ones(Int64,n)
  for i in 2:(n-2)
    for j in (i+1):n-1
      Mat[i,j]=rand(Binomial())
    end
  end

  for i in 2:n-2
    for j in 1:(i-1)
      if Mat[n-i,n-j]==1
        for k in (n-j+1):(n-1)
          Mat[n-i,k]=max(Mat[n-i,k],Mat[n-j,k])
        end
      end
    end
  end
  return Mat
end

#for some matrices C and P check if Markov Process is Mobius down monotone

function is_mobius_down(C::Array{Int64,2},P::Array{Float64,2})
  R=convert_to_min_Int(P)
  return minimum(inv(C)*R*C .>= 0)
end

#for some matrices C and P check if Markov Process is Mobius up monotone
function is_mobius_up(C::Array{Int64,2},P::Array{Float64,2})
  R=convert_to_min_Int(P)
  return minimum(inv(C')*R*C' .>= 0)
end

#for some matrix C which represent poset give predecessors matrix
function predecessors_matrix(C::Array{Int64,2})
  C1=deepcopy(C)
  d=spzeros(1, size(C1)[1])
  for i=1:size(C1)[1]
    for k=i+1:size(C1)[1]
      if C1[i,k]==1
        C1[i,1:size(C1)[1]]-=C1[k,1:size(C1)[1]]
        d=vcat(d,sparse([1,1,1],[i,k,size(C1)[1]],[1,-1,0]))
      end
    end
  end
  return dense(d)[2:size(d)[1],[1:size(d)[2]]]
end

#for some matrices C and P check if Markov Process is weak down monotone
function is_weak_down(C::Array{Int64,2},P::Array{Float64,2})
  R=convert_to_min_Int(P)
  return minimum(predecessors_matrix(C)*R*C .>= 0)
end

#for some matrices C and P check if Markov Process is weak up monotone
function is_weak_up(C::Array{Int64,2},P::Array{Float64,2})
  R=convert_to_min_Int(P)
  return minimum(predecessors_matrix(C)*R*C' .<= 0)
end

#extended matrix C of up sets with minimal element to matrix of all up sets
function up_sets(C::Array{Int64,2})
  d=zeros(1, size(C)[1])
  for i=1:size(C)[1]
    D=deepcopy(C[i,1:size(C)[1]])
    for k=(i+1):size(C)[1]
      if D[1,k]==0
        for j=1:size(D)[1]
          if D[j,k]==0
            D=vcat(D,D[j,1:size(C)[1]]+C[k,1:size(C)[1]])
          end
        end
      end
    end
    d=vcat(d,D)
  end
  return sign(d[2:size(d)[1],[1:size(d)[2]]])'
end

#for some matrices C and P check if Markov Process is stochastically monotone
function is_stochastic(C::Array{Int64,2},P::Array{Float64,2})
  C1=deepcopy(C)
  R=convert_to_min_Int(P)
  return minimum(predecessors_matrix(C1)*R*up_sets(C) .<= 0)
end

#multiplying elements of matrix by lcm of their denominators geting matrix with (small) integral elements
function convert_to_min_Int(P::Array{Float64,2})
  return convert(Array{Int64,2},lcm(map(den,convert(Array{Rational,2},P)))*P)
end

#from transition matrix P give any realization (in rows we have transitions on some range to some state)
function convert_to_min_Realization(P::Array{Float64,2})
  R=convert_to_min_Int(P)
  F=convert(Array{Int64,2},zeros(1, 1+sum(R[1,1:size(R)[2]])))
  for i=1:size(R)[1]
    f=[0]
    for j=1:size(R)[2]
      for k=1:R[i,j]
        f=hcat(f,j)
      end
    end
    F=vcat(F,f)
  end
  return F[2:size(F)[1],2:size(F)[2]]
end

#all posible realizations from transition matrix P (=all permutations on each state of any realization)
function all_Functions(P::Array{Float64,2})
  F=convert_to_min_Realization(P)

  # Perm=unique(permutations([x::Int64 for x in F[1,1:size(F)[2]]]))
# for i=2:size(F)[1]
#   Perm=vcat(Perm,unique(permutations([x::Int64 for x in F[i,1:size(F)[2]]])),(1,4))
  #end

  Perm = [unique(permutations([x::Int64 for x in F[1,i:size(F)[2]]])) for i = 1:size(F)[1]]
  return Perm
end

#checking if vector is not less that other considering ordering of matrix C
function stoch_not_less(vec1::Array{Int64,2},vec2::Array{Int64,2},C::Array{Int64,2})
  for i=1:length(vec1)
    if C[vec1[i],vec2[i]]!=1
      return true
    end
  end
  return false
end

#checking if Realization_function[j] is posible function on state 'j' (=is bigger then other, i<j, functions)
function check_vector_j(Realization_function::Array{Int64,2},j::Int64,C::Array{Int64,2})
  for i=1:j-1
    if stoch_not_less(Realization_function[i,1:size(Realization_function)[2]],Realization_function[j,1:size(Realization_function)[2]],C)
      return false
    end
  end
  return true
end

#searching for posible function on state 'j' when we know previous, i<j, functions
function Realization_j(j::Int64,Realization_function::Array{Int64,2},C::Array{Int64,2},aF)
  if j>size(aF)[1]
    return true
  end
  state=start(aF[j])
  while !done(aF[j], state)
    (i, state) = next(aF[j], state)
    Realization_function[j,1:length(i)]=i'
    if check_vector_j(Realization_function,j,C)
      if Realization_j(j+1,Realization_function,C,aF)==true
        return true
      end
    end
  end
  return false
end

#for some matrices C and P check if there exist (some kind of) realization of Markov Process
function is_realizable(C::Array{Int64,2},P::Array{Float64,2})
  aF=all_Functions(P)
  Realization_function=convert_to_min_Realization(P)
  if Realization_j(2,Realization_function,C,aF)==true
    return true
  end
  return false
end

#for some matrices C and P check all definitions of monotonicity
function which_is(C::Array{Int64,2},P::Array{Float64,2})
  return ["mobius_down" "mobius_up" "weak_down" "weak_up" "stochastic" "realizable"; is_mobius_down(C,P) is_mobius_up(C,P) is_weak_down(C,P) is_weak_up(C,P) is_stochastic(C,P) is_realizable(C,P)]'
end

#for some matrices C and P check all definitions of monotonicity except realizable monotonicity
function which_is_not_real(C::Array{Int64,2},P::Array{Float64,2})
  return ["mobius_down" "mobius_up" "weak_down" "weak_up" "stochastic"; is_mobius_down(C,P) is_mobius_up(C,P) is_weak_down(C,P) is_weak_up(C,P) is_stochastic(C,P)]'
end

function check_stoch_mob_not_real(C,P)
  if is_mobius_up(C,P) if is_mobius_down(C,P) if is_stochastic(C,P) return !(is_realizable(C,P))
      end end end
  return false
end

function searching_for_stoch_mob_not_real(N::Int64, n:Int64)
  for i in 1:N
    m=rand(7:12)
    C=RandomGraphMatrix(n)
    for j in 1:1000
      P=RandomTransitionMatrix(n,m)
      if check_stoch_mob_not_real(C,P)
        return (C,P)
      end
    end
  end
  return false
end

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
## USER FRIENDLY INSPECTION OF COMPOSITE MACHINES

try_scalarize(v) = length(v) == 1 ? v[1] : v

function machines_given_model_name(mach::Machine{M}) where  M<:Composite
    network_model_names = mach.cache.network_model_names
    names = unique(filter(name->!(name === nothing), network_model_names))
    network_models = MLJBase.models(glb(mach))
    network_machines = MLJBase.machines(glb(mach))
    ret = LittleDict{Symbol,Any}()
    for name in names
        mask = map(==(name), network_model_names)
        _models = network_models[mask]
        _machines = filter(mach->mach.model in _models, network_machines)
        ret[name] = _machines
    end
    return ret
end

function tuple_keyed_on_model_names(item_given_machine, mach)
    dict = MLJBase.machines_given_model_name(mach)
    names = tuple(keys(dict)...)
    named_tuple_values = map(names) do name
        [item_given_machine[m] for m in dict[name]] |> try_scalarize
    end
    return NamedTuple{names}(named_tuple_values)
end

function report(mach::Machine{<:Composite})
    machines = mach.report.machines
    dict = mach.report.report_given_machine
    return merge(tuple_keyed_on_model_names(dict, mach), mach.report)
end

function fitted_params(mach::Machine{<:Composite})
    fp = fitted_params(mach.model, mach.fitresult)
    _machines = fp.machines
    dict = fp.fitted_params_given_machine
    return merge(MLJBase.tuple_keyed_on_model_names(dict, mach), fp)
end


text:
using Pkg
Pkg.add(PackageSpec(url="https://github.com/sdwfrost/Gillespie.jl", rev="master"))

using Gillespie
using Random
using StatsPlots
using BenchmarkTools

function sir_rates(x,parms)
  (S,I,R) = x
  (β,c,γ) = parms
  N = S+I+R
  infection = β*c*I/N*S
  recovery = γ*I
  [infection,recovery]
end;

sir_transitions = [[-1 1 0];[0 -1 1]];

tmax = 40.0;

u0 = [990,10,0]; # S,I,R

p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

sol_jump = ssa(u0,sir_rates,sir_transitions,p,tmax);

df_jump = ssa_data(sol_jump);

@df df_jump plot(:time,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark ssa(u0,sir_rates,sir_transitions,p,tmax)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


text:
function filter_insertions(seed::Dict)

    tmpname = tempname() #PlmDCA takes a file as input
    #tmpname = "tmp.fasta"
    tmpfile = open(tmpname, "w")
    for name in keys(seed)
        @printf(tmpfile, ">%s\n", name)
	seq = seed[name]
        s = ""
        for i in 1:length(seq)
	    if isuppercase(seq[i]) || seq[i] == '-'
                s = s * seq[i]
            end
        end
        @printf(tmpfile, "%s\n", s)
    end
    flush(tmpfile)
    return tmpname

end

function print_results(filename_ins::String, l_o::Array{Float64,1}, l_e::Array{Float64,1}, filename_par::String, PlmData::Any, ctype::Symbol, L::Int64, filename_gap::String, muext::Float64, muint::Float64)

    println("### Printing insertions penalties in ", filename_ins, " ###")
    file_ins = open(filename_ins, "w")
    for i in 1:L
        @printf(file_ins, "%.4f\t%.4f\n", l_o[i], l_e[i])
    end
    flush(file_ins)
    filep = open(filename_par, "w")
    if ctype == :amino
        q = 21
    elseif ctype == :nbase
        q = 5
    end

    println("### Printing DCA model in ", filename_par, " ###")
    for i in 1:L, j in i+1:L, a in 1:q, b in 1:q
    	@printf(filep, "J %d %d %d %d %f\n", a,b,i,j,PlmData.Jtensor[a,b,i,j])
    end
    flush(filep)
    for i in 1:L, a in 1:q
    	@printf(filep, "h %d %d %f\n", a,i,PlmData.htensor[a,i])
    end
    flush(filep)

    println("### Printing gap penalties in ", filename_gap, " ###")
    fileg = open(filename_gap, "w")
    @printf(fileg, "%f %f\n", muext,muint)
    flush(fileg)

end

"""
function readfull(filefull; ctype::Symbol="amino", pos::Bool=true)
    ffull = FastaIO.FastaReader(filefull)
    dheader = Dict{String,String}()
    for (_name,seq) in ffull
        name=String(split(_name)[1])
        if occursin('|', name)
            (aux,name) = split(name, "|")
        end
        if occursin('/', name) && !pos
            (name,aux) = split(name,"/")
        end
        if !haskey(dheader,name)
            dheader[name] = seq
        end
    end
    close(ffull)
    dheader
end
"""

function extract_full_seq(full::Dict, seed::Dict, Mtest::Int64; ctype=:amino)

    dheader = Dict{String, String}()

    for (_name, seq) in seed
        name = String(split(_name)[1])
        if occursin('/', name) && ctype == :amino
            (name,aux) = split(name, "/")
        end
        try
            seqfull = full[name]
            if !haskey(dheader, name)
                dheader[name] = seqfull
            end
        catch
            println("Error: ", name, " not found in full length sequences")
        end
    end

    #println("tot full seq", length(full), "tot align seq", length(seed), "tot new dict", length(dheader))
    tmpname = tempname()
    fulltmp = open(tmpname, "w")

    if Mtest == length(seed)
        for (name, seq) in dheader
            @printf(fulltmp, ">%s\n", name)
            @printf(fulltmp, "%s\n", seq)
        end
    else
        idx = sample(1:length(dheader), Mtest, replace= false)
        names = collect(keys(dheader))
        for i in 1:Mtest
            a = idx[i]
            name = names[a]
            seq = dheader[name]
            @printf(fulltmp, ">%s\n", name)
            @printf(fulltmp, "%s\n", seq)
        end
    end
    flush(fulltmp)

    return tmpname
end

function compute_average_dist(trueseed::Dict, tmpseed::Dict)

    i = 1
    d = zeros(length(trueseed))
    hd, tmp1, tmp2, tmp3 = 0, 0, 0, 0
    for (_name,seq) in trueseed
        L = length(seq)
        if occursin('/',_name)
            (name, hit) = split(_name, "/")
            (b,en) = split(hit, "-")
            b = parse(Int64,b)
            en = parse(Int64, en)
        else
            b = 0
            en = L
            name = _name
        end
        th = 0.10*en
        for (tmpname, tmpseq) in tmpseed
            if occursin('/', tmpname)
                (tmpname, hit) = split(tmpname, "/")
                (tmpb, tmpen) = split(hit, "-")
                tmpb = parse(Int64, tmpb)
                tmpen = parse(Int64, tmpen)
            else
                tmpb = 0
                tmpen = L
            end
            if occursin(tmpname, name)
                if abs(b-tmpb) < th && abs(en - tmpen) < th
                    hd, tmp1, tmp2, tmp3 = DCAlign.hammingdist(seq, tmpseq)
                    continue
                end
            end
        end
        d[i] = hd / (L*1.0)
        i = i + 1
    end

    return mean(d)

end


text:
abstract type HyperPlaneType end

abstract type OptimalityCut <: HyperPlaneType end
abstract type FeasibilityCut <: HyperPlaneType end
abstract type LinearConstraint <: HyperPlaneType end
abstract type Infeasible <: HyperPlaneType end
abstract type Unbounded <: HyperPlaneType end

struct HyperPlane{H <: HyperPlaneType, T <: Real, A <: AbstractVector}
    δQ::A
    q::T
    id::Int

    function (::Type{HyperPlane})(δQ::AbstractVector, q::Real, id::Int, ::Type{H}) where H <: HyperPlaneType
        T = promote_type(eltype(δQ), Float32)
        δQ_ = convert(AbstractVector{T}, δQ)
        new{H, T, typeof(δQ_)}(δQ_, q, id)
    end
end
OptimalityCut(δQ::AbstractVector, q::Real, id::Int) = HyperPlane(δQ, q, id, OptimalityCut)
FeasibilityCut(δQ::AbstractVector, q::Real, id::Int) = HyperPlane(δQ, q, id, FeasibilityCut)
LinearConstraint(δQ::AbstractVector, q::Real, id::Int) = HyperPlane(δQ, q, id, LinearConstraint)
Infeasible(id::Int) = HyperPlane(sparsevec(Float64[]), 1e10, id, Infeasible)
Unbounded(id::Int) = HyperPlane(sparsevec(Float64[]), 1e10, id, Unbounded)

SparseHyperPlane{T <: Real} = HyperPlane{<:HyperPlaneType, T, SparseVector{T,Int64}}
SparseOptimalityCut{T <: Real} = HyperPlane{OptimalityCut, T, SparseVector{T,Int64}}
SparseFeasibilityCut{T <: Real} = HyperPlane{FeasibilityCut, T, SparseVector{T,Int64}}
SparseLinearConstraint{T <: Real} = HyperPlane{LinearConstraint, T, SparseVector{T,Int64}}

function (hyperplane::HyperPlane{FeasibilityCut})(x::AbstractVector)
    return Inf
end
function (cut::HyperPlane{OptimalityCut})(x::AbstractVector)
    if length(cut.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match", length(cut.δQ), length(x))))
    end
    return cut.q-cut.δQ⋅x
end
function (hyperplane::HyperPlane{Infeasible})(x::AbstractVector)
    return Inf
end
function (hyperplane::HyperPlane{Unbounded})(x::AbstractVector)
    return -Inf
end

infeasible(hyperplane::HyperPlane) = true
infeasible(hyperplane::HyperPlane{Infeasible}) = true
bounded(hyperplane::HyperPlane) = true
bounded(hyperplane::HyperPlane{Unbounded}) = false
function optimal(cut::HyperPlane{OptimalityCut},x::AbstractVector, θ::Real, τ::Real)
    Q = cut(x)
    return θ > -Inf && abs(θ-Q) <= τ*(1+abs(Q))
end
function active(hyperplane::HyperPlane, x::AbstractVector, τ::Real)
    return abs(gap(hyperplane,x)) <= τ
end
function satisfied(hyperplane::HyperPlane, x::AbstractVector, τ::Real)
    return gap(hyperplane,x) >= -τ
end
function satisfied(cut::HyperPlane{OptimalityCut}, x::AbstractVector, θ::Real, τ::Real)
    Q = cut(x)
    return θ > -Inf && θ >= Q - τ
end
function gap(hyperplane::HyperPlane,x::AbstractVector)
    if length(hyperplane.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match", length(hyperplane.δQ), length(x))))
    end
    return hyperplane.δQ⋅x-hyperplane.q
end
function gap(cut::HyperPlane{OptimalityCut}, x::AbstractVector, θ::Real)
    if θ > -Inf
        return θ-cut(x)
    else
        return Inf
    end
end
function lowlevel(hyperplane::HyperPlane{H,T,SparseVector{T,Int}}) where {H <: HyperPlaneType, T <: Real}
    return hyperplane.δQ.nzind, hyperplane.δQ.nzval, hyperplane.q, Inf
end
function lowlevel(cut::HyperPlane{OptimalityCut,T,SparseVector{T,Int}}) where T <: Real
    nzind = copy(cut.δQ.nzind)
    nzval = copy(cut.δQ.nzval)
    push!(nzind,length(cut.δQ)+cut.id)
    push!(nzval,1.0)
    return nzind, nzval, cut.q, Inf
end

# Constructors #
# ======================================================================== #
function OptimalityCut(subproblem::SubProblem)
    λ = getduals(subproblem.solver)
    π = subproblem.π
    cols = zeros(length(subproblem.masterterms))
    vals = zeros(length(subproblem.masterterms))
    for (s,(i,j,coeff)) in enumerate(subproblem.masterterms)
        cols[s] = j
        vals[s] = -π*λ[i]*coeff
    end
    δQ = sparsevec(cols,vals, length(subproblem.x))
    q = π*getobjval(subproblem.solver)+δQ⋅subproblem.x

    return OptimalityCut(δQ, q, subproblem.id)
end
ArtificialCut(val::Real,dim::Int,id::Int) = OptimalityCut(sparsevec(zeros(dim)), val, id)

function FeasibilityCut(subproblem::SubProblem)
    λ = getduals(subproblem.feasibility_solver)
    cols = zeros(length(subproblem.masterterms))
    vals = zeros(length(subproblem.masterterms))
    for (s, (i,j,coeff)) in enumerate(subproblem.masterterms)
        cols[s] = j
        vals[s] = -λ[i]*coeff
    end
    G = sparsevec(cols,vals,length(subproblem.x))
    g = getobjval(subproblem.feasibility_solver)+G⋅subproblem.x

    return FeasibilityCut(G, g, subproblem.id)
end

function LinearConstraint(constraint::JuMP.LinearConstraint, i::Integer)
    sense = JuMP.sense(constraint)
    if sense == :range
        throw(ArgumentError("Cannot handle range constraints"))
    end
    cols = map(v->v.col, constraint.terms.vars)
    vals = constraint.terms.coeffs * (sense == :(>=) ? 1 : -1)
    G = sparsevec(cols, vals, constraint.terms.vars[1].m.numCols)
    g = JuMP.rhs(constraint) * (sense == :(>=) ? 1 : -1)

    return LinearConstraint(G, g, i)
end

function linearconstraints(m::JuMP.Model)
    constraints = Vector{HyperPlane{LinearConstraint}}(length(m.linconstr))
    for (i, c) in enumerate(m.linconstr)
        constraints[i] = LinearConstraint(c, i)
    end
    return constraints
end

Infeasible(subprob::SubProblem) = Infeasible(subprob.id)
Unbounded(subprob::SubProblem) = Unbounded(subprob.id)

function zero(h::HyperPlane{H,T,A}) where {H <: HyperPlaneType, T <: Real, A <: AbstractVector}
    return HyperPlane(zero(h.δQ), zero(T), h.id, H)
end
function +(h1::HyperPlane{H,T,A},h2::HyperPlane{H,T,A}) where {H <: HyperPlaneType, T <: Real, A <: AbstractVector}
    return HyperPlane(h1.δQ + h2.δQ, h1.q + h2.q, h1.id, H)
end

mutable struct CutBundle{T <: Real}
    cuts::Vector{SparseOptimalityCut{T}}
    q::T

    function (::Type{CutBundle})(::Type{T}) where T <: Real
        new{T}(Vector{SparseOptimalityCut{T}}(), zero(T))
    end
end
length(bundle::CutBundle) = length(bundle.cuts)
function aggregate!(bundle::CutBundle)
    aggregated_cut = zero(bundle.cuts[1])
    for j = 1:length(bundle)
        aggregated_cut += pop!(bundle.cuts)
    end
    return aggregated_cut
end

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


text:
using CommonRLInterface
using Test

mutable struct LQREnv <: AbstractEnv
    s::Float64
end

@testset "from README" begin

    function CommonRLInterface.reset!(m::LQREnv)
        m.s = 0.0
    end

    CommonRLInterface.actions(m::LQREnv) = [-1.0, 0.0, 1.0]
    CommonRLInterface.observe(m::LQREnv) = m.s
    CommonRLInterface.terminated(m::LQREnv) = false

    function CommonRLInterface.act!(m::LQREnv, a)
        r = -m.s^2 - a^2
        m.s = m.s + a + randn()
        return r
    end

    env = LQREnv(0.0)
    reset!(env)
    acts = actions(env)
    rsum = 0.0
    step = 1
    while !terminated(env) && step <= 10
        rsum += act!(env, rand(acts))
        step += 1
    end
    @show rsum
end

# a reference MDP Env
mutable struct MyEnv <: AbstractEnv
    state::Int
end
MyEnv() = MyEnv(1)
function CommonRLInterface.reset!(env::MyEnv)
    env.state = 1
end
CommonRLInterface.actions(env::MyEnv) = [-1, 0, 1]
CommonRLInterface.observe(env::MyEnv) = env.state
CommonRLInterface.terminated(env::MyEnv) = false
function CommonRLInterface.act!(env::MyEnv, a)
    env.state = clamp(env.state + a, 1, 10)
    return -o^2
end
env = MyEnv(1)

# a reference Game Env
mutable struct MyGame <: AbstractEnv
    state::Int
end
MyGame() = MyGame(1)
function CommonRLInterface.reset!(env::MyGame)
    env.state = 1
end
CommonRLInterface.actions(env::MyGame) = [-1, 1]
CommonRLInterface.observe(env::MyGame) = env.state
CommonRLInterface.terminated(env::MyGame) = false
function CommonRLInterface.act!(env::MyGame, a)
    env.state = clamp(env.state + a, 1, 10)
    return -o^2
end
@provide CommonRLInterface.player(env::MyGame) = 1 + iseven(env.state)
game = MyGame()

function f end

# h needs to be out here for the inference to work for some reason
function h(x)
    if provided(f, x)
        return f(x)
    else
        return sin(x)
    end
end

@testset "providing" begin
    global f, h

    @test_throws MethodError f(2)
    @test provided(f, 2) == false

    @provide f(x::Int) = x^2

    @test provided(f, 2) == true
    @test provided(f, Tuple{Int}) == true
    @test f(2) == 4

    @provide function f(s::String)
        s*"^2"
    end
    @test provided(f, "2") == true
    @test provided(f, Tuple{String}) == true
    @test f("2") == "2^2"

    @test provided(f, 2.0) == false
    @test provided(f, Tuple{Float64}) == false

    @provide f(x::AbstractArray{N}) where {N<:Number} = x.^2

    @test provided(f, Tuple{Vector{Float64}})
    @test f([1,2]) == [1, 4]

    @test @inferred(h(2)) == f(2)::Int

    g() = nothing

    @provide g(a, b::AbstractVector{<:Number}) = a.*b

    @test provided(g, 1, [1.0]) == true
    @test provided(g, typeof((1, [1.0]))) == true
    @test g(1, [1.0]) == [1.0]
    @test provided(g, 1, ["one"]) == false

    @test provided(reset!, MyEnv())
    @test !provided(clone, MyEnv())
    @test !provided(player, MyEnv())
    @test !provided(actions, MyEnv(), 1)

    @test provided(player, MyGame())
end

@testset "environment" begin
    Base.:(==)(a::MyEnv, b::MyEnv) = a.state == b.state
    Base.hash(x::MyEnv, h=zero(UInt)) = hash(x.state, h)
    @provide CommonRLInterface.clone(env::MyEnv) = MyEnv(env.state)
    @test provided(clone, env)
    @test clone(env) == MyEnv(1)

    @provide CommonRLInterface.render(env::MyEnv) = "MyEnv with state $(env.state)"
    @test provided(render, MyEnv(1))
    @test render(MyEnv(1)) == "MyEnv with state 1"

    @provide CommonRLInterface.state(env::MyEnv) = env.state
    @provide function CommonRLInterface.setstate!(env::MyEnv, s) env.state = s end
    @test provided(state, MyEnv(1))
    @test state(MyEnv(1)) == 1
    env1 = MyEnv(0)
    setstate!(env1, 1)
    @test state(env1) == 1
end

@testset "spaces" begin
    @provide CommonRLInterface.valid_actions(env::MyEnv) = [0, 1]
    @test provided(valid_actions, env)
    @test issubset(valid_actions(env), actions(env))

    @provide CommonRLInterface.valid_action_mask(env::MyEnv) = [false, true, true]
    @test provided(valid_action_mask, env)
    va = valid_actions(env)
    vam = actions(env)[valid_action_mask(env)]
    @test issubset(va, vam)
    @test issubset(vam, va)

    @provide CommonRLInterface.observations(env::MyEnv) = 1:10
    @test provided(observations, env)
    @test observations(env) == 1:10
end

include("examples/gridworld.jl")
include("examples/tictactoe.jl")
include("examples/rock_paper_scissors.jl")

include("wrappers.jl")


text:
# This file was generated, do not modify it. # hide
corestrict(string.(1:10), f, 2)


text:
module ImageComparisons
    using Winston, Compat, Colors
    using Printf, Random
    include("examples.jl")
    include("issues.jl")
    include("plot.jl")
end

function read_png_data(fn::String)
    surface = Cairo.read_from_png(fn)
    w = Cairo.width(surface)
    h = Cairo.height(surface)
    a = Array{UInt32}(undef, convert(Int, w), convert(Int, h))

    p = ccall((:cairo_image_surface_get_data,Cairo._jl_libcairo),
              Ptr{UInt32}, (Ptr{Cvoid},), surface.ptr)

    for i = 1:length(a)
        a[i] = unsafe_load(p, i)
    end
    return a
end

function img_dist(img1::Matrix{UInt32}, img2::Matrix{UInt32})
    @assert size(img1) == size(img2)
    s = 0.
    for i = 1:length(img1)
        a = convert(RGB, reinterpret(RGB24, img1[i]))
        b = convert(RGB, reinterpret(RGB24, img2[i]))
        s += colordiff(a, b)
    end
    s
end

@testset "Image comparison" begin
    root = dirname(task_local_storage()[:SOURCE_PATH])
    dir1, dir2 = joinpath(root, "_baseline"), joinpath(root, "_current")
    isdir(dir1) || mkdir(dir1)
    isdir(dir2) || mkdir(dir2)

    for name in sort!(names(ImageComparisons))
        name == :ImageComparisons && continue

        @testset "$name" begin
            println("Testing $name")
            figure()
            func = eval(:(ImageComparisons.$name))
            global p = func()

            fn1 = joinpath(dir1, "$name.png")
            fn2 = joinpath(dir2, "$name.png")
            isfile(fn1) || savefig(p, fn1)

            savefig(p, fn2)

            img1 = read_png_data(fn1)
            img2 = read_png_data(fn2)
            @test size(img1) == size(img2)

            @test img_dist(img1, img2) < 0.1
        end
    end
end


text:
include("datatypes.jl")
include("utils.jl")
include("consts.jl")
include("integral.jl")
include("rules.jl")