text:
# adaptive step
# with noise
@everywhere using Distributions
@everywhere using StatsBase
@everywhere using DifferentialEquations
@everywhere using Interpolations
include("../src/SMCABC_RKHS_RJ.jl")
include("Models_rj.jl")
cd(@__DIR__)

# data generation for test---------
@everywhere function generate_dataset()
    timepoints = linspace(0,5, 500)
    μ = 0
    σ = 10
    kernel = Normal(μ, σ)
    x=rand(kernel, length(timepoints));

    itp = interpolate(x, BSpline(Linear()), OnGrid());
    sitp = Interpolations.scale(itp, timepoints);
    model = model_generator(0, sitp)

    y0 = [0.0,0.0]
    tspan = (timepoints[1], timepoints[end])
    params = [1., 0.05, 50., 1e+3]

    alg = CVODE_BDF(max_hnil_warns = 0, max_error_test_failures = 0)
    prob = ODEProblem(model, y0, tspan, params)
    res = solve(prob, alg, saveat=timepoints, maxiters=1e7, abstol = 1e-8, reltol = 1e-8)
    return [hcat(res.u...)']
end
@everywhere dataset = generate_dataset()
#--------------------------------
srand(9357478685)

@everywhere N_model = 2
N_particle = 5000
ϵ₀ = Inf
Δϵ_tol = 1e-8
α = 0.98
target = 1-α
N_sample = 100
Δϵ_window = 20
save_interval = 50
parallel = true
output = joinpath(@__DIR__, "result_RKHS_rj2")
@everywhere μ = 0
@everywhere σ = 10
@everywhere timepoints = linspace(0,5, 500)
@everywhere lbs_set = [[0.1, 0.005, 5., 10^2],[0.1, 0.005, 5., 10^2, 10^4]]
@everywhere ubs_set = [[10, 0.5, 500., 10^4],[10, 0.5, 500., 10^4, 10^6]]
@everywhere kernel = Normal(μ, σ)
@everywhere x=rand(kernel, length(timepoints));
@everywhere itp = interpolate(x, BSpline(Linear()), OnGrid());
@everywhere sitp = Interpolations.scale(itp, timepoints);
@everywhere models = [model_generator(i-1, sitp) for i in 1:N_model]
@everywhere obj_funcs = [f(param) = obj_func(dataset, models[i], param, lbs_set[i], ubs_set[i]) for i in 1:N_model]

run_rjsmc(N_particle, lbs_set, ubs_set, obj_funcs, ϵ₀, Δϵ_tol, α, Δϵ_window, output;
          save_interval=save_interval, target=target, N_sample=N_sample, parallel=parallel)

# plot -----------------------------------------------
# est_params = readcsv("result_rj/est_param_M2_655.csv")
#
# nrow=5
# ncol=1
# figsize=(5,10)
# fig, axes = subplots(nrow, ncol, figsize=figsize);
#
# for i in 1:5
#     axes[i][:hist](est_params[:,i], bins=100, normed=true, color="b");
#     axes[i][:axvline](params[i], color="r");
# end
# fig[:savefig](joinpath(output, "param_dist.png"), dpi=300);
# PyPlot.plt[:close]()
#
#
# y0 = [0.0,0.0]
# tspan = (timepoints[1], timepoints[end])
# params = [1., 0.05, 50., 1e+3, 1e+5]
#
# trajectories = zeros(size(est_params)[1], length(timepoints), 2)
# for i in 1:size(est_params)[1]
#     param = est_params[i,:]
#     prob = ODEProblem(model, y0, tspan, param)
#     res_temp = solve(prob, alg, saveat=timepoints, maxiters=1e7,
#                 abstol = 1e-8, reltol = 1e-8)
#     trajectories[i,:,:] = hcat(res_temp.u...)'
# end
#
# ave = mean(trajectories, Dims(1))
# sd = std(trajectories, Dims(1))
# y1=ave.-2sd
# y2=ave.+2sd
#
# prob = ODEProblem(model, y0, tspan, params)
# res = solve(prob, alg, saveat=timepoints, maxiters=1e7,
#             abstol = 1e-8, reltol = 1e-8)
#
# plot(timepoints, vec(ave[:,:,1]), color="red", ls="--", lw=1);
# fill_between(timepoints, vec(y1[:,:,1]), vec(y2[:,:,1]));
#
# plot(timepoints, res[1, :], lw=1, color="k", ls="--");
#
# savefig(joinpath(output, "predict.png"), dpi=300);
# PyPlot.plt[:close]()


text:
module Optimise

using LinearAlgebra
using ProximalOperators

abstract type AbstractOptimiser end;
abstract type AbstractSubspaceOptimiser end

include("./strridge.jl")
include("./admm.jl")
include("./sr3.jl")

#Nullspace for implicit sindy
include("./adm.jl")

export init, init!, fit!, set_threshold!, set_threshold
export STRRidge, ADMM, SR3
export ADM

end

ext: jl


text:
@testset "Day5.jl" begin
    # Write your tests here.
    @test AdventOfCode2020.row("FFFFFFF") == 0
    @test AdventOfCode2020.row("BBBBBBB") == 127
    @test AdventOfCode2020.row("FFFFFFB") == 1
    @test AdventOfCode2020.row("FBFBBFF") == 44

    @test AdventOfCode2020.col("LLL") == 0
    @test AdventOfCode2020.col("RRR") == 7
    @test AdventOfCode2020.col("LLR") == 1
    @test AdventOfCode2020.col("RLR") == 5

    @test AdventOfCode2020.toseat("FBFBBFFRLR") == AdventOfCode2020.SeatLoc(44,5)
    @test AdventOfCode2020.toseat("FBFBBFFRLR").id == 357
    @test AdventOfCode2020.toseat("BFFFBBFRRR").id == 567
    @test AdventOfCode2020.toseat("FFFBBBFRRR").id == 119
    @test AdventOfCode2020.toseat("BBFFBBFRLL").id == 820
end


text:
abstract type AbstractTestClass end
function func()::String
    return "test"
end

mutable struct TestClass <: AbstractTestClass
end
function func(self::TestClass)::String
    return "test2"
end

function test()
    num::Int64 = 2
    teststr::String = "test"
    function inner_test()::String
        return repeat(teststr, num)
    end

    function inner_test_2(num)::String
        return num * teststr
    end

    @assert(inner_test() == "testtest")
    @assert(inner_test_2(4) == "testtesttesttest")
end

if abspath(PROGRAM_FILE) == @__FILE__
    @assert(func() == "test")
    testClass = TestClass()
    @assert(func(testClass) == "test2")
    test()
end


text:
"""
Training functions.
"""

"""
Custom training loop to save and plot losses.
"""
function logging_train!(loss, ps, data, opt, losses)
    for (i,d) in enumerate(data)
        # back is a method that computes the product of the gradient so far with its argument.
        train_loss, back = Zygote.pullback(() -> loss(d...), ps)
        # Insert whatever code you want here that needs training_loss, e.g. logging.
        # Apply back() to the correct type of 1.0 to get the gradient of loss.
        gs = back(one(train_loss))
        push!(losses, train_loss)
        # Insert what ever code you want here that needs gradient.
        # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
        Flux.update!(opt, ps, gs)
        # Here you might like to check validation set accuracy, and break out to do early stopping.
    end
end

"""
Handles the entire training session. Saves best model to "dir/models/ONTOLOGY/model_name"
"""
function Train(model, dir,train_data, test_data, opt, loss, accuracy; epochs=100::Integer)
    ps = Flux.params(model)
    best_acc = 0.0
    # best_loss = 9999999
    last_improvement = 0
    test_x = test_data[1]
    test_y = test_data[2]
    model_name = ""
    temp_dir = "$dir/models/temp"
    loss_log = []
    @info("Beginning training loop...")
    for epoch_idx in 1:epochs
        losses = []
        # Train for a single epoch
        # Flux.train!(loss, ps, train_data, opt)
        logging_train!(loss, ps, train_data, opt, losses)
        # Calculate loss
        last_loss = losses[end]
        append!(loss_log, last_loss)
        # Calculate accuracy:
        acc = accuracy(test_x, test_y)
        @info(@sprintf("[%d]: Validation loss: %.6f", epoch_idx, last_loss))
        @info(@sprintf("[%d]: Validation accuracy: %f", epoch_idx, acc))

        # If our accuracy is good enough, quit out.
        if acc >= 0.999
            @info(" -> Early-exiting: We reached our target accuracy of 99.9%")
            break
        end

        if acc > best_acc
            model_name = "$ONTOLOGY-$(now())_acc_$(acc*100).bson"
            @info(" -> New best loss! Saving model out to $temp_dir")
            if !ispath(temp_dir)
                mkpath(temp_dir)
            else
                rm(temp_dir, recursive=true)
                mkpath(temp_dir)
            end
            BSON.@save "$temp_dir/$model_name" model epoch_idx opt acc
            last_improvement = epoch_idx
            best_acc = acc
        end

        # If we haven't seen improvement in 20 epochs, drop our learning rate:
        if epoch_idx - last_improvement >= 5 && opt.eta > 1e-6
            opt.eta /= 10.0
            @warn(" -> Haven't improved in a while, dropping learning rate to $(opt.eta)!")

            # After dropping learning rate, give it a few epochs to improve
            last_improvement = epoch_idx
        end

        if epoch_idx - last_improvement >= 10
            @warn("Accuracy no longer increasing. Exiting early.")
            BSON.@load "$temp_dir/$model_name" model_to_save epoch_idx opt acc
            BSON.@save "$dir/models/$ONTOLOGY/$model_name" model_to_save epoch_idx opt acc
            rm(temp_dir, recursive=true)
            break
        end
    end
    ## Export loss log for plotting
    export_csv(losses, "$dir/models/$ONTOLOGY/logs", "$model_name.csv")
end


text:
using Random

mutable struct SocietyMember <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    mass::Float64
    days_infected::Int  # number of days since is infected
    status::Symbol  # :S, :I or :R
    β::Float64
end

const steps_per_day = 24

"""
```julia
social_distancing(;
    infection_period = 30 * steps_per_day,
    detection_time = 14 * steps_per_day,
    reinfection_probability = 0.05,
    isolated = 0.5, # in percentage
    interaction_radius = 0.012,
    dt = 1.0,
    speed = 0.002,
    death_rate = 0.044, # from website of WHO
    N = 1000,
    initial_infected = 5,
    seed = 42,
    βmin = 0.4,
    βmax = 0.8,
)
```
Same as in [Continuous space social distancing for COVID-19](@ref).
"""
function social_distancing(;
        infection_period = 30 * steps_per_day,
        detection_time = 14 * steps_per_day,
        reinfection_probability = 0.05,
        isolated = 0.5, # in percentage
        interaction_radius = 0.012,
        dt = 1.0,
        speed = 0.002,
        death_rate = 0.044, # from website of WHO
        N = 1000,
        initial_infected = 5,
        seed = 42,
        βmin = 0.4,
        βmax = 0.8,
        spacing = 0.02
    )

    properties = Dict(
        :infection_period => infection_period,
        :reinfection_probability => reinfection_probability,
        :detection_time => detection_time,
        :death_rate => death_rate,
        :interaction_radius => interaction_radius,
        :dt => dt,
    )
    space = ContinuousSpace((1, 1), spacing)
    model = ABM(SocietyMember, space, properties = properties)

    Random.seed!(seed)
    for ind in 1:N
        pos = Tuple(rand(2))
        status = ind ≤ N - initial_infected ? :S : :I
        isisolated = ind ≤ isolated * N
        mass = isisolated ? Inf : 1.0
        vel = isisolated ? (0.0, 0.0) : sincos(2π * rand()) .* speed
        β = (βmax - βmin) * rand() + βmin
        add_agent!(pos, model, vel, mass, 0, status, β)
    end

    return model, social_distancing_agent_step!, social_distancing_model_step!
end

function sd_transmit!(a1, a2, rp)
    count(a.status == :I for a in (a1, a2)) ≠ 1 && return
    infected, healthy = a1.status == :I ? (a1, a2) : (a2, a1)

    rand() > infected.β && return

    if healthy.status == :R
        rand() > rp && return
    end
    healthy.status = :I
end

function social_distancing_model_step!(model)
    r = model.interaction_radius
    for (a1, a2) in interacting_pairs(model, r, :nearest)
        sd_transmit!(a1, a2, model.reinfection_probability)
        elastic_collision!(a1, a2, :mass)
    end
end

function social_distancing_agent_step!(agent, model)
    move_agent!(agent, model, model.dt)
    sd_update!(agent)
    social_distancing_recover_or_die!(agent, model)
end

sd_update!(agent) = agent.status == :I && (agent.days_infected += 1)

function social_distancing_recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand() ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end


text:
using Indicators
using Base.Test

# Read sample data and calculations
const tol = 1e-8
data = X = readcsv("AAPL-2015.csv", header=true)[1]
ma = readcsv("ma.csv", header=true)[1]
o = data[:,1]
h = data[:,2]
l = data[:,3]
c = x = data[:,4]
v = data[:,5]
hlc = data[:,2:4]
ohlc = data[:,1:4]

# Compare against `Indicators` package function results
@test maximum(sma(x)[10:end]-ma[10:end,1]) < tol
@test maximum(ema(x)[10:end]-ma[10:end,2]) < tol
@test maximum(wma(x)[10:end]-ma[10:end,3]) < tol
@test maximum(hma(x)[23:end]-ma[23:end,4]) < tol


text:
"""
    Chain(layers...)

Chain multiple layers / functions together, so that they are called in sequence
on a given input.

`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.

# Examples
```jldoctest
julia> m = Chain(x -> x^2, x -> x+1);

julia> m(5) == 26
true

julia> m = Chain(Dense(10, 5), Dense(5, 2));

julia> x = rand(10);

julia> m(x) == m[2](m[1](x))
true
```
"""
struct Chain{T<:Tuple}
  layers::T
  Chain(xs...) = new{typeof(xs)}(xs)
end

@forward Chain.layers Base.getindex, Base.length, Base.first, Base.last,
  Base.iterate, Base.lastindex

functor(::Type{<:Chain}, c) = c.layers, ls -> Chain(ls...)

applychain(::Tuple{}, x) = x
applychain(fs::Tuple, x) = applychain(tail(fs), first(fs)(x))

(c::Chain)(x) = applychain(c.layers, x)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)

