text:
# Abstraction of size information
abstract type AbstractSize end

struct CountSize <: AbstractSize
    val::Int  # objects
end

struct MemorySize <: AbstractSize
    val::Int  # bytes
end

# show methods
show(io::IO, sz::CountSize) =
    sz.val == 1 ? print(io, "1 object") : print(io, "$(sz.val) objects")

show(io::IO, sz::MemorySize) = print(io, "$(sz.val/1024) KiB")

# Object size
object_size(dd::Dict, ::Type{CountSize}) = length(dd)

object_size(dd::Dict, ::Type{MemorySize}) =
    isempty(dd) ? 0 : mapreduce(x->summarysize(x), +, values(dd))

object_size(object, ::Type{CountSize}) = 1

object_size(object, ::Type{MemorySize}) = summarysize(object)

# Cache struct
abstract type AbstractCache end

mutable struct Cache{T<:Function, O, S<:AbstractSize} <: AbstractCache
    func::T                                 # function being cached
    name::String                            # name of the object
    filename::String                        # file name
    func_def::Union{Nothing, String}        # function code (if any)
    cache::Dict{UInt, O}                    # cache dictionary
    offsets::Dict{UInt, Tuple{UInt, UInt}}  # existing hash - to - file positions
    history::Deque{UInt}                    # order in which the calls were executed (by hash)
    max_size::S                             # maximum size of the cache
end

# Overload constructor
function Cache(f::T;
               name::String = string(f),
               filename::String = generate_cache_filename(),
               func_def = nothing,
               output_type::Type=Any,
               max_size::S=CountSize(MAX_CACHE_SIZE)
              ) where {T<:Function, S<:AbstractSize}
    # The main constructor call (separating comment ;)
    return Cache(f,
                 name,
                 abspath(filename),
                 func_def,
                 Dict{UInt, output_type}(),
                 Dict{UInt, Tuple{UInt, UInt}}(),
                 Deque{UInt}(),
                 max_size)
end

# Show method
show(io::IO, cache::Cache) = begin
     memory_size = length(cache.cache)
     disk_size = length(cache.offsets)
     total_size = length(symdiff(keys(cache.cache), keys(cache.offsets))) +
         length(intersect(keys(cache.cache), keys(cache.offsets)))
     _en = ifelse(total_size == 1, "entry", "entries")
     print(io, "$(cache.name) (cache with $total_size $_en, ",
               "$memory_size in memory $disk_size on disk)")
end

# Other useful functions
length(cache::Cache) = length(cache.cache)

max_cache_size(cache::Cache) = cache.max_size.val

object_size(cache::Cache{T,O,S}) where {T<:Function, O, S<:AbstractSize} =
    object_size(cache.cache, S)

# Call method (caches only to memory, caching to disk has to be explicit)
function (cache::Cache{T,O,S})(args...; kwargs...) where {T<:Function, O, S<:AbstractSize}
    # ~~ Caclculate hash ~~
    __hash__ = arghash(args...; kwargs...)
    # ~~~
    _, decompressor = get_transcoders(cache.filename)
    if haskey(cache.cache, __hash__)
        # Move hash from oldest to most recent
        # so that the next entry does not remove it;
        # return the cached value
        if max_cache_size(cache) <= object_size(cache) && front(cache.history) == __hash__
            push!(cache.history, __hash__)
            popfirst!(cache.history)
        end
        return cache.cache[__hash__]
    elseif haskey(cache.offsets, __hash__)
        # Entries found only on disk do not update the history,
        # this is just a load operation; only an explicit synchonization
        # will load into memory
        @debug "Memory cache miss, loading hash=0x$(string(__hash__, base=16))..."
        open(cache.filename, "r") do fid
            startpos, endpos = cache.offsets[__hash__]
            datum = load_disk_cache_entry(fid, startpos, endpos, decompressor)
            if !(typeof(datum) <: O)
                throw(AssertionError("Output type is not a subtype $O"))
            end
            return datum::O
        end
    else
        # A hash miss: a new value must be cached. This requires a check
        # of the size of the returned object and dynamic removal of existing
        # entries if the maximum size is reached.
        @debug "Full cache miss, adding hash=0x$(string(__hash__, base=16))..."
        # Check that the output type matches
        if !(return_type(cache.func, typeof.(args)) <: O)
            throw(AssertionError("Output type is not a subtype $O"))
        end
        # Run function
        out::O = cache.func(args...; kwargs...)
        # Check if cache is full when limit set in KiB
        if S<:MemorySize && object_size(out,S) > max_cache_size(cache)  # check for memory
            @warn "Cannot cache result of size $(S(object_size(out,S))) " *
            "(maximum cache size is $(cache.max_size))."
            return out
        end
        # Delete objects and cache
        while object_size(cache) + object_size(out, S) > max_cache_size(cache)
            delete!(cache.cache, popfirst!(cache.history))
        end
        cache.cache[__hash__] = out
        push!(cache.history, __hash__)
        return out
    end
end

# @cache macro
macro cache(expression,
            filename::String=generate_cache_filename(),
            max_size::Number=MAX_CACHE_SIZE)
    # Check size, no need to check the filename
    @assert max_size > 0 "The maximum size has to be > 0 (objects or KiB)."
    if max_size isa Int
        _max_size = CountSize(max_size)
    elseif max_size isa Real
        _max_size = MemorySize(max_size*1024)
    else
        @error "The maximum size has to be an Int or a Float64."
    end

    # Parse macro input
    if expression isa Symbol
        #######################
        # julia> @cache foo   #-> foo remains a function, new Cache object returned
        #######################
        _name = String(expression)
        ex = quote
                Base.@__doc__ Cache($expression;
                                    name=$_name,
                                    filename=$filename,
                                    func_def=nothing,
                                    max_size=$_max_size)
        end
    elseif expression isa Expr && expression.head == :(::) && length(expression.args) == 2
        ############################
        # julia> @cache foo::Int   #-> foo remains a function, new Cache object returned
        ############################
        _symb = expression.args[1]
        _typesymbol = expression.args[2]
        _type = eval(_typesymbol)
        _name = String(_symb)
        @assert _type isa Type "The right-hand argument of `::` is not a type."
        ex = quote
                Base.@__doc__ Cache($_symb;
                                    name=$_name,
                                    filename=$filename,
                                    func_def=nothing,
                                    output_type=$_type,
                                    max_size=$_max_size)
        end
    elseif expression isa Expr && expression.head == :(=)
        #############################
        # julia> @cache foo=x->x+1  #-> foo becomes a Caching.Cache object
        #############################
        f_output_type = Any
        if @capture(expression, f_name_ = arg_::input_type_->body_)
            if input_type == nothing
                input_type = Any
            end
        elseif @capture(expression, f_name_ = arg_->body_)
            input_type = Any
        else
            @error "Only one input argument lambdas are supported."
        end
        new_expression = :(($arg::$input_type)->$body)
        new_expression_str = string(new_expression)
        _func = eval(new_expression)
        _t = eval(input_type)
        ex = quote
            Base.@__doc__ $f_name = Caching.Cache(eval($new_expression);
                                                  name=$(string(f_name)),
                                                  filename=$filename,
                                                  func_def=$new_expression_str,
                                                  output_type=Core.Compiler.return_type($_func,($_t,)),
                                                  max_size=$_max_size)
        end

    else
        #######################################
        # julia> @cache function foo(x::Int)  #
        #                   x+1               #-> foo becomes a Caching.Cache object
        #               end                   #
        #######################################
        expression_str = string(expression)
        f_parsed = splitdef(expression)
        f_name = f_parsed[:name]
        f_output_type = get(f_parsed, :rtype, :Any)
        random_name = Symbol(:f_, randstring(20))
        new_definition = copy(f_parsed)
        new_f_name = random_name
        new_definition[:name] = new_f_name
        ex = quote
                $(MacroTools.combinedef(new_definition))  # reconstruct function def and run it
                Base.@__doc__ $f_name = Caching.Cache($random_name;
                                                      name=$(string(f_name)),
                                                      filename=$filename,
                                                      func_def=$expression_str,
                                                      output_type=$f_output_type,
                                                      max_size=$_max_size)
        end
    end
    return esc(ex)
