text:
# compare-hyperparameters-gmm-posterior
#
# This script compares the quality of samples generated from different
# step size hyperparameter settings of a Stochastic Gradient Langevin Dynamics
# sampler for a Gaussian mixture model posterior target distribution.

include("experiment_utils.jl")

# set the discrepancy
@parsestringcli discrepancytype "k" "discrepancytype" "inversemultiquadric"
# set the sampler
@parsestringcli sampler "q" "sampler" "sgld"
# the number of likelihood terms to sample (0 means all)
@parseintcli likelihood_n "l" "likelihoodn" 0
# thinning posterior samples to save computation
@parseintcli thinningn "t" "thinningn" 1
# how many likelihood terms
@parseintcli L "L" "likelihoodterms" 100
# how many samples
@parseintcli n "n" "numsamples" 1000
# should we use multiple cores?
@parseintcli numcores "n" "numcores" 1
# parallelize it
if (numcores > 1)
    if (numcores > CPU_CORES)
        error(@sprintf("Requested %d cores, only %d available.", numcores, CPU_CORES))
    end
    addprocs(numcores - 1)
    @assert (numcores == nprocs())
end

using Iterators: product
import SteinDistributions
using SteinDiscrepancy:
    SteinInverseMultiquadricKernel,
    ksd,
    gsd
using SteinSamplers:
    runsgld,
    runapproxslice

# Run one complete experiment per random number generator seed
seeds = 1:50
# Select a solver for graph Stein discrepancy optimization problem
#solver = "clp"
solver = "gurobi"
## Sampler settings for SGLD
distname = "gmm-posterior"
# Generate random dataset from GMM
truex = [0.0, 1.0] # Parameter vector used to generate dataset
data_generation_seed = 1
sigma2y = 2.0
y = @setseed data_generation_seed SteinDistributions.randgmm(truex, L; sigma2y=sigma2y)
# save the y values for later
save_json(
    Dict{Any, Any}("y" => y);
    dir="compare-hyperparameters-$(distname)-y",
    numsamples=L,
    x=string(truex)
)
# Create posterior distribution based on dataset
@eval @everywhere target = SteinDistributions.SteinGMMPosterior($y)
@eval @everywhere likelihoodn = $likelihood_n
# define gradlogp
if likelihood_n == 0
    @everywhere function gradlogp(x::Array{Float64,1})
        SteinDistributions.gradlogdensity(target, x)
    end
else
    @everywhere function gradlogp(x::Array{Float64,1})
        L = SteinDistributions.numdatapoints(target)
        SteinDistributions.gradlogdensity(target, x; idx=randperm(L)[1:likelihoodn])
    end
end
# set the kernel
kernel = nothing
if discrepancytype == "inversemultiquadric"
    kernel = SteinInverseMultiquadricKernel(0.5)
end
# the max number of likelihood evals to run for
# (if seed=7, epsilon=0.1, sampler=approxslice, L=100, n=1000, maxeval=147695)
# (if seed=7, epsilon=0.2, sampler=approxslice, L=100, n=2000, maxeval=152875)
maxevals = 147695
# the SGLD epsilons to run the experiemnt
sgld_epsilons = [5e-2, 1e-2, 5e-3, 1e-3, 5e-4, 1e-4, 5e-5, 1e-5]
# the approx slice epsilons to run the experiemnt
approxslice_epsilons = [0.0, 1e-2, 1e-1]
# Fix mini-batch size
batchsize = 5

if sampler == "sgld"
    # the number of sweeps to make through the data
    numsweeps = round(Int, ceil(n / floor(SteinDistributions.numdatapoints(target)/batchsize)))
    epsilons = sgld_epsilons
    runsampler(x0,epsilon) = runsgld(target, x0;
                                     epsilonfunc=t -> epsilon,
                                     numsweeps=numsweeps,
                                     batchsize=batchsize)
elseif sampler == "approxslice"
    epsilons = approxslice_epsilons
    runsampler(x0,epsilon) = runapproxslice(target, x0;
                                            epsilon=epsilon,
                                            numlikelihood=maxevals,
                                            batchsize=batchsize)
end

## Draw and evaluate samples for different hyperparameter and seed combinations
d = SteinDistributions.numdimensions(target)
for (epsilon, seed) in product(epsilons, seeds)
    # Generate starting sample value x0
    x0 = @setseed seed randn(d)

    println("Generating sample for seed=$(seed), epsilon=$(epsilon)")
    # Run sampler
    points = []; numgrad = NaN
    @trycatchcontinue(
        begin
            (points, numgrad) = @setseed seed runsampler(x0,epsilon)
            if any(isnan(points))
                println("[$(distname):seed=$(seed),epsilon=$(epsilon)] NaNs found. Skipping.")
                continue
            end
        end,
        println("[$(distname):seed=$(seed),epsilon=$(epsilon)]:")
    )

    # thin the sample if necessary
    thinpoints = points[thinningn:thinningn:size(points,1),:]

    # Compute Stein discrepancy for first n points in sample
    # (with equal point weights)
    println("Computing Stein discrepancy for n=$(n), d=$(d), dist=$(distname)")
    objectivevalue = nothing
    edgetime = nothing
    solvetime = nothing

    if discrepancytype == "graph"
        res = nothing
        @trycatchcontinue(
            begin
                res = gsd(points=thinpoints, gradlogdensity=gradlogp, solver=solver)
            end,
            println("[n=$(n)|sampler=$(distname)]:")
        )

        objectivevalue = res.objectivevalue
        edgetime = res.edgetime
        solvetime = res.solvetime
    elseif kernel != nothing
        res = ksd(points=thinpoints, gradlogdensity=gradlogp, kernel=kernel)
        objectivevalue = sqrt(res.discrepancy2)
        solvetime = res.solvetime
    end

    println("\tn = $(n), objective = $(objectivevalue)")
    # Package results
    instance_data = Dict{Any, Any}(
        "distname" => distname,
        "n" => size(points, 1),
        "sampler" => sampler,
        "solver" => solver,
        "d" => d,
        "seed" => seed,
        "epsilon" => epsilon,
        "X" => points,
        "q" => res.weights,
        "objectivevalue" => objectivevalue,
        "likelihoodn" => likelihood_n,
        "thinningn" => thinningn,
        "numgrad" => numgrad,
        "edgetime" => edgetime,
        "solvetime" => solvetime,
        "ncores" => nprocs(),
    )

    # Save results
    save_json(
        instance_data;
        dir="compare-hyperparameters-$(distname)",
        distname=distname,
        discrepancytype=discrepancytype,
        sampler=sampler,
        n=n,
        solver=solver,
        epsilon=epsilon,
        ncores=nprocs(),
        likelihoodn=likelihood_n,
        thinningn=thinningn,
        L=L,
        d=d,
        seed=seed,
    )
end


text:
export
    Ising

import Base: *

"""
HOBO classes. Current HOBO object should not contain zero(T) values, which may result
in unnecessary memory blow.
"""
struct Ising{T} <: BinaryOptimizationModel{T}
    vars::Variables
    values::Dict{Vector{Vector{Int}},T}
    function Ising{T}(vars::Variables,
                      values::Dict{Vector{Vector{Int}},T}) where T<:Real
        final_values = Dict{Vector{Vector{Int}},T}()

        for (key, val) = values
            final_key = clean_key(Ising, key)
            @assert all(length(i) == dim(vars) for i = key)
            @assert all(all(0 .< ind .<= indexcapacities(vars)) for ind=final_key)
            _dict_add!(final_values, final_key, val)
        end
        new{T}(vars, final_values)
    end
end

function Ising(::Type{T}, vars::Variables) where T
    Ising{T}(vars, Dict{Vector{Vector{Int}},T}())
end

Ising(vars::Variables) = Ising(Float64, vars::Variables)
Ising{T}(vars::Variables) where T<:Real = Ising(T, vars::Variables)

function Ising(num::T, vars::Variables) where T<:Real
    Ising{T}(vars, Dict(Vector{Int}[] => num))
end

function clean_key(::Type{Ising}, key::Vector{Vector{Int}})
    unique_keys = unique(key)
    result = Vector{Int}[]
    for k in unique_keys
        if isodd(count(x -> x == k, key))
            push!(result, k)
        end
    end
    sort(result)
end

clean_key(::Type{Ising{T}}, key::Vector{Vector{Int}}) where T<:Real = clean_key(HOBO, key)

## may be improved - smart tree walker
function (ising::Ising{T})(keys_list::Set{Vector{Int}}) where T<:Real
    result = zero(T)
    for (key, val) = values(ising)
        if isodd(length(setdiff(Set(key), Set(keys_list))))
            result -= val
        else
            result += val
        end
    end
    result
end

#(hobo::HOBO{<:Real})(keys_list) = (hobo::HOBO)(Set(keys_list...))
#(hobo::HOBO{<:Real})(keys_list::Vector{Vector{Int}}) = (hobo::HOBO)([keys_list])

ext: jl


text:
## Lecture-3-Least-squares-CVX

## Here are the packages we'll use
using Convex
using Plots
plotlyjs() # need to use plotlyjs on my machine for plots.

## Step 1: Create a set of data with a linear model
m=40
n=2
A = randn(m,n)
xex = [5;1]      # "b" = 5 in ax+b and a=1
pts = -10+20*rand(m,1)
A = [ones(m,1) pts]
b = A*xex + .5*randn(m,1)

## Show the linear model
scatter(pts,b)
xlabel!("x")
ylabel!("y")

##
x = Variable(n)
problem = minimize(sumsquares(b - A*x))
solve!(problem)
xls = x.value
@show x
## Show the least squares fit
scatter(pts,b;label="data")
plot!([-11; 11], [1 -11; 1 11]*xls;label="fit")
#xaxis!([-11 11])
title!("Least-square fit")
xlabel!("x")
ylabel!("y")

## Now we add outliers
outliers = [-9.5; 9]
outvals = [20; -15]
A = [A; ones(length(outliers),1) outliers]
b = [b; outvals]
m = size(A,1)
pts = [pts;outliers]

## Show the new data
scatter(pts,b)
xlabel!("x")
ylabel!("y")

## Look at the LS fit
x = Variable(n)
problem = minimize(sumsquares(b - A*x))
solve!(problem)
xls = x.value
plot!([-11; 11], [1 -11; 1 11]*xls;label="fit")

## Solve the Huber problem and look at the fit
x = Variable(n)
problem = minimize(sum(huber(b - A*x)))
solve!(problem)
xr = x.value
plot!([-11; 11], [1 -11; 1 11]*xr;label="fit_huber")


text:
module Tests

using AnonymousTypes
using Base.Test

a = @type x = 1 y = 2.0

@test a.x == 1
@test a.y == 2.0

@test fieldtype(typeof(a), :x) == Int
@test fieldtype(typeof(a), :y) == Float64

@test (a.x += 1) == 2

b = @immutable x = 1 y = Int[] z = ""

@test b.x == 1
@test b.y == Int[]
@test b.z == ""

@test fieldtype(typeof(b), :x) == Int
@test fieldtype(typeof(b), :y) == Vector{Int}
@test fieldtype(typeof(b), :z) == ASCIIString

@test_throws ErrorException (b.x += 1)

c = @type a b

@test c.a == a
@test c.b == b

@test_throws ErrorException AnonymousTypes.buildcall(:m, 1)
@test_throws ErrorException AnonymousTypes.buildcall(:m, :(x, y))

f(:: @Anon(a, b)) = 1
f(:: @Anon(x :: Int, y :: AbstractString)) = 2
f(:: @Anon(:: Integer, :: AbstractFloat, :: Vector)) = 3

g(:: @Type(a, b)) = 1
g(:: @Type(x :: Int, y :: AbstractString)) = 2
g(:: @Type(:: Integer, :: AbstractFloat, :: Vector)) = 3

g(:: @Immutable(a, b)) = 4
g(:: @Immutable(x :: Int, y :: AbstractString)) = 5
g(:: @Immutable(:: Integer, :: AbstractFloat, :: Vector)) = 6

x_1 = @type      a = 1 b = 2
x_2 = @immutable a = 2 b = 3

y_1 = @type      x = 1 y = ""
y_2 = @immutable x = 1 y = ""

z_1 = @type      a = 1 b = 1.0 c = [1, 2, 3]
z_2 = @immutable a = 1 b = 1.0 c = [1, 2, 3]

@test f(x_1) == 1
@test g(x_1) == 1
@test f(x_2) == 1
@test g(x_2) == 4

@test f(y_1) == 2
@test g(y_1) == 2
@test f(y_2) == 2
@test g(y_2) == 5

@test f(z_1) == 3
@test g(z_1) == 3
@test f(z_2) == 3
@test g(z_2) == 6

@test_throws ErrorException AnonymousTypes.buildsig(true, 1)

end

max_stars_repo_name: MichaelHatherly/AnonymousTypes

max_issues_repo_name: MichaelHatherly/AnonymousTypes

max_forks_repo_name: MichaelHatherly/AnonymousTypes


text:
"""
    SeisRadonFreqFor(m, nt; <keyword arguments>)

Transform a tau-p gather to a time-offset gather using a frequency domain
forward parabolic or linear Radon operator.

# Arguments
* `m::Array{T<:Real,2}`: 2D Radon panel, `m[1:ntau,1:np]`, where `ntau` is the
number of intercept times and `np` the number of curvatures or ray parameters.
* `nt::Int`: number of time samples in the data domain.

# Keyword arguments
* `order::String="parab"`: `"parab"` for parabolic transform, `"linear"`
for linear transform.
* `dt::Real=0.004`: sampling interval in seconds.
* `h::Vector{Real}=collect(0.0:20.0:1000.0)`: offset vector; `h[1:nh]`.
* `href::Real=0.0`: reference offset for parabolic Radon Transform. If the
defautl value `href=0.0` is given, `href` is set to `max(abs(h))`.
* `p::Vector{Real}=collect(-0.05:0.01:2.2)`: `p[1:np]`. If `order="parab"`, `p`
is a vector of residual moveout ("curvatures") at reference offset `href` in
seconds; if `order=linear`, `p` is a vector of ray parameters in s/m.
* `flow::Real=0.0`: minimum frequency in the data in Hz.
* `fhigh::Real=125.0`: maximum frequency in the data in Hz.

# Output
* `d`: data synthetized via forward Radon modeling, `d[1:nt, 1:nh]`.

# References
*  Hampson, D., 1986, Inverse velocity stacking for multiple elimination:
Canadian Journal of Exploration Geophysics, 22, 44-55.
*  Sacchi, M.D. and Ulrych, T.J., 1995, High-resolution velocity gathers and
offset space reconstruction: Geophysics, 60, 1169-1177.
"""