function Base.show(io::IO, c::Chain)
  print(io, "Chain(")
  join(io, c.layers, ", ")
  print(io, ")")
end

# This is a temporary and naive implementation
# it might be replaced in the future for better performance
# see issue https://github.com/FluxML/Flux.jl/issues/702
# Johnny Chen -- @johnnychen94
# only slightly changed to better handle interaction with Zygote @dsweber2
"""
    activations(c::Chain, input)

Calculate the forward results of each layers in Chain `c` with `input` as model input.
"""
function activations(c::Chain, input)
    extraChain(c.layers, input)
end

function extraChain(fs::Tuple, x)
    res = first(fs)(x)
    return (res, extraChain(Base.tail(fs), res)...)
end

extraChain(::Tuple{}, x) = ()

"""
    Dense(in, out, σ=identity; initW=glorot_uniform, initb=zeros, bias=true)
    Dense(W, b, σ=identity)

Create a traditional `Dense` layer with in×out weight matrix `W` and
bias vector  `b` of length `out`. The forward pass is given by:

    y = σ.(W * x .+ b)

The input `x` must be a vector of length `in`, a batch of vectors represented
as an `in × N` matrix, or a higher order tensor where all dimensions
after the first one will be treated as batch dimensions.

The out `y` will be a vector  of length `out` or a batch whose first
dimension is `out` and the remaining dimensions are the same as in the input.

Setting `bias` to `false` will switch the bias  off for the layer.

`initW` and `initb` are callables used to initialize weights and biases respectively,
through the calls `initW(out, in)` and `initb(out)`.

# Examples

```julia-repl
julia> d = Dense(5, 2)
Dense(5, 2)

julia> d(rand(Float32, 5))
2-element Array{Float32,1}:
 -0.16210233
  0.123119034

julia> d = Dense(5, 2; bias=false)
Dense(5, 2)
```
"""
struct Dense{F,S<:AbstractArray,T<:Union{Zeros, AbstractVector}}
  W::S
  b::T
  σ::F
end

Dense(W, b) = Dense(W, b, identity)

function Dense(in::Integer, out::Integer, σ = identity;
               initW = glorot_uniform, initb = zeros, bias=true)
  return Dense(initW(out, in), create_bias(bias, initb, out), σ)
end

@functor Dense

function (a::Dense)(x::AbstractArray)
  W, b, σ = a.W, a.b, a.σ
  sz = size(x)
  x = reshape(x, sz[1], :) # reshape to handle dims > 1 as batch dimensions
  x = σ.(W*x .+ b)
  return reshape(x, :, sz[2:end]...)
end

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

"""
    Diagonal(α, β)
    Diagonal(sz::Integer...; initα=ones, initβ=zeros)

Create an element-wise linear layer with learnable
arrays `α` and `β` of size `sz`. The layer performs

    y = α .* x .+ β

The input `x` must have size broadcast-compatible with `α` and `β`.
The parameters will be created with the calls
`α = initα(sz)` and `β = initβ(sz)`.
"""
struct Diagonal{T}
  α::T
  β::T
end

function Diagonal(sz::Integer...;
      initα = i -> ones(Float32, i),
      initβ = i -> zeros(Float32, i))
  Diagonal(initα(sz), initβ(sz))
end

@functor Diagonal

(a::Diagonal)(x) = a.α .* x .+ a.β

function Base.show(io::IO, l::Diagonal)
  print(io, "Diagonal(", size(l.α), ")")
end

"""
    Maxout(over)

The [Maxout](https://arxiv.org/abs/1302.4389) layer has a number of
internal layers which all receive the same input. It returns the elementwise
maximum of the internal layers' outputs.

Maxout over linear dense layers satisfies the univeral approximation theorem.
"""
struct Maxout{FS<:Tuple}
    over::FS
end

"""
    Maxout(f, n_alts)

Construct a Maxout layer over `n_alts` instances of the layer given by `f`.
The function takes no arguments and should return some callable layer.
Conventionally, this is a linear dense layer.

# Examples

This constructs a `Maxout` layer over 4 internal dense linear layers, each
identical in structure (784 inputs, 128 outputs):
```jldoctest
julia> insize = 784;

julia> outsize = 128;

julia> Maxout(()->Dense(insize, outsize), 4);
```
"""
function Maxout(f, n_alts)
  over = Tuple(f() for _ in 1:n_alts)
  return Maxout(over)
end

@functor Maxout

function (mo::Maxout)(input::AbstractArray)
    mapreduce(f -> f(input), (acc, out) -> max.(acc, out), mo.over)
end

"""
    SkipConnection(layer, connection)

Create a skip connection which consists of a layer or `Chain` of consecutive
layers and a shortcut connection linking the block's input to the output
through a user-supplied 2-argument callable. The first argument to the callable
will be propagated through the given `layer` while the second is the unchanged,
"skipped" input.

The simplest "ResNet"-type connection is just `SkipConnection(layer, +)`.
Here is a more complicated example:
```jldoctest
julia> m = Conv((3,3), 4 => 7, pad=(1,1));

julia> x = ones(Float32, 5, 5, 4, 10);

julia> size(m(x)) == (5, 5, 7, 10)
true

julia> sm = SkipConnection(m, (mx, x) -> cat(mx, x, dims=3));

julia> size(sm(x)) == (5, 5, 11, 10)
true
```
"""
struct SkipConnection{T,F}
  layers::T
  connection::F  #user can pass arbitrary connections here, such as (a,b) -> a + b
end

@functor SkipConnection

function (skip::SkipConnection)(input)
  skip.connection(skip.layers(input), input)
end

function Base.show(io::IO, b::SkipConnection)
  print(io, "SkipConnection(", b.layers, ", ", b.connection, ")")
end

"""
    Bilinear(in1, in2, out)

Creates a Bilinear layer, which operates on two inputs at the same time.
It has parameters `W` and `b`, and its output given vectors `x`, `y` is of the form

    z[i] = σ.(x' * W[i,:,:] * y .+ b[i])

If `x` and `y` are matrices, then each column of the output `z = B(x, y)` is of this form,
given that `B` is a Bilinear layer of appropriate size.

If `y` is not given, it is taken to be equal to `x`, i.e. `B(x) == B(x, x)`
The two inputs may also be provided as a tuple, `B((x, y)) == B(x, y)`,
which is accepted as the input to a `Chain`.

```julia
# using Bilinear to generate interactions, on one input
x = randn(Float32, 11, 7)
B = Bilinear(11, 11, 3)
size(B(x)) == (3, 7)

# using Bilinear on two data streams at once, as a tuple
x = randn(Float32, 10, 9)
y = randn(Float32, 2, 9)
m = Chain(Bilinear(10, 2, 3), Dense(3, 1))
size(m((x, y))) == (1, 9)

# using Bilinear as the recombinator in a SkipConnection
x = randn(Float32, 10, 9)
sc = SkipConnection(Dense(10, 10), Bilinear(10, 10, 5))
size(sc(x)) == (5, 9)
```
"""
struct Bilinear{A,B,S}
  W::A
  b::B
  σ::S
end

@functor Bilinear

Bilinear(W, b) = Bilinear(W, b, identity)

function Bilinear(in1::Integer, in2::Integer, out::Integer, σ = identity;
  initW = glorot_uniform, initb = zeros)
  return Bilinear(initW(out, in1, in2), initb(out), σ)
end

function (a::Bilinear)(x::AbstractMatrix, y::AbstractMatrix)
  W, b, σ = a.W, a.b, a.σ

  d_z, d_x, d_y = size(W)
  d_x == size(x,1) && d_y == size(y,1) || throw(DimensionMismatch("number of rows in data must match W"))
  size(x,2) == size(y,2) || throw(DimensionMismatch("Data inputs must agree on number of columns, got $(size(x,2)) and $(size(y,2))"))

  # @einsum Wy[o,i,s] := W[o,i,j] * y[j,s]
  Wy = reshape(reshape(W, (:, d_y)) * y, (d_z, d_x, :))

  # @einsum Z[o,s] := Wy[o,i,s] * x[i,s]
  Wyx = batched_mul(Wy, reshape(x, (d_x, 1, :)))
  Z = reshape(Wyx, (d_z, :))

  # @einsum out[o,s] := σ(Z[o,i] + b[o])
  σ.(Z .+ b)
end

(a::Bilinear)(x::AbstractVecOrMat) = a(x, x)
(a::Bilinear)(x::AbstractVector, y::AbstractVector) = vec(a(reshape(x, :,1), reshape(y, :,1)))
(a::Bilinear)(x::NTuple{2, AbstractArray}) = a(x[1], x[2])

function Base.show(io::IO, l::Bilinear)
  print(io, "Bilinear(", size(l.W, 2), ", ", size(l.W, 3), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

"""
Parallel(connection, layers...)

Create a 'Parallel' layer that passes an input array to each path in
`layers`, reducing the output with `connection`.

Called with one input `x`, this is equivalent to `reduce(connection, [l(x) for l in layers])`.
If called with multiple inputs, they are `zip`ped with the layers, thus `Parallel(+, f, g)(x, y) = f(x) + g(y)`.

# Examples

```jldoctest
julia> model = Chain(Dense(3, 5),
                     Parallel(vcat, Dense(5, 4), Chain(Dense(5, 7), Dense(7, 4))),
                     Dense(8, 17));

julia> size(model(rand(3)))
(17,)

julia> model = Parallel(+, Dense(10, 2), Dense(5, 2))
Parallel(+, Dense(10, 2), Dense(5, 2))

julia> size(model(rand(10), rand(5)))
(2,)
```
"""
struct Parallel{F, T}
  connection::F
  layers::T
end

Parallel(connection, layers...) = Parallel(connection, layers)

@functor Parallel

(m::Parallel)(x::AbstractArray) = mapreduce(f -> f(x), m.connection, m.layers)
(m::Parallel)(xs::Vararg{<:AbstractArray}) = mapreduce((f, x) -> f(x), m.connection, m.layers, xs)
(m::Parallel)(xs::Tuple) = m(xs...)

Base.getindex(m::Parallel, i::Integer) = m.layers[i]
Base.getindex(m::Parallel, i::AbstractVector) = Parallel(m.connection, m.layers[i]...)

function Base.show(io::IO, m::Parallel)
  print(io, "Parallel(", m.connection, ", ")
  join(io, m.layers, ", ")
  print(io, ")")
end


text:
using Distributed
using Random
using Revise
# using ProgressMeter
# using ArgParse

const IN_SLURM = "SLURM_JOBID" in keys(ENV)
IN_SLURM && using ClusterManagers

ARGS_KWARGS = [
    Dict(:ARGS=>(0.5/8, 0.1, 1.0, 32, 8),   :KWARGS=>Dict([:n_episodes=>1000])),
    Dict(:ARGS=>(0.5/4, 0.1, 1.0, 32, 4),   :KWARGS=>Dict([:n_episodes=>1000])),
    Dict(:ARGS=>(0.5/2, 0.1, 1.0, 32, 2),   :KWARGS=>Dict([:n_episodes=>1000])),
    Dict(:ARGS=>(0.5/16, 0.1, 1.0, 32, 16), :KWARGS=>Dict([:n_episodes=>1000]))
]

function parallel_experiment_runs_args_episodic(experiment_file, args_kwargs_list, n_runs=10; exp_func_name=:run, num_workers=5, make_results=nothing)

    pids = nothing
    if IN_SLURM
        pids = addprocs(SlurmManager(parse(Int, ENV["SLURM_NTASKS"])))
        print("\n")
    else
        println(num_workers, " ", nworkers())
        if nworkers() == 1
            pids = addprocs(num_workers)
        elseif nworkers() < num_workers
            pids = addprocs((num_workers) - nworkers())
        end
    end

    println(nworkers(), " ", pids)

    @everywhere global exp_file=$experiment_file
    @everywhere include(exp_file)
    exp_func = getfield(Main, Symbol(exp_func_name))
    experiment(r, ak) = exp_func(ak[:ARGS]...; seed=1002349+r, ak[:KWARGS]...)

    res_dict = Dict{Int64, Any}()
    futures = Array{Array{Future,1}, 1}()

    for (args_kwargs_idx, args_kwargs) in enumerate(args_kwargs_list)

        in_futures = Array{Future,1}()

        for run in 1:n_runs
            append!(in_futures, [@spawn experiment(run, args_kwargs)])
        end
        if make_results == nothing
            results = zeros(n_runs, args_kwargs[:KWARGS][:n_episodes])
        else
            results = make_results(n_runs, args_kwargs)
        end

        append!(futures, [in_futures])

        res_dict[args_kwargs_idx] = Dict([
            "args"=>args_kwargs,
            "results"=>results
        ])
    end

    for (args_kwargs_idx, args_kwargs) in enumerate(args_kwargs_list)
        for run in 1:n_runs
            res_dict[args_kwargs_idx]["results"][run,:] .= fetch(futures[args_kwargs_idx][run])
        end
    end

    return res_dict
end

function parallel_experiment_runs(experiment_file, args_kwargs_list, id, n_runs=10)

end

function experiment(experiment_file, arg_file, id, r=1)

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
using CUDA
using PaddedViews
function k_elem_mul_3d!(A, B)
    tx = (blockIdx().x-1)*blockDim().x + threadIdx().x
    ty = (blockIdx().y-1)*blockDim().y + threadIdx().y
    tz = (blockIdx().z-1)*blockDim().z + threadIdx().z
    sx,sy,sz = size(A)
    if tx > sx || ty > sy || tz > sz
        return nothing
    end
    @inbounds a = A[tx,ty,tz]
    @inbounds b = B[tx,ty,tz]
    @inbounds A[tx,ty,tz] = a*b
    return nothing
end

function elem_mul_3d!(A,B)
    sx,sy,sz = size(A)
    tx = 16; ty = 8; tz = 8

    bx = 1 + sx÷tx
    by = 1 + sy÷ty
    bz = 1 + sz÷tz
    @cuda blocks=(bx,by,bz) threads=(tx,ty,tz) k_elem_mul_3d!(A,B)
    nothing
end

function fft_stencil_3d(data, template, t_steps=1)
    # template is a 3d square array
    #r = (size(template,1)-1)÷2

    r2 = size(template,1)
    r = 32
    while r < r2
        r += 32
    end
    r = r + 1
    dx,dy,dz = size(data)
    padded_data = PaddedView(0, data, (1:dx+r-1, 1:dy+r-1, 1:dz+r-1), (1:dx,1:dy,1:dz))
    padded_template = PaddedView(0, template, (1:dx+r-1, 1:dy+r-1, 1:dz+r-1), (1:r2,1:r2,1:r2))
    ## Initial Upload
    fd = CuArray{CUDA.cuFloatComplex}(padded_data)
    ft = CuArray{CUDA.cuFloatComplex}(padded_template)
    temp = similar(ft)
    CUFFT.fft!(fd)
    CUFFT.fft!(ft)
    #@time m = fd.*ft
    #@time elem_mul_3d!(fd,ft)
    ## Time Loop
    CUFFT.ifft!(fd)
    r2 = r2÷2 + 1
    rn = r2:(r2+dx-1)
    return real(Array(fd))[rn,rn,rn]
end

# function stride_bench(data, rsize)
#     template = rand(Float32, rsize, rsize, rsize)
#     @info "running for $rsize"
#     fft_stencil_3d(data, template, 1)
# end
#
# d = [1<<7 for i = 1:3]
# data = rand(Float32, d...)
# for i in [3 5 9 17 33]
#     global data
#     stride_bench(data, i)
# end
#
# CUDA.cuProfilerStart()
# for i in [3 5 9 17 33]
#     NVTX.@range "template$(i)" begin
#         global data
#         stride_bench(data, i)
#     end
# end

# function locFFTstencil(st_inst::StencilInstance, org_data, t_steps::Integer, multi=1)
#     @show z = size(org_data, 3)
#     template_single = matrix_from_expression(st_inst.stencil_sym, st_inst.max_radius)
#     template_multi = nothing
#     if multi > 1
#         template_multi = DSP.conv(template_single, template_single)
#         for i = 2:multi
#             template_multi = DSP.conv(template_multi, template_single)
#         end
#     end
include("../../module/CustomStencil.jl")
include("../../misc/misc.jl")
include("../../misc/cpu_stencils.jl")
coefs = [0.75;0.5;0.25]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1,2, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst1 = NewStencilInstance(st_def, m_step=false)

radius = 4
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)
gpu_out = ApplyStencil(st_inst1, data, 1)