end  # macro


text:
using MyExampleILIC
using Test

@testset "MyExampleILIC.jl" begin
    @test my_f(2,1) == 5
    @test my_f(2, 2) == 6
    @test greetings() == "hello"
end

ext: jl


text:
module Agents

include("Agents/BaseAgent.jl")
include("Agents/RandomChan.jl")
include("Agents/CooperativeQ.jl")
include("Agents/IndividualQ.jl")
include("Agents/OptHighestSNR.jl")
include("Agents/ContextAware.jl")

end


text:
module EconFunctions

using FastGaussQuadrature

include("quadrature.jl")
include("utility.jl")

end # module


text:
function cartesian_from_polar(r, θ)
    s, c = sincos(θ)
    x = r * c
    y = r * s
    return SVector(x, y)
end

# TODO: tighter enclosures for wide intervals
polar_from_cartesian(xy::AbstractVector) = (hypot(xy[1], xy[2]), atan(xy[2], xy[1]))

function polar_from_cartesian(xy::AbstractVector{arb_series})
    r = sqrt(xy[1]^2 + xy[2]^2)

    if contains_zero(xy[2][0]) && !iszero(xy[2][0]) && !isnonnegative(xy[1][0])
        θ = atan(xy[2], -xy[1]) + base_ring(xy[1].poly)(π)
    else
        θ = atan(xy[2], xy[1])
    end

    return r, θ
end

function polar_from_cartesian(xy::AbstractVector{arb})
    r = hypot(xy[1], xy[2])

    if contains_zero(xy[2]) && !iszero(xy[2]) && !isnonnegative(xy[1])
        θ = atan(xy[2], -xy[1]) + parent(xy[1])(π)
    else
        θ = atan(xy[2], xy[1])
    end

    return r, θ
end

function cartesian(θ, ϕ)
    sθ, cθ = sincos(θ)
    sϕ, cϕ = sincos(ϕ)
    SVector(sθ * cϕ, sθ * sϕ, cθ)
end

cartesian((θ, ϕ)) = cartesian(θ, ϕ)

function spherical(xyz)
    (θ = acos(xyz[3]), ϕ = atan(xyz[2], xyz[1]))
end

function Base.Float64((θ, ϕ)::NamedTuple{(:θ, :ϕ),Tuple{arb,arb}})
    (θ = Float64(θ), ϕ = Float64(ϕ))
end

function LinearAlgebra.dot(xyz1::SVector{3,arb}, xyz2::SVector{3,arb})
    sum(xyz1 .* xyz2)
end

LinearAlgebra.norm(x::SVector{N,arb}) where {N} = sqrt(sum(x .^ 2))

LinearAlgebra.normalize(x::SVector{N,arb}) where {N} = x ./ sqrt(sum(x .^ 2))

LinearAlgebra.normalize(x::SVector{N,arb_series}) where {N} = x ./ sqrt(sum(x .^ 2))


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

module EnzymeMFABC
include("../example.jl")
end

using FileIO, JLD2
using Statistics

Random.seed!(111)
MFABC_output = EnzymeMFABC.do_mfabc.([40_000, 80_000, 160_000, 320_000, 640_000];
    δ = 1e3,
    N_batch = 1,
    GD_batch = 1,
)

c = map(MFABC_output) do (S, μ)
	sum(S.c_lo .+ S.Σc_hi)
end
v = map(MFABC_output) do (S, μ)
    Ḡ = mean(EnzymeMFABC.G, S)
	w_mf = S.L_lo .+ (S.ΣL_mf./S.μ)
	sum((w_mf .* (EnzymeMFABC.G.(S.θ) .- Ḡ)).^2) / sum(w_mf)^2
end

S, μ, μ_t, ν_t = MFABC_output[end]

interesting = findall(p->(!iszero(p.ΣL_mf2)), S)
w_mf = S.L_lo .+ (S.ΣL_mf ./ S.μ)
w_mf_sparse = hcat(interesting, w_mf[interesting])

save("MFABC_output.jld2", Dict(
    "c" => c,
    "v" => v,
    "w_mf" => w_mf_sparse,
    "μ" => μ,
    "μ_t" => μ_t[100:100:end, :],
    "ν_t" => ν_t[100:100:end, :],
))

max_stars_repo_name: tpprescott/mf-lf

max_issues_repo_name: tpprescott/mf-lf

max_forks_repo_name: tpprescott/mf-lf


text:
module Stars

using ..Ahorn, Maple

placements = Maple.Stars

function Ahorn.canFgBg(effect::Maple.Stars)
    return true, true
end

end


text:
#JuliaRL/src/conf.jl

# Configuration values for the DQN algorithms are defined here
module Conf

# BasicDQN, DQN, PDQN

# duration in steps of the training
duration = 150_000
# size of mini-batches
batch_size = 64
# number of transitions that should be experienced before updating the approximator
min_replay_history = 100
# decay_steps for EpsilonGreedyExplorer
decay_steps = 3000
# capacity (in steps) of the experience buffer
capacity = duration
# frequency at which the agent is saved during training
save_freq = div(duration,2)

# DQN, PDQN

# the frequency of updating the approximator
update_freq = 4
# the frequency of updating the target
target_update_freq = 100

end

max_stars_repo_name: FeLusiani/JuliaRL

max_issues_repo_name: FeLusiani/JuliaRL

max_forks_repo_name: FeLusiani/JuliaRL


text:
module _XCondEn
export XCondEn
using StatsBase: fit, Histogram
using Statistics: mean, std
    """
        XCond, SEw, SEz = XCondEn(Sig)

    Returns the corrected cross-conditional entropy estimates (`XCond`) and the
    corresponding Shannon entropies (m: SEw, m+1: SEz) for m = [1,2]
    estimated for the data sequences contained in `Sig` using the default
    parameters:  embedding dimension = 2, time delay = 1, number of symbols = 6,
    logarithm = natural
    ** Note: XCondEn is direction-dependent. Therefore, the order of the
    data sequences in `Sig` matters. If the first row/column of `Sig` is the
    sequence 'y', and the second row/column is the sequence 'u', the `XCond` is
    the amount of information carried by y(i) when the pattern u(i) is found.**

        XCond, SEw, SEz = XCondEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1, c::Int=6, Logx::Real=exp(1), Norm::Bool=false)

    Returns the corrected cross-conditional entropy estimates (`XCond`) for
    the data sequences contained in `Sig` using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, an integer > 1        [default: 2]  \n
    `tau`   - Time Delay, a positive integer             [default: 1]  \n
    `c`     - Number of symbols, an integer > 1          [default: 6]  \n
    `Logx`  - Logarithm base, a positive scalar          [default: natural] \n
    `Norm`  - Normalisation of `XCond` values:
                [false]  no normalisation                  [default]\n
                [true]   normalises w.r.t cross-Shannon entropy.  \n

    # See also `XFuzzEn`, `XSampEn`, `XApEn`, `XPermEn`, `CondEn`, `XMSEn`

    # References:
        [1] Alberto Porta, et al.,
            "Conditional entropy approach for the evaluation of the
            coupling strength."
            Biological cybernetics
            81.2 (1999): 119-129.

    """
    function XCondEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1,
        c::Int=6, Logx::Real=exp(1), Norm::Bool=false)

    (size(Sig,2) > size(Sig,1)) ? Sig = transpose(Sig) : nothing

    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 1) ? nothing : error("m:     must be an integer > 1")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing : error("Logx:     must be a positive number > 0")
    (c>1) ? nothing :  error("c:     must be an integer > 1")

    S1 = (Sig[:,1] .- mean(Sig[:,1]))/std(Sig[:,1],corrected=false)
    S2 = (Sig[:,2] .- mean(Sig[:,2]))/std(Sig[:,2],corrected=false)
    Edges = range(minimum(S1),maximum(S1),length=c+1)
    Sx1 = map(x -> sum(Edges[1:c].<=x), S1)
    Edges = range(minimum(S2),maximum(S2),length=c+1)
    Sx2 = map(x -> sum(Edges[1:c].<=x), S2)

    SEw = zeros(m-1)
    SEz = zeros(m-1)
    Prcm = zeros(m-1)
    Xi = zeros(Int,N,m)
    for k = 1:m-1
        Nx = N-(k-1)*tau
        Xi[1:Nx,m-(k-1)] = Sx1[(k-1)*tau+1:N]
        #Wi = transpose(c.^(k-1:-1:0))*transpose(Xi[1:Nx,m-k+1:m])
        Wi = Xi[1:Nx,m-k+1:m]*(c.^(k-1:-1:0))
        #Zi = (c^k)*transpose(Sx2[(k-1)*tau+1:N]) .+ Wi
        Zi = (c^k)*Sx2[(k-1)*tau+1:N]  + Wi
        Pw = fit(Histogram, Wi[:], minimum(Wi)-.5:maximum(Wi)+.5).weights
        Pz = fit(Histogram, Zi[:], minimum(Zi)-.5:maximum(Zi)+.5).weights
        Prcm[k] = sum(Pz.==1)/Nx

        (sum(Pw)!= Nx || sum(Pz)!= Nx) ? @warn("Potential error estimating probabilities.") : nothing

        Pw = Pw[Pw.!=0];     Pw = Pw/N
        Pz = Pz[Pz.!=0];     Pz = Pz/N
        SEw[k] = -transpose(Pw)*log.(Logx, Pw)
        SEz[k] = -transpose(Pz)*log.(Logx, Pz)
    end

    Temp = fit(Histogram,Sx2,nbins=c).weights/N
    Temp = Temp[Temp.!=0]
    Sy = -transpose(Temp)*log.(Logx, Temp)
    XC = SEz - SEw + Prcm*Sy
    XC = vcat(Sy, XC)
    Norm ? XC = XC/Sy : nothing

    return XC, SEw, SEz
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using UniqueVectors
using Test
using InteractiveUtils: @which