function SeisRadonFreqFor{Tm<:Real, Th<:Real, Tp<:Real
                          }(m::Array{Tm,2}, nt::Int; order::String="parab",
                            dt::Real=0.004, href::Real=0.0,
                            h::Vector{Th}=collect(0.0:20.0:1000.0),
                            p::Vector{Tp}=collect(-0.05:0.01:2.2),
                            flow::Real=0.0, fhigh::Real=125.0)
    if order=="parab"
        I = 2
        href == 0 && (href = maximum(abs(h)))
    elseif order=="linear"
        I = 1
        href = 1.0
    else
        error("Order should be equal to \"parab\" or \"linear\"")
    end
    ntau = size(m, 1)
    np = length(p)
    np == size(m, 2) || error("lenght(p) must be equal to size(m, 2)")
    nh = length(h)
    nw = nextpow2(ntau)
    m = cat(1, m, zeros(Tm, nw-ntau, np))
    M = fft(m, 1)
    iw_low = round(Int, flow*dt*nw+1)
    iw_high = round(Int, fhigh*dt*nw+1)
    D = zeros(Complex{Tm}, nw, nh)
    for iw = iw_low:iw_high
        w  = 2.0*pi*(iw-1)/(nw*dt)
        L = zeros(Complex{Tm}, nh, np)
        for ip = 1:np
            for ih = 1:nh
                phi = w*p[ip]*(h[ih]/href)^I
      	        L[ih, ip] = exp(-im*phi)
            end
        end
        D[iw, :] = L*M[iw,:]
    end
    for iw = round(Int, nw/2)+2:nw
        D[iw, :] = conj(D[nw-iw+2, :])
    end
    d = real(ifft(D, 1))
    d = d[1:nt, :]
    return d
end


text:
@testset "Estimators" begin
  # floating point tolerance
  tol = 10eps()

  # create some data
  dim = 3; nobs = 10
  X = rand(dim, nobs)
  z = rand(nobs)

  # some target location
  xₒ = rand(dim)

  γ = GaussianVariogram(sill=1., range=1., nugget=0.)
  simkrig = SimpleKriging(X, z, γ, mean(z))
  ordkrig = OrdinaryKriging(X, z, γ)
  unikrig = UniversalKriging(X, z, γ, 1)
  driftkrig = ExternalDriftKriging(X, z, γ, [x->1.])

  # Kriging is an interpolator
  for j=1:nobs
    SKestimate, SKvar = predict(simkrig, X[:,j])
    OKestimate, OKvar = predict(ordkrig, X[:,j])
    UKestimate, UKvar = predict(unikrig, X[:,j])
    DKestimate, DKvar = predict(driftkrig, X[:,j])

    # estimate checks
    @test SKestimate ≈ z[j]
    @test OKestimate ≈ z[j]
    @test UKestimate ≈ z[j]
    @test DKestimate ≈ z[j]

    # variance checks
    @test SKvar + tol ≥ 0
    @test OKvar + tol ≥ 0
    @test UKvar + tol ≥ 0
    @test DKvar + tol ≥ 0
    @test SKvar ≤ OKvar + tol
  end

  # save results on a particular location xₒ
  SKestimate, SKvar = predict(simkrig, xₒ)
  OKestimate, OKvar = predict(ordkrig, xₒ)
  UKestimate, UKvar = predict(unikrig, xₒ)
  DKestimate, DKvar = predict(driftkrig, xₒ)

  # Kriging is translation-invariant
  h = rand(dim)
  simkrig_h = SimpleKriging(X .+ h, z, γ, mean(z))
  ordkrig_h = OrdinaryKriging(X .+ h, z, γ)
  unikrig_h = UniversalKriging(X .+ h, z, γ, 1)
  driftkrig_h = ExternalDriftKriging(X .+ h, z, γ, [x->1.])
  SKestimate_h, SKvar_h = predict(simkrig_h, xₒ .+ h)
  OKestimate_h, OKvar_h = predict(ordkrig_h, xₒ .+ h)
  UKestimate_h, UKvar_h = predict(unikrig_h, xₒ .+ h)
  DKestimate_h, DKvar_h = predict(driftkrig_h, xₒ .+ h)
  @test SKestimate_h ≈ SKestimate
  @test SKvar_h ≈ SKvar
  @test OKestimate_h ≈ OKestimate
  @test OKvar_h ≈ OKvar
  @test UKestimate_h ≈ UKestimate
  @test UKvar_h ≈ UKvar
  @test DKestimate_h ≈ DKestimate
  @test DKvar_h ≈ DKvar

  # Kriging estimate is invariant under covariance scaling
  # Kriging variance is multiplied by the same factor
  α = rand()
  γ_α = GaussianVariogram(sill=α, range=1., nugget=0.)
  simkrig_α = SimpleKriging(X, z, γ_α, mean(z))
  ordkrig_α = OrdinaryKriging(X, z, γ_α)
  unikrig_α = UniversalKriging(X, z, γ_α, 1)
  driftkrig_α = ExternalDriftKriging(X, z, γ_α, [x->1.])
  SKestimate_α, SKvar_α = predict(simkrig_α, xₒ)
  OKestimate_α, OKvar_α = predict(ordkrig_α, xₒ)
  UKestimate_α, UKvar_α = predict(unikrig_α, xₒ)
  DKestimate_α, DKvar_α = predict(driftkrig_α, xₒ)
  @test SKestimate_α ≈ SKestimate
  @test SKvar_α ≈ α*SKvar
  @test OKestimate_α ≈ OKestimate
  @test OKvar_α ≈ α*OKvar
  @test UKestimate_α ≈ UKestimate
  @test UKvar_α ≈ α*UKvar
  @test DKestimate_α ≈ DKestimate
  @test DKvar_α ≈ α*DKvar

  # Kriging variance is a function of data configuration, not data values
  δ = rand(nobs)
  z_δ = z .+ δ
  simkrig_δ = SimpleKriging(X, z_δ, γ, mean(z_δ))
  ordkrig_δ = OrdinaryKriging(X, z_δ, γ)
  unikrig_δ = UniversalKriging(X, z_δ, γ, 1)
  driftkrig_δ = ExternalDriftKriging(X, z_δ, γ, [x->1.])
  SKestimate_δ, SKvar_δ = predict(simkrig_δ, xₒ)
  OKestimate_δ, OKvar_δ = predict(ordkrig_δ, xₒ)
  UKestimate_δ, UKvar_δ = predict(unikrig_δ, xₒ)
  DKestimate_δ, DKvar_δ = predict(driftkrig_δ, xₒ)
  @test SKvar_δ ≈ SKvar
  @test OKvar_δ ≈ OKvar
  @test UKvar_δ ≈ UKvar
  @test DKvar_δ ≈ DKvar

  # Ordinary Kriging ≡ Universal Kriging with 0th degree drift
  unikrig_0th = UniversalKriging(X, z, γ, 0)
  OKestimate, OKvar = predict(ordkrig, xₒ)
  UKestimate, UKvar = predict(unikrig_0th, xₒ)
  @test OKestimate ≈ UKestimate
  @test OKvar ≈ UKvar

  # Ordinary Kriging ≡ Kriging with constant external drift
  driftkrig_const = ExternalDriftKriging(X, z, γ, [x->1.])
  OKestimate, OKvar = predict(ordkrig, xₒ)
  DKestimate, DKvar = predict(driftkrig_const, xₒ)
  @test OKestimate ≈ DKestimate
  @test OKvar ≈ DKvar

  # Non-stationary variograms are allowed
  γ_ns = PowerVariogram()
  ordkrig_ns = OrdinaryKriging(X, z, γ_ns)
  unikrig_ns = UniversalKriging(X, z, γ_ns, 1)
  driftkrig_ns = ExternalDriftKriging(X, z, γ_ns, [x->1.])
  for j=1:nobs
    OKestimate, OKvar = predict(ordkrig_ns, X[:,j])
    UKestimate, UKvar = predict(unikrig_ns, X[:,j])
    DKestimate, DKvar = predict(driftkrig_ns, X[:,j])

    # estimate checks
    @test OKestimate ≈ z[j]
    @test UKestimate ≈ z[j]
    @test DKestimate ≈ z[j]

    # variance checks
    @test OKvar + tol ≥ 0
    @test UKvar + tol ≥ 0
    @test DKvar + tol ≥ 0
  end

  # Floating point precision checks
  X_f  = rand(Float32, dim, nobs)
  z_f  = rand(Float32, nobs)
  xₒ_f = rand(Float32, dim)
  γ_f = GaussianVariogram(sill=1f0, range=1f0, nugget=0f0)
  X_d  = Float64.(X_f)
  z_d  = Float64.(z_f)
  xₒ_d = Float64.(xₒ_f)
  γ_d = GaussianVariogram(sill=1., range=1., nugget=0.)
  simkrig_f = SimpleKriging(X_f, z_f, γ_f, mean(z_f))
  ordkrig_f = OrdinaryKriging(X_f, z_f, γ_f)
  unikrig_f = UniversalKriging(X_f, z_f, γ_f, 1)
  driftkrig_f = ExternalDriftKriging(X_f, z_f, γ_f, [x->1f0])
  simkrig_d = SimpleKriging(X_d, z_d, γ_d, mean(z_d))
  ordkrig_d = OrdinaryKriging(X_d, z_d, γ_d)
  unikrig_d = UniversalKriging(X_d, z_d, γ_d, 1)
  driftkrig_d = ExternalDriftKriging(X_d, z_d, γ_d, [x->1.])
  SKestimate_f, SKvar_f = predict(simkrig_f, xₒ_f)
  OKestimate_f, OKvar_f = predict(ordkrig_f, xₒ_f)
  UKestimate_f, UKvar_f = predict(unikrig_f, xₒ_f)
  DKestimate_f, DKvar_f = predict(driftkrig_f, xₒ_f)
  SKestimate_d, SKvar_d = predict(simkrig_d, xₒ_d)
  OKestimate_d, OKvar_d = predict(ordkrig_d, xₒ_d)
  UKestimate_d, UKvar_d = predict(unikrig_d, xₒ_d)
  DKestimate_d, DKvar_d = predict(driftkrig_d, xₒ_d)
  @test isapprox(SKestimate_f, SKestimate_d, atol=1e-4)
  @test isapprox(SKvar_f, SKvar_d, atol=1e-4)
  @test isapprox(OKestimate_f, OKestimate_d, atol=1e-4)
  @test isapprox(OKvar_f, OKvar_d, atol=1e-4)
  @test isapprox(UKestimate_f, UKestimate_d, atol=1e-4)
  @test isapprox(UKvar_f, UKvar_d, atol=1e-4)
  @test isapprox(DKestimate_f, DKestimate_d, atol=1e-4)
  @test isapprox(DKvar_f, DKvar_d, atol=1e-4)
end


text:
"""
`openit(path::AbstractString)`
open the path of the file/folder.
"""
function openit(path::AbstractString)
    if isempty(path)
        path = pwd();
    end

    if !ispath(path)
        @warn("You have entered an invalid file system entity.");
        return nothing
    end

    path = abspath(path); # this is required otherwise in windows the cmd command to run won't work properly.

    if Sys.iswindows()
        # if isfile(path)
        #     run(`cmd /c start $path`, wait=false);
        # else
            run(`cmd /c start %windir%\\explorer.exe $path`, wait=false);
        # end
    elseif Sys.islinux()
        run(xdg-open $path, wait=false);
    else
        @warn("Couldn't open $path");
    end
    return nothing
end


text:
module FFT
using FFTW
using ..Stats

export fft

function fft_frequency_label(index, min_index_2, min_index_4, min_index_6, min_index_8, min_index_10, min_index_12, length)
  if index == min_index_2
    "1/2"
  elseif index == min_index_4
    "1/4"
  elseif index == min_index_6
    "1/6"
  elseif index == min_index_8
    "1/8"
  elseif index == min_index_10
    "1/10"
  elseif index == min_index_12
    "1/12"
  elseif index == 1
    "1/$length"
  else
    ""
  end
end

function fft_frequency_labels(array)
  min_2, index_2 = findmin((x -> abs(x - 1/2)).(array))
  min_4, index_4 = findmin((x -> abs(x - 1/4)).(array))
  min_6, index_6 = findmin((x -> abs(x - 1/6)).(array))
  min_8, index_8 = findmin((x -> abs(x - 1/8)).(array))
  min_10, index_10 = findmin((x -> abs(x - 1/10)).(array))
  min_12, index_12 = findmin((x -> abs(x - 1/12)).(array))
  min, index = findmin((x -> abs(x - 1/size(array, 1))).(array))

  fft_frequency_label.(1:size(array, 1), index_2, index_4, index_6, index_8, index_10, index_12, size(array, 1))
end

function fft(values)
  fft_arr = (FFTW.fft(values, 1))
  length = convert(Int64, ceil(size(fft_arr, 1) / 2))
  fft_results = abs.(fft_arr[1:length])
  frequencies = fft_frequency_labels(FFTW.fftfreq(size(values, 1), 1))[1:length]
  vcat(hcat(frequencies...), hcat((fft_results ./ Stats.max(fft_results))...))
end
end


text:
using Pkg

const project = Pkg.project()

const modulefile = joinpath(dirname(project.path), "src", "$(project.name).jl")

if project.ispackage && !isfile(modulefile)
    @info "Project is a package, but has no src/\$JULIA_PROJECT.jl, creating an empty one."
    mkpath(dirname(modulefile))
    open(modulefile, "w") do io
        println(io, "module $(project.name)")
        println(io, "end")
    end
end

# Skip generating a system image when there are no dependencies,
# or all deps are stdlibs.
if isempty(project.dependencies) || values(project.dependencies) ⊆ keys(Pkg.Types.stdlibs())
    exit(0)
end

Pkg.add(name="PackageCompiler", version="1")

using PackageCompiler, UUIDs