t = matrix_from_expression(st_inst1.stencil_sym, st_inst1.max_radius)
a = fft_stencil_3d(data, t, 1)

b = a
b[.*(a .< 0.00001)] .= 0
c = gpu_out[10:20,10:20,16]
d = b[10:20,10:20,18]

using Plots

heatmap(a[:,:,1])
heatmap(gpu_out[:,:,4])

sum(a .- gpu_out)
sum(b .- gpu_out)
max((a.-gpu_out)...)
max((b.-gpu_out)...)

sum(b .< 0.00001)


text:
#####################################################
#Helpers

function getsparsity(input::Array{Float64, 1})
	length(findall(x -> (x == 0),input))/length(input)
end

function getsmallimg()
    global patternindex = rand(1:size(smallimgs)[2])
    smallimgs[:, patternindex]
end

function getsample()
    global patternindex = rand(1:n_samples)
    smallimgs[:, patternindex]
end

function getlabel(x)
    [labels[patternindex] == i for i in 0:9]
end

function gethiddenreps(sae)
		global patternindex = rand(1:size(smallimgs)[2])
		#sae.reps[n_ae][1][:, patternindex]
		sae.reps[n_ae-1][1][:, patternindex]
end

function gethiddenrepstest(sae)
		global patternindex = rand(1:size(smallimgstest)[2])
		#sae.reps[n_ae][1][:, patternindex]
		sae.reps[n_ae-1][2][:, patternindex]
end

function _generatehiddenreps!(sae,n_ae::Int64, smallimgs, smallimgstest, conf, func_forwardprop::Function;
		no_recurrence = false, Földiak_model = false, rand_patch_pos = false, colorimage = false)
	if typeof(conf) == configAE
		if n_ae == 1
			@showprogress for i in 1:size(smallimgs)[2]
				sae.aes[n_ae].x[1] = smallimgs[:, i]
				func_forwardprop(sae.aes[n_ae]; nonlinearity = conf.nonlinearity)
				sae.reps[n_ae][1][:,i] = deepcopy(sae.aes[n_ae].x[Int(mean([1,sae.aes[n_ae].nl+1]))])
			end
			@showprogress for i in 1:size(smallimgstest)[2]
				sae.aes[n_ae].x[1] = smallimgstest[:, i]
				func_forwardprop(sae.aes[n_ae]; nonlinearity = conf.nonlinearity)
				sae.reps[n_ae][2][:,i] = deepcopy(sae.aes[n_ae].x[Int(mean([1,sae.aes[n_ae].nl+1]))])
			end
		else
			@showprogress for i in 1:size(smallimgs)[2]
				sae.aes[n_ae].x[1] = deepcopy(sae.reps[n_ae-1][1][:,i])
				func_forwardprop(sae.aes[n_ae]; nonlinearity = conf.nonlinearity)
				sae.reps[n_ae][1][:,i] = deepcopy(sae.aes[n_ae].x[Int(mean([1,sae.aes[n_ae].nl+1]))])
			end
			@showprogress for i in 1:size(smallimgstest)[2]
				sae.aes[n_ae].x[1] = deepcopy(sae.reps[n_ae-1][2][:,i])
				func_forwardprop(sae.aes[n_ae]; nonlinearity = conf.nonlinearity)
				sae.reps[n_ae][2][:,i] = deepcopy(sae.aes[n_ae].x[Int(mean([1,sae.aes[n_ae].nl+1]))])
			end
		end
	elseif typeof(conf) == configAE_sparse
		if n_ae == 1
			for i in 1:size(smallimgs)[2]
				sae.aes[n_ae].x[1] = smallimgs[:, i]
				func_forwardprop(sae.aes[n_ae], conf.lambda; nonlinearity = conf.nonlinearity,
					no_recurrence = no_recurrence, Földiak_model = Földiak_model,
					rand_patch_pos = rand_patch_pos, colorimage = colorimage)
				sae.reps[n_ae][1][:,i] = deepcopy(sae.aes[n_ae].x[end])
			end
			for i in 1:size(smallimgstest)[2]
				sae.aes[n_ae].x[1] = smallimgstest[:, i]
				func_forwardprop(sae.aes[n_ae], conf.lambda; nonlinearity = conf.nonlinearity,
					no_recurrence = no_recurrence, Földiak_model = Földiak_model,
					rand_patch_pos = rand_patch_pos, colorimage = colorimage)
				sae.reps[n_ae][2][:,i] = deepcopy(sae.aes[n_ae].x[end])
			end
		else
			for i in 1:size(smallimgs)[2]
				sae.aes[n_ae].x[1] = deepcopy(sae.reps[n_ae-1][1][:,i])
				func_forwardprop(sae.aes[n_ae], conf.lambda; nonlinearity = conf.nonlinearity,
					no_recurrence = no_recurrence, Földiak_model = Földiak_model,
					rand_patch_pos = rand_patch_pos, colorimage = colorimage)
				sae.reps[n_ae][1][:,i] = deepcopy(sae.aes[n_ae].x[end])
			end
			for i in 1:size(smallimgstest)[2]
				sae.aes[n_ae].x[1] = deepcopy(sae.reps[n_ae-1][2][:,i])
				func_forwardprop(sae.aes[n_ae], conf.lambda; nonlinearity = conf.nonlinearity,
					no_recurrence = no_recurrence, Földiak_model = Földiak_model,
					rand_patch_pos = rand_patch_pos, colorimage = colorimage)
				sae.reps[n_ae][2][:,i] = deepcopy(sae.aes[n_ae].x[end])
			end
		end
	end
end

#for sparse coding (only 2-layer AEs!)
function generatehiddenreps!(sae::SAE_sparse,n_ae::Int64, smallimgs, smallimgstest, conf::configAE_sparse)
	print(string("compute sparse hidden representations (for train and test sets) of AE nr.",n_ae,"...\n"))
	_generatehiddenreps!(sae,n_ae, smallimgs, smallimgstest, conf, forwardprop!)
end

#for sparse coding à la Brito (only 2-layer AEs!)
function generatehiddenreps_Brito!(sae::SAE_sparse,n_ae::Int64, smallimgs, smallimgstest, conf::configAE_sparse;
		no_recurrence = false, Földiak_model = false, rand_patch_pos = false, colorimage = false)
	print(string("compute sparse hidden representations (for train and test sets) of AE nr.",n_ae,"...\n"))
	_generatehiddenreps!(sae,n_ae, smallimgs, smallimgstest, conf, forwardprop_Brito_full!;
		no_recurrence = no_recurrence, Földiak_model = Földiak_model, rand_patch_pos = rand_patch_pos, colorimage = colorimage) #forwardprop_Brito!
end
function generatehiddenreps_Brito!(sae::SAE_sparse_patchy,n_ae::Int64, smallimgs, smallimgstest, conf::configAE_sparse;
		no_recurrence = false, Földiak_model = false, rand_patch_pos = false, colorimage = false)
	print(string("compute sparse hidden representations (for train and test sets) of AE nr.",n_ae,"...\n"))
	_generatehiddenreps!(sae,n_ae, smallimgs, smallimgstest, conf, forwardprop_Brito_full!;
		no_recurrence = no_recurrence, Földiak_model = Földiak_model, rand_patch_pos = rand_patch_pos, colorimage = colorimage) #forwardprop_Brito!
end

function geterrors(net, imgs, labels; nonlinearity = nonlinearity)
	error = 0
	noftest = size(imgs)[2]
	for i in 1:noftest
		net.x[1] = imgs[:, i]
		forwardprop!(net, nonlinearity = nonlinearity)
		error += findmax(net.x[end])[2] != labels[i] + 1
	end
	error/noftest
end
function printerrors(net, smallimgs, labels, smallimgstest, labelstest, conf)
	error_train = geterrors(net, smallimgs, labels, nonlinearity = conf.nonlinearity)
	error_test = geterrors(net, smallimgstest, labelstest, nonlinearity = conf.nonlinearity)
	print(string(100 * error_train, " % on training set\n"))
	print(string(100 * error_test, " % on test set\n"))

	return error_train, error_test
end

function _rescale_weights!(net, i, input_dim, patch_size, mode)
	if mode == "randn"
		net.w[1][i,:] .*= 3 * input_dim / (patch_size) #rescale weights for same input strength as full connectivity
	elseif mode == "spherical"
		net.w[1][i,:] .*= 1. / norm(net.w[1][i,:]) # normalize random patches: since sampled from randn they lie on the surface of a hypersphere now! norm grows as sqrt(eff. in-dim)
		net.w[1][i,:] .*= 3.#sqrt(patch_size) # heuristic rescaling for best performance
	end
end
function set_connectivity!(net,patch_size; mode = "randn", colorimage = false)
	if colorimage
		input_size = Int(length(net.x[1])/3)
	else
		input_size = length(net.x[1])
	end
	input_dim = Int(sqrt(input_size))
	nr_hidden_neurons = length(net.x[2])
	mask = zeros(input_dim,input_dim) #initialize mask
	mask[1:patch_size,1:patch_size] = 1.
	if colorimage
		for i in 1:nr_hidden_neurons
			shifts = rand(0:input_dim-patch_size,2)#draw random shifts of mask
			for j in 1:3
				net.w[1][i,(j-1)*input_size+1:j*input_size] = reshape(reshape(net.w[1][i,(j-1)*input_size+1:j*input_size],input_dim,input_dim).*circshift(mask,shifts),input_size)
			end
			_rescale_weights!(net, i, input_dim, patch_size, mode)
		end
	else
		for i in 1:nr_hidden_neurons
			shifts = rand(0:input_dim-patch_size,2)#draw random shifts of mask
			net.w[1][i,:] = reshape(reshape(net.w[1][i,:],input_dim,input_dim).*circshift(mask,shifts),input_size)
			_rescale_weights!(net, i, input_dim, patch_size, mode)
		end
	end
	w_connect = clamp(ceil(abs(net.w[1])),0,1)
	return w_connect
end


text: include("openml.jl")


text:
module MiscUtilities

#include("baseimports.jl")

include("inspections.jl")
include("randsubset.jl")
include("methodlists.jl")
include("objectof.jl")
include("summation.jl")
include("ast.jl")
include("sparsegen.jl")
include("threads.jl")
include("exfiltrate.jl")
include("copyelementtypes.jl")
include("ilog2.jl")

end # module


text:
############
# Starring #
############