@test length(UniqueVector([1,5,6,3])) == 4
@test_throws ArgumentError UniqueVector([1,3,5,6,3])

uv = UniqueVector{String}()

@test isempty(uv)
@test allunique(uv)
@test_throws ArgumentError pop!(uv)
@test findfirst(isequal("cat"), uv) == nothing
@test findfirst(isequal("cat"), uv) == nothing
@test findfirst!(isequal("cat"), uv) == 1
@test !isempty(uv)
@test "cat" in uv
@test "dog" ∉ uv
@test count(isequal("cat"), uv) == 1
@test count(isequal("dog"), uv) == 0
@test findall(isequal("cat"), uv) == [1]
@test findall(isequal("dog"), uv) == Int[]
@test findfirst!(isequal("dog"), uv) == 2
@test findfirst!(isequal("cat"), uv) == 1
@test findfirst!(isequal("mouse"), uv) == 3
@test findfirst!(isequal("dog"), uv) == 2
@test findfirst(isequal("cat"), uv) == 1
@test findlast(isequal("cat"), uv) == 1
@test findfirst(isequal("dog"), uv) == 2
@test findfirst(isequal("mouse"), uv) == 3
@test uv[1] == "cat"
@test uv[2] == "dog"
@test uv[3] == "mouse"
@test uv[:] == ["cat", "dog", "mouse"]
@test size(uv) == (3,)
@test length(uv) == 3
@test lastindex(uv) == 3
@test allunique(uv)
@test unique!(uv) === uv
@test unique(uv) == uv
@test unique(uv) == uv.items
@test unique(uv) !== uv.items
empty!(unique(uv))
@test length(uv.items) == 3

uv2 = UniqueVector([1, 2, 3])
@test eltype(uv2) == Int
@test findfirst(isequal(0x02), uv2) == 2
@test findfirst!(isequal(0x02), uv2) == 2
@test uv2 == UniqueVector(collect(i for i in 1:3))
for elt in [3,2,1]
    @test pop!(uv2) == elt
end
@test isempty(uv2)

uv2 = copy(uv)
@test uv2 == uv

@test empty!(uv) === uv
@test isempty(uv)
@test findfirst(isequal("cat"), uv) == nothing
@test findfirst(isequal("cat"), uv) == nothing
@test findfirst!(isequal("horse"), uv) == 1
@test_throws ArgumentError push!(uv, "horse")
@test length(uv) == 1
@test push!(uv, "human") === uv
@test findfirst(isequal("human"), uv) == 2
@test pop!(uv) == "human"
@test length(uv) == 1
@test uv[:] == ["horse"]
@test findfirst(isequal("human"), uv) == nothing
@test findlast(isequal("human"), uv) == nothing

@test uv2[:] == ["cat", "dog", "mouse"]
@test findfirst(isequal("cat"), uv2) == 1

# Test swap!
let uv = UniqueVector(["cat", "dog", "mouse", "human"]), original = copy(uv)

    uv = swap!(uv, 2, 2)

    @test uv == original
    @test_throws BoundsError swap!(uv, 5, 5)

    swap!(uv, 2, 3)

    @test findfirst(isequal("mouse"), uv) == 2
    @test findfirst(isequal("mouse"), original) == 3

    @test findfirst(isequal("dog"), uv) == 3
    @test findfirst(isequal("dog"), original) == 2
end

# Test permute! and invpermute!
let uv = UniqueVector(["cat", "dog", "mouse", "human"]), original = copy(uv)
    perm = sortperm(rand(4))
    permute!(uv, perm)
    @test uv == original[perm]
    @test uv[findfirst(isequal("dog"), uv)] == "dog"

    invpermute!(uv, perm)
    @test uv == original

    invpermute!(uv, perm)
    @test uv == original[invperm(perm)]
    @test uv[findfirst(isequal("dog"), uv)] == "dog"
end

@test UniqueVector([1,2,3,4]) == UniqueVector(1:4)

# Test it works with `Any` datatype
let uv3 = UniqueVector([1,"cat",2,"dog"])
    @test eltype(uv3) == Any
    @test findfirst(isequal(1), uv3) == 1
    @test findlast(isequal(1), uv3) == 1
    @test findall(isequal(1), uv3) == [1]
    @test count(isequal(1), uv3) == 1
    @test findfirst!(isequal("dog"), uv3) == 4
    @test findfirst!(isequal("horse"), uv3) == 5
end

# Test setindex!
uv4 = UniqueVector(["cat", "dog", "mouse"])
@test_throws BoundsError uv4[4] = "horse"
uv4[2] = "horse"
uv4[3] = "dog"
@test uv4[:] == ["cat", "horse", "dog"]
uv4[1] = "cat"
@test_throws ArgumentError uv4[2] = "dog"
push!(uv4, "mouse")
@test uv4[:] == ["cat", "horse", "dog", "mouse"]
@test uv4[1:2] == ["cat", "horse"]

# Test view
@test view(uv4, 2:3) isa SubArray
@test findfirst(isequal("dog"), view(uv4, 2:3)) == 2
@test @which(findfirst(isequal("dog"), view(uv4, 2:3))).module == UniqueVectors

uv5 = UniqueVector{Float64}()
push!(uv5, 3)
@test uv5[:] == [3.0]
uv5[1] = 4
@test uv5[:] == [4.0]
@test 4 ∈ uv5
@test findfirst(isequal(4), uv5) == 1
@test findlast(isequal(4), uv5) == 1
@test findall(isequal(4), uv5) == [1]

# Test indexin and findall(in)
@test indexin([1,2,34,0,5,56], UniqueVector([34,56,35,1,5,0])) == [4,nothing,1,6,5,2]
@test indexin([1,2,34,0,5,56], UniqueVector([34,56,35,1,5,0])) == [4,nothing,1,6,5,2]
@test findall(in(UniqueVector([34,56,35,1,5,0])), [1,2,34,0,5,56]) == [1,3,4,5,6]
@test findall(in(UniqueVector([34,56,35,1,5,0])), (1,2,34,0,5,56)) == [1,3,4,5,6]
@test findall(in(UniqueVector([5,7,9])), [5 6; 7 8]) == findall(in([5,7,9]), [5 6; 7 8])