function main()
    packages = [Symbol(k) for k in keys(project.dependencies)]
    # exclude some packages that can make sysimage creation or usage fail.
    exclude = [:PackageCompiler]
    println("excluding $exclude")
    packages = setdiff(packages, exclude)
    println("Creating sys image with deps being tracked from a registry:")
    println(packages)
    create_sysimage(packages; replace_default=true, cpu_target = "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)")
end

main()
Pkg.rm("PackageCompiler")


text:
#==============================================================================#
# mime.jl
#
# Simple MIME Multipart encoder.
#
# Copyright OC Technology Pty Ltd 2015 - All rights reserved
#==============================================================================#

const Part = Tuple{String,String,String}

"""
    mime_multipart([header,] parts)

Encode `parts` as a
[MIME Multipart](https://www.w3.org/Protocols/rfc1341/7_2_Multipart.html)
message.

`parts` is a Vector of `(filename, content_type, content)` Tuples.
"""
mime_multipart(parts::Array) = mime_multipart("", parts::Vector{Part})

function mime_multipart(header::AbstractString, parts::Vector{Part})

    boundary = "=PRZLn8Nm1I82df0Dtj4ZvJi="

    mime =
    """
    MIME-Version: 1.0
    Content-Type: multipart/mixed; boundary="$boundary"
    $header

    """

    for (filename, content_type, content) in parts

        mime *= "--$boundary\n"

        if filename != ""
            mime *= "Content-Type: $content_type;\n    name=$filename\n"
            mime *= "Content-Disposition: attachment;\n    filename=$filename\n"
        else
            mime *= "Content-Type: $content_type\n"
        end

        if isa(content, AbstractString)
            mime *= "Content-Transfer-Encoding: binary\n"
        else
            mime *= "Content-Transfer-Encoding: base64\n"
            b64 = base64encode(content)
            b64 = [b64[i:min(length(b64),i+75)] for i in 1:76:length(b64)]
            content = join(b64, "\n")
        end

        mime *= "\n$content\n"
    end

    mime *= "--$boundary--\n"
    return mime
end

#==============================================================================#
# End of file
#==============================================================================#


text:
using GNSSTools
using Statistics
using Random
using StatsBase
using LinearAlgebra
using ProgressMeter
using JLD
using PyPlot
pygui(true)

struct ConstellationHistogram{T1,T2}
    nbins::Int
    dopplers::Array{Float64,1}
    doppler_rates::Array{Float64,1}
    histogram::T1
    weights::T2
    fd_rate::Vector{Float64}
    f_d::Vector{Float64}
end

function initiate_hist(dopplers, doppler_rates; nbins=1000)
    histogram = normalize(fit(Histogram, (doppler_rates, dopplers),
                              nbins=nbins), mode=:probability)
    weights = collect(Iterators.flatten(histogram.weights))
    weights = Weights(weights)
    fd_rate, f_d = meshgrid(Array(histogram.edges[1][1:end-1]),
                            Array(histogram.edges[2][1:end-1]))
    fd_rate = collect(Iterators.flatten(fd_rate))
    f_d = collect(Iterators.flatten(f_d))
    return ConstellationHistogram(nbins, dopplers, doppler_rates,
                                  histogram, weights, fd_rate, f_d)
end

function sample_distribution(histogram::ConstellationHistogram)
    f_d = sample(histogram.f_d, histogram.weights)
    fd_rate = sample(histogram.fd_rate, histogram.weights)
    return (f_d, fd_rate)
end

# "doppler_gps",
# "doppler_iridium",
# "doppler_starlink",
# "doppler_oneweb",
# "doppler_rate_gps",
# "doppler_rate_iridium",
# "doppler_rate_starlink",
# "doppler_rate_oneweb",
# "doppler_bounds_gps",
# "doppler_bounds_iridium",
# "doppler_bounds_starlink",
# "doppler_bounds_oneweb",
# "doppler_raw_gps",
# "doppler_raw_iridium",
# "doppler_raw_starlink",
# "doppler_raw_oneweb",
# "doppler_rate_raw_gps",
# "doppler_rate_raw_iridium",
# "doppler_rate_raw_starlink",
# "doppler_rate_raw_oneweb",
# "freq",
# "user_lla",
# "t_range",
# "min_elevation")