function stargazers(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/stargazers"; options...)
    return map(Owner, results), page_data
end

function starred(user; options...)
    results, page_data = gh_get_paged_json("/users/$(name(user))/starred"; options...)
    return map(Repo, results), page_data
end

star(repo; options...) = gh_put("/user/starred/$(name(repo))"; options...)

unstar(repo; options...) = gh_delete("/user/starred/$(name(repo))"; options...)

############
# Watching #
############

function watchers(repo; options...)
    results, page_data = gh_get_paged_json("/repos/$(name(repo))/subscribers"; options...)
    return map(Owner, results), page_data
end

function watched(owner; options...)
    results, page_data = gh_get_paged_json("/users/$(name(owner))/subscriptions"; options...)
    return map(Repo, results), page_data
end

watch(repo; options...) = gh_put("/repos/$(name(repo))/subscription"; options...)

unwatch(repo; options...) = gh_delete("/repos/$(name(repo))/subscription"; options...)


text:
using Pkg
using LinearAlgebra
using DataFrames
import Statistics.mean
import Statistics.std
using Printf
using CSV

# load functions
include(pwd()*"/src/SDEMEM OU neuron data/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU neuron data/mcmc.jl")

seed = 100

M_subjects = 1

y, prior_parameters_η, prior_parameters_σ_ϵ = set_up(M_subjects, seed)

startval_ϕ = ones(M_subjects,3)

for j = 1:3

    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_ϕ[:,j,1] = μ_0_j*ones(M_subjects)

end

startval_log_σ_ϵ = log(0.5)

# set time step
dt = diff(y[1].mSec)[1]

# posterior mean ests from kalman filter
log_σ_ϵ = -3.60203
ϕ  = [-3.26986 -0.888105 -0.780767]

nbr_particles = 100*ones(Int64, M_subjects)
ρ = 0.0

function run_bridge_pf(particels)

    nbr_particles = particels*ones(Int64, M_subjects)
    ρ = 0.0

    # set number of obs for each data set
    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y[i].mV); end

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    run_time = @elapsed ll1 =  bridge_pf(y, exp(log_σ_ϵ), ϕ, dt, u_prop_old, u_resample_old, nbr_particles, false)

    return ll1, run_time

end

################################################################################
# run at one specific # particles
################################################################################

#=

nbr_particels = 1000

ll1, run_time = @time run_bridge_pf(nbr_particels)

@printf "Runtime Bridge filter: %.4f\n" run_time

nbr_pf_eval = 50
ll_vec = zeros(nbr_pf_eval)
run_time_vec = zeros(nbr_pf_eval)
ll_all = zeros(M_subjects, nbr_pf_eval)

for i = 1:nbr_pf_eval

    println(i)

    ll1, run_time = run_bridge_pf(nbr_particels)
    ll_all[:,i] = ll1
    ll_vec[i] = sum(ll1)
    run_time_vec[i] = run_time

end

@printf "----------------\n"
@printf "Bridge filter\n"
@printf "Nbr particles: %.2f\n" nbr_particels
@printf "Runtime: %.4f\n" mean(run_time_vec)
@printf "Loglik (mean): %.4f\n" mean(ll_vec)
@printf "Loglik (var): %.4f\n" var(ll_vec)
@printf "Loglik (std): %.4f\n" std(ll_vec)
@printf "Loglik (var) ind ll est:\n"
println(var(ll_all, dims = 2))

=#

################################################################################
# run over multiple number of # particles
################################################################################

function loglik_est(nbr_particles_pf, nbr_pf_eval)

    ll = zeros(nbr_pf_eval)

    for i in 1:nbr_pf_eval
        ll_part, ~ = run_bridge_pf(nbr_particles_pf)
        ll[i] = sum(ll_part)
    end

    return ll

end

nbr_pf_eval = 100
nbr_particles_pf = [1, 2, 5, 10, 20, 50, 100]

loglik_matrix = zeros(nbr_pf_eval, length(nbr_particles_pf))
run_times = zeros(length(nbr_particles_pf))

i = 1
run_times[i]  = @elapsed loglik_matrix[:,i] = loglik_est(nbr_particles_pf[i], nbr_pf_eval)

Random.seed!(1)

for i in 1:length(nbr_particles_pf)
    @printf "Current nbr of particels: %.2f\n" nbr_particles_pf[i]
    run_times[i]  = @elapsed loglik_matrix[:,i]  = loglik_est(nbr_particles_pf[i], nbr_pf_eval)
end

using PyPlot

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, mean(loglik_matrix, dims = 1)[:], "--*")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Avg. loglik est. val.")
#PyPlot.savefig("figures/neuronal_avg_loglik_val_vs_nbr_particels_pf.png")

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, var(loglik_matrix, dims = 1)[:], "--*")
PyPlot.plot(nbr_particles_pf, 2*ones(length(nbr_particles_pf),1), "k")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Var of loglik est.")
#PyPlot.savefig("figures/neuronal_loglik_var_vs_nbr_particles_pf.png")

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, run_times/100, "--*")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Run time (sec)")
#PyPlot.savefig("figures/neuronal_run_time_vs_nbr_particls_pf.png")

max_stars_repo_path: src/SDEMEM OU neuron data/tune_nbr_particles_bridge_pf_one_subject.jl


text:
function step_change_first!(integrator)
  integrator.p[2].dyn_injections[1].P_ref += integrator.p[1][1]
end

function step_change_second!(integrator)
    integrator.p[2].dyn_injections[2].P_ref += integrator.p[1][1]
end


text:
module DietRecommender

using ProgressMeter
using Serialization
using DataFrames
using CSV
using Query
using JuMP
using Ipopt
using Juniper
using StatsBase: sample
include("get_data.jl")
include("optimize.jl")
include("return_diet.jl")

end # module


text:
rf001,cf001=simpleprover("data/dataf/extf001.cnf",5,5)
rf002,cf002=simpleprover("data/dataf/extf002.cnf",5,5)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
module L1Homotopy

export L1_homotopy

"""
    L1_homotopy(A, y, λₑ)

Solve the optimization problem

  minimize  J(x) = ‖A*x - y‖₂² + λₑ‖x‖₁
     x

  subject to  xᵢ > 0 ∀ i

using the homotopy algorithm.  Here we largely follow the notation from the
paper "Fast solution to L1-norm minimization problems when the solution may be
sparse" by David Donoho and Yaakov Tsaig.  The formulation above is labelled
(D_λ) in the paper, but there are several other equivalent formulations,
including two constrained formulations labelled (P_1) and (L_q) in the paper.
L_q is popular for statistical model selection and is called Lasso in the
statistical literature.

# Devdocs

The algorithm works by following the optimal solution from λ=∞ down to λ=λₑ.
It turns out this path is a series of linear segments x0 -> x1 -> ... -> xf
through the space of x's, with `x0=zeros(N)`.
"""
function L1_homotopy(A, y, λₑ)
    # TODO: ATA may be precached, or computed on the fly.
    ATA = Symmetric(A'*A)
    N = size(A,2)
    # Residual correlations of signal with system model,
    #    c(x) ≡ A'*(y - A*x)
    # initialized where x=0
    c = A'*y
    # Initial active set
    λ = maximum(c)
    Iₐ = find(c .== λ)
    ET = eltype(A)
    T = typeof((zero(ET)*one(ET) + zero(ET)) / one(ET))
    xₐ = zeros(Iₐ, T)
    iter = 0
    while λ > λₑ
        # Cache inactive indices
        Iᵢ = setdiff(1:length(y), Iₐ)
        # Update direction, restricted to the active set
        dₐ = Symmetric(ATA[Iₐ,Iₐ]) \ sign.(c[Iₐ])
        # Compute update magnitude, γ. There's two ways this can occur
        #
        # 1) An inactive set correlation residual increases to the current λ,
        #    causing an inactive element to enter the active set
        γ⁺ = Inf
        new_active = 0
        if !isempty(Iᵢ)
            dcᵢ_dγ = ATA[Iᵢ,Iₐ]*dₐ
            γs = (λ .- c[Iᵢ]) ./ (one(T) .- dcᵢ_dγ)
            for (i, γ) ∈ zip(Iᵢ, γs)
                if γ > 0 && γ <= γ⁺
                    # TODO: Multiple additions to active set
                    γ⁺ = γ
                    new_active = i
                end
            end
        end
        # 2) An active set element changes sign (ie, decreases to zero in the
        #    positively constrained formulation)
        # Active element to leave the active set
        γ⁻ = Inf
        old_active_ind_ind = 0
        if !isempty(Iₐ)
            γs = -xₐ ./ dₐ
            for (i, γ) ∈ enumerate(γs)
                if γ > 0 && γ <= γ⁻
                    # TODO: Mutliple removals from active set
                    γ⁻ = γ
                    old_active_ind_ind = i
                end
            end
        end
        γ = min(γ⁺, γ⁻)
        # Incremental update to  c ≡ A'*(y - A*x)
        c .-= γ .* ATA[:,Iₐ]*dₐ
        xₐ .+= γ .* dₐ
        λ = maximum(c)
        if γ⁺ < γ⁻
            push!(Iₐ, new_active)
            push!(xₐ, zero(T))
        else
            splice!(Iₐ, old_active_ind_ind)
            splice!(xₐ, old_active_ind_ind)
        end
        iter += 1
    end
    x = zeros(N)
    x[Iₐ] = xₐ
    x
end

end


text:
################################################################################
######          Particle filter and help functions for pf                  #####
################################################################################

doc"""
    pf(y::Array{Float64}, theta::Array{Float64},theta_known::Float64,N::Int64,
plotflag::Bool=false, return_weigths_and_particles::Bool=false)

pf runs the bootstrap particel filter for the Ricker model.
"""
function pf(y::Array{Float64}, theta::Array{Float64},theta_known::Float64,N::Int64,
plotflag::Bool=false, return_weigths_and_particles::Bool=false)

  # set parameter values
  r = exp(theta[1])
  phi = exp(theta[2])
  sigma = exp(theta[3])

  # set startvalue for loglik
  loglik = 0.

  # set length pf data
  T = length(y)

  # pre-allocate matriceis
  x = zeros(N,T) # particels
  w = zeros(N,T) # weigts
  x_anc = zeros(N,T+1) # ancestral particels

  # set start values
  xint = rand(Uniform(1,30),N,1)
  #xint = rand(Uniform(1,200),N,1)

  x_anc[:,1] = xint # set anc particels for t = 1

  # set gaussian noise
  e = rand(Normal(0,sigma), N,T)

  for t = 1:T

  if t == 1 # first iteration

    # propagate particels
    x[:,1] = r*xint.*exp(-xint .+ e[:,1]);

    # calc weigths and update loglik
    (w[:,t], loglik) = calc_weigths(y[t],x[:,t],phi,loglik,N)

  else

    # resample particels
    ind = stratresample(w[:,t-1], N)
    x_resample = x[ind,t-1]

    x_anc[:,t+1] = x_resample # store ancestral particels

    # propagate particels
    x[:,t] = r*x_resample.*exp(-x_resample .+ e[:,t])

    # calc weigths and update loglik
    (w[:,t], loglik) = calc_weigths(y[t],x[:,t],phi,loglik,N)
  end

  end

  if plotflag # plot ESS at last iteration
  @printf "ESS: %.4f\n" 1/sum(w[:,end].^2)
  end

  if return_weigths_and_particles
  # return loglik, weigths and particels
  loglik, w, x
  else
  # return loglik
  return loglik
  end
end

# help functions for particle filter

doc"""
    calc_weigths(y::Array{Float64},x::Array{Float64},phi::Float64,loglik::Float64,N::Int64)

Calculates the weigths in the particel filter and the estiamtes the loglikelihood value.
"""
function calc_weigths(y::Float64,x::Array{Float64},phi::Float64,loglik::Float64,N::Int64)
  logweigths = y*log(x.*phi)  .- x*phi # compute logweigths
  constant = maximum(logweigths) # find largets wegith
  weigths = exp(logweigths - constant) # subtract largets weigth and compute weigths
  loglik =  loglik + constant + log(sum(weigths)) - log(N) # update loglik
  return weigths/sum(weigths), loglik
end

doc"""
    stratresample(p , N)

Stratified resampling.

Sample N times with repetitions from a distribution on [1:length(p)] with probabilities p. See [link](http://www.cornebise.com/julien/publis/isba2012-slides.pdf).
"""
function stratresample(p , N)

  p = p/sum(p)  # normalize, just in case...

  cdf_approx = cumsum(p)
  cdf_approx[end] = 1
  #I = zeros(N,1)
  indx = zeros(Int64, N)
  U = rand(N,1)
  U = U/N + (0:(N - 1))/N
  index = 1
  for k = 1:N
    while (U[k] > cdf_approx[index])
      index = index + 1
    end
    indx[k] = index
  end

  return indx

end

max_stars_repo_path: Ricker model/pf.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: Ricker model/pf.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: Ricker model/pf.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
using DashBootstrapComponents

badge = dbc_button(
    [
        "Notifications",
        dbc_badge(
            "99+",
            color = "danger",
            pill = true,
            text_color = "white",
            className = "position-absolute top-0 start-100 translate-middle",
        ),
    ],
    color = "primary",
    className = "position-relative",
)

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ContainerState <: SwaggerModel
    running::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ContainerStateRunning } # spec name: running
    terminated::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ContainerStateTerminated } # spec name: terminated
    waiting::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ContainerStateWaiting } # spec name: waiting

    function IoK8sApiCoreV1ContainerState(;running=nothing, terminated=nothing, waiting=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ContainerState, Symbol("running"), running)
        setfield!(o, Symbol("running"), running)
        validate_property(IoK8sApiCoreV1ContainerState, Symbol("terminated"), terminated)
        setfield!(o, Symbol("terminated"), terminated)
        validate_property(IoK8sApiCoreV1ContainerState, Symbol("waiting"), waiting)
        setfield!(o, Symbol("waiting"), waiting)
        o
    end
end # type IoK8sApiCoreV1ContainerState