# Test findnext and findprev
@test findnext(isequal(7), UniqueVector([3,5,7,9]), 1) == 3
@test findnext(isequal(7), UniqueVector([3,5,7,9]), 2) == 3
@test findnext(isequal(7), UniqueVector([3,5,7,9]), 3) == 3
@test findnext(isequal(7), UniqueVector([3,5,7,9]), 4) == nothing
@test findnext(isequal(7), UniqueVector([3,5,7,9]), 4) == nothing
@test findprev(isequal(7), UniqueVector([3,5,7,9]), 1) == nothing
@test findprev(isequal(7), UniqueVector([3,5,7,9]), 1) == nothing
@test findprev(isequal(7), UniqueVector([3,5,7,9]), 2) == nothing
@test findprev(isequal(7), UniqueVector([3,5,7,9]), 2) == nothing
@test findprev(isequal(7), UniqueVector([3,5,7,9]), 3) == 3
@test findprev(isequal(7), UniqueVector([3,5,7,9]), 4) == 3

# Test constructor behavior (#7)
let items = [1,2,3]
    @test UniqueVector(items).items !== items
    @test UniqueVector{Int}(items).items !== items
end


text:
using SVM
using Test

@testset "SVM.jl" begin

end


text:
"""
    CWindow(N::Vector{Float64}, NCut::Int64)

This function evaluates the smoothing window function as defined
[in Eq. (C.1) of this paper](https://arxiv.org/abs/1603.04826).
"""
function CWindow(N::Vector{Float64}, NCut::Int64)
    NRight = last(N) - NCut
    NR = filter(x->x>=NRight, N)
    ThetaRight = (last(N).-NR) ./ (last(N) - NRight - 1)
    W = ones(length(N))
    W[findall(x->x>=NRight, N)] = ThetaRight .- 1 ./ (2*π) .* sin.(2 .* π .*
	ThetaRight)
    return W
end

function GetCM!(FFTLog::FFTLog)
    FFTLog.CM = FFTW.rfft(FFTLog.FXArray .* FFTLog.XArray .^ (-FFTLog.ν))
    FFTLog.M = Array(0:length(FFTLog.XArray)/2)
    FFTLog.CM = FFTLog.CM .* CWindow(FFTLog.M, floor(Int,
	FFTLog.CWindowWidth*FFTLog.N/2))
end

function EvaluateηM!(FFTLog::FFTLog)
    FFTLog.ηM = 2 .* π ./ (FFTLog.N .* FFTLog.DLnX) .* FFTLog.M
    return FFTLog.ηM
end

function GMVals(Mu::Float64, Q::Array{Complex{Float64},1})
    if(Mu+1+real(Q[1]) == 0)
        println("Γ(0) encountered. Please change another ν value! Try ν=1.1 .")
    end
    ImagQ = imag(Q)
    GM = zeros(ComplexF64, length(Q))
    cut = 200
    AsymQ = filter(x->abs(imag(x)) + abs(Mu) >cut, Q)
    AsymPlus = (Mu .+ 1 .+ AsymQ) ./2
    AsymMinus = (Mu .+ 1 .- AsymQ) ./2

    QGood = filter(x->abs.(imag(x)) .+ abs(Mu) .<=cut && x != Mu + 1 ,Q)
    AlphaPlus  = (Mu .+1 .+ QGood) ./2
    AlphaMinus = (Mu .+1 .- QGood) ./2
    GM[findall(x->abs.(imag(x)) .+ abs(Mu) .<=cut && x != Mu + 1 , Q)] .=
	SpecialFunctions.gamma.(AlphaPlus) ./ SpecialFunctions.gamma.(AlphaMinus)
    GM[findall(x->abs.(imag(x))+abs(Mu) .> cut && x != Mu + 1 , Q)] = exp.(
	(AsymPlus .- 0.5) .* log.(AsymPlus) .- (AsymMinus .- 0.5) .*
	log.(AsymMinus) .-
	AsymQ .+ 1/12 .* (1 ./AsymPlus .- 1 ./ AsymMinus) .+ 1/360 .* (1 ./
	AsymMinus .^3 .- 1 ./AsymPlus .^3) +1/1260 .* (1 ./AsymPlus .^ 5  .- 1 ./
	AsymMinus .^ 5) )
    return GM
end

function GL(Ell::Float64, ZArray::Vector{ComplexF64})
    GL = 2 .^ZArray .* GMVals(Ell+0.5, ZArray .- 1.5)
    return GL
end

function GL1(Ell::Float64, ZArray::Vector{ComplexF64})
    GL1 = -2 .^(ZArray .- 1) * (z_array -1) .* GMVals(Ell+0.5, ZArray .- 2.5)
    return GL1
end

function GL2(Ell::Float64, ZArray::Vector{ComplexF64})
    GL2 = 2 .^ (ZArray .- 2) .* (ZArray .- 2) .* (ZArray .- 1) .*
	GMVals(Ell+0.5, ZArray .- 3.5)
    return GL
end

function LogExtrap(X::Vector{Float64}, NExtrapLow::Int64, NExtrapHigh::Int64)
    DLnXLow = log(X[2]/X[1])
    DLnXHigh= log(reverse(X)[1]/reverse(X)[2])
    if NExtrapLow != 0
        LowX = X[1] .* exp.(DLnXLow .* Array(-NExtrapLow:-1))
        X = vcat(LowX, X)
    end
    if NExtrapHigh != 0
        HighX = last(X) .* exp.(DLnXLow .* Array(1:NExtrapHigh))
        X = vcat(X,HighX)
    end
    return X
end

function ZeroPad!(FFTLog::FFTLog)
    ZeroArray = zeros(FFTLog.NPad)
    FFTLog.XArray = LogExtrap(FFTLog.XArray, FFTLog.NPad, FFTLog.NPad)
    FFTLog.FXArray = vcat(ZeroArray, FFTLog.FXArray, ZeroArray)
end

function CheckNumberElements!(FFTLog::FFTLog)
    if iseven(FFTLog.N+1)
        deleteat!(FFTLog.XArray, length(FFTLog.XArray))
        deleteat!(FFTLog.FXArray, length(FFTLog.FXArray))
        FFTLog.N -= 1
        if (FFTLog.NExtrapHigh != 0)
            FFTLog.NExtrapHigh -= 1
        end
    end
end

function EvaluateFFTLog(FFTLog::FFTLog, Ell::Vector{Float64})
    FFTLog.XArray = LogExtrap(FFTLog.XArray, FFTLog.NExtrapLow,
	FFTLog.NExtrapHigh)
    FFTLog.FXArray = LogExtrap(FFTLog.FXArray, FFTLog.NExtrapLow,
	FFTLog.NExtrapHigh)
    ZeroPad!(FFTLog)
    GetCM!(FFTLog)
    EvaluateηM!(FFTLog)
    X0 = FFTLog.XArray[1]
    ZAr = FFTLog.ν .+ im .* FFTLog.ηM
    YArray = zeros(length(Ell), length(FFTLog.XArray))
    HMArray = zeros(ComplexF64, length(Ell), length(FFTLog.CM))
    FYArray = zeros(length(Ell), length(FFTLog.XArray))
    @inbounds for myl in 1:length(Ell)
        YArray[myl,:] = (Ell[myl] + 1) ./ reverse(FFTLog.XArray)
        HMArray[myl,:]  = FFTLog.CM .* (FFTLog.XArray[1] .* YArray[myl,1] ) .^
		(-im .*FFTLog.ηM) .* GL(Ell[myl], ZAr)
        FYArray[myl,:] = FFTW.irfft(conj(HMArray[myl,:]),
		length(FFTLog.XArray)) .* YArray[myl,:] .^ (-FFTLog.ν) .* sqrt(π) ./4
    end
    return YArray[:,FFTLog.NExtrapLow+FFTLog.NPad+1:FFTLog.NExtrapLow+
	FFTLog.NPad+FFTLog.OriginalLenght], FYArray[:,FFTLog.NExtrapLow+FFTLog.NPad+
	1:FFTLog.NExtrapLow+FFTLog.NPad+FFTLog.OriginalLenght]