function benchmark_tracking(sigtype::GNSSTools.SignalType,
                            histogram::ConstellationHistogram,
                            Ts, CN0s, fd_rates, f_s; iterations=100,
                            include_phase_noise=true,
                            include_thermal_noise=true,
                            include_adc=true,
                            include_carrier=true,
                            h_parms=[1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2,
                            code_start_idx=missing, phi=missing,
                            nADC=8, fd_range=5000, f_if=1.25e6, prn=1, M=20,
                            fd_center=0, Pd_threshold=0.9, p_fa=1e-7,
                            acquisition_T=1e-3, fine_acq_T=10e-3,
                            fine_acq_method=:fft, use_fine_acq=true,
                            h₀=h_pamrs[3], h₋₂=h_parms[1], dll_b=0.1,
                            state_num=2, q_a=1, t_length=1, q_mult=1,
                            cov_mult=1, R_mult=1, dynamickf=true)
    if sigtype.include_I && sigtype.include_Q
        B = max(sigtype.B_I, sigtype.B_Q)
    elseif sigtype.include_I && !sigtype.include_Q
        B = sigtype.B_I
    elseif !sigtype.include_I && sigtype.include_Q
        B = sigtype.B_Q
    else
        error("No CodeType defined in SigType struct.")
    end
    signal = definesignal(sigtype, f_s, t_length;
                          prn=prn, f_if=f_if, nADC=nADC,
                          include_phase_noise=include_phase_noise,
                          include_thermal_noise=include_thermal_noise,
                          include_adc=include_adc,
                          include_carrier=include_carrier,
                          receiver_h_parms=h_parms)
    Pds = Array{Float64}(undef, length(Ts), length(CN0s))
    p = Progress(length(Ts)*length(CN0s)*iterations, 1, "Processing...")
    for i in 1:length(Ts)
        tracking_T = Ts[i]
        N_over_2 = floor(Int, t_length/tracking_T/2)
        for j in 1:length(CN0s)
            CN0 = CN0s[j]

            for iteration in 1:iterations
                if ismissing(phi)
                    ϕ₀_truth = rand(0:0.0001:2π)
                else
                    ϕ₀_truth = phi
                end
                if ismissing(code_start_idx)
                    n0_truth = rand(1:1:replica.sample_num)
                else
                    n0_truth = code_start_idx
                end
                if include_phase_noise
                    new_phase_noise = true
                else
                    new_phase_noise = false
                end
                if include_thermal_noise
                    new_thermal_noise = true
                else
                    new_thermal_noise = false
                end
                f_d, fd_rate = sample_distribution(histogram)
                definesignal!(signal; CN0=CN0, phi=ϕ₀_truth, f_d=f_d,
                              code_start_idx=n0_truth, fd_rate=fd_rate,
                              new_phase_noise=new_phase_noise,
                              new_thermal_noise=new_thermal_noise)
                generatesignal!(signal)
                acqresults, trackresults, corr_result, SNR_est,
                P_d, above_threshold = process(signal,
                                               sigtype,
                                               prn, "I"; fd_range=fd_range,
                                               fine_acq_method=fine_acq_method,
                                               return_corrresult=true,
                                               return_Pd=true,
                                               h₀=h₀,
                                               h₋₂=h₋₂,
                                               q_mult=q_mult,
                                               cov_mult=cov_mult,
                                               R_mult=R_mult,
                                               q_a=q_a,
                                               dynamickf=dynamickf,
                                               state_num=state_num,
                                               dll_b=dll_b,
                                               acquisition_T=acquisition_T,
                                               fine_acq_T=fine_acq_T,
                                               tracking_T=tracking_T,
                                               show_plot=false,
                                               M=M);
                if P_d > Pd_threshold
                    dphi_meas_std = std(trackresults.dphi_meas[N_over_2:end])
                    code_err_meas_std = std(trackresults.code_err_meas[N_over_2:end])
                    P = mean.(trackresults.P[N_over_2:end])
                end
                next!(p)
            end
        end
    end
end

file_name = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/ch1_doppler_curves.jld"
doppler_raw_gps, doppler_rate_raw_gps, sig_freq = load(file_name,
                                                       "doppler_raw_gps",
                                                       "doppler_rate_raw_gps",
                                                       "freq")

gps_hist = initiate_hist(doppler_raw_gps, doppler_rate_raw_gps; nbins=1000)
# f_d, fd_rate = sample_distribution(gps_hist)
f_s = 2.5e6  # Hz
Ts = [1, 2, 5, 10] .* 1e-3  # seconds
CN0s = [25, 30, 35, 40, 45, 50]  # dB
fd_rates_gps = unique(gps_hist.fd_rate)[1:30:end]  # Hz/second
sigtype = define_l1ca_code_type(maximum(Ts))


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, ColorSchemes
let
    xs = 0:0.01:1
    p = -10:1:10
    p= filter(x->x != 0,collect(p))
    ys = zeros(length(xs), length(p))
    for (indx,i) in enumerate(p)
        if i <= -1
            ys[:,indx]= xs .^(1/abs(i))
        elseif i>= 1
            ys[:,indx]= xs .^i
        end
    end
    cbarPal = :thermal
    cmap = get(colorschemes[cbarPal], LinRange(0,1,length(p)))

    fig = Figure(resolution = (700, 450), font =:sans, fontsize = 18)
    ax = Axis(fig, aspect = 1, xlabel = "x", ylabel = "xᵖ")
    lins = [lines!(xs, ys[:,v], color = cmap[v]) for v in 1:length(p)]
    leg = Legend(fig, lins, string.(p), "p", nbanks = 2, labelsize = 12,
     valign = :center)
    fig[1, 1] = ax
    fig[1, 2] = leg
    fig
    #save("MulticoloredLinesLegOut.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "MulticoloredLinesLegOut.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "ColorSchemes"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
println("Create README.md")

import FrameFun
using Conda: PYTHONDIR
jupyter = PYTHONDIR * "/jupyter"
FRAMEFUNSRC = pathof(FrameFun)
FRAMEFUNPATH = splitdir(splitdir(FRAMEFUNSRC)[1])[1]
run(`$jupyter nbconvert --execute --to markdown --output $FRAMEFUNPATH/README.md $FRAMEFUNPATH/README.ipynb --ExecutePreprocessor.timeout=180`)


text:
function digit_product(s::String)
	ans=1

	for c in s
		ans *= c-'0'

	end
	return ans
end

function compute()
	ans=0
	for i=1:length(NUMBER)-ADJACENT
		ans=max(digit_product(NUMBER[i:1:i+ADJACENT-1]) ,ans)
	end

	return ans

end

NUMBER = "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"
ADJACENT = 13

println(compute())


text:
module LearningToRank

using CatBoost
using DataFrames

train, test = load_dataset(:msrank_10k)
x_train = train.drop([0, 1]; axis=1).values
y_train = train[1].values
queries_train = train[2].values

x_test = test.drop([0, 1]; axis=1).values
y_test = test[1].values
queries_test = test[2].values

# Important dims.
num_documents, num_features = size(x_train)
num_queries = size(unique(queries_train))[1]
println("Data dims: $((num_documents, num_features))")
println("Num queries: $(num_queries)")

# Get relevance labels.
max_relevance = maximum(y_train)
y_train /= max_relevance
y_test /= max_relevance

# Create pools.
train = Pool(; data=x_train, label=y_train, group_id=queries_train)
test = Pool(; data=x_test, label=y_test, group_id=queries_test)

# small number of iterations to not slow down CI too much
default_parameters = Dict("iterations" => 10, "loss_function" => "RMSE",
                          "custom_metric" => ["MAP:top=10", "PrecisionAt:top=10",
                                              "RecallAt:top=10"], "verbose" => false,
                          "random_seed" => 314159)

function fit_model(params, train_pool, test_pool)
    model = catboost.CatBoost(params)
    model.fit(train_pool; eval_set=test_pool, plot=false)
    return model
end

println("Training with:")
display(default_parameters)
model = fit_model(default_parameters, train, test)

end # module


text:
function make_gaussian(grid::Grid)
    u0 = zeros(grid.ny, grid.nx)
    x = grid.xRange
    y = grid.yRange

    x₀ = 0.0
    y₀ = 0.0
    var₀ = 0.9
    T₀ = 1.0
    for j = 1 : grid.ny, i = 1 : grid.nx
        u0[j,i] = T₀*exp(-((x[i]-x₀)^2 + (y[j]-y₀)^2)/(2*var₀))
    end
    return u0
end
function solution_gaussian(grid::Grid, t, κ)
    uana = zeros(grid.ny, grid.nx)
    x = grid.xRange
    y = grid.yRange

    x₀ = 0.0
    y₀ = 0.0
    var₀ = 0.9
    T₀ = 1.0
    #t₀ = -0.1

    for j = 1 : grid.ny, i = 1 : grid.nx
        uana[j,i] = T₀/sqrt(1+4*κ*t/var₀)*
                exp(-((x[i]-x₀)^2 + (y[j]-y₀)^2)/(2*var₀+4*κ*t))
    end
    return uana
end


text:
__precompile__

mutable struct Settings
    # grid settings
    # number spatial interfaces
    Nx::Int64;
    # start and end point
    a::Float64;
    b::Float64;
    # grid cell width
    dx::Float64

    # time settings
    # end time
    tEnd::Float64;
    # time increment
    dt::Float64;
    # cfl number
    cfl::Float64;

    # number of quadrature points
    Nq::Int64;
    # definition of quadrature type
    quadratureType::String;

    # maximal polynomial degree
    N::Int64;

    # DLR rank
    r::Int64

    # solution values
    uL::Float64;
    uR::Float64;

    # initial shock position at xi = 0
    x0::Float64;

    # grid
    x

    sigma::Float64;

    # spatial limiter
    limiterType::String;
    # time update type
    rkType::String;
    # number of RK stages
    rkStages::Int;
    # filter
    filterType::String;
    lambda::Float64;
    filterOrder::Int64;

    useStabilizingTermsS::Bool;
    useStabilizingTermsL::Bool;
    stabilization::Int;

    # conservation settings
    NCons::Int; # number of conserved basis functions
    iCons; #::Array{Int64,2}; # indices of conserved basis funtions

    # initial condition
    IC::Function;
    solutionExact::Function;

    function Settings(Nx=600,N=10,r=9)
        # define spatial grid
        a = 0.0;
        b = 1.0;
        x = range(a,stop = b,length = Nx)
        dx = (b-a)/(Nx-1.0);

        # define time settings
        tEnd = 0.01;#0.01;#0.01115;
        cfl = 0.9;#0.5;

        # define test case settings
        x0 = 0.3; # 0.5
        x1 = 0.6;
        uL = 12.0;
        uM = 6.0;
        uR = 1.0;
        sigma0 = 0.2;
        sigma1 = 5.0;#5.0;

        quadratureType = "Gauss"; # Possibilities are Gauss and ClenshawCurtis

        # compute time step size
        dt = cfl*dx/uL;

        # determine quadrature points
        Nq = ceil(1.5*N+1);

        # filter parameters
        filterType = "L2" # L2, EXP
        lambda = 0.0;#0.00001
        filterOrder = 1;

        # limiter type
        limiterType = "None";#"Minmod"

        # Runge-Kutta type Euler, Heun, SSP
        rkType = "Euler"
        rkStages = 1;

        useStabilizingTermsS = true;
        useStabilizingTermsL = true;

        # stabilization method for PS: 0 - standard, 1 - Lax Wendroff in L and S, 2 - stabilization in all steps
        stabilization = 0;

        # conservation settings
        NCons = 2;#3; # number of conserved basis functions
        iCons = [1 1; 1 2];#[1 1;1 2;1 3]; # indices of conserved basis funtions

        #NCons = 0;#3; # number of conserved basis functions
        #iCons = 0;

        # build class
        new(Nx,a,b,dx,tEnd,dt,cfl,Nq,quadratureType,N,r,uL,uR,x0,x,sigma0,limiterType,rkType,rkStages,filterType,lambda,filterOrder,useStabilizingTermsS,useStabilizingTermsL,stabilization,NCons,iCons,
            #(xV,xi,eta)->IC3(xV,xi,eta,sigma0,sigma1,uL,uM,uR,x0,x1),
            #(t,xV,xi,eta)->IC3Exact(t,xV,xi,eta,sigma0,sigma1,uL,uM,uR,x0,x1))
            (xV,xi,eta)->IC5(xV,xi,eta,sigma0,sigma1,uL,uR,x0),
            (t,xV,xi,eta)->IC5Exact(t,xV,xi,eta,sigma0,sigma1,uL,uR,x0))
    end

end

function IC1(x,xi,sigma,uL,uR,x0,x1)
    y = zeros(size(xi));
    for j = 1:length(y);
        if x < x0+sigma*xi[j]
            y[j] = uL;
        elseif x < x1+sigma*xi[j]
            y[j] = uL + (uR - uL)*(x-sigma*xi[j]-x0)/(x1-x0);
        else
            y[j] = uR;
        end
    end
    return y;
end

function IC1Exact(t::Float64,x,xi::Float64,sigma::Float64,uL::Float64,uR::Float64,x0::Float64,x1::Float64)
    y = zeros(length(x));

    if t >= (x1-x0)/(uL-uR);
        tS = (x1-x0)/(uL-uR);
        x0BeforeShock = x0+sigma*xi + tS*uL;
        x1BeforeShock = x1+sigma*xi + tS*uR;
        x0 = x0BeforeShock + (t-tS)*(uL+uR)*0.5;
        x1 = x0 - 1.0;
    else
        x0 = x0+sigma*xi + t*uL;
        x1 = x1+sigma*xi + t*uR;
    end

    for j = 1:length(y);
        if x[j] < x0
            y[j] = uL;
        elseif x[j] < x1
            y[j] = uL + (uR - uL)*(x[j]-x0)/(x1-x0);
        else
            y[j] = uR;
        end
    end

    return y;
end

function IC2(x,xi,sigma,uL,uR,x0,x1)
    y = zeros(size(xi));
    K0 = 12;
    K1 = 1;
    div = 1/(x0^3+3*x0*x1^2-x1^3-3*x1*x0^2);
    a = -2*(K0-K1)*div;
    b = 3*(K0-K1)*(x0+x1)*div;
    c = -6*(K0-K1)*x0*x1*div;
    d = (-K0*x1^3+3*K0*x0*x1^2+K1*x0^3-3*K1*x1*x0^2)*div;

    for k = 1:length(xi)
        if x-sigma*xi[k] < x0
            y[k] = K0;
        elseif x-sigma*xi[k] < x1
            y[k] = a*(x-sigma*xi[k])^3+b*(x-sigma*xi[k])^2+c*(x-sigma*xi[k])+d;
        else
            y[k] = K1;
        end
    end
    return y;
end

function IC2Exact(t::Float64,x::Float64,xi,sigma,uL,uR,x0,x1,problem,advectionSpeed)
    y = zeros(size(xi));
    K0 = 12;
    K1 = 1;
    div = 1/(x0^3+3*x0*x1^2-x1^3-3*x1*x0^2);
    a = -2*(K0-K1)*div;
    b = 3*(K0-K1)*(x0+x1)*div;
    c = -6*(K0-K1)*x0*x1*div;
    d = (-K0*x1^3+3*K0*x0*x1^2+K1*x0^3-3*K1*x1*x0^2)*div;

    for k = 1:length(xi)
        if x-sigma*xi[k] < x0+t*advectionSpeed
            y[k] = K0;
        elseif x-sigma*xi[k] < x1+t*advectionSpeed
            y[k] = a*(x-sigma*xi[k]-t*advectionSpeed)^3+b*(x-sigma*xi[k]-t*advectionSpeed)^2+c*(x-sigma*xi[k]-t*advectionSpeed)+d;
        else
            y[k] = K1;
        end
    end
    return y;
end

function IC3Exact(t::Float64,x,xi,eta,sigma0,sigma1,uL,uM,uR,x0,x1)
    y = 0.0
    println(x);
    for j = 1:length(y);
        tStar = 2.0*(x0-x1)/(uR-uL-sigma0*xi[j]);
        if t < tStar
            v0 = 0.5*(uL+uM+sigma0*xi[j]+sigma1*eta[j]);
            v1 = 0.5*(uR+uM+sigma1*eta[j]);
            x0T = x0+v0*t;
            x1T = x1+v1*t
            if x < x0T
                y[j] = uL+ sigma0*xi[j];
            elseif x < x1T
                y[j] = uM+ sigma1*eta[j];
            else
                y[j] = uR;
            end
        end
        if t >= tStar
            xStar = (x0-x1)*(uM-sigma1*eta+uR)/(uR-uL-sigma0*xi)+x1;
            vStar = 0.5*(uL + sigma0*xi +uR);
            if x < xStar + (t-tStar)*vStar
                y[j] = uL + sigma0*xi[j];
            else
                y[j] = uR;
            end
        end
    end
    return y;
end

function IC3(x,xi,eta,sigma0,sigma1,uL,uM,uR,x0,x1)
    y = zeros(length(xi));
    for j = 1:length(y);
        if x < x0
            y[j] = uL + sigma0*xi[j];
        elseif x < x1
            y[j] = uM + sigma1*eta[j];
        else
            y[j] = uR;
        end
    end
    return y;
end

function IC4(x,xi,eta,sigma0,sigma1,uL,uM,uR,x0,x1)
    y = zeros(length(xi));
    for j = 1:length(y);
        if x < x0+ sigma0*xi[j]
            y[j] = uL + sigma1*eta[j];
        elseif x < x1
            y[j] = uM;;
        else
            y[j] = uR;
        end
    end
    return y;
end

function IC4Exact(t::Float64,x,xi::Float64,eta::Float64,sigma0::Float64,sigma1::Float64,uL::Float64,uM::Float64,uR::Float64,x0::Float64,x1::Float64)
    y = zeros(length(x));
    uLXi = uL+eta*sigma1;
    tS = (x0+sigma0*xi-x1)/(uM-uLXi);
    #println("tS = ",tS)
    #println("x0 = ",x0+sigma0*xi + tS*uLXi+ (t-tS)*(uLXi+uR)*0.5)

    if t >= tS;
        x0BeforeShock = x0+sigma0*xi + tS*uLXi;
        x0 = x0BeforeShock + (t-tS)*(uLXi+uR)*0.5;
        x1 = x0 - 1.0;
    else
        x0 = x0+sigma0*xi + t*uLXi;
        x1 = x1 + t*uM;
    end

    for j = 1:length(y);
        if x[j] < x0
            y[j] = uLXi;
        elseif x[j] < x1
            y[j] = uM;
        else
            y[j] = uR;
        end
    end

    return y;
end

function IC5(x,xi,eta,sigma0::Float64,sigma1::Float64,uL::Float64,uR::Float64,x0::Float64)
    y = zeros(length(xi));
    for j = 1:length(y);
        if x < x0+ sigma0*xi[j]
            y[j] = uL;
        else
            y[j] = uR+ sigma1*(eta[j]+1)*0.5;
        end
    end
    return y;
end

function IC5Exact(t::Float64,x,xi::Float64,eta::Float64,sigma0::Float64,sigma1::Float64,uL::Float64,uR::Float64,x0::Float64)

    y = zeros(length(x));
    uRXi = uR+ sigma1*(eta+1)*0.5;

    x0 = x0+sigma0*xi + t*(uL+uRXi)/2;

    for j = 1:length(y);
        if x[j] < x0
            y[j] = uL;
        else
            y[j] = uRXi;
        end
    end

    return y;
end


text:
"""Module for a Polya-Gamma Sampler"""
module PGSampler

using Distributions
const __TRUNC = 0.64;
const __TRUNC_RECIP = 1.0 / __TRUNC;
export PolyaGammaDist, draw
"""Sampler object"""
mutable struct PolyaGammaDist{T}
  # For sum of Gammas.
  trunc::Int64;
  bvec::Vector{T}
  # Helper.
  #Constructor
  function PolyaGammaDist{T}(trunc::Int) where {T<:Real}
    this = new{T}()
	set_trunc(this,trunc)
	return this
  end
end

function PolyaGammaDist(;trunc::Int=200)
	PolyaGammaDist{Float64}(trunc)
end

## Utility

function set_trunc(pg::PolyaGammaDist{T},trunc::Int64) where {T<:Real}

  if trunc < 1
    @warn "PolyaGamma(trunc::Int): trunc < 1. Set trunc=1."
    trunc = 1;
  end

  pg.trunc = trunc;
  pg.bvec = zeros(T,pg.trunc);

  for k in 1:pg.trunc
    d = k - 0.5;
    pg.bvec[k] = 4* pi^2 * d * d;
  end
end # set_trunc

function a(n::Int64, x::T) where {T<:Real}
  K = (n + 0.5) * pi;
  y = 0;
  if x > __TRUNC
    y = K * exp( -0.5 * K*K * x );
  elseif x > 0
    expnt = -1.5 * (log(0.5 * pi)  + log(x)) + log(K) - 2.0 * (n+0.5)*(n+0.5) / x;
    y = exp(expnt);
  end
  return y;
end

function mass_texpon(Z::T) where {T<:Real}
  t = __TRUNC;

  fz = 0.125 * pi*pi + 0.5 * Z^2;
  b = sqrt(1.0 / t) * (t * Z - 1);
  a = sqrt(1.0 / t) * (t * Z + 1) * -1.0;

  x0 = log(fz) + fz * t;
  xb = x0 - Z + logcdf(Distributions.Normal(),b);
  xa = x0 + Z + logcdf(Distributions.Normal(),a);

  qdivp = 4 / pi * ( exp(xb) + exp(xa) );

  return 1.0 / (1.0 + qdivp);
end # mass_texpon

function rtigauss(Z::T) where {T<:Real}
  Z = abs(Z);
  t = __TRUNC;
  X = t + 1.0;
  Unif = Uniform(); # between 0 and 1
  Normal = Distributions.Normal(0.0,1.0);
  if __TRUNC_RECIP > Z
    alpha = 0.0;
    rate = 1.0;
    Expo = Exponential(1.0/rate);
    while (rand(Unif) > alpha)
      E1 = rand(Expo); E2 = rand(Expo);
      while  E1*E1 > 2 * E2 / t
	       E1 = rand(Expo); E2 = rand(Expo);
      end
      X = 1 + E1 * t;
      X = t / (X * X);
      alpha = exp(-0.5 * Z*Z * X);
    end
  else
    mu = 1.0 / Z;
    while (X > t)
      Y = rand(Normal); Y = Y^2;
      half_mu = 0.5 * mu;
      mu_Y    = mu  * Y;
      X = mu + half_mu * mu_Y - half_mu * sqrt(4 * mu_Y + mu_Y * mu_Y);
      if rand(Unif) > mu / (mu + X)
	       X = mu*mu / X;
       end
    end
  end
  return X;
end

# ////////////////////////////////////////////////////////////////////////////////
# 				  // Sample //
# ////////////////////////////////////////////////////////////////////////////////

function draw(pg::PolyaGammaDist{T},n::Real, z::Real) where {T<:Real}
  if n == 0
    # @warn "draw(PolyaGamma): n < 1.  Set n = 1." maxlog=1
    return 0.0
  end
  sum = 0.0;
  for i in 1:n
    sum = sum + draw_like_devroye(z);
  end
  return sum;
end # draw

function draw_like_devroye(Z::T) where {T<:Real}
  # Change the parameter.
  Z = abs(Z) * 0.5;

  # Now sample 0.25 * J^*(1, Z := Z/2).
  fz = 0.125 * pi*pi + 0.5 * Z*Z;
  # ... Problems with large Z?  Try using q_over_p.
  # double p  = 0.5 * __PI * exp(-1.0 * fz * __TRUNC) / fz;
  # double q  = 2 * exp(-1.0 * Z) * pigauss(__TRUNC, Z);

  X = 0.0;
  S = 1.0;
  Y = 0.0;
  # int iter = 0; If you want to keep track of iterations.
  Unif = Uniform()
  Expo = Exponential(1.0/1.0)
  while (true)

    if rand(Unif) < mass_texpon(Z)
      X = __TRUNC + rand(Expo) / fz;
    else
      X = rtigauss(Z);
    end
    S = a(0, X);
    Y = rand(Unif) * S;
    n = 0;
    go = true;

    # Cap the number of iterations?
    while (go)
      n = n+1;
      if n%2==1
	       S = S - a(n, X);
	       if Y<=S; return 0.25 * X;end;
      else
	       S = S + a(n, X);
	       if Y>S; go = false; end;
      end

    end
    # Need Y <= S in event that Y = S, e.g. when X = 0.

  end
end # draw_like_devroye

end #module PGSampler


text:
# # [Custom objective terms](@id custom_terms)
#
# ProximalAlgorithms relies on the first-order primitives defined in [ProximalCore](https://github.com/JuliaFirstOrder/ProximalCore.jl).
# While a rich library of function types, implementing such primitives, is provided by [ProximalOperators](https://github.com/JuliaFirstOrder/ProximalOperators.jl),
# one may need to formulate problems using custom objective terms.
# When that is the case, one only needs to implement the right first-order primitive,
# ``\nabla f`` or ``\operatorname{prox}_{\gamma f}`` or both, for algorithms to be able
# to work with ``f``.
#
# Defining the proximal mapping for a custom function type requires adding a method for [`ProximalCore.prox!`](@ref).
#
# To compute gradients, ProximalAlgorithms provides a fallback definition for [`ProximalCore.gradient!`](@ref),
# relying on [Zygote](https://github.com/FluxML/Zygote.jl) to use automatic differentiation.
# Therefore, you can provide any (differentiable) Julia function wherever gradients need to be taken,
# and everything will work out of the box.
#
# If however one would like to provide their own gradient implementation (e.g. for efficiency reasons),
# they can simply implement a method for [`ProximalCore.gradient!`](@ref).
#
# ```@docs
# ProximalCore.prox
# ProximalCore.prox!
# ProximalCore.gradient
# ProximalCore.gradient!
# ```
#
# ## Example: constrained Rosenbrock
#
# Let's try to minimize the celebrated Rosenbrock function, but constrained to the unit norm ball. The cost function is

rosenbrock2D(x) = 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2

# To enforce the constraint, we define the indicator of the unit ball, together with its proximal mapping:
# this is simply projection onto the unit norm ball, so it is sufficient to normalize any given point that lies
# outside of the set.

using LinearAlgebra
using ProximalCore

struct IndUnitBall end

(::IndUnitBall)(x) = norm(x) > 1 ? eltype(x)(Inf) : eltype(x)(0)

function ProximalCore.prox!(y, ::IndUnitBall, x, gamma)
    if norm(x) > 1
        y .= x ./ norm(x)
    else
        y .= x
    end
    return zero(eltype(x))
end

# We can now minimize the function, for which we will use `PANOC`, which is a Newton-type method:

using ProximalAlgorithms

panoc = ProximalAlgorithms.PANOC()
solution, iterations = panoc(x0=-ones(2), f=rosenbrock2D, g=IndUnitBall())

# Plotting the solution against the cost function contour and constraint, gives an idea of its correctness.

using Plots

contour(-2:0.1:2, -2:0.1:2, (x,y) -> rosenbrock2D([x, y]), fill=true, framestyle=:none, background=nothing)
plot!(Shape(cos.(0:0.01:2*pi), sin.(0:0.01:2*pi)), opacity=.5, label="feasible set")
scatter!([solution[1]], [solution[2]], color=:red, markershape=:star5, label="computed solution")

# ## Example: counting operations
#
# It is often interesting to measure how many operations (gradient- or prox-evaluation) an algorithm is taking.
# In fact, in algorithms involving backtracking or some other line-search logic, the iteration count may not be entirely
# representative of the amount of operations are being performed; or maybe some specific implementations require
# additional operations to be performed when checking stopping conditions. All of this makes it difficult to quantify the exact
# iteration complexity.
#
# We can achieve this by wrapping functions in a dedicated `Counting` type:

mutable struct Counting{T}
    f::T
    gradient_count::Int
    prox_count::Int
end

Counting(f::T) where T = Counting{T}(f, 0, 0)

# Now we only need to intercept any call to `gradient!` and `prox!` and increase counters there:

function ProximalCore.gradient!(y, f::Counting, x)
    f.gradient_count += 1
    return ProximalCore.gradient!(y, f.f, x)
end

function ProximalCore.prox!(y, f::Counting, x, gamma)
    f.prox_count += 1
    return ProximalCore.prox!(y, f.f, x, gamma)
end

# We can run again the previous example, this time wrapping the objective terms within `Counting`:

f = Counting(rosenbrock2D)
g = Counting(IndUnitBall())

solution, iterations = panoc(x0=-ones(2), f=f, g=g)

# and check how many operations where actually performed:

println(f.gradient_count)
println(g.prox_count)


text:
# the maximum length in bytes of a varint32.
const global K_MAX32 = 5

"""
    parse32(buf::Vector{UInt8}, offset::UInt)

Reads and returns a varint32 from `buf`. Raises an error if no varint
can be parsed. Returns the value and the index just past the last
byte of the varint32.
"""
@inline function parse32(buf::Vector{UInt8}, offset::Integer)
    offset > length(buf) && @goto failure
    local b::UInt32 = convert(UInt32, buf[offset])
    local result::UInt32 = b & 0x7f

    (b < 0x80) && return (result, offset+1)
    (offset+=1) > length(buf) && @goto failure

    b = buf[offset]; result |= (b & 0x7f) << 7
    (b < 0x80) && return (result, offset+1)
    (offset+=1) > length(buf) && @goto failure

    b = buf[offset]; result |= (b & 0x7f) << 14
    (b < 0x80) && return (result, offset+1)
    (offset+=1) > length(buf) && @goto failure

    b = buf[offset]; result |= (b & 0x7f) << 21
    (b < 0x80) && return (result, offset+1)
    (offset+=1) > length(buf) && @goto failure

    b = buf[offset]; result |= (b & 0x7f) << 28
    (b < 0x10) && return (result, offset+1)

    @label failure
    error("Could not decode varint32.")
end

"""
    encode32!(buf::Vector{UInt8}, offset::Int, value::UInt32)

Encodes `value` into `buf` starting at `offset`. Returns the index
just past the last byte of the varint32.
"""
@inline function encode32!(buf::Vector{UInt8}, offset::Integer, value::UInt32)
    if value < (1 << 7)
        buf[offset] = value % UInt8
        elseif value < (1 << 14)
            buf[offset] = (value | 128) % UInt8
            buf[offset+=1] = (value >> 7) % UInt8
        elseif value < (1 << 21)
            buf[offset] = (value | 128) % UInt8
            buf[offset+=1] = (value >> 7 | 128) % UInt8
            buf[offset+=1] = (value >> 14) % UInt8
        elseif value < (1 << 28)
            buf[offset] = (value | 128) % UInt8
            buf[offset+=1] = (value >> 7 | 128) % UInt8
            buf[offset+=1] = (value >> 14 | 128) % UInt8
            buf[offset+=1] = (value >> 21) % UInt8
        else
            buf[offset] = (value | 128) % UInt8
            buf[offset+=1] = (value >> 7 | 128) % UInt8
            buf[offset+=1] = (value >> 14 | 128) % UInt8
            buf[offset+=1] = (value >> 21 | 128) % UInt8
            buf[offset+=1] = (value >> 28) % UInt8
        end
    return offset + 1
end


text:
module AdvancedHMC

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_AHMC", "0")))

using Statistics: mean, var, middle
using LinearAlgebra: Symmetric, UpperTriangular, mul!, ldiv!, dot, I, diag, cholesky
using StatsFuns: logaddexp
using LazyArrays: BroadcastArray
using Random: GLOBAL_RNG, AbstractRNG
using ProgressMeter
using Parameters: @unpack, reconstruct
using ArgCheck: @argcheck

import StatsBase: sample

# Notations
# ℓπ: log density of the target distribution
# ∂ℓπ∂θ: gradient of the log density of the target distribution
# θ: position variables / model parameters
# r: momentum variables

include("adaptation/Adaptation.jl")
export UnitEuclideanMetric, DiagEuclideanMetric, DenseEuclideanMetric
export NesterovDualAveraging, Preconditioner, NaiveHMCAdaptor, StanHMCAdaptor
using .Adaptation
import .Adaptation: adapt!, NesterovDualAveraging

include("hamiltonian.jl")
export Hamiltonian
include("integrator.jl")
export Leapfrog, JitteredLeapfrog, TemperedLeapfrog
include("trajectory.jl")
export StaticTrajectory, find_good_eps, HMCDA, NUTS, SliceTS, MultinomialTS, ClassicNoUTurn, GeneralisedNoUTurn

include("diagnosis.jl")
include("sampler.jl")
export sample

end # module


text:
"""
    symbol, idxnext = tryparseelement(s, idx, elements)

Extract the element symbol starting at position `idx` in `s`. It must be one of the items in `elements`.

`symbol` is the element symbol (a string), or `nothing` if no valid element can be found.
`idxnext` is the first position in `s` after the parsed element.
"""
function tryparseelement(s::AbstractString, idx::Int, elements)
    len = length(s)
    idx > len && return nothing, idx
    c = s[idx]
    isletter(c) || return nothing, idx
    idxnext = nextind(s, idx)
    cnext = idxnext <= len ? s[idxnext] : '1'
    if islowercase(cnext)
        # Try a 2-letter element
        element = c*cnext
        element ∈ elements && return element, nextind(s, idxnext)
    end
    # Try a 1-letter element
    element = string(c)
    element ∈ elements && return element, idxnext
    return nothing, idx
end

"""
    i, idxnext = tryparseint16(s, idx)

Parse an `Int16` starting at position `idx` in `s`.

`i` is the `Int16` value, or `nothing` if `s` does not contain an integer starting at `idx`.
`idxnext` is the first position in `s` after the parsed integer.
"""
function tryparseint16(s::AbstractString, idx::Int)
    len = length(s)
    lastidx = prevind(s, idx)
    while lastidx < len
        i = nextind(s, lastidx)
        isdigit(s[i]) || break
        lastidx = i
    end
    return lastidx >= idx ? (parse(Int16, SubString(s, idx, lastidx)), nextind(s, lastidx)) : (nothing, idx)
end

"""
    z, idxnext = tryparseint16(s, idx)

Parse charge state (as an `Int`) starting at position `idx` in `s`.

`z` is the charge, or `nothing` if `s` does not contain a charge starting at `idx`.
`idxnext` is the first position in `s` after the parsed charge.
"""
function tryparsecharge(s::AbstractString, idx::Int)
    len = length(s)
    idx > len && return nothing, idx
    c = s[idx]
    ispm(c) || return nothing, idx
    sgn = c == '+' ? 1 : -1
    idxnext = nextind(s, idx)
    count = 1
    while idxnext <= len
        cnext = s[idxnext]
        if cnext == c
            count += 1
            @warn("++ and -- are deprecated in the OpenSMILES specification, please use + or - followed by multiplicity.", maxlog=1)
            idxnext = nextind(s, idxnext)
        elseif ispm(cnext)
            error("+- and -+ are not allowed for charge in SMILES")
        elseif isdigit(cnext)
            count == 1 || error("cannot use ++ or -- in charge specification followed by multiplicity")
            break
        else
            return sgn*count, idxnext
        end
    end
    idxnext > len && return sgn*count, idxnext
    mult, idxnext = tryparseint16(s, idxnext)
    mult === nothing && return sgn, idxnext
    @assert count == 1
    return sgn*mult, idxnext
end

function tryparsechirality(s::AbstractString, idx::Int)
    m = match(r"(@TH[12]|@AL[12]|@SP[1-3]|@TB\d\d?|@OH\d\d?|@@?)", s, idx)
    (m === nothing || m.offset > idx) && return nothing, idx
    chirality = m.captures[1]::AbstractString
    return chirality, m.offset + lastindex(chirality)
end


text:
module EquivalentCircuits

    export circuitevolution
    export parameteroptimisation
    export loadpopulation
    export Circuit, EquivalentCircuit
    using Random, Combinatorics, GeneralizedGenerated, DelimitedFiles, Distributions, Optim
    import Base: isless, length

    include("Circuits.jl")
    include("CircuitFunction.jl")
    include("EvolutionOperators.jl")
    include("ObjectiveFunction.jl")
    include("OptimizeParameters.jl")
    include("SimulateImpedance.jl")
    include("CircuitSimplification.jl")
    include("CircuitLibrary.jl")
    include("RedundancyTesting.jl")
    include("CircuitEvolution.jl")

end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Abstract base type for time integration schemes of storage class `2N`
abstract type SimpleAlgorithm2N end

"""
    CarpenterKennedy2N54()

The following structures and methods provide a minimal implementation of
the low-storage explicit Runge-Kutta method of

    Carpenter, Kennedy (1994) Fourth order 2N storage RK schemes, Solution 3

using the same interface as OrdinaryDiffEq.jl.
"""
struct CarpenterKennedy2N54 <: SimpleAlgorithm2N
  a::SVector{5, Float64}
  b::SVector{5, Float64}
  c::SVector{5, Float64}

  function CarpenterKennedy2N54()
    a = SVector(0.0, 567301805773.0 / 1357537059087.0,2404267990393.0 / 2016746695238.0,
                3550918686646.0 / 2091501179385.0, 1275806237668.0 / 842570457699.0)
    b = SVector(1432997174477.0 / 9575080441755.0, 5161836677717.0 / 13612068292357.0,
                1720146321549.0 / 2090206949498.0, 3134564353537.0 / 4481467310338.0,
                2277821191437.0 / 14882151754819.0)
    c = SVector(0.0, 1432997174477.0 / 9575080441755.0, 2526269341429.0 / 6820363962896.0,
                2006345519317.0 / 3224310063776.0, 2802321613138.0 / 2924317926251.0)

    new(a, b, c)
  end
end

"""
      CarpenterKennedy2N43()

Carpenter, Kennedy (1994) Third order 2N storage RK schemes with error control
"""
struct CarpenterKennedy2N43 <: SimpleAlgorithm2N
  a::SVector{4, Float64}
  b::SVector{4, Float64}
  c::SVector{4, Float64}

  function CarpenterKennedy2N43()
    a = SVector(0, 756391 / 934407, 36441873 / 15625000, 1953125 / 1085297)
    b = SVector(8 / 141, 6627 / 2000, 609375 / 1085297, 198961 / 526383)
    c = SVector(0, 8 / 141, 86 / 125, 1)

    new(a, b, c)
  end
end

# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1
mutable struct SimpleIntegrator2NOptions{Callback}
  callback::Callback # callbacks; used in Trixi
  adaptive::Bool # whether the algorithm is adaptive; ignored
  dtmax::Float64 # ignored
  maxiters::Int # maximal numer of time steps
  tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored
end

function SimpleIntegrator2NOptions(callback, tspan; maxiters=typemax(Int), kwargs...)
  SimpleIntegrator2NOptions{typeof(callback)}(
    callback, false, Inf, maxiters, [last(tspan)])
end

# This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77
# This implements the interface components described at
# https://diffeq.sciml.ai/v6.8/basics/integrator/#Handing-Integrators-1
# which are used in Trixi.
mutable struct SimpleIntegrator2N{RealT<:Real, uType, Params, Sol, Alg, SimpleIntegrator2NOptions}
  u::uType #
  du::uType
  u_tmp::uType
  t::RealT
  dt::RealT # current time step
  dtcache::RealT # ignored
  iter::Int # current number of time step (iteration)
  p::Params # will be the semidiscretization from Trixi
  sol::Sol # faked
  alg::Alg
  opts::SimpleIntegrator2NOptions
  finalstep::Bool # added for convenience
end

# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1
function solve(ode::ODEProblem, alg::T;
               dt, callback=nothing, kwargs...) where {T<:SimpleAlgorithm2N}
  u = copy(ode.u0)
  du = similar(u)
  u_tmp = similar(u)
  t = first(ode.tspan)
  iter = 0
  integrator = SimpleIntegrator2N(u, du, u_tmp, t, dt, zero(dt), iter, ode.p,
                  (prob=ode,), alg,
                  SimpleIntegrator2NOptions(callback, ode.tspan; kwargs...), false)

  # initialize callbacks
  if callback isa CallbackSet
    for cb in callback.continuous_callbacks
      error("unsupported")
    end
    for cb in callback.discrete_callbacks
      cb.initialize(cb, integrator.u, integrator.t, integrator)
    end
  elseif !isnothing(callback)
    error("unsupported")
  end

  solve!(integrator)
end

function solve!(integrator::SimpleIntegrator2N)
  @unpack prob = integrator.sol
  @unpack alg = integrator
  t_end = last(prob.tspan)
  callbacks = integrator.opts.callback

  integrator.finalstep = false
  @trixi_timeit timer() "main loop" while !integrator.finalstep
    if isnan(integrator.dt)
      error("time step size `dt` is NaN")
    end

    # if the next iteration would push the simulation beyond the end time, set dt accordingly
    if integrator.t + integrator.dt > t_end || isapprox(integrator.t + integrator.dt, t_end)
      integrator.dt = t_end - integrator.t
      terminate!(integrator)
    end

    # one time step
    integrator.u_tmp .= 0
    for stage in eachindex(alg.c)
      t_stage = integrator.t + integrator.dt * alg.c[stage]
      prob.f(integrator.du, integrator.u, prob.p, t_stage)

      a_stage    = alg.a[stage]
      b_stage_dt = alg.b[stage] * integrator.dt
      @trixi_timeit timer() "Runge-Kutta step" begin
        @threaded for i in eachindex(integrator.u)
          integrator.u_tmp[i] = integrator.du[i] - integrator.u_tmp[i] * a_stage
          integrator.u[i] += integrator.u_tmp[i] * b_stage_dt
        end
      end
    end
    integrator.iter += 1
    integrator.t += integrator.dt

    # handle callbacks
    if callbacks isa CallbackSet
      for cb in callbacks.discrete_callbacks
        if cb.condition(integrator.u, integrator.t, integrator)
          cb.affect!(integrator)
        end
      end
    end

    # respect maximum number of iterations
    if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep
      @warn "Interrupted. Larger maxiters is needed."
      terminate!(integrator)
    end
  end

  return TimeIntegratorSolution((first(prob.tspan), integrator.t),
                                (prob.u0, integrator.u),
                                integrator.sol.prob)
end

# get a cache where the RHS can be stored
get_du(integrator::SimpleIntegrator2N) = integrator.du
get_tmp_cache(integrator::SimpleIntegrator2N) = (integrator.u_tmp,)

# some algorithms from DiffEq like FSAL-ones need to be informed when a callback has modified u
u_modified!(integrator::SimpleIntegrator2N, ::Bool) = false

# used by adaptive timestepping algorithms in DiffEq
function set_proposed_dt!(integrator::SimpleIntegrator2N, dt)
  integrator.dt = dt
end

# stop the time integration
function terminate!(integrator::SimpleIntegrator2N)
  integrator.finalstep = true
  empty!(integrator.opts.tstops)
end

# used for AMR
function Base.resize!(integrator::SimpleIntegrator2N, new_size)
  resize!(integrator.u, new_size)
  resize!(integrator.du, new_size)
  resize!(integrator.u_tmp, new_size)
end

end # @muladd


text:
@json struct BranchChecks
    enforcement_level::String
    contexts::Vector{String}
end

@json struct BranchProtection
    enabled::Bool
    required_status_checks::BranchChecks
end

@json struct BranchLinks
    html::String
    self::String
end

@json struct Branch
    name::String
    commit::Commit
    _links => links::BranchLinks
    protected::Bool
    protection::BranchProtection
    protection_url::String
end

endpoint(::GitHubAPI, ::typeof(get_branch), owner::AStr, repo::AStr, branch::AStr) =
    Endpoint(:GET, "/repos/$owner/$repo/branches/$branch")
into(::GitHubAPI, ::typeof(get_branch)) = Branch


text:
################################################################################
#general functions
# work around https://github.com/JuliaStats/Distributions.jl/issues/846
rand(::AbstractRNG, sd::StandardizedDistribution) where {T} = rand(sd)

#loop invariant part of the kernel
@inline kernelinvariants(::Type{<:StandardizedDistribution}, coefs) = ()
################################################################################
#standardized
"""
    Standardized{D<:ContinuousUnivariateDistribution, T}  <: StandardizedDistribution{T}
A wrapper type for standardizing a distribution from Distributions.jl.
"""
struct Standardized{D<:ContinuousUnivariateDistribution, T<:AbstractFloat}  <: StandardizedDistribution{T}
    coefs::Vector{T}
end
Standardized{D}(coefs::T...) where {D, T} = Standardized{D, T}([coefs...])
Standardized{D}(coefs::Vector{T}) where {D, T} = Standardized{D, T}([coefs...])
rand(s::Standardized{D, T}) where {D, T} = (rand(D(s.coefs...))-mean(D(s.coefs...)))./std(D(s.coefs...))
@inline logkernel(S::Type{<:Standardized{D, T1} where T1}, x, coefs::Vector{T}) where {D, T} = (try sig=std(D(coefs...)); logpdf(D(coefs...), mean(D(coefs...)) + sig*x)+log(sig); catch; T(-Inf); end)
@inline logconst(S::Type{<:Standardized{D, T1} where T1}, coefs::Vector{T}) where {D, T} = zero(T)
nparams(S::Type{<:Standardized{D, T} where T}) where {D} = length(fieldnames(D))
coefnames(S::Type{<:Standardized{D, T}}) where {D, T} = [string.(fieldnames(D))...]
distname(S::Type{<:Standardized{D, T}}) where {D, T} = (io = IOBuffer(); sprint(io->Base.show_type_name(io, Base.typename(TDist{Float64}))))

function quantile(s::Standardized{D, T}, q::Real) where {D, T}
    (quantile(D(s.coefs...), q)-mean(D(s.coefs...)))./std(D(s.coefs...))
end

function constraints(S::Type{<:Standardized{D, T1} where T1}, ::Type{T})  where {D, T}
    lower = Vector{T}(undef, nparams(S))
    upper = Vector{T}(undef, nparams(S))
    fill!(lower, T(-Inf))
    fill!(upper, T(Inf))
    lower, upper
end

function startingvals(S::Type{<:Standardized{D, T1} where {T1}}, data::Vector{T})  where {T, D}
    svals = Vector{T}(undef, nparams(S))
    fill!(svals, eps(T))
    return svals
end

#for rand to work
Base.eltype(::StandardizedDistribution{T}) where {T} = T

"""
    fit(::Type{SD}, data; algorithm=BFGS(), kwargs...)

Fit a standardized distribution to the data, using the MLE. Keyword arguments
are passed on to the optimizer.
"""
function fit(::Type{SD}, data::Vector{T};
             algorithm=BFGS(), kwargs...
             ) where {SD<:StandardizedDistribution, T<:AbstractFloat}
    nparams(SD) == 0 && return SD{T}()
    obj = x -> -loglik(SD, data, x)
    lower, upper = constraints(SD, T)
    x0 = startingvals(SD, data)
    res = optimize(obj, lower, upper, x0, Fminbox(algorithm); kwargs...)
    coefs = res.minimizer
    return SD(coefs)
end

function loglik(::Type{SD}, data::Vector{<:AbstractFloat},
                coefs::Vector{T2}
                ) where {SD<:StandardizedDistribution, T2}
    T = length(data)
    length(coefs) == nparams(SD) || throw(NumParamError(nparams(SD), length(coefs)))
    @inbounds begin
        LL = zero(T2)
        iv = kernelinvariants(SD, coefs)
        @fastmath for t = 1:T
            LL += logkernel(SD, data[t], coefs, iv...)
        end#for
    end#inbounds
    LL += T*logconst(SD, coefs)
end#function

################################################################################
#StdNormal
"""
    StdNormal{T} <: StandardizedDistribution{T}

The standard Normal distribution.
"""
struct StdNormal{T} <: StandardizedDistribution{T}
    coefs::Vector{T}
    function StdNormal{T}(coefs::Vector) where {T}
        length(coefs) == 0 || throw(NumParamError(0, length(coefs)))
        new{T}(coefs)
    end
end
"""
    StdNormal(T::Type=Float64)
    StdNormal(v::Vector)
    StdNormal{T}()

Construct an instance of StdNormal.
"""
StdNormal(T::Type{<:AbstractFloat}=Float64) = StdNormal(T[])
StdNormal{T}() where {T<:AbstractFloat} = StdNormal(T[])
StdNormal(v::Vector{T}) where {T} = StdNormal{T}(v)
rand(::StdNormal{T}) where {T} = randn(T)

@inline logkernel(::Type{<:StdNormal}, x, coefs) = -abs2(x)/2
@inline logconst(::Type{<:StdNormal}, coefs::Vector{T}) where {T} =  -T(log2π)/2
nparams(::Type{<:StdNormal}) = 0
coefnames(::Type{<:StdNormal}) = String[]
distname(::Type{<:StdNormal}) = "Gaussian"

function constraints(::Type{<:StdNormal}, ::Type{T})  where {T<:AbstractFloat}
    lower = T[]
    upper = T[]
    return lower, upper
end

function startingvals(::Type{<:StdNormal}, data::Vector{T})  where {T<:AbstractFloat}
    return T[]
end

function quantile(::StdNormal, q::Real)
    norminvcdf(q)
end

################################################################################
#StdT

"""
    StdT{T} <: StandardizedDistribution{T}

The standardized (mean zero, variance one) Student's t distribution.
"""
struct StdT{T} <: StandardizedDistribution{T}
    coefs::Vector{T}
    function StdT{T}(coefs::Vector) where {T}
        length(coefs) == 1 || throw(NumParamError(1, length(coefs)))
        new{T}(coefs)
    end
end

"""
    StdT(ν)

Create a standardized t distribution with `ν` degrees of freedom. `ν`` can be passed
as a scalar or vector.
"""
StdT(ν) = StdT([ν])
StdT(ν::Integer) = StdT(float(ν))
StdT(ν::Vector{T}) where {T} = StdT{T}(ν)
(rand(d::StdT{T})::T) where {T}  =  (ν=d.coefs[1]; tdistrand(ν)*sqrt((ν-2)/ν))
@inline kernelinvariants(::Type{<:StdT}, coefs) = (1/ (coefs[1]-2),)
@inline logkernel(::Type{<:StdT}, x, coefs, iv) = (-(coefs[1] + 1) / 2) * log1p(abs2(x) *iv)
@inline logconst(::Type{<:StdT}, coefs)  = (lgamma((coefs[1] + 1) / 2)
                                               - log((coefs[1]-2) * pi) / 2
                                               - lgamma(coefs[1] / 2)
                                               )
nparams(::Type{<:StdT}) = 1
coefnames(::Type{<:StdT}) = ["ν"]
distname(::Type{<:StdT}) = "Student's t"

function constraints(::Type{<:StdT}, ::Type{T}) where {T}
    lower = T[20/10]
    upper = T[Inf]
    return lower, upper
end

function startingvals(::Type{<:StdT}, data::Array{T}) where {T}
    #mean of abs(t)
    eabst(ν)=2*sqrt(ν-2)/(ν-1)/beta(ν/2, 1/2)
    ##alteratively, could use mean of log(abs(t)):
    #elogabst(ν)=log(ν-2)/2-digamma(ν/2)/2+digamma(1/2)/2
    ht = T[]
    lht = T[]
    zt = T[]
    at = T[]
    loglik!(ht, lht, zt, at,  GARCH{1, 1}, StdNormal, Intercept(0.), data, vcat(startingvals(GARCH{1, 1}, data), startingvals(Intercept(0.), data)))
    lower = convert(T, 2)
    upper = convert(T, 30)
    z = mean(abs.(data.-mean(data))./sqrt.(ht))
    z > eabst(upper) ? [upper] : [find_zero(x -> z-eabst(x), (lower, upper))]
end

function quantile(dist::StdT, q::Real)
    ν = dist.coefs[1]
    tdistinvcdf(ν, q)*sqrt((ν-2)/ν)
end

################################################################################
#StdGED

"""
    StdGED{T} <: StandardizedDistribution{T}

The standardized (mean zero, variance one) generalized error distribution.
"""
struct StdGED{T} <: StandardizedDistribution{T}
    coefs::Vector{T}
    function StdGED{T}(coefs::Vector) where {T}
        length(coefs) == 1 || throw(NumParamError(1, length(coefs)))
        new{T}(coefs)
    end
end

"""
    StdGED(p)

Create a standardized generalized error distribution parameter `p`. `p` can be passed
as a scalar or vector.
"""
StdGED(p) = StdGED([p])
StdGED(p::Integer) = StdGED(float(p))
StdGED(v::Vector{T}) where {T} = StdGED{T}(v)

(rand(d::StdGED{T})::T) where {T} = (p = d.coefs[1]; ip=1/p;  (2*rand()-1)*gammarand(1+ip, 1)^ip * sqrt(gamma(ip) / gamma(3*ip)) )

@inline logconst(::Type{<:StdGED}, coefs)  = (p = coefs[1]; ip = 1/p; lgamma(3*ip)/2 - lgamma(ip)*3/2 - logtwo  - log(ip))
@inline logkernel(::Type{<:StdGED}, x, coefs, s) = (p = coefs[1]; -abs(x*s)^p)
@inline kernelinvariants(::Type{<:StdGED}, coefs) = (p = coefs[1]; ip = 1/p; (sqrt(gamma(3*ip) / gamma(ip)),))

nparams(::Type{<:StdGED}) = 1
coefnames(::Type{<:StdGED}) = ["p"]
distname(::Type{<:StdGED}) = "GED"

function constraints(::Type{<:StdGED}, ::Type{T}) where {T}
    lower = [zero(T)]
    upper = T[Inf]
    return lower, upper
end

function startingvals(::Type{<:StdGED}, data::Array{T}) where {T}
    ht = T[]
    lht = T[]
    zt = T[]
    at = T[]
    loglik!(ht, lht, zt, at, GARCH{1, 1}, StdNormal, Intercept(0.), data, vcat(startingvals(GARCH{1, 1}, data), startingvals(Intercept(0.), data)))
    z = mean((abs.(data.-mean(data))./sqrt.(ht)).^4)
    lower = T(0.05)
    upper = T(25.)
    f(r) = z-gamma(5/r)*gamma(1/r)/gamma(3/r)^2
    f(lower)>0 && return [lower]
    f(upper)<0 && return [upper]
    return T[find_zero(f, (lower, upper))]

end

function quantile(dist::StdGED, q::Real)
    p = dist.coefs[1]
    ip = 1/p
    qq = 2*q-1
    return sign(qq) * (gammainvcdf(ip, 1., abs(qq)))^ip/kernelinvariants(StdGED, [p])[1]
end

################################################################################
#Hansen's SKT-Distribution

"""
    StdSkewT{T} <: StandardizedDistribution{T}

Hansen's standardized (mean zero, variance one) Skewed Student's t distribution.
"""
struct StdSkewT{T} <: StandardizedDistribution{T}
    coefs::Vector{T}
    function StdSkewT{T}(coefs::Vector) where {T}
        length(coefs) == 2 || throw(NumParamError(2, length(coefs)))
        new{T}(coefs)
    end
end

"""
    StdSkewT(v,λ)

Create a standardized skewed t distribution with `v` degrees of freedom and `λ` shape parameter. `ν,λ`` can be passed
as scalars or vectors.
"""
StdSkewT(ν,λ) = StdSkewT([float(ν), float(λ)])
StdSkewT(coefs::Vector{T}) where {T} = StdSkewT{T}(coefs)

(rand(d::StdSkewT{T}) where {T} = (quantile(d, rand())))

@inline a(d::Type{<:StdSkewT}, coefs) =  (ν=coefs[1];λ=coefs[2]; 4λ*c(d,coefs) * ((ν-2)/(ν-1)))
@inline b(d::Type{<:StdSkewT}, coefs) =  (ν=coefs[1];λ=coefs[2]; sqrt(1+3λ^2-a(d,coefs)^2))
@inline c(d::Type{<:StdSkewT}, coefs) =  (ν=coefs[1];λ=coefs[2]; gamma((ν+1)/2) / (sqrt(π*(ν-2)) * gamma(ν/2)))

@inline kernelinvariants(::Type{<:StdSkewT}, coefs) = (1/ (coefs[1]-2),)
@inline function logkernel(d::Type{<:StdSkewT}, x, coefs, iv)
    ν=coefs[1]
    λ=coefs[2]
    c = gamma((ν+1)/2) / (sqrt(π*(ν-2)) * gamma(ν/2))
    a = 4λ * c * ((ν-2)/(ν-1))
    b = sqrt(1 + 3λ^2 -a^2)
    λsign = x < (-a/b) ? -1 : 1
    (-(ν + 1) / 2) * log1p(1/abs2(1+λ*λsign) * abs2(b*x + a) *iv)
end
@inline logconst(d::Type{<:StdSkewT}, coefs)  = (log(b(d,coefs))+(log(c(d,coefs))))

nparams(::Type{<:StdSkewT}) = 2
coefnames(::Type{<:StdSkewT}) = ["ν", "λ"]
distname(::Type{<:StdSkewT}) = "Hansen's Skewed t"

function constraints(::Type{<:StdSkewT}, ::Type{T}) where {T}
    lower = T[20/10, -one(T)]
    upper = T[Inf,one(T)]
    return lower, upper
end

startingvals(::Type{<:StdSkewT}, data::Array{T}) where {T<:AbstractFloat} = [startingvals(StdT, data)..., zero(T)]

function quantile(d::StdSkewT{T}, q::T) where T
    ν = d.coefs[1]
    λ = d.coefs[2]
    a_val = a(typeof(d),d.coefs)
    b_val = b(typeof(d),d.coefs)
    λconst = q < (1 - λ)/2 ? (1 - λ) : (1 + λ)
    quant_numer = q < (1 - λ)/2 ? q : (q + λ)
    1/b_val * ((λconst) * sqrt((ν-2)/ν) * tdistinvcdf(ν, quant_numer/λconst) - a_val)
end


text:
function molecule(::Molecule"OH")
  return """
O            0.000000000000     0.000000000000    -0.056372602727
H            0.000000000000     0.000000000000     0.894674113719
  """
end


text:
function soft_threshold_operator(βi::Float64, λ::Float64)
    if βi > λ
        return βi - λ
    elseif βi < -λ
        return βi + λ
    else
        return 0.0
    end
end

function pathwise_coordinate_optimization!(
        X::Vector{Vector{Float64}}, r::Vector{Float64}, β::Vector{Float64},
        λ::Float64, update_function::Function=identity!,
        args::Tuple=(); tol::Float64=1e-7, max_iter::Int=10000, intercept::Bool=true)
    # Get dimensions
    N = length(r)
    p = length(β)
    # Update coeficients iterativelly
    list_index = collect(1:p)
    # Main loop
    df = 0
    k = 0
    active_set = falses(p)
    full_iteration = true
    fullfil_tolerance = false
    change_on_active_set = false
    while true
        change_on_active_set = false
        fullfil_tolerance = true
        df = 0
        for i in list_index
            if ~full_iteration && ~active_set[i]
                continue
            end
            update_function(X, i, N, r, args...)
            βi = β[i]
            x = X[i]
            if βi != 0.0
                Base.LinAlg.BLAS.axpy!(βi, x, r)
            end
            β[i] = 1/dot(x, x)*soft_threshold_operator(dot(x, r), N*λ)
            if β[i] != 0.0
                Base.LinAlg.BLAS.axpy!(-β[i], x, r)
                df += 1
                if ~active_set[i]
                    change_on_active_set = true
                end
                active_set[i] = true
            end
            if abs(β[i] - βi) > tol
                fullfil_tolerance = false
            end
        end
        k += 1
        # Next iteration logic
        if full_iteration
            if change_on_active_set
                full_iteration = false
            elseif fullfil_tolerance
                break
            end
        elseif fullfil_tolerance
            full_iteration = true
        end
        if k > max_iter
            break
        end
    end
    return k, df
end

function pathwise_coordinate_optimization(
        X::Matrix{Float64}, y::Vector{Float64}, λ::Real,
        update_function::Function=identity!, args::Tuple=(); kwargs...)
        y_copy = copy(y)
        X_views, N, p = matrix2subviews(X)
        β = zeros(p)
        df, k = pathwise_coordinate_optimization!(X_views, y_copy, β, λ,
                                          update_function, args; kwargs...)
        return β, λ, df, k
end

function pathwise_coordinate_optimization(
        X::Matrix{Float64}, y::Vector{Float64}, λ_list::Vector{Float64},
        update_function::Function=identity!, args::Tuple=(); max_df::Int=1000,
        kwargs...)
    n = length(λ_list)
    X_views, N, p = matrix2subviews(X)
    r = copy(y)
    β_aux = zeros(p)
    β_matrix = zeros(p, n)
    df_matrix = zeros(Int, n)
    total_iter = 0
    for (i, λ) in enumerate(λ_list)
        k, df = pathwise_coordinate_optimization!(X_views, r,
                                                  β_aux, λ,
                                                  update_function, args;
                                                  kwargs...)
        β_matrix[:, i] .= β_aux
        df_matrix[i] = df
        total_iter += k
        if df >= max_df
            return β_matrix[:, 1:i], λ_list[1:i], df_matrix[1:i], total_iter
        end
    end
    return β_matrix, λ_list, df_matrix, total_iter
end

function pathwise_coordinate_optimization(
        X::Matrix{Float64}, y::Vector{Float64},
        update_function::Function=identity!, args::Tuple=();
        λ_min_ratio::Float64=1e-4, nλ::Int=100, kwargs...)
    N, p = size(X)
    β = zeros(p)
    λ_max = 1/N*maximum(abs.(X'*y))
    λ_min = λ_max*λ_min_ratio
    λ_list = logspace(log10(λ_min), log10(λ_max), nλ)[end:-1:1]
    pathwise_coordinate_optimization(X, y, λ_list, update_function, args;
                                     kwargs...)
end


text:
using Random

@testset "RL.jl" begin

    include("representation/representation.jl")
    include("nn_structures/nn_structures.jl")
    include("utils.jl")

end


text:
"""
    assigndist_normal(μ, lower, upper; nσ = 1,
        trunc_lower = -Inf, trunc_upper = Inf,
        tolerance = 1e-2μ)

Assign parameters to a Normal distribution based on location
parameter `μ`, together with `lower` and `upper` uncertainty
bounds. `nσ` indicates how many standard deviations away
from `μ` `lower`/`upper` are (they must be equidistant
from `μ`). `trunc_lower` and `trunc_upper` truncated the
distribution if specified (defaults to `-Inf` and `Inf`).
"""
function assigndist_normal(μ, lower, upper; nσ = 1,
        trunc_lower = -Inf, trunc_upper = Inf,
        tolerance = 1e-7)

    dist_from_μ_upper = upper - μ
    dist_from_μ_lower = μ - lower

    if abs((upper - μ) - (μ - lower)) > tolerance
        throw(DomainError("(μ - lower, upper - μ) = ($dist_from_μ_lower, $dist_from_μ_upper): lower and upper bounds are not equidistant-ish from μ. Cannot create normal distribution."))
    end
    σ_estim = dist_from_μ_upper/nσ
    Truncated(Normal(μ, σ_estim), trunc_lower, trunc_upper)
end

function assigndist_normal(μ, σ; nσ = 1, trunc_lower = -Inf, trunc_upper = Inf,
        tolerance = 1e-7)

    σ_estim = σ/nσ
    Truncated(Normal(μ, σ_estim), trunc_lower, trunc_upper)
end

export assigndist_normal


text:
##
# Assign Job Ids. A one up number
global __next_job_id = Int64(0)::JobId
function _next_job_id()
    global __next_job_id
    __next_job_id += 1
end

##
# Worker task type. One of: map, file_info, reduce(in future)
type WorkerTaskMap <: WorkerTask
    jid::JobId
    blk_url::AbstractString
end

type WorkerTaskFileInfo <: WorkerTask
    jid::JobId
    source::MRHdfsFileInput
end

type WorkerTaskInitJob <: WorkerTask
    jid::JobId
    inp_typ::Any
    fn_find_rec::Function
    fn_map::Function
    fn_collect::Function
    fn_reduce::FuncNone
end

type WorkerTaskUnloadJob <: WorkerTask
    jid::JobId
end

type WorkerTaskFetchCollected <: WorkerTask
    jid::JobId
end

##
# Iterator for MapInputReader using the find_rec function
type MapInputIterator
    r::MapInputReader
    fn_find_rec::Function
    is_done::Bool
    rec::Any

    function MapInputIterator(r::MapInputReader, url::AbstractString, fn_find_rec::Function)
        reset_pos(r, url)
        new(r, fn_find_rec, false, nothing)
    end
end

function start(iter::MapInputIterator)
    iter.rec, iter.is_done, state = iter.fn_find_rec(iter.r, nothing)
    state
end
function next(iter::MapInputIterator, state)
    ret_rec = iter.rec
    iter.rec, iter.is_done, state = iter.fn_find_rec(iter.r, state)
    (ret_rec, state)
end
done(iter::MapInputIterator, state) = iter.is_done

function _worker_task(t::WorkerTaskInitJob)
    #logmsg("WorkerTaskInitJob jid:$(t.jid)")
    js = ((myid() == 1) ? _def_wrkr_job_store : _job_store)
    js[t.jid] = HdfsJobCtx(t.jid, t.inp_typ, t.fn_find_rec, t.fn_map, t.fn_collect, t.fn_reduce)
    #logmsg("WorkerTaskInitJob jid:$(t.jid) initialized")
    t.jid
end

function _worker_task(t::WorkerTaskUnloadJob)
    d = ((myid() == 1) ? _def_wrkr_job_store : _job_store)
    j = d[t.jid]
    delete!(d, t.jid)
    isa(j.info, HdfsJobRunInfo) && isa(j.info.rdr, MapStreamInputReader) && close(j.info.rdr)
    t.jid
end

function _worker_task(t::WorkerTaskFetchCollected)
    ((myid() == 1) ? _def_wrkr_job_store : _job_store)[t.jid].info.results
end

function _worker_task(t::WorkerTaskMap)
    j = ((myid() == 1) ? _def_wrkr_job_store : _job_store)[t.jid]
    jinfo = j.info
    recs = 0
    fn_map = j.fn_map
    fn_collect = j.fn_collect

    results = jinfo.results
    for rec in MapInputIterator(jinfo.rdr, t.blk_url, j.fn_find_rec)
        #logmsg("MapInputIterator got: rec")
        (nothing == rec) && continue
        for mapped_rec in fn_map(rec)
            results = fn_collect(results, mapped_rec)
            recs += 1
        end
    end
    jinfo.results = results
    #logmsg((recs == 0) ? "no usable record in block" : "no usable record at end of block")
    recs
end

function _worker_task(t::WorkerTaskFileInfo)
    expand_file_inputs(t.source)
end

##
# scheduler function
# TODO: queued task priorities to be linked to the chain
function _sched()
    set_priorities((_1,qt,_2)->Float64(qt.wtask.jid))
end

function _callback(t::WorkerTaskInitJob, ret)
    jid = t.jid
    #logmsg("CB> WorkerTaskInitJob jid:$jid/$ret")
    !haskey(_job_store, jid) && return # if job is deleted because of a past error
    j = _job_store[jid]

    (ret != jid) && return _set_status(j, STATE_ERROR, ret)
    if(0 == (j.info.num_pending_inits -= 1))
        #logmsg("distributing blocks for job $jid")
        @async begin
            try
                _distribute(jid, j.info.source_spec)
            catch ex
                #logmsg("error distributing job $jid")
                _set_status(j, STATE_ERROR, ex)
            end
        end
    end
end

function _callback(t::WorkerTaskFetchCollected, ret)
    jid = t.jid
    !haskey(_job_store, jid) && return # if job is deleted because of a past error
    j = _job_store[jid]

    if(isa(ret,Exception))
        _set_status(j, STATE_ERROR, ret)
        dequeue_worker_task(t)
        return
    end
    ji = j.info
    try
        ji.red = j.fn_reduce(ji.red, ret)
        ji.num_reduces_done += 1
        (ji.num_reduces == ji.num_reduces_done) && _set_status(j, STATE_COMPLETE)
    catch ex
        #logmsg("error reducing job $jid")
        _set_status(j, STATE_ERROR, ex)
        return
    end
end

function _callback(t::WorkerTaskFileInfo, ret)
    jid = t.jid
    j = _job_store[jid]
    try
        #logmsg("WorkerTaskFileInfo callback got $(typeof(ret))")
        !isa(ret, MRHdfsFileInput) && throw(ret)

        #logmsg("WorkerTaskFileInfo callback got $(ret.file_list)")
        nparts = 0
        for (idx,fname) in enumerate(ret.file_list)
            #logmsg("WorkerTaskFileInfo idx: $idx fname: $fname")
            blk_dist = ret.file_blocks[idx]
            #logmsg("WorkerTaskFileInfo blk_dist: $blk_dist")
            nparts += length(blk_dist)
            for (blk_id,macs) in enumerate(blk_dist)
                #logmsg("WorkerTaskFileInfo blk_id:$blk_id macs:$macs")
                #TODO: use actual offset instead of block id
                maclist = macs[2]
                wt = WorkerTaskMap(jid, string(fname, '#', blk_id))
                #logmsg("WorkerTaskFileInfo wt: $wt")
                qt = QueuedWorkerTask(wt, _worker_task, _callback, maclist)
                #logmsg("WorkerTaskFileInfo qt: $qt")
                queue_worker_task(qt)
                #logmsg("WorkerTaskFileInfo queued")
            end
        end
        _start_running(jid, nparts)
    catch ex
        logmsg("error starting job $(jid): $ex")
        _set_status(j, STATE_ERROR, ex)
    end
end

function _callback(t::WorkerTaskMap, ret)
    jid = t.jid
    j = _job_store[jid]
    (j.info.state == STATE_ERROR) && return
    !isa(ret, Integer) && (return _set_status(j, STATE_ERROR, ret))

    logmsg("num records mapped: $(ret)")
    j.info.num_parts_done += 1

    if(j.info.num_parts == j.info.num_parts_done)
        (j.fn_reduce == nothing) && (return _set_status(j, STATE_COMPLETE))
        j.info.num_reduces = num_remotes()
        queue_worker_task(QueuedWorkerTask(WorkerTaskFetchCollected(jid), _worker_task, _callback, :wrkr_all))
    end
end

##
# Job store and scheduler
type HdfsJobSchedInfo
    source_spec::MRInput                                    # job_id or url (optionally with wildcard)
    num_parts::Int
    num_parts_done::Int
    num_pending_inits::Int
    num_reduces::Int
    num_reduces_done::Int
    state::Int
    state_info::Any
    red::Any
    trigger::RemoteRef
    begin_time::Float64                                     # submitted
    sched_time::Float64                                     # scheduled (after dependent jobs finish)
    end_time::Float64                                       # completed

    function HdfsJobSchedInfo(source_spec::MRInput)
        new(source_spec, 0, 0, 0, 0, 0, STATE_STARTING, nothing, nothing, RemoteRef(), time(), 0.0, 0.0)
    end
end

type HdfsJobRunInfo
    rdr::MapInputReader
    results::Any
end

type HdfsJobCtx
    jid::JobId
    fn_find_rec::Function
    fn_map::Function
    fn_collect::Function
    fn_reduce::FuncNone

    info::Union{HdfsJobRunInfo, HdfsJobSchedInfo}           # holds sched info at master node and run info at worker nodes

    # constructor at the scheduler end
    function HdfsJobCtx(jid::JobId, source::MRInput, fn_map::Function, fn_collect::Function, fn_reduce::FuncNone)
        new(jid, source.reader_fn, fn_map, fn_collect, fn_reduce, HdfsJobSchedInfo(source))
    end

    # constructor at the worker end
    function HdfsJobCtx(jid::JobId, inp_typ, fn_find_rec::Function, fn_map::Function, fn_collect::Function, fn_reduce::FuncNone)
        rdr = get_input_reader(inp_typ...)
        new(jid, fn_find_rec, fn_map, fn_collect, fn_reduce, HdfsJobRunInfo(rdr,nothing))
    end
end

const _def_wrkr_job_store = Dict{JobId, HdfsJobCtx}()
const _job_store = Dict{JobId, HdfsJobCtx}()

const STATE_ERROR = 0
const STATE_STARTING = 1
const STATE_RUNNING = 2
const STATE_COMPLETE = 3

function _start_running(jid::JobId, nparts)
    j = _job_store[jid]
    j.info.num_parts = nparts
    if(j.info.state != STATE_ERROR)
        #logmsg("distributed blocks for job $jid")
        _set_status(j, STATE_RUNNING)
        _sched()
        #logmsg("scheduled blocks for job $jid")
    end
end

function _distribute(jid::JobId, source::MRMapInput)
    qtarr = map(x->QueuedWorkerTask(WorkerTaskMap(jid, string(x)), _worker_task, _callback, :wrkr_all), source.job_list)
    for src_jid in source.job_list
        if(STATE_COMPLETE != wait(src_jid))
            j = _job_store[jid]
            _set_status(j, STATE_ERROR, "source job id $(src_jid) has errors")
            return
        end
    end
    nparts = length(qtarr) * num_remotes()
    for qt in qtarr queue_worker_task(qt) end
    _start_running(jid, nparts)
end
function _distribute(jid::JobId, source::MRHdfsFileInput)
    queue_worker_task(QueuedWorkerTask(WorkerTaskFileInfo(jid,source), _worker_task, _callback, :wrkr_any))
end
function _distribute(jid::JobId, source::MRFsFileInput)
    j = _job_store[jid]
    try
        expand_file_inputs(source)
        nparts = 0
        for (idx,fname) in enumerate(source.file_list)
            nfileblks = round(Int, ceil(source.file_info[idx].size/source.block_sz))
            nparts += nfileblks
            for blk_id in 1:nfileblks
                queue_worker_task(QueuedWorkerTask(WorkerTaskMap(jid, string(fname, '#', blk_id)), _worker_task, _callback, :wrkr_any))
            end
        end
        _start_running(jid, nparts)
    catch ex
        logmsg("error starting job $(jid)")
        _set_status(j, STATE_ERROR, ex)
    end
end

dmap(source::MRInput, fn_map::Function, fn_collect::Function) = dmapreduce(source, fn_map, fn_collect, nothing)
function dmapreduce(source::MRInput, fn_map::Function, fn_collect::Function, fn_reduce::FuncNone)
    prep_remotes()

    jid = _next_job_id()
    j = HdfsJobCtx(jid, source, fn_map, fn_collect, fn_reduce)
    _job_store[jid] = j

    j.info.num_pending_inits = num_remotes()
    queue_worker_task(QueuedWorkerTask(WorkerTaskInitJob(jid, input_reader_type(source), source.reader_fn, fn_map, fn_collect, fn_reduce), _worker_task, _callback, :wrkr_all))
    jid
end

wait(jid::JobId) = wait(_job_store[jid])
wait(j::HdfsJobCtx) = fetch(j.info.trigger)

status(jid::JobId, desc::Bool=false) = status(_job_store[jid], desc)
function status(j::HdfsJobCtx, desc::Bool=false)
    ji = j.info
    function status_str()
        (ji.state == STATE_RUNNING) && (return "running")
        (ji.state == STATE_COMPLETE) && (return "complete")
        (ji.state == STATE_STARTING) && (return "starting")
        (ji.state == STATE_ERROR) && (return "error")
        error("job in unknown state")
    end

    !desc && return status_str()

    descinfo = ji.state_info
    (ji.state == STATE_RUNNING) && (descinfo = (0 == ji.num_parts) ? 0 : round(Int, ji.num_parts_done*100/ji.num_parts))
    (status_str(), descinfo)
end

function _set_status(j::HdfsJobCtx, state::Int, desc=nothing)
    j.info.state = state
    (nothing != desc) && (j.info.state_info = desc)
    (STATE_ERROR == state) && dequeue_worker_task((x)->(!isa(x.wtask, WorkerTaskMap) || (j.jid != x.wtask.jid)))
    (STATE_RUNNING == state) && (j.info.sched_time = time())
    ((STATE_COMPLETE == state) || (STATE_ERROR == state)) && (j.info.end_time = time(); !isready(j.info.trigger) && put!(j.info.trigger, state))
    state
end

results(jid::JobId) = results(_job_store[jid])
results(j::HdfsJobCtx) = (status(j), j.info.red)

function unload(jid::JobId)
    j = _job_store[jid]
    (j.info.state == STATE_RUNNING) && error("can't unload running job")
    queue_worker_task(QueuedWorkerTask(WorkerTaskUnloadJob(jid), _worker_task, nothing, :wrkr_all))
    delete!(_job_store, jid)
    jid
end

function times(jid::JobId)
    ji = _job_store[jid].info
    t_total = (ji.end_time - ji.begin_time)
    t_sched = max(0, (ji.sched_time - ji.begin_time))
    t_run = (t_sched > 0) ? (ji.end_time - ji.sched_time) : 0
    (t_total, t_sched, t_run)
end

max_stars_repo_name: tanmaykm/HadoopBlocks

max_issues_repo_name: tanmaykm/HadoopBlocks

max_forks_repo_name: tanmaykm/HadoopBlocks


text:
using Slurm.Directives

export Sbatch,
    prefix,
    ARRAY,
    ACCOUNT,
    COMMENT,
    CPUS_PER_GPU,
    CPUS_PER_TASK,
    DISTRIBUTION,
    INPUT,
    JOB_NAME,
    MAIL_TYPE,
    MAIL_USER,
    MEM,
    MEM_PER_CPU,
    NODES,
    NTASKS,
    NTASKS_PER_CORE,
    NTASKS_PER_NODE,
    OUTPUT,
    PARTITION,
    PRIORITY,
    QUIET,
    RESERVATION,
    CORE_SPEC,
    TIME,
    UID,
    VERSION,
    VERBOSE,
    WAIT,
    EXCLUDE

struct Sbatch{A <: AbstractVector{<: Directive}, B <: AbstractString}
    directives::A
    script::B
end

prefix(::Sbatch) = "SBATCH"

ARRAY(v) = Directive("array", 'a', v)
ACCOUNT(v) = Directive("account", 'A', v)
COMMENT(v) = Directive("comment", nothing, v)
CPUS_PER_GPU(v) = Directive("cpus-per-gpu", nothing, v)
CPUS_PER_TASK(v) = Directive("cpus-per-task", 'c', v)
DISTRIBUTION(v) = Directive("distribution", 'm', v)
INPUT(v::AbstractString) = Directive("input", 'i', v)
JOB_NAME(v::AbstractString) = Directive("job-name", 'J', v)
MAIL_TYPE(v) = Directive("mail-type", nothing, v)
MAIL_USER(v::AbstractString) = Directive("mail-user", nothing, v)
MEM(v) = Directive("mem", nothing, v)
MEM_PER_CPU(v) = Directive("mem-per-cpu", nothing, v)
NODES(v) = Directive("nodes", 'N', v)
NTASKS(v) = Directive("ntasks", 'n', v)
NTASKS_PER_CORE(v::Integer) = Directive("ntasks-per-core", nothing, v)
NTASKS_PER_NODE(v::Integer) = Directive("ntasks-per-node", nothing, v)
OUTPUT(v::AbstractString) = Directive("output", 'o', v)
PARTITION(v) = Directive("partition", 'p', v)
PRIORITY(v) = Directive("priority", nothing, v)
QUIET(v::Bool) = Directive("quiet", 'Q', v)
RESERVATION(v) = Directive("reservation", nothing, v)
CORE_SPEC(v) = Directive("core-spec", 'S', v)
TIME(v) = Directive("time", 't', v)
UID(v) = Directive("uid", nothing, v)
VERSION(v) = Directive("version", 'V', v)
VERBOSE(v) = Directive("verbose", 'v', v)
WAIT(v) = Directive("wait", 'W', v)
EXCLUDE(v) = Directive("exclude", 'x', v)


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    varshape(d::Distributions.Distribution)::AbstractValueShape

Get the value shape of the variates of distribution `d`.
"""
varshape(d::UnivariateDistribution) = ScalarShape{Real}()
varshape(d::MultivariateDistribution) = ArrayShape{Real}(size(d)...)
varshape(d::MatrixDistribution) = ArrayShape{Real}(size(d)...)

@deprecate valshape(d::Distribution) varshape(d)

"""
    vardof(d::Distributions.Distribution)

Get the number of degrees of freedom of the variates of distribution `d`.
Equivalent to `totalndof(varshape(d))`.
"""
vardof(d::Distribution) = totalndof(varshape(d))