const _property_map_IoK8sApiCoreV1ContainerState = Dict{Symbol,Symbol}(Symbol("running")=>Symbol("running"), Symbol("terminated")=>Symbol("terminated"), Symbol("waiting")=>Symbol("waiting"))
const _property_types_IoK8sApiCoreV1ContainerState = Dict{Symbol,String}(Symbol("running")=>"IoK8sApiCoreV1ContainerStateRunning", Symbol("terminated")=>"IoK8sApiCoreV1ContainerStateTerminated", Symbol("waiting")=>"IoK8sApiCoreV1ContainerStateWaiting")
Base.propertynames(::Type{ IoK8sApiCoreV1ContainerState }) = collect(keys(_property_map_IoK8sApiCoreV1ContainerState))
Swagger.property_type(::Type{ IoK8sApiCoreV1ContainerState }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ContainerState[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ContainerState }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ContainerState[property_name]

function check_required(o::IoK8sApiCoreV1ContainerState)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ContainerState }, name::Symbol, val)
end


text:
include("MQuery.jl")
@testset "MQuery" begin
using Base.Enums

@enum TypeChecking Dynamic Static
df = Dict(
        Symbol("Type checking") => [
            Dynamic, Static, Static, Dynamic, Static, Dynamic, Dynamic, Static
        ],
        :name => [
            "Julia", "C#", "F#", "Ruby", "Java", "JavaScript", "Python", "Haskell"
        ],
        :year => [
            2012, 2000, 2005, 1995, 1995, 1995, 1990, 1990
        ]
)

res = df |>
@where !startswith(_.name, "Java"),
@groupby _."Type checking" => TC, endswith(_.name, "#") => is_sharp,
@having TC === Dynamic || is_sharp,
@select join(_.name, " and ") => result, _.TC => TC

@info res
@test res[:result][map(==(Dynamic), res[:TC]) ] == ["Julia and Ruby and Python"]
@test res[:result][map(==(Static), res[:TC])  ] == ["C# and F#"]

res = df |>
@select _.(!startswith("Type"))
@test Set(keys(res)) == Set([:name, :year])

end


text:
## Exercise 2-1
## To check your understanding, type the following statements in the Julia REPL and see what they do:
## in repl
# 5  # >> 5
# x = 5  # >> 5
# x + 1  # >> 6

# Now put the same statements in a script and run it. What is the output? Modify the script by transforming each expression into a print statement and then run it again.
println("Ans: ")
5
x = 5
x + 1
# no output, lets try using print.

println(5)
# println(x = 5)  # error, LoadError: MethodError
x = 5
println(x)
println(x + 1)

println("End.")


text:
module HelmholtzOperators2DSymmetric

using ..VectorFields
using ..Simulations
using ..Dispersions
using SparseArrays

import ..Symmetric, ..Unsymmetric
import LinearAlgebra: I
import ..Helmholtz
import ..helmholtz_susceptability!

function Helmholtz(
            sim::TSIM;
            m::Int = 1) where TSIM<:Simulation{2,Symmetric}

    ka, kb, kc = 0, 0, 0
    laplacian = sim.laplacian.l0

    Σs = _compute_Σs(sim)
    f = sim.self_energy.f

    D² = laplacian - Σs[1]*f[1](0,ka) - Σs[2]*f[2](0,ka) - Σs[3]*f[3](0,ka) - Σs[4]*f[4](0,ka) # Left + Right + Bottom + Top

    αε = sim.laplacian.complex_scaling*spdiagm(0=>sim.ε)
    αεpFχ = copy(αε)
    Fχ = spdiagm(0=>zeros(ComplexF64,size(αε,1)))

    Fχs = Vector{SparseMatrixCSC{ComplexF64,Int}}(undef,m)
    dFχdψr = zeros(ComplexF64,length(sim),m,m)
    dFχdψi = copy(dFχdψr)
    dFχdω = copy(dFχdψi)
    dFχdϕ = copy(dFχdω)
    for μ ∈ eachindex(Fχs) Fχs[μ] = spdiagm(0=>zeros(ComplexF64,size(αε,1))) end
    Helmholtz{2,Symmetric,typeof(Σs),TSIM}(D²+I,D²,αεpFχ,laplacian,Σs,αε,Fχ,Fχs,dFχdψr,dFχdψi,dFχdω,dFχdϕ,sim,ka,kb,kc)
end

@inline function (h::Helmholtz{2,Symmetric})(ω::Number,args...)
    f = h.sim.self_energy.f
    ka = 0
    fL = f[1](ω,ka)
    fR = f[2](ω,ka)
    fB = f[3](ω,ka)
    fT = f[4](ω,ka)

    rows = rowvals(h.D²)
    vals = nonzeros(h.D²)
    _, n = size(h.D²)
    @inbounds for i ∈ 1:n
        col = i
        @fastmath @inbounds @simd for j ∈ nzrange(h.D², i)
            row = rows[j]
            vals[j] = h.laplacian[row,col] - h.Σs[1][row,col]*fL - h.Σs[2][row,col]*fR - h.Σs[3][row,col]*fB - h.Σs[4][row,col]*fT
        end
    end

    helmholtz_susceptability!(h,ω,args...)
    ω² = ω^2

    rows = rowvals(h.A)
    vals = nonzeros(h.A)
    @inbounds for col ∈ 1:n
        @fastmath @inbounds @simd for j ∈ nzrange(h.A, col)
            row = rows[j]
            vals[j] = h.D²[row,col] + ω²*h.αεpFχ[row,col]
        end
    end
    return h.A
end

################################################################################
function _compute_Σs(sim::Simulation{2,Symmetric})
    Σ0L = sim.self_energy.Σ0[1]
    Σ0R = sim.self_energy.Σ0[2]
    Σ0B = sim.self_energy.Σ0[3]
    Σ0T = sim.self_energy.Σ0[4]
    return (Σ0L, Σ0R, Σ0B, Σ0T)
end

end

using .HelmholtzOperators2DSymmetric


text:
# Requires Julia 1.3 (in beta)
# Run with JULIA_NUM_THREADS=36 benchmarks/partr_fib.jl

import Base.Threads.@spawn

function fib(n::Int)
    if n < 2
        return n
    end
    x = @spawn fib(n - 1)
    y = fib(n - 2)
    return fetch(x) + y
end

function main()
  n = parse(Int, ARGS[1])
  println(fib(n))
end

main()

max_stars_repo_name: Vindaar/weave

max_issues_repo_name: Vindaar/weave

max_forks_repo_name: Vindaar/weave


text:
#Código de Charly

using FileIO
using Base.Threads
using BenchmarkTools

const order = 1
const L = 5
const N = 80000
#const N = 20000
const J = 50
const κ = 2π/L # for Fourier Transform
const dx = L/J

include("aux_functions.jl")

println("threads = $(nthreads())")

function get_current_threads_soa(r, v, p)
  L, N, J, κ, dx, order = p

  TS = zeros(J, nthreads())
  @threads for i in 1:N
    @inbounds j, y = get_index_and_y(r[i], J, L)
    for l in (-order):-j
      @inbounds TS[J + j + l, threadid()] += W(order, -y + l) * v[i] / dx;
    end
    for l in max(-order,-j+1):min(order,J-j)
      @inbounds TS[j + l, threadid()] += W(order, -y + l) * v[i] / dx;
    end
    for l in J-j+1:order
      @inbounds TS[j - J + l, threadid()] += W(order, -y + l) * v[i] / dx;
    end
    # for l in (-order):order
    #   @inbounds TS[mod1(j + l, J), threadid()] += W(order, -y + l) * v[i] / dx;
    # end
  end

  S = zeros(J)
  @threads for i in 1:J
    for t in 1:nthreads()
      @inbounds S[i] += TS[i, t]
    end
  end
  S
end

function get_current_threads(u, p)
  L, N, J, κ, dx, order = p

  TS = zeros(J, nthreads())
  @threads for i in 1:N
    @inbounds j, y = get_index_and_y(u[i], J, L)
    for l in (-order):-j
      @inbounds TS[J + j + l, threadid()] += W(order, -y + l) * u[N+i] / dx;
    end
    for l in max(-order,-j+1):min(order,J-j)
      @inbounds TS[j + l, threadid()] += W(order, -y + l) * u[N+i] / dx;
    end
    for l in J-j+1:order
      @inbounds TS[j - J + l, threadid()] += W(order, -y + l) * u[N+i] / dx;
    end
    # for l in (-order):order
    #   @inbounds TS[mod1(j + l, J), threadid()] += W(order, -y + l) * v[i] / dx;
    # end
  end

  S = zeros(J)
  @threads for i in 1:J
    for t in 1:nthreads()
      @inbounds S[i] += TS[i, t]
    end
  end
  S
end

function get_current_threads_v2(u, p)
  L, N, J, κ, dx, order = p

  TS = zeros(J, nthreads())
  j = fill(Int64(1),nthreads())
  y = zeros(nthreads())
  @threads for i in 1:N
    @inbounds j[threadid()], y[threadid()] = get_index_and_y(u[i], J, L)
    for l in (-order):-j[threadid()]
      @inbounds TS[J + j[threadid()] + l, threadid()] += W(order, -y[threadid()] + l) * u[N+i] / dx;
    end
    for l in max(-order,-j[threadid()]+1):min(order,J-j[threadid()])
      @inbounds TS[j[threadid()] + l, threadid()] += W(order, -y[threadid()] + l) * u[N+i] / dx;
    end
    for l in J-j[threadid()]+1:order
      @inbounds TS[j[threadid()] - J + l, threadid()] += W(order, -y[threadid()] + l) * u[N+i] / dx;
    end
    # for l in (-order):order
    #   @inbounds TS[mod1(j + l, J), threadid()] += W(order, -y + l) * v[i] / dx;
    # end
  end

  S = zeros(J)
  @threads for i in 1:J
    for t in 1:nthreads()
      @inbounds S[i] += TS[i, t]
    end
  end
  S
end

function get_current_threads_aos(rv, p)
  L, N, J, κ, dx, order = p

  TS = zeros(J, nthreads())
  @threads for i in 1:N
    @inbounds j, y = get_index_and_y(rv[1,i], J, L)
    for l in (-order):-j
      @inbounds TS[J + j + l, threadid()] += W(order, -y + l) * rv[2,i] / dx;
    end
    for l in max(-order,-j+1):min(order,J-j)
      @inbounds TS[j + l, threadid()] += W(order, -y + l) * rv[2,i] / dx;
    end
    for l in J-j+1:order
      @inbounds TS[j - J + l, threadid()] += W(order, -y + l) * rv[2,i] / dx;
    end
    # for l in (-order):order
    #   @inbounds S[mod1(j + l, J), threadid()] += W(order, -y + l) * rv[2,i] / dx;
    # end
  end

  S = zeros(J)
  @threads for i in 1:J
    for t in 1:nthreads()
      @inbounds S[i] += TS[i, t]
    end
  end
  S
end

"""
reorder particle vector so that position and velocity are contiguous
"""
function split_particles_soa!(u, r, v)
  r[:] = u[1:N]
  v[:] = u[N+1:2N]
end

function split_particles_aos!(u, rv)
  for i in 1:N
    rv[1, i] = u[i]
    rv[2, i] = u[N+i]
  end
end

function load_soa()
  pars = zeros(2N)
  #pars = dones((N,), workers(), [nworkers()] );
  par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta001.jld2", "par_dis");

  r = zeros(N)
  v = zeros(N)
  split_particles_soa!(par_dis, r, v)
  (r, v)
end

function load_aos()
  pars = zeros(2N)
  #pars = dones((N,), workers(), [nworkers()] );
  par_dis = load("Initial_Distributions/par_dis_L5_N8_5_theta001.jld2", "par_dis");

  rv = zeros(2,N)
  split_particles_aos!(par_dis, rv)
  rv
end

x = [dx*(i-1) for i in 1:J]
p = (L, N, J, κ, dx, order)

#rv = load_aos()
#@benchmark get_current_threads_aos(rv, p)

(r,v) = load_soa()
#@benchmark get_current_threads_soa(r, v, p)

u = [r;v]
#@benchmark get_current_threads(u, p)

if nthreads() == 1
  S = zeros(J)
  get_current!(u,S,p)
else
  S = zeros(J)
  S = get_current_threads_v2(u,p)
end

save("current_th$(nthreads()).jld2", Dict("s" => S))


text:
instructions = []
open("day08/day08.in") do io
    global instructions
    for line in readlines(io)
        line = split(line, " ")
        push!(instructions, [line[1], parse(Int64, line[2])])
    end
end

executed = zeros(Bool, length(instructions))
exec_pos = 1
accumulator = 0
while !executed[exec_pos]
    global accumulator, exec_pos
    executed[exec_pos] = true
    instruction = instructions[exec_pos]
    exec_pos += 1
    if instruction[1] == "acc"
        accumulator += instruction[2]
    elseif instruction[1] == "jmp"
        exec_pos += instruction[2]-1
    end
end
println(accumulator)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
```
Central difference

```
function central_diff(
    y::AbstractArray{<:AbstractFloat,1},
    x::AbstractArray{<:AbstractFloat,1},
)

    dy = zeros(eltype(y), axes(y))

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    dy[i0] = (y[i0+1] - y[i0]) / (x[i0+1] - x[i0])
    dy[i1] = (y[i1] - y[i1-1]) / (x[i1] - x[i1-1])
    for i = i0+1:i1-1
        dy[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1])
    end

    return dy

end

function central_diff(y::AbstractArray{<:AbstractFloat,1}, dx::Real)
    x = ones(eltype(y), axes(y)) .* dx
    dy = central_diff(y, x)

    return dy
end

function central_diff!(
    dy::AbstractArray{<:AbstractFloat,1},
    y::AbstractArray{<:AbstractFloat,1},
    x::AbstractArray{<:AbstractFloat,1},
)

    @assert axes(dy) == axes(y) == axes(x)

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    dy[i0] = (y[i0+1] - y[i0]) / (x[i0+1] - x[i0])
    dy[i1] = (y[i1] - y[i1-1]) / (x[i1] - x[i1-1])
    for i = i0+1:i1-1
        dy[i] = (y[i+1] - y[i-1]) / (x[i+1] - x[i-1])
    end

end