end

function EvaluateFFTLogDJ(FFTLog::FFTLog, Ell::Vector{Float64})
    FFTLog.XArray = LogExtrap(FFTLog.XArray, FFTLog.NExtrapLow,
	FFTLog.NExtrapHigh)
    FFTLog.FXArray = LogExtrap(FFTLog.FXArray, FFTLog.NExtrapLow,
	FFTLog.NExtrapHigh)
    ZeroPad!(FFTLog)
    GetCM!(FFTLog)
    EvaluateηM!(FFTLog)
    X0 = FFTLog.XArray[1]
    ZAr = FFTLog.ν .+ im .* FFTLog.ηM
    YArray = zeros(length(Ell), length(FFTLog.XArray))
    HMArray = zeros(ComplexF64, length(Ell), length(FFTLog.CM))
    FYArray = zeros(length(Ell), length(FFTLog.XArray))
    @inbounds for myl in 1:length(Ell)
        YArray[myl,:] = (Ell[myl] + 1) ./ reverse(FFTLog.XArray)
        HMArray[myl,:]  = FFTLog.CM .* (FFTLog.XArray[1] .* YArray[myl,1] ) .^
		(-im .*FFTLog.ηM) .* GL1(Ell[myl], ZAr)
        FYArray[myl,:] = FFTW.irfft(conj(HMArray[myl,:]),
		length(FFTLog.XArray)) .* YArray[myl,:] .^ (-FFTLog.ν) .* sqrt(π) ./4
    end
    return YArray[:,FFTLog.NExtrapLow+FFTLog.NPad+1:FFTLog.NExtrapLow+
	FFTLog.NPad+FFTLog.OriginalLenght], FYArray[:,FFTLog.NExtrapLow+FFTLog.NPad+
	1:FFTLog.NExtrapLow+FFTLog.NPad+FFTLog.OriginalLenght]
end

function EvaluateFFTLogDDJ(FFTLog::FFTLog, Ell::Vector{Float64})
    FFTLog.XArray = LogExtrap(FFTLog.XArray, FFTLog.NExtrapLow,
	FFTLog.NExtrapHigh)
    FFTLog.FXArray = LogExtrap(FFTLog.FXArray, FFTLog.NExtrapLow,
	FFTLog.NExtrapHigh)
    ZeroPad!(FFTLog)
    GetCM!(FFTLog)
    EvaluateηM!(FFTLog)
    X0 = FFTLog.XArray[1]
    ZAr = FFTLog.ν .+ im .* FFTLog.ηM
    YArray = zeros(length(Ell), length(FFTLog.XArray))
    HMArray = zeros(ComplexF64, length(Ell), length(FFTLog.CM))
    FYArray = zeros(length(Ell), length(FFTLog.XArray))
    @inbounds for myl in 1:length(Ell)
        YArray[myl,:] = (Ell[myl] + 1) ./ reverse(FFTLog.XArray)
        HMArray[myl,:]  = FFTLog.CM .* (FFTLog.XArray[1] .* YArray[myl,1] ) .^
		(-im .*FFTLog.ηM) .* GL2(Ell[myl], ZAr)
        FYArray[myl,:] = FFTW.irfft(conj(HMArray[myl,:]),
		length(FFTLog.XArray)) .* YArray[myl,:] .^ (-FFTLog.ν) .* sqrt(π) ./4
    end
    return YArray[:,FFTLog.NExtrapLow+FFTLog.NPad+1:FFTLog.NExtrapLow+
	FFTLog.NPad+FFTLog.OriginalLenght], FYArray[:,FFTLog.NExtrapLow+FFTLog.NPad+
	1:FFTLog.NExtrapLow+FFTLog.NPad+FFTLog.OriginalLenght]
end


text:
@deprecate
function hamiltonian!(ds, s, p::ArrayPartition, t)
    J = p.x[1]
    h = p.x[2]
    anisotropy = p.x[3]
    hamiltonian!(ds, s, J, h, anisotropy, t)
end

function DTWAProblem(hamiltonian!, Ψ, saveat, p, )
    J = p.x[1]
    prob = ODEProblem(hamiltonian!, permutedims(dtwa_state(Ψ)), (saveat[1], saveat[end]), p,
        dense=false,
        dt=1/maximum(maximum(J, dims=1))/4,
        force_dtmin=true,
        adaptive=true,
        saveat=saveat
    )

    EnsembleProblem(
        prob,
        prob_func = (prob, i, repeat) -> dtwa_prob_func(prob, Ψ),
        output_func = (sol, i) -> (vcat(mean.(sol.u, dims=1)...), false),
        reduction = reduction,
        u_init = zeros(Float64, (length(saveat), 3)),
    )
end

DTWAProblem(hamiltonian!, p_0, tspan, p, s::Sequence) = DTWAProblem(
    hamiltonian!, p_0, tspan, p, s.durations
)

function solve_dtwa(ensemble_prob, trajectories; kwargs...)
    tspan = ensemble_prob.prob.tspan
    nsteps = size(ensemble_prob.u_init, 1)
    solve(
        ensemble_prob,
        Tsit5(),
        # AutoTsit5(Rosenbrock23()),
        # EnsembleSerial(),
        EnsembleDistributed(),
        trajectories=trajectories,
        batch=8;
        #saveat=collect(LinRange(tspan[1], tspan[2], nsteps));
        kwargs...
    )
end

function solve_dtwa(ensemble_prob, tracectories, s::Sequence; kwargs...)
    cb = IterativeCallback(s)
    solve_dtwa(ensemble_prob, tracectories; callback=cb, kwargs...)
end

max_stars_repo_name: titusfranz/DTWA

max_issues_repo_name: titusfranz/DTWA

max_forks_repo_name: titusfranz/DTWA


text:
mutable struct pbsrtfProblem{T<:Real}
    y       :: Vector{T}              # response
    x       :: Vector{T}              # design
    k       :: Int
    restriction :: String             # 'increasing','convex','decreasing','concave'
                                      # 'inc-convex','inc-concave','dec-convex','dec-concave'
    lb      :: Union{Vector{T},Nothing} # lower bound vector for β
    ub      :: Union{Vector{T},Nothing} # upper bound vector for β
    λ       :: Union{T,Nothing}       # Moreau-Yosida envelope parameter
    s       :: T                      # IG(s, r) prior for σ²
    r       :: T
    μ       :: Union{T,Nothing}       # exp(μ) prior for α
    m       :: Int                    # number of observations
    n       :: Int                    # number of grid locations
    xgrid   :: Vector{T}              # strictly increasing grid locations (length n)
    w       :: Vector{T}              # number of observations at each location (length n)
    ybar    :: Vector{T}              # the average of y at each location (length n)
    sse     :: T                      # sum of squared errors
    D¹      :: SparseMatrixCSC{T,Int}
    D²      :: SparseMatrixCSC{T,Int}
    D       :: SparseMatrixCSC{T,Int}
    GRB_ENV :: Any                    # gurobi environment
end

function pbsrtfProblem(y::Vector{T}, x::Vector{T}, k::Int, restriction::String;
                       lb::Union{Vector{T},Nothing}=nothing,
                       ub::Union{Vector{T},Nothing}=nothing,
                       thinning::Bool=false,nbins::Int=100,λ::Union{T,Nothing}=nothing,
                       s::T=1e-3,r::T=1e-3,μ::Union{T,Nothing}=nothing) where T<: Real
    if (k!=1) & (k!=2)
        throw(ArgumentError("Currenly only k=1 and k=2 are accepted for pbsrtf."))
    end
    m = length(x)
    xgrid = sort(unique(x))
    n = length(xgrid)
    w = zeros(T,n)
    ybar = zeros(T,n)
    sse = 0.0
    for i in 1:n
        w[i] = length(x[x.==xgrid[i]])
        ybar[i] = mean(y[x.==xgrid[i]])
        sse += sum(abs2, y[x.==xgrid[i]].-ybar[i])
    end
    if thinning
        xgrid,w,ybar,sse = thin(xgrid,w,ybar,sse,nbins)
        n = nbins
    end
    if λ == nothing
        λ = var(y)*0.25*1e-3
    end
    if μ == nothing
       μ = 3.0
    end
    D¹ = getD(0,n,xgrid)
    D² = getD(1,n,xgrid)
    D  = getD(k,n,xgrid)
    GRB_ENV = Env()
    GRBsetdblparam(GRB_ENV,"FeasibilityTol",1e-6)
    GRBsetdblparam(GRB_ENV,"OptimalityTol",1e-6)
    pbsrtfProblem{T}(y,x,k,restriction,lb,ub,λ,s,r,μ,m,n,xgrid,w,ybar,sse,D¹,D²,D,GRB_ENV)
end

dimension(problem::pbsrtfProblem) = problem.n + 2
capabilities(::Type{<:pbsrtfProblem}) = LogDensityOrder{1}()

function logdensity_and_gradient(problem::pbsrtfProblem, z)
    @unpack y,x,k,restriction,lb,ub,λ,s,r,μ,m,n,xgrid,w,ybar,sse,D¹,D²,D,GRB_ENV=problem
    β        = view(z,1:n)
    logσ²    = z[n+1]
    σ²       = exp(logσ²)
    logα     = z[n+2]
    α        = exp(logα)
    η = Variable(n)
    t = Variable()
    p_gurobi = minimize(sumsquares(β-η)+sumsquares(α-t),norm(D*η,1)<=t)
    if restriction=="increasing"
    elseif restriction=="decreasing"
        p_gurobi.constraints += D¹*η<=0
    elseif restriction=="convex"
        p_gurobi.constraints += D²*η>=0
    elseif restriction=="concave"
        p_gurobi.constraints += D²*η<=0
    elseif restriction=="inc-convex"
        p_gurobi.constraints += [D¹*η>=0,D²*η>=0]
    elseif restriction=="dec-convex"
        p_gurobi.constraints += [D¹*η<=0,D²*η>=0]
    elseif restriction=="inc-concave"
        p_gurobi.constraints += [D¹*η>=0,D²*η<=0]
    elseif restriction=="dec-concave"
        p_gurobi.constraints += [D¹*η<=0,D²*η<=0]
    end
    if !isnothing(lb)
        p_gurobi.constraints += η>=lb
    end
    if !isnothing(ub)
        p_gurobi.constraints += η<=ub
    end
    solve!(p_gurobi,Optimizer(GRB_ENV),silent_solver=true,warmstart=true,verbose=false)
    #while p_gurobi.status!=Convex.MOI.OPTIMAL
        #β = β .+ 1e-3 .*randn(n)
        #p_gurobi.objective = sumsquares(β-η)+sumsquares(α-t)
        #solve!(p_gurobi,Gurobi.Optimizer(GRB_ENV),silent_solver=true,warmstart=true)
        #println("Warning, solution failed!")
    #end
    β_proj = evaluate(η)
    α_proj = evaluate(t)
    res      = ybar .- β
    ssm      = sum(abs2.(res).*w)
    qf       = (ssm+sse+2*r)/(2*σ²)
    logl     = -qf-(m/2+s)*logσ²+logα-μ*α
    dist²    = norm(β.-β_proj)^2+(α-α_proj)^2
    logl    -= dist²/(2*λ)
    ∇        = similar(z)
    ∇[1:n]   = (w.*res)./σ².-((β.-β_proj)./λ)
    ∇[n+1]   = qf-(m/2+s)
    ∇[n+2]   = 1-α/λ*(α-α_proj)-μ*α
    logl, ∇
end


text:
using .Plots
function corner(
    table,
    labels = map(string, Tables.columnnames(table)),
    units = ["" for _ in labels],
    ;
    title="",
    plotcontours=true,
    plotscatter=true,
    plotpercentiles=[15,50,84],
    histfunc=prepare_hist,
    filterscatter=true,
    hist_kwargs=(;),
    hist2d_kwargs=(;),
    contour_kwargs=(;),
    scatter_kwargs=(;),
    percentiles_kwargs=(;),
    appearance=(;),
    lens=nothing,
    lens_kwargs=(;),
    bonusplot=nothing,

    # Format of 1D histogram titles
    titlefmt="\$\\mathrm{%s = %.2f^{+%.2f}_{-%.2f} \\; %s}\$",

    # By how much should we oversize the limits? Use this option with caution.
    lim_factor = 1.00,

    # Plot layout options (subject to change)
    # Overall plot scale factor
    sf = 1,
    pad_top = 20sf,
    pad_left = nothing,
    w=40sf,
    h=40sf,
    p=nothing, # 1sf noramlly, 10sf if 3D
    pad_right = 0sf,
    pad_bottom = 10sf,
    pad_bonus = nothing,
    # Remaining kwargs are for the overall plot group
    kwargs...
)

if !Tables.istable(table)
    error("You must supply data in a format consistent with the definition in Tables.jl, e.g. a Named Tuple of vectors.")
end

if any(!isascii, labels) ||  any(!isascii, units)
    @warn "Non-ascii labels or units detected. Some plotting backends require passing these using LaTeX escapes, e.g. \\alpha instead of α"
end

columns = Tables.columnnames(table)
n = length(columns)

# filter out empty (0 variance) variables
ii_empty = findall(==(0),  var.(collect(Tables.columns(table))))
if length(ii_empty) > 0
    @warn "Skipping entries with zero variance" skip=columns[ii_empty]
end

if isnothing(pad_left)
    pad_left = 20sf + 1sf*length(columns)
end

# Merge keywords for different subplots & series
hist_kwargs=merge((; nbins=20, color=:black, yticks=[], minorgrid=false, minorticks=false, ylims=(0,NaN)), hist_kwargs)
hist2d_kwargs=merge((; nbins=32, color=:Greys, colorbar=:none), hist2d_kwargs)
contour_kwargs=merge((; color=:black, linewidth=1.5), contour_kwargs)
scatter_kwargs=merge((; color=:black, alpha=0.5, markersize=0.6, markerstrokewidth=0), scatter_kwargs)
percentiles_kwargs=merge((;linestyle=:dash, color=:black), percentiles_kwargs)
appearance = merge((;
    framestyle=:box,
    xrotation = 45,
    yrotation = 45,
    # margin=-6Measures.mm,
    grid=:none,
    xformatter=:plain,
    yformatter=:plain,
    zformatter=t->"",
    titlefontsize=10,
    tick_direction=:out,
    label="",
), appearance)

threeD = get(hist2d_kwargs, :seriestype, nothing) == :wireframe
# Increase padding for 3D plots
if isnothing(p)
    p = threeD ? 10sf : 1sf
end

if isnothing(pad_bonus)
    pad_bonus = iseven(n) ? 20sf : 5sf
end

# Calculate a default reasonable size unless overriden
# s = 250*n
# kwargs = merge((;size = (s,s)), kwargs)

# Rather than messing with linking axes, let's just take over setting the plot limits ourselves
ex = extrema.(collect(Tables.columns(table)))

mids = [(ex[1]+ex[2])/2 for ex in ex]
widths = [(ex[2]-ex[1])/2 for ex in ex]
mins = [mid-width*lim_factor for (mid,width) in zip(mids,widths)]
maxs = [mid+width*lim_factor for (mid,width) in zip(mids,widths)]
# mins = [ex[1] for ex in ex]
# maxs = [ex[2] for ex in ex]

full_width = (pad_left+(w+p)*n+pad_right)
full_height = (pad_top+(h+p)*n+pad_bottom)
parent = RecipesBase.plot(
    # size=((w+p)*n, ((w+p)*n+pad_top)),
    size=(full_width*4,full_height*4),
    framestyle=:none,
    title=title;
    # left_margin=-10Measures.mm,
    background_color_inside=:transparent,
    kwargs...
)

fmt(label, unit) = "\$\\mathrm{$label \\; $unit}\$"