function central_diff!(
    dy::AbstractArray{<:AbstractFloat,1},
    y::AbstractArray{<:AbstractFloat,1},
    dx::Real,
)
    x = ones(eltype(y), axes(y)) .* dx
    central_diff!(dy, y, x)
end

```
Upwind difference

```
function upwind_diff(
    y::AbstractArray{<:Real,1},
    x::AbstractArray{<:Real,1};
    stream = :right::Symbol,
)

    dy = zeros(eltype(y), axes(y))

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    if stream == :right
        dy[i0] = 0.0
        for i = i0+1:i1
            dy[i] = (y[i] - y[i-1]) / (x[i] - x[i-1])
        end
    elseif stream == :left
        dy[i1] = 0.0
        for i = i0:i1-1
            dy[i] = (y[i+1] - y[i]) / (x[i+1] - x[i])
        end
    else
        throw("streaming direction should be :left or :right")
    end

    return dy

end

function upwind_diff(y::AbstractArray{<:AbstractFloat,1}, dx::Real; stream = :right::Symbol)
    x = ones(eltype(y), axes(y)) .* dx
    dy = upwind_diff(y, x, stream = stream)

    return dy
end

function upwind_diff!(
    dy::AbstractArray{<:AbstractFloat,1},
    y::AbstractArray{<:AbstractFloat,1},
    x::AbstractArray{<:AbstractFloat,1};
    stream = :right::Symbol,
)

    @assert axes(dy) == axes(y) == axes(x)

    idx = eachindex(y) |> collect
    i0 = idx[1]
    i1 = idx[end]

    if stream == :right
        dy[i0] = 0.0
        for i = i0+1:i1
            dy[i] = (y[i] - y[i-1]) / (x[i] - x[i-1])
        end
    elseif stream == :left
        dy[i1] = 0.0
        for i = i0:i1-1
            dy[i] = (y[i+1] - y[i]) / (x[i+1] - x[i])
        end
    else
        throw("streaming direction should be :left or :right")
    end

    return dy

end

function upwind_diff!(
    dy::AbstractArray{<:AbstractFloat,1},
    y::AbstractArray{<:AbstractFloat,1},
    dx::Real;
    stream = :right::Symbol,
)
    x = ones(eltype(y), axes(y)) .* dx
    upwind_diff!(dy, y, x, stream = stream)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Ccache_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Ccache")
JLLWrappers.@generate_main_file("Ccache", UUID("c6bc53e4-e6d8-51d9-af6d-d7969aa08f60"))
end  # module Ccache_jll


text:
using WeibullCount

function test_weibull_count_pdf(x, λ, c; k = 200, transformed=true)
    for T1 ∈ [Int64, Int128, BigInt]
        for T2 ∈ [Float64, BigFloat]
            WeibullCount.cache_clear!()

            print("$T1 and $T2: ")
            println(@time weibull_count_pdf(T1(x), T2(λ), T2(c), T2(1.0); k = T1(k), transform=transformed))
        end
    end
end

# Ensure `weibull_count_pdf_approx` works with different types.
@testset "Parametric types" begin
    x = 2
    c = 0.36787944117144233
    λ = 20.08553692

    for T1 ∈ [Int64, Int128, BigInt]
        for T2 ∈ [Float64, BigFloat]
            WeibullCount.cache_clear!()

            # print("$T1 and $T2: ")
            res, converged = WeibullCount.weibull_count_pdf_approx(T1(x), T2(λ), T2(c), T2(1.0))
            @test isa(res, T2)
        end
    end
end

# Failure of convergence
@testset "DivergedError" begin
    x = 2
    c = 0.36787944117144233
    λ = 20.08553692

    WeibullCount.cache_clear!()
    @test_throws WeibullCount.DivergedError WeibullCount.weibull_count_pdf(x, λ, c)

    WeibullCount.cache_clear!()
    @test WeibullCount.weibull_count_pdf(BigInt(x), BigFloat(λ), BigFloat(c)) ≥ 0
end

# Difficult parameters
@testset "Run-times" begin
    x = 2
    c = 0.36787944117144233
    λ = 20.08553692

    WeibullCount.cache_clear!()
    t1 = @elapsed WeibullCount.weibull_count_pdf_old(BigInt(x), BigFloat(λ), BigFloat(c), BigFloat(1.0); k = BigInt(200))

    WeibullCount.cache_clear!()
    t2 = @elapsed weibull_count_pdf(BigInt(x), BigFloat(λ), BigFloat(c), BigFloat(1.0); k = BigInt(200), increment = BigInt(10), tol = BigFloat(0.001))

    @test t2 < t1
end


text:
#

# Low storage Runge-Kutta coefficients
rk4a = [ 0.0 ...
         -567301805773.0/1357537059087.0 ...
         -2404267990393.0/2016746695238.0 ...
         -3550918686646.0/2091501179385.0  ...
         -1275806237668.0/842570457699.0]
rk4b = [ 1432997174477.0/9575080441755.0 ...
         5161836677717.0/13612068292357.0 ...
         1720146321549.0/2090206949498.0  ...
         3134564353537.0/4481467310338.0  ...
         2277821191437.0/14882151754819.0]
rk4c = [ 0.0  ...
         1432997174477.0/9575080441755.0 ...
         2526269341429.0/6820363962896.0 ...
         2006345519317.0/3224310063776.0 ...
         2802321613138.0/2924317926251.0]

"
Integrates the RHS function rhsfunc
from t0 to t1, starting with solution u

This version modifies u
"
function lserk!(rhsfunc, u, t0, t1, dt; args=())
    Nsteps = ceil((t1 - t0)/dt)
    dt = (t1 - t0)/Nsteps

    # Runge-Kutta residual storage
    resu = similar(u)

    time = t0
    # outer time step loop
    for tstep=1:Nsteps
        for INTRK = 1:5
            timelocal = time + rk4c[INTRK]*dt
            rhsu = rhsfunc(u, timelocal, args...)
            resu = rk4a[INTRK]*resu + dt*rhsu
            u += rk4b[INTRK]*resu
        end
        # Increment time
        time += dt
    end

    return u
end

"
Integrates the RHS function rhsfunc
from t0 to t1, starting with solution u0
"
function lserk(rhsfunc, u0, t0, t1, dt; args=())
    u = copy(u0)
    return lserk!(rhsfunc, u, t0, t1, dt; args=args)
end


text:
abstract AbstractField

# Tensor fields (right now 6 components)
abstract AbstractTensor <: AbstractField
get_ncomponents{T <: AbstractTensor}(::Type{T}) = 6

# Vector fields
abstract AbstractVector <: AbstractField
get_ncomponents{T <: AbstractVector}(::Type{T}) = 3

# Scalar fields
abstract AbstractScalar <: AbstractField
get_ncomponents{T <: AbstractScalar}(::Type{T}) = 1
#get_ncomponents(::Type{AbstractTensor}) = 1 # Not needed for scalars?

abstract AbstractFullTensor <: AbstractField
get_ncomponents{T <: AbstractFullTensor}(::Type{T}) = 9

type KAlpha <: AbstractField end
get_ncomponents(::Type{KAlpha}) = 2

type Kappa <: AbstractField end
get_ncomponents(::Type{Kappa}) = 2

type τ <: AbstractField end
get_ncomponents(::Type{τ}) = 2
# TODO: Maybe we need a distinction between primary and secondary fields?
# Exporting a primary field as cell_data makes no sense? Or does it?

# Tensors
type Strain <: AbstractTensor end
type Stress <: AbstractTensor end

type VonMises <: AbstractScalar end

# Vectors
type Displacement <: AbstractVector end

# Scalars
type Pressure <: AbstractScalar end

type FullStress <: AbstractFullTensor end
type FullStrain <: AbstractFullTensor end
type InvFp <: AbstractFullTensor end
type FullStrain <: AbstractFullTensor end


text:
abstract type AbstractProgressiveHedgingSolver <: AbstractStructuredModel end

nscenarios(ph::AbstractProgressiveHedgingSolver) = ph.nscenarios

# Initialization #
# ======================================================================== #
function init!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver)
    # Initialize progress meter
    ph.progress.thresh = ph.parameters.τ
    # Finish initialization based on solver traits
    init_subproblems!(ph, subsolver)
    # Initialize penalty parameter (if applicable)
    ph.data.δ₁ = 1.0
    init_penalty!(ph)
end
# ======================================================================== #

# Functions #
# ======================================================================== #
function set_params!(ph::AbstractProgressiveHedgingSolver; kwargs...)
    for (k,v) in kwargs
        setfield!(ph.parameters, k, v)
    end
end

function current_objective_value(ph::AbstractProgressiveHedgingSolver, Qs::AbstractVector)
    return sum(Qs)
end
current_objective_value(ph) = current_objective_value(ph, ph.subobjectives)

function get_objective_value(ph::AbstractProgressiveHedgingSolver)
    if !isempty(ph.Q_history)
        return ph.Q_history[end]
    else
        return calculate_objective_value(ph)
    end
end

function iterate_nominal!(ph::AbstractProgressiveHedgingSolver)
    # Resolve all subproblems at the current optimal solution
    Q = resolve_subproblems!(ph)
    if Q == Inf
        return :Infeasible
    elseif Q == -Inf
        return :Unbounded
    end
    ph.data.Q = Q
    # Update iterate
    update_iterate!(ph)
    # Update subproblems
    update_subproblems!(ph)
    # Get dual gap
    update_dual_gap!(ph)
    # Update penalty (if applicable)
    update_penalty!(ph)
    # Update progress
    @unpack δ₁, δ₂ = ph.data
    ph.data.δ = sqrt(δ₁ + δ₂)/(1e-10+norm(ph.ξ,2))
    # Log progress
    log!(ph)
    # Check optimality
    if check_optimality(ph)
        # Optimal
        return :Optimal
    end
    # Just return a valid status for this iteration
    return :Valid
end

function log!(ph::AbstractProgressiveHedgingSolver)
    @unpack Q, δ, δ₂ = ph.data
    @unpack keep, offset, indent = ph.parameters
    push!(ph.Q_history, Q)
    push!(ph.dual_gaps, δ₂)
    ph.data.iterations += 1
    if ph.parameters.log
        ProgressMeter.update!(ph.progress,δ,
                              showvalues = [
                                  ("$(indentstr(indent))Objective",Q),
                                  ("$(indentstr(indent))δ",δ)
                              ], keep = keep, offset = offset)
    end
end

function indentstr(n::Integer)
    return repeat(" ", n)
end

function check_optimality(ph::AbstractProgressiveHedgingSolver)
    @unpack τ = ph.parameters
    @unpack δ = ph.data
    return δ <= τ
end
# ======================================================================== #
function show(io::IO, ph::AbstractProgressiveHedgingSolver)
    println(io, typeof(ph).name.name)
    println(io, "State:")
    show(io, ph.data)
    println(io, "Parameters:")
    show(io, ph.parameters)
end

function show(io::IO, ::MIME"text/plain", ph::AbstractProgressiveHedgingSolver)
    show(io, ph)
end


text:
using Test
using Unitful: m, s, cm
using UnitfulRecipes: recipe!, UnitFormatter
import RecipesBase

Attributes = Dict{Symbol, Any}
@testset "One Array" begin
    attr = Attributes()
    ys_val = [1, 2.3]
    ys = ys_val * m
    ys_ret = recipe!(attr, ys)
    @test ys_ret ≈ ys_val
    @test attr[:yformatter] == UnitFormatter(m)

    attr = Attributes(:yunit => cm)
    ys_ret = recipe!(attr, ys)
    @test ys_ret ≈ ys_val * 100
    @test !haskey(attr, :yunit)

    attr = Attributes(:ylims => (100cm, 2m))
    ys_ret = recipe!(attr, ys)
    @test ys_ret ≈ ys_val
    @test attr[:ylims] == (1, 2)
end

@testset "Multi Array" begin
    attr = Attributes()
    xs_val = randn(3)
    ys_val = randn(3)
    xu = s
    yu = m/s
    xs = xs_val * xu
    ys = ys_val * yu
    xs_ret, ys_ret = recipe!(attr, xs_val, ys)
    @test xs_ret ≈ xs_val
    @test ys_ret ≈ ys_val
    @test !haskey(attr, :xformatter)
    @test haskey(attr, :yformatter)

    xs_ret, ys_ret = recipe!(attr, xs, ys)
    @test xs_ret ≈ xs_val
    @test ys_ret ≈ ys_val
    @test haskey(attr, :xformatter)
    @test haskey(attr, :yformatter)

    zs_val = randn(3)
    xs_ret, ys_ret, zs_ret = recipe!(attr, xs, ys, zs_val)
    @test xs_ret ≈ xs_val
    @test ys_ret ≈ ys_val
    @test zs_ret ≈ zs_val
    @test haskey(attr, :xformatter)
    @test haskey(attr, :yformatter)
    @test !haskey(attr, :zformatter)
end

@testset "format" begin
    @test UnitFormatter(cm)(170.0) == "170.0cm"
    @test_broken UnitFormatter(cm)(170.000000000001) == "170.0cm"
end


text:
function f(x::Array{Float64,2})
  x[1]^2 + x[2]^2
end
function gradient(x::Array{Float64,2})
  [2*x[1] 2*x[2]]
end

x = [1.0 3.0]
dd = DifferentiableDensity(f,gradient)
@assert dd.f(x) == 10.0
@assert dd.gradient(x) == [2.0 6.0]


text:
module IOLA
using Statistics

export analyze, Codec

module Codec
import MDCT, DSP, ...IOLA

module Transform
@enum TransformType begin
    MDCT = 1
end
end

struct CodecParams
    transform::Transform.TransformType
    window::Function
    length::Integer
    hop::Integer
end

@enum CodecType begin
    MP3 = 1
    AAC = 2
    AC3 = 3
    OGG = 4
    WMA = 5
end