# Build grid of nxn plots
k = 1
for row in 1:n, col in 1:n

    if row < col
        # push!(subplots, emptyplot())
        continue
    end
    k+=1

    # fmt(label) = "\$\\mathrm{$label}\$"
    kw = (;
        title = row == col ? labels[row] : "",
        xguide = threeD || row == n ? fmt(labels[col], units[col]) : "",
        yguide = col == 1  && row > 1 ? fmt(labels[row], units[row]) : "",
        xformatter = threeD || row == n ? appearance.xformatter : t->"",
        yformatter = threeD || col == 1 && row > 1 ? appearance.yformatter : t->"",

        # For crudlely tweaking aspect ratios of the plots to account for
        # top_margin = row == 1 ? 10Measures.mm : :match,
        # left_margin = col == 1 ? 8Measures.mm : :match,
        # right_margin = col == n ? 6Measures.mm : :match,

        xlims=(mins[col], maxs[col]),
        label="",
        subplot=k,
    )
    x=(col-1)*(w+p)+pad_left
    y=(row-1)*(h+p)+pad_top
    inset = bbox(x*Measures.mm, y*Measures.mm, w*Measures.mm, h*Measures.mm, :left)

    if row != col
        kw = (;kw..., ylims=(mins[row], maxs[row]))
    end

    table_keys = Tables.columnnames(table)
    if row == col
        # 1D histogram
        unit = units[row]
        hist(Tables.getcolumn(table, table_keys[row]), histfunc, merge(appearance, kw, hist_kwargs, (;inset)), plotpercentiles, merge(kw, percentiles_kwargs), titlefmt, unit)
    else row > col
        # 2D histogram
        hist(Tables.getcolumn(table, table_keys[row]), Tables.getcolumn(table, table_keys[col]), histfunc, merge(appearance, kw, hist2d_kwargs, (;inset)), merge(kw,contour_kwargs), merge(kw,scatter_kwargs), plotcontours, plotscatter, filterscatter)
    end
    # push!(subplots, subplot)
end

if !isnothing(lens)
    k += 1
    nspan = floor(mean(1:n))
    x=nspan*(w+p)+pad_bonus+pad_left
    y=pad_top
    bonuswidth = full_width-x
    bonusheight = full_height-pad_bottom-pad_top-nspan*(h+p) - pad_bonus
    inset = bbox(x*Measures.mm, y*Measures.mm, bonuswidth*Measures.mm, bonusheight*Measures.mm, :left)
    if lens ∈ columns
        kw = (;label="", subplot=k, title=labels[findfirst(==(lens), columns)])
        hist(Tables.getcolumn(table, lens), histfunc, merge(appearance, kw, hist_kwargs, (;inset), get(lens_kwargs,:hist_kwargs,(;))), plotpercentiles, merge(kw, percentiles_kwargs, get(lens_kwargs,:hist_kwargs,(;))))
    elseif lens[1] ∈ columns && lens[2] ∈ columns
        kw = (;label="", subplot=k, yguide=fmt(labels[findfirst(==(lens[1]), columns)]), xguide=fmt(labels[findfirst(==(lens[2]), columns)]))
        hist(Tables.getcolumn(table, lens[1]), Tables.getcolumn(table, lens[2]), histfunc, merge(appearance, kw, hist2d_kwargs, (;inset), lens_kwargs, get(lens_kwargs,:hist2d_kwargs,(;))), merge(kw,contour_kwargs,get(lens_kwargs,:contour_kwargs,(;))), merge(kw,scatter_kwargs,get(lens_kwargs,:scatter_kwargs,(;))), get(lens_kwargs,:plotcontours,plotcontours), get(lens_kwargs,:plotscatter,plotscatter), get(lens_kwargs,:filterscatter,filterscatter))
    end
elseif !isnothing(bonusplot)
    k += 1
    nspan = floor(mean(1:n))
    x=nspan*(w+p)+pad_bonus+pad_left
    y=pad_top
    bonuswidth = full_width-x-pad_bonus
    bonusheight = full_height-pad_bottom-pad_top-nspan*(h+p) - pad_bonus
    inset = bbox(x*Measures.mm, y*Measures.mm, bonuswidth*Measures.mm, bonusheight*Measures.mm, :left)
    kw = (;inset,subplot=k)
    bonusplot(kw)
end

return parent
end

max_stars_repo_name: sefffal/CornerPlots

max_issues_repo_name: sefffal/CornerPlots


text:
using DataConvenience
using DataFrames
using CSV
using Test
using Random: randstring

@testset "read csv in chunks" begin
    filepath = "tmp.csv"

    M = 1000
    str_base = [randstring(8) for i in 1:1_000]
    @time df = DataFrame(int = rand(Int32, M), float=rand(M), str = rand(str_base, M))

    @time CSV.write(filepath, df)

    # read the file 100 bytes at a time
    chunks = CsvChunkIterator(filepath, 100)

    dfs = [DataFrame(chunk) for chunk in chunks]

    made = reduce(vcat, dfs)
    actual = CSV.read(filepath, DataFrame)
    @test nrow(made) == nrow(actual)
    @test ncol(made) == ncol(actual)

    # read the file 500 bytes
    chunks = CsvChunkIterator(filepath, 500)
    dfs = [DataFrame(chunk) for chunk in chunks]

    made = reduce(vcat, dfs)
    @test nrow(made) == nrow(actual)
    @test ncol(made) == ncol(actual)

    # read the file 7 rows at a time
    # where the file is of size 32 rows
    chunks = CsvChunkIterator(filepath)
    dfs = [DataFrame(chunk) for chunk in chunks]

    made = reduce(vcat, dfs)
    @test nrow(made) == nrow(actual)
    @test ncol(made) == ncol(actual)

    collect(CsvChunkIterator(filepath; header=0))
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
# using Images, BenchmarkTools

# const c = complex(-.79,.15)

# function calc(z::ComplexF64)::RGB
#     for i in 1:10000
#         z = z*z + c
#         if abs(z) > 2
#             return RGB(i > 255 ? 1. : i/255, 0., 0.)
#         end
#     end
#     return RGB(0.,0.,0.)
# end

# function main(height::Int,width::Int)::Matrix{RGB}
#     img::Matrix{RGB} = zeros(RGB,height,width)
#     for x in 1:width
#         a = x * 3.75 / width - 1.9
#         Threads.@threads for y in 1:height
#             b = y * 2.5 / height - 1.25
#             img[y,x] = calc(complex(a,b))
#         end
#     end
#     return img
# end

# @benchmark main(1000,1500)
# #@time save("julia.png", main(1000,1500))
using JLD
points = open("./sets/")

max_stars_repo_name: BoardTM/fractals

max_issues_repo_name: BoardTM/fractals

max_forks_repo_name: BoardTM/fractals


text:
using SpikingNetworks
parseParameters"""
τ=10
τᵢ=5
σ=0.015
σᵢ=20"""
parseEquations"""
dv/dt=(0.9+0.5*I-v+σ*ξ(t))/τ
dI/dt=(-I+σᵢ*η(t))/τᵢ"""
net=init"""
N=25
v=rand(N)
I=rand()*ones(N)"""
spike"""
v.>1
5"""
reset"v[i]=0"
net.dt=0.1
run!(net,500/net.dt)

using MatlabPlot
t,neuron=rasterPlot(net)
mplot(t,neuron,".k")


text:
module cumulantexpansion

using ..bases
using ..states
using ..operators
using ..operators_lazy
using ..ode_dopri

import Base: *, full
import ..operators

type ProductDensityOperator <: Operator
    basis_l::CompositeBasis
    basis_r::CompositeBasis
    operators::Vector{DenseOperator}
    function ProductDensityOperator(operators::DenseOperator...)
        basis_l = tensor([op.basis_l for op in operators]...)
        basis_r = tensor([op.basis_r for op in operators]...)
        new(basis_l, basis_r, collect(operators))
    end
end

ProductDensityOperator() = error("ProductDensityOperator needs at least one operator.")
ProductDensityOperator(states::Ket...) = ProductDensityOperator(DenseOperator[tensor(state, dagger(state)) for state in states]...)