function getparams(type::CodecType)
    return if type == MP3
        CodecParams(Transform.MDCT,
                    DSP.Windows.cosine,
                    1152,
                    div(1152, 2))
    elseif type == AAC
        CodecParams(Transform.MDCT,
                    (n) -> IOLA.Windows.KBD(n, 4),
                    2048,
                    div(2048, 2))
    elseif type == OGG
        CodecParams(Transform.MDCT,
                    IOLA.Windows.slope,
                    2048,
                    div(2048, 2))
    elseif type == AC3
        CodecParams(Transform.MDCT,
                    (n) -> IOLA.Windows.KBD(n, 5),
                    512,
                    div(512, 2))
    elseif type == WMA
        CodecParams(Transform.MDCT,
                    DSP.Windows.cosine,
                    4096,
                    div(4096, 2))

    else
        error("Codec $type not yet implemented");
    end
end

function gettransform(params::CodecParams)
    if params.transform == Transform.MDCT
        win = params.window(params.length)
        mdct_fn = let mdct_plan = MDCT.plan_mdct(win)
            function(X::AbstractArray{T}) where T<:Number
                mdct_plan*X
            end
        end
        return IOLA.Utils.shorttermify(mdct_fn, win, params.hop)
    else
        error("Transform $(params.transform) not yet implemented!")
    end
end

end

module Windows
import DSP

function KBD(n::Integer, α::Real)
    if isodd(n)
        throw(ArgumentError("KBD window length must be even"))
    end
    N = div(n, 2)
    kais = DSP.Windows.kaiser(N+1, α)
    out = zeros(n)
    out[n] = (out[1] = kais[1])
    @inbounds for i = 2:N
        out[n+1-i] = (out[i] = kais[i] + out[i-1])
    end
    return sqrt.(out./out[N])
end

function slope(n::Integer)
    return sin.(0.5π.*sin.(π/n.*((1:n).-0.5)).^2)
end

end

module Utils

using Statistics
export findperiod, radiusofmean

function fast_sum_abs_log10_abs!(array::AbstractArray{SignalType}) where SignalType <: Number
    N = length(array)
    while N > 1
        for it = 1:2:N
            a1 = abs(array[it])
            if a1 == 0; return Inf; end
            it2 = div(it+1,2)
            if it+1 <= N
                a2 = abs(array[it+1])
                if a2 == 0; return Inf; end
                if (a1 >= 1 && a2 >= 1) || (a1 < 1 && a2 < 1)
                    a1 *= a2
                else
                    a1 /= a2
                end
                if a1 == 0 || a1 == Inf
                    return sum(abs.(log10.(view(array,1:it2-1)))) .+
                        sum(abs.(log10.(view(array,it:N))))
                end
            end
            array[it2] = a1
        end
        N = div(N+1,2)
    end
    return abs(log10(array[1]))
end

function shorttermify(fun::Function, window::AbstractVector{SignalType},
                      hop::Integer) where SignalType <: Number
    return let winlen = length(window), outlen = length(fun(window))
        function(signal::AbstractVector{SignalType})
            siglen = length(signal)
            N = div(siglen, hop)-div(winlen, hop)+1
            pad = div(siglen - (N-1)*hop - winlen, 2)
            signal_padded = zeros(SignalType, 2*pad + siglen + 1)
            signal_padded[pad+1:siglen+pad] = signal
            out = zeros(outlen, N)

            for i = 1:N
                sind = (i-1)*hop+1
                @views out[:,i] = fun(signal_padded[sind:(sind+winlen-1)].*window)
            end
            return out
        end
    end
end

function findperiod(vector::AbstractVector{T},
                    domain::AbstractVector) where T<:Number
    mod_val = typemax(T)
    p_min = 0
    for p in domain
        per = mod.(vector, p)
        val = std(per; corrected=false)/mean(per)
        if val <= mod_val
            mod_val = val
            p_min = p
        end
    end
    return p_min
end

function radiusofmean(radiuses::AbstractVector,
                      angles::AbstractVector,
                      na_rm = true)
    M = length(radiuses)
    idxs = na_rm ? isfinite.(radiuses) : Colon()
    sqrt(sum(radiuses[idxs].*sin.(angles[idxs]))^2 + sum(radiuses[idxs].*cos.(angles[idxs]))^2)/M
end

end

function energydiff(signal::AbstractVector{SignalType}, transform::Function,
                    segment_length::Integer, start_index::Integer,
                    end_index::Integer) where SignalType <: Number
    N = end_index-start_index+1
    energies = zeros(N)
    out = zeros(N-1)
    for i = start_index:end_index
        arind = i-start_index+1
        segment = transform(view(signal, i:(i+segment_length-1)))
        energies[arind] = Utils.fast_sum_abs_log10_abs!(segment)./segment_length
        if i > start_index
            out[arind-1] = abs(energies[arind-1] - energies[arind])
        end
    end
    return out
end

function analyze(signal::AbstractVector{SignalType}, transform::Function,
                 segment_length::Integer, overlap::Integer,
                 normalize::Bool = false) where SignalType <: Number
    siglen = length(signal)
    return let N = div(siglen-overlap, segment_length), out = zeros(N, 2)
        Threads.@threads for i = 1:N
            sind = (i-1)*segment_length+1
            segment_diffs = energydiff(signal, transform, segment_length, sind, sind+overlap)
            if normalize
                segment_diffs = (segment_diffs.-mean(segment_diffs))./std(segment_diffs)
            end
            out[i, :] = [findmax(segment_diffs)...]
        end
        out
    end
end

end


text:
function flux_Gxp_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)

    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]

    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S1 = ut*CUDAnative.sqrt(beta)
    B1 = 0.5*CUDAnative.exp(-S1*S1)/CUDAnative.sqrt(pi*beta)
    A1pos = 0.5*(1 + CUDAnative.erf(S1))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    shared[thread_idx] = op((rho*(ut*A1pos + B1)), shared[thread_idx])
    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1pos + ut*B1
    shared[thread_idx + block_dim] = op((rho*temp2), shared[thread_idx + block_dim])
    temp1 = ut*un*A1pos + un*B1
    shared[thread_idx + block_dim * 2] = op((rho*temp1), shared[thread_idx + block_dim * 2])

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1pos
    temp1 = (6*pr_by_rho) + u_sqr
    shared[thread_idx + block_dim * 3] = op((rho*(temp2 + 0.5*temp1*B1)), shared[thread_idx + block_dim * 3])

    return nothing
end

function flux_Gxn_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)

    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]

    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S1 = ut*CUDAnative.sqrt(beta)
    B1 = 0.5*CUDAnative.exp(-S1*S1)/CUDAnative.sqrt(pi*beta)
    A1neg = 0.5*(1 - CUDAnative.erf(S1))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un
    shared[thread_idx] = op((rho*(ut*A1neg - B1)), shared[thread_idx])

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1neg - ut*B1
    shared[thread_idx + block_dim] = op((rho*temp2), shared[thread_idx + block_dim])
    temp1 = ut*un*A1neg - un*B1
    shared[thread_idx + block_dim * 2] = op((rho*temp1), shared[thread_idx + block_dim * 2])
    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1neg
    temp1 = (6*pr_by_rho) + u_sqr
    shared[thread_idx + block_dim * 3] = op((rho*(temp2 - 0.5*temp1*B1)), shared[thread_idx + block_dim * 3])
    return nothing
end

function flux_Gyp_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)
    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]

    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S2 = un*CUDAnative.sqrt(beta)
    B2 = 0.5*CUDAnative.exp(-S2*S2)/CUDAnative.sqrt(pi*beta)
    A2pos = 0.5*(1 + CUDAnative.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    shared[thread_idx] = op((rho*(un*A2pos + B2)), shared[thread_idx])
    temp1 = pr_by_rho + un*un
    temp2 = temp1*A2pos + un*B2
    temp1 = ut*un*A2pos + ut*B2
    #TODO - Verify this is correct
    shared[thread_idx + block_dim] = op((rho*temp1), shared[thread_idx + block_dim])
    shared[thread_idx + block_dim * 2] = op((rho*temp2), shared[thread_idx + block_dim * 2])

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*un*temp1*A2pos
    temp1 = (6*pr_by_rho) + u_sqr
    shared[thread_idx + block_dim * 3] = op((rho*(temp2 + 0.5*temp1*B2)), shared[thread_idx + block_dim * 3])
    return nothing
end

function flux_Gyn_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)
    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]

    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    S2 = un*CUDAnative.sqrt(beta)
    B2 = 0.5*CUDAnative.exp(-S2*S2)/CUDAnative.sqrt(pi*beta)
    A2neg = 0.5*(1 - CUDAnative.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    shared[thread_idx] = op((rho*(un*A2neg - B2)), shared[thread_idx])
    temp1 = pr_by_rho + un*un
    temp2 = temp1*A2neg - un*B2
    temp1 = ut*un*A2neg - ut*B2
    #TODO - Verify this is correct
    shared[thread_idx + block_dim] = op((rho*temp1), shared[thread_idx + block_dim])
    shared[thread_idx + block_dim * 2] = op((rho*temp2), shared[thread_idx + block_dim * 2])

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*un*temp1*A2neg
    temp1 = (6*pr_by_rho) + u_sqr
    shared[thread_idx + block_dim * 3] = op((rho*(temp2 - 0.5*temp1*B2)), shared[thread_idx + block_dim * 3])

    return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
"""
**GibbsSampling**

Draw samples from the true posterior via Gibbs Sampling.

```julia
GibbsSampling(;ϵ::T=1e-5,nBurnin::Int=100,samplefrequency::Int=10)
```
**Keywords arguments**

    - `ϵ::T` : convergence criteria
    - `nBurnin::Int` : Number of samples discarded before starting to save samples
    - `samplefrequency::Int` : Frequency of sampling
"""
mutable struct GibbsSampling{T<:Real} <: Inference{T}
    nBurnin::Integer
    samplefrequency::Integer
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 # Number of samples
    nSamplesUsed::Int64 #Size of mini-batches
    MBIndices::AbstractVector #Indices of the minibatch
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    sample_store::AbstractVector{AbstractVector{AbstractVector{T}}}
    x::SubArray{T,2,Matrix{T}}#,Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true}
    y::LatentArray{SubArray}
    function GibbsSampling{T}(nBurnin::Int,samplefrequency::Int,ϵ::T,nIter::Integer,Stochastic::Bool,nSamples::Integer,nSamplesUsed::Integer,MBIndices::AbstractVector,ρ::T,flag::Bool) where T
        return new{T}(nBurnin,samplefrequency,ϵ,nIter,Stochastic,nSamples,nSamplesUsed,MBIndices,ρ,flag)
    end
end

function GibbsSampling(;ϵ::T=1e-5,nBurnin::Int=100,samplefrequency::Int=10) where {T<:Real}
    GibbsSampling{Float64}(nBurnin,samplefrequency,ϵ,0,false,1,1,[1],1.0,true)
end

function Base.show(io::IO,inference::GibbsSampling{T}) where {T<:Real}
    print(io,"Gibbs Sampler")
end

function init_inference(inference::GibbsSampling{T},nLatent::Integer,nFeatures::Integer,nSamples::Integer,nSamplesUsed::Integer) where {T<:Real}
    inference.nSamples = nSamples
    inference.nSamplesUsed = nSamples
    inference.MBIndices = collect(1:nSamples)
    inference.sample_store = [Vector{Vector{T}}() for _ in 1:nLatent]
    return inference
end

function variational_updates!(model::VGP{T,L,GibbsSampling{T}}) where {T,L}
    sample_local!(model)
    sample_global!(model)
    if model.inference.nIter > model.inference.nBurnin && (model.inference.nIter-model.inference.nBurnin)%model.inference.samplefrequency==0
        for k in 1:model.nLatent
            push!(model.inference.sample_store[k],model.μ[k])
        end
    end
end

function ELBO(model::AbstractGP{T,<:Likelihood,<:GibbsSampling}) where {T}
    return NaN
end

function sample_global!(model::VGP{T,<:Likelihood,<:GibbsSampling}) where {T}
    model.Σ .= inv.(Symmetric.(Diagonal.(2.0.*∇E_Σ(model)).+model.invKnn))
    model.μ .= rand.(MvNormal.(model.Σ.*(∇E_μ(model).+model.invKnn.*model.μ₀),model.Σ))
    return nothing
end

function post_process!(model::AbstractGP{T,<:Likelihood,<:GibbsSampling}) where {T}
    for k in 1:model.nLatent
        model.μ[k] = vec(mean(hcat(model.inference.sample_store[k]...),dims=2))
        model.Σ[k] = Symmetric(cov(hcat(model.inference.sample_store[k]...),dims=2))
    end
    nothing
end


text:
@testset verbose=true "Inequalities                " begin

	n = 16
	m = 5

	# no bounds, lower bound, upper bound, both bounds
	divsize = n ÷ 4
	xl = [-Inf.*ones(divsize); randn(divsize); -Inf.*ones(divsize); randn(divsize)]
	xu = [Inf.*ones(divsize); Inf.*ones(divsize); randn(divsize); randn(divsize)]

	x = [randn(n÷4);
		xl[divsize+1:2*divsize] .+ rand(0:2, divsize);
		xu[2*divsize+1:3*divsize] .- rand(0:2, divsize);
		xl[3*divsize+1:4*divsize] .+ rand(divsize).*(xu[3*divsize+1:4*divsize] .- xl[3*divsize+1:4*divsize])]

	xaug = zeros(2*n)
	xaug[1:n] .= x

	idata = LFPSQP.InequalityData(xl, xu)

	# test aspects of Inequality data
	@testset "InequalityData" begin
		@test all(idata.q .≈ [zeros(3*divsize); ones(divsize)])
		@test all(idata.r .≈ [zeros(divsize);
								xl[divsize+1:2*divsize];
								xu[2*divsize+1:3*divsize];
								xl[3*divsize+1:4*divsize]./2 .+ xu[3*divsize+1:4*divsize]./2])
		@test all(idata.s .≈ [zeros(divsize); -ones(divsize); ones(2*divsize)])
		@test all(idata.t .≈ [zeros(divsize);
								xl[divsize+1:2*divsize];
								xu[2*divsize+1:3*divsize];
								(xu[3*divsize+1:4*divsize] .- xl[3*divsize+1:4*divsize]).^2 ./ 4])

		@test all(idata.isline .== vcat(trues(divsize), falses(3*divsize)))
		@test all(idata.isparabola .== vcat(falses(divsize), trues(2*divsize), falses(divsize)))
	end

	# test y generation satisfies inequalities
	@testset "Initial y" begin
		LFPSQP.generate_initial_y!(xaug, idata)

		xview = view(xaug, 1:n)
		yview = view(xaug, n+1:2*n)

		@test all(isapprox.(idata.q.*(xview .- idata.r).^2 .+ (1.0 .- idata.q.^2).*xview .+
			idata.s.*(yview .- idata.r).^2 .- (1.0 .- idata.s.^2).*yview .- idata.t, 0.0, atol=1e-15))

		cvalaug = zeros(n+m)
		LFPSQP.calculate_h!(cvalaug, xaug, idata)

		@test all(isapprox.(cvalaug, 0.0, atol=1e-15))
	end

	# generate random decomposition
	Jct = randn(n, m)

	∇hx = 2.0*idata.q.*(x .- idata.r) .+ (1.0 .- idata.q.^2)
	∇hy = 2.0*idata.s.*(xaug[n+1:end] .- idata.r) .- (1.0 .- idata.s.^2)

	S = sqrt.(∇hx.^2 .+ ∇hy.^2)
	Dx = ∇hx ./ S
	Dy = ∇hy ./ S

	R = Dx .* Jct

	PJct = zeros(2*n, m)
	PJct[1:n, :] .= (1.0 .- Dx.*Dx) .* Jct
	PJct[n+1:2*n, :] .= -1.0 .* Dy.*Dx .* Jct

	U, Σ, V = svd(PJct)
	Vt = Matrix(V')

	bigA = [diagm(∇hx) Jct; diagm(∇hy) zeros(n, m)]
	bigQ = hcat([diagm(Dx); diagm(Dy)], U)
	bigR = [diagm(S) R; zeros(m, n) diagm(Σ)*Vt]

	idecomp = LFPSQP.InequalityDecomp(U, Σ, Vt, [Array{Float64}(undef, n) for i in 1:3]..., Jct, m)

	# test InequalityDecomposition
	@testset "InequalityDecomposition" begin

		LFPSQP.inequality_gradient!(idecomp, xaug, idata)

		@test all(isapprox.(idecomp.Dx, Dx, atol=1e-15))
		@test all(isapprox.(idecomp.Dy, Dy, atol=1e-15))
		@test all(isapprox.(idecomp.S, S, atol=1e-15))

		@test all(isapprox.(bigA, bigQ*bigR, atol=1e-14))					# decomposition is correct
		@test all(isapprox.(bigQ'*bigQ, Matrix(I, n+m, n+m), atol=1e-14))	# orthogonality
	end

	@testset "Multiplication" begin
		v = randn(n+m)
		w = randn(2*n)

		destv = zeros(2*n)
		destw = zeros(n+m)

		# projections
		mul!(destv, LFPSQP.InequalityDecompProject(idecomp), v)
		@test all(isapprox.(destv, bigQ*v, atol=1e-15))

		mul!(destw, LFPSQP.InequalityDecompProject(idecomp)', w)
		@test all(isapprox.(destw, bigQ'*w, atol=1e-15))

		destv .= 1.0
		mul!(destv, LFPSQP.InequalityDecompProject(idecomp), v, 2.0, 3.0)
		@test all(isapprox.(destv, 2*bigQ*v .+ 3, atol=1e-15))

		# Jacobians
		mul!(destv, idecomp, v)
		@test all(isapprox.(destv, bigA*v, atol=1e-15))

		mul!(destw, idecomp', w)
		@test all(isapprox.(destw, bigA'*w, atol=1e-15))

		destv .= 1.0
		mul!(destv, idecomp, v, 2.0, 3.0)
		@test all(isapprox.(destv, 2*bigA*v .+ 3, atol=3e-15))

		# test lower rank
		idecomp.rank = m - 2

		v = randn(n+m-2)
		destv = zeros(2*n)
		destw = zeros(n+m-2)

		mul!(destv, LFPSQP.InequalityDecompProject(idecomp), v)
		@test all(isapprox.(destv, bigQ[:, 1:n+m-2]*v, atol=1e-15))

		mul!(destw, LFPSQP.InequalityDecompProject(idecomp)', w)
		@test all(isapprox.(destw, bigQ[:, 1:n+m-2]'*w, atol=1e-15))

		destv .= 1.0
		mul!(destv, LFPSQP.InequalityDecompProject(idecomp), v, 2.0, 3.0)
		@test all(isapprox.(destv, 2*bigQ[:, 1:n+m-2]*v .+ 3, atol=1e-15))

		idecomp.rank = m
	end

	@testset "Lagrange multipliers" begin
		d = randn(2*n)
		Qt∇f = zeros(n+m)
		mul!(Qt∇f, LFPSQP.InequalityDecompProject(idecomp)', d)

		λ = zeros(n+m)
		λ_kkt = view(λ, n+1:n+m)
		λy_kkt = view(λ, 1:n)

		LFPSQP.calculate_λ_kkt!(λ_kkt, λy_kkt, Qt∇f, idecomp)

		@test all(isapprox.(λ, bigA \ d, atol=1e-14))
	end

	@testset "Hessian action" begin

		A = randn(n, n)
		C = randn(n, n, m)
		λ_kkt = randn(m)
		λy_kkt = randn(n)

		dest = zeros(2*n)
		v = randn(2*n)

		function hess_lag_vec!(dest, src, x, λ)
			dest .= A*src + sum(λ[i]*C[:,:,i]*src for i in 1:m)
		end

		LFPSQP.augmented_hess_lag_vec!(dest, v, hess_lag_vec!, xaug, λ_kkt, λy_kkt, idata)

		H = A + sum(λ_kkt[i]*C[:,:,i] for i in 1:m)
		bigH = [H+2*diagm(λy_kkt.*idata.q) zeros(n, n); zeros(n,n) 2*diagm(λy_kkt.*idata.s)]

		@test all(isapprox.(dest, bigH*v, atol=1e-14))
	end

	@testset "y retraction" begin
		# take step in tangent plane
		d = randn(2*n)
		tmp_m = zeros(n+m)
		ineqproject = LFPSQP.InequalityDecompProject(idecomp)

		mul!(tmp_m, ineqproject', d)				# tmp_m = Q'd
		mul!(d, ineqproject, tmp_m, -1.0, 1.0)		# d = d - Q*tmp_m

		# test really long step length
		xnewaug = xaug .+ d
		xaugcopy = copy(xaug)
		LFPSQP.y_retract!(xnewaug, xaug, idata)

		# ensure y retraction worked and that inequality constraints are 0
		cvalaug = zeros(n+m)
		LFPSQP.calculate_h!(cvalaug, xnewaug, idata)

		@test all(isapprox.(cvalaug, 0.0, atol=1e-13))
		@test all(xaug .== xaugcopy)					# ensure xaug not overwritten
	end
end


text:
"""
    mutable struct DeterministicSingleTimeSeries <: AbstractDeterministic
        single_time_series::SingleTimeSeries
        initial_timestamp::Dates.DateTime
        interval::Dates.Period
        count::Int
        horizon::Int
    end

A deterministic forecast for a particular data field in a Component that wraps a SingleTimeSeries.

# Arguments
- `single_time_series::SingleTimeSeries`: wrapped SingleTimeSeries object
- `initial_timestamp::Dates.DateTime`: time series availability time
- `interval::Dates.Period`: time step between forecast windows
- `count::Int`: number of forecast windows
- `horizon::Int`: length of this time series
"""
mutable struct DeterministicSingleTimeSeries <: AbstractDeterministic
    "wrapped SingleTimeSeries object"
    single_time_series::SingleTimeSeries
    "time series availability time"
    initial_timestamp::Dates.DateTime
    "time step between forecast windows"
    interval::Dates.Period
    "number of forecast windows"
    count::Int
    "length of this time series"
    horizon::Int
end

function DeterministicSingleTimeSeries(;
    single_time_series,
    initial_timestamp,
    interval,
    count,
    horizon,
)
    return DeterministicSingleTimeSeries(
        single_time_series,
        initial_timestamp,
        interval,
        count,
        horizon,
    )
end

get_name(value::DeterministicSingleTimeSeries) = get_name(value.single_time_series)
"""Get [`DeterministicSingleTimeSeries`](@ref) `single_time_series`."""
get_single_time_series(value::DeterministicSingleTimeSeries) = value.single_time_series
"""Get [`DeterministicSingleTimeSeries`](@ref) `initial_timestamp`."""
get_initial_timestamp(value::DeterministicSingleTimeSeries) = value.initial_timestamp
"""Get [`DeterministicSingleTimeSeries`](@ref) `interval`."""
get_interval(value::DeterministicSingleTimeSeries) = value.interval
"""Get [`DeterministicSingleTimeSeries`](@ref) `count`."""
get_count(value::DeterministicSingleTimeSeries) = value.count
"""Get [`DeterministicSingleTimeSeries`](@ref) `horizon`."""
get_horizon(value::DeterministicSingleTimeSeries) = value.horizon

"""Set [`DeterministicSingleTimeSeries`](@ref) `single_time_series`."""
set_single_time_series!(value::DeterministicSingleTimeSeries, val) =
    value.single_time_series = val
"""Set [`DeterministicSingleTimeSeries`](@ref) `initial_timestamp`."""
set_initial_timestamp!(value::DeterministicSingleTimeSeries, val) =
    value.initial_timestamp = val
"""Set [`DeterministicSingleTimeSeries`](@ref) `interval`."""
set_interval!(value::DeterministicSingleTimeSeries, val) = value.interval = val
"""Set [`DeterministicSingleTimeSeries`](@ref) `count`."""
set_count!(value::DeterministicSingleTimeSeries, val) = value.count = val
"""Set [`DeterministicSingleTimeSeries`](@ref) `horizon`."""
set_horizon!(value::DeterministicSingleTimeSeries, val) = value.horizon = val

eltype_data(ts::DeterministicSingleTimeSeries) = eltype_data(ts.single_time_series)
get_scaling_factor_multiplier(ts::DeterministicSingleTimeSeries) =
    get_scaling_factor_multiplier(ts.single_time_series)

function get_array_for_hdf(forecast::DeterministicSingleTimeSeries)
    return get_array_for_hdf(forecast.single_time_series)
end

get_resolution(val::DeterministicSingleTimeSeries) = get_resolution(val.single_time_series)

function get_window(
    forecast::DeterministicSingleTimeSeries,
    initial_time::Dates.DateTime;
    len::Union{Nothing, Int} = nothing,
)
    tdiff = Dates.Millisecond(initial_time - forecast.initial_timestamp)
    interval_ms = Dates.Millisecond(forecast.interval)
    if interval_ms != Dates.Millisecond(0) && tdiff % interval_ms != Dates.Millisecond(0)
        throw(ArgumentError("initial_time=$initial_time is not on a window boundary"))
    end

    if len === nothing
        len = forecast.horizon
    end

    ta = get_data(forecast.single_time_series)
    resolution = get_resolution(forecast)
    end_time = initial_time + (len - 1) * resolution
    timestamps = TimeSeries.timestamp(ta)
    for timestamp in (initial_time, end_time)
        @assert timestamp >= first(timestamps) && timestamp <= last(timestamps) "invalid " *
                                                                                "timestamp=$timestamp is not within $(first(timestamps)) - $(last(timestamps))"
    end

    return ta[initial_time:resolution:end_time]
end

function iterate_windows(forecast::DeterministicSingleTimeSeries)
    if get_count(forecast) == 1
        return (get_window(forecast, get_initial_timestamp(forecast)),)
    end

    initial_times =
        range(forecast.initial_timestamp; step = forecast.interval, length = forecast.count)
    return (get_window(forecast, it) for it in initial_times)
end

function deserialize_deterministic_from_single_time_series(
    storage::TimeSeriesStorage,
    ts_metadata::DeterministicMetadata,
    rows,
    columns,
    last_index,
)
    @debug "deserializing a SingleTimeSeries" _group = LOG_GROUP_TIME_SERIES
    horizon = get_horizon(ts_metadata)
    interval = get_interval(ts_metadata)
    resolution = get_resolution(ts_metadata)
    if length(rows) != horizon
        throw(
            ArgumentError(
                "Transforming SingleTimeSeries to Deterministic requires a full horizon: $rows",
            ),
        )
    end

    sts_rows =
        _translate_deterministic_offsets(horizon, interval, resolution, columns, last_index)
    sts = deserialize_time_series(
        SingleTimeSeries,
        storage,
        SingleTimeSeriesMetadata(ts_metadata),
        sts_rows,
        UnitRange(1, 1),
    )
    initial_timestamp =
        get_initial_timestamp(ts_metadata) + (columns.start - 1) * get_interval(ts_metadata)
    return DeterministicSingleTimeSeries(
        sts,
        initial_timestamp,
        interval,
        length(columns),
        horizon,
    )
end

function _translate_deterministic_offsets(
    horizon,
    interval,
    resolution,
    columns,
    last_index,
)
    interval = Dates.Millisecond(interval)
    interval_offset = Int(interval / resolution)
    s_index = (columns.start - 1) * interval_offset + 1
    e_index = (columns.stop - 1) * interval_offset + horizon
    @debug "translated offsets" _group = LOG_GROUP_TIME_SERIES horizon columns s_index e_index last_index
    @assert_op s_index <= last_index
    @assert_op e_index <= last_index
    return UnitRange(s_index, e_index)
end