*(a::ProductDensityOperator, b::ProductDensityOperator) = (check_multiplicable(a.basis_r, b.basis_l); ProductDensityOperator(a.basis_l, b.basis_r, DenseOperator[a_i*b_i for (a_i, b_i) in zip(a.operators, b.operators)]))
function *(a::LazyTensor, b::ProductDensityOperator)
    check_multiplicable(a.basis_r, b.basis_l);
    operators = DenseOperator[]
    for (alpha, b_alpha) in enumerate(b.operators)
        if alpha in keys(a.operators)
            push!(operators, a.operators[alpha]*b_alpha)
        else
            push!(operators, deepcopy(b_alpha))
        end
    end
    if a.factor != 1.
        operators[0] *= a.factor
    end
    ProductDensityOperator(operators...)
end
function *(a::ProductDensityOperator, b::LazyTensor)
    check_multiplicable(a.basis_r, b.basis_l);
    operators = DenseOperator[]
    for (alpha, a_alpha) in enumerate(a.operators)
        if alpha in keys(b.operators)
            push!(operators, a_alpha*b.operators[alpha])
        else
            push!(operators, deepcopy(a_alpha))
        end
    end
    if b.factor != 1.
        operators[0] *= b.factor
    end
    ProductDensityOperator(operators...)
end

dims(bl::CompositeBasis, br::CompositeBasis) = [length(bl_i)*length(br_i) for (bl_i, br_i) in zip(bl.bases, br.bases)]
dims(x::ProductDensityOperator) = dims(x.basis_l, x.basis_r)

operators.expect(op::Operator, rho::ProductDensityOperator) = trace(op*rho)
operators.expect(op::LazySum, rho::ProductDensityOperator) = sum([expect(x, rho) for (f, x) in zip(op.factors, op.operators)])

traces(x::ProductDensityOperator) = [trace(op) for op in x.operators]
operators.trace(x::ProductDensityOperator) = prod(traces(x))

function Base.full(x::ProductDensityOperator)
    tensor(x.operators...)
end

function fill!(x::ProductDensityOperator, alpha::Number)
    for op in x.operators
        Base.fill!(op.data, complex(alpha))
    end
end

# Ignores the factor in the LazyTensor
function operators.gemm!(alpha, a::LazyTensor, b::ProductDensityOperator, beta, result::ProductDensityOperator)
    @assert abs(beta)==0.
    for (k, a_k) in a.operators
        operators.gemm!(complex(1.), a_k, b.operators[k], complex(0.), result.operators[k])
    end
end

function dmaster(rho0::ProductDensityOperator, H::LazySum,
                 J::Vector{LazyTensor}, Jdagger::Vector{LazyTensor},
                 drho::ProductDensityOperator, tmp::ProductDensityOperator, tmp2::ProductDensityOperator)
    fill!(drho, 0.)
    for h_k in H.operators
        operators.gemm!(1., h_k, rho0, 0., tmp)
        subtraces = traces(tmp)
        for (alpha, h_k_alpha) in h_k.operators
            factor = h_k.factor
            for gamma in keys(h_k.operators)
                if alpha!=gamma
                    factor *= subtraces[gamma]
                end
            end
            operators.gemm!(factor*complex(0,-1.), h_k_alpha, rho0.operators[alpha], complex(1.), drho.operators[alpha])
            operators.gemm!(factor*complex(0,1.), rho0.operators[alpha], h_k_alpha, complex(1.), drho.operators[alpha])
        end
    end
    for k=1:length(J)
        operators.gemm!(1., J[k], rho0, 0., tmp)
        operators.gemm!(1., Jdagger[k], tmp, 0., tmp2)
        subtraces = traces(tmp)
        subindices = keys(J[k].operators)
        for alpha in subindices
            factor = Jdagger[k].factor * J[k].factor
            for gamma in subindices
                if alpha!=gamma
                    factor *= subtraces[gamma]
                end
            end
            operators.gemm!(complex(factor), J[k].operators[alpha], rho0.operators[alpha], complex(0.), tmp.operators[alpha])
            operators.gemm!(complex(1.), tmp.operators[alpha], Jdagger[k].operators[alpha], complex(1.), drho.operators[alpha])

            operators.gemm!(complex(-0.5), Jdagger[k].operators[alpha], tmp.operators[alpha], complex(1.), drho.operators[alpha])

            operators.gemm!(complex(-0.5*factor), rho0.operators[alpha], Jdagger[k].operators[alpha], complex(0.), tmp.operators[alpha])
            operators.gemm!(complex(1.), tmp.operators[alpha], J[k].operators[alpha], complex(1.), drho.operators[alpha])
        end
    end
end

function dmaster_timedependent(t::Float64, rho0::ProductDensityOperator, f::Function,
                 drho::ProductDensityOperator, tmp::ProductDensityOperator, tmp2::ProductDensityOperator)
    H, J, Jdagger = f(t, rho0)
    dmaster(rho0, H, J, Jdagger, drho, tmp, tmp2)
end

dims(rho::ProductDensityOperator) = [length(op.basis_l)*length(op.basis_r) for op in rho.operators]

function as_vector(rho::ProductDensityOperator, x::Vector{Complex128})
    @assert length(x) == sum(dims(rho))
    i = 0
    for op in rho.operators
        N = length(op.basis_l)*length(op.basis_r)
        x[i+1:i+N] = reshape(op.data, N)
        i += N
    end
    x
end

function as_operator(x::Vector{Complex128}, rho::ProductDensityOperator)
    @assert length(x) == sum(dims(rho))
    i = 0
    for op in rho.operators
        N = length(op.basis_l)*length(op.basis_r)
        reshape(op.data, N)[:] = x[i+1:i+N]
        i += N
    end
    rho
end

function integrate_master(dmaster::Function, tspan, rho0::ProductDensityOperator;
                fout::Union{Function,Void}=nothing, kwargs...)
    x0 = as_vector(rho0, zeros(Complex128, sum(dims(rho0))))
    f = (x->x)
    if fout==nothing
        tout = Float64[]
        xout = ProductDensityOperator[]
        function fout_(t, rho::ProductDensityOperator)
            push!(tout, t)
            push!(xout, deepcopy(rho))
        end
        f = fout_
    else
        f = fout
    end
    tmp = deepcopy(rho0)
    f_(t, x::Vector{Complex128}) = f(t, as_operator(x, tmp))
    ode(dmaster, float(tspan), x0, f_; kwargs...)
    return fout==nothing ? (tout, xout) : nothing
end

function master(tspan, rho0::ProductDensityOperator, H::LazySum, J::Vector{LazyTensor};
                fout::Union{Function,Void}=nothing,
                kwargs...)

    Jdagger = LazyTensor[dagger(j) for j=J]
    rho = deepcopy(rho0)
    drho = deepcopy(rho0)
    tmp = deepcopy(rho0)
    tmp2 = deepcopy(rho0)
    function dmaster_(t, x::Vector{Complex128}, dx::Vector{Complex128})
        dmaster(as_operator(x, rho), H, J, Jdagger, drho, tmp, tmp2)
        as_vector(drho, dx)
    end
    integrate_master(dmaster_, tspan, rho0; fout=fout, kwargs...)
end

master(tspan, rho0, H::LazyTensor, J; kwargs...) = master(tspan, rho0, LazySum(H), J; kwargs...)

"""
f
    Function f(t, rho_t) -> (H_t, J_t, Jdagger_t)
"""
function master_timedependent(tspan, rho0::ProductDensityOperator, f::Function;
                fout::Union{Function,Void}=nothing,
                kwargs...)
    rho = deepcopy(rho0)
    drho = deepcopy(rho0)
    tmp = deepcopy(rho0)
    tmp2 = deepcopy(rho0)
    function dmaster_(t, x::Vector{Complex128}, dx::Vector{Complex128})
        dmaster_timedependent(t, as_operator(x, rho), f, drho, tmp, tmp2)
        as_vector(drho, dx)
    end
    integrate_master(dmaster_, tspan, rho0; fout=fout, kwargs...)
end

end # module