text:
function my_s3_get_range(aws::AWSConfig, bucket, path, rstart::Int, rend::Int; version="",
                                              retry=true,
                                              raw=false)

    hdr=Dict("Range"=>"bytes=$rstart-$rend")

    @repeat 4 try

        return AWSS3.s3(aws, "GET", bucket; path = path,
                                      version = version,
                                      return_raw = raw,
                                      headers = hdr)

    catch e
        @delay_retry if retry && ecode(e) in ["NoSuchBucket", "NoSuchKey"] end
    end
end

my_s3_get_range(a...; b...) = my_s3_get_range(default_aws_config(), a...; b...)


text:
using Clp
using Compat
using Compat.Test
import Compat.Pkg

@testset "MathProgBase" begin
    include("mathprog.jl")
end

@testset "MathOptInterface" begin
    include("MOIWrapper.jl")
end

ext: jl


text:
using ..Tagging, ..Traits, ..Var, ..Space, ..Basis
# @inline hasintervene(ω) = hastag(ω, Val{:intervene})

@inline tagintervene(::trait(Intervene), ω, intervention) =
  let i = mergeinterventions(intervention, ω.tags.intervene.intervention)
    mergetag(ω, (intervene = (intervention = i, intctx = ω.tags.intervene),))
  end

@inline tagintervene(trait, ω, intervention) =
  mergetag(ω, (intervene = (intervention = intervention, intctx = NoIntervention()),))

@inline tagintervene(ω, intervention) =
  tagintervene(traits(ω), ω, intervention)
  # tag(ω, (intervene = (intervention = intervention, intctx = ictx(ω),),), mergetags)

@inline Var.recurse(x::Intervened, ω) = x.x(tagintervene(ω, x.i))
# @inline (x::Intervened)(ω) = x.x(tagintervene(ω, x.i))

@inline Var.recurse(x::Intervened{X, <: HigherIntervention}, ω) where X =
  x.x(tagintervene(ω, x.i(ω)))

replaceintervene(ω, intctx) = mergetag(ω, (intervene = intctx,))
rmintervention(ω) = Basis.rmtag(ω, Val{:intervene})

@inline applyintervention(i::ValueIntervention, ω, intctx) = i.v
@inline applyintervention(i::Intervention, ω, intctx) = i.v(replaceintervene(ω, intctx))
@inline applyintervention(i::Intervention, ω, intctx::NoIntervention) = i.v(rmintervention(ω))

@inline function passintervene(traits, i::Intervention{X, V}, intctx, x::X, ω) where {X, V}
  # If the variable ` x` to be applied to ω is the variable to be replaced
  # then replace it, and apply its replacement to ω instead.
  # Othewise proceed as normally
  if i.x == x
    applyintervention(i, ω, intctx)
    # i.v(replacetags(ω, intctx))
  else
    prepostapply(traits, x, ω)
  end
end

@inline function passintervene(traits, i::ValueIntervention{X, V}, intctx, x::X, ω) where {X, V}
  if i.x == x
    i.v
  else
    prepostapply(traits, x, ω)
  end
end

# @inline function passintervene(traits, i::ValueIntervention{X, V}, intctx, x::X2, ω) where {X, X2, V}
#   @show X, X2
#   @assert false
# end

function passintervene(traits,
  i::Union{MultiIntervention{<:Tuple{SlightlyLessAbstractIntervention{X1, V1},
                                   SlightlyLessAbstractIntervention{X2, V2}}},
           MultiIntervention{<:Tuple{SlightlyLessAbstractIntervention{X1, V1},
                                   SlightlyLessAbstractIntervention{X2, V2},
                                   SlightlyLessAbstractIntervention{X3, V3}}},
           MultiIntervention{<:Tuple{SlightlyLessAbstractIntervention{X1, V1},
                                   SlightlyLessAbstractIntervention{X2, V2},
                                   SlightlyLessAbstractIntervention{X3, V3},
                                   SlightlyLessAbstractIntervention{X4, V4}}},
           MultiIntervention{<:Tuple{SlightlyLessAbstractIntervention{X1, V1},
                                   SlightlyLessAbstractIntervention{X2, V2},
                                   SlightlyLessAbstractIntervention{X3, V3},
                                   SlightlyLessAbstractIntervention{X4, V4},
                                   SlightlyLessAbstractIntervention{X5, V5}}}},
  intctx,
  x::Union{X1, X2, X3, X4, X5},
                       ω) where {X1, X2, X3, X4, X5, V1, V2, V3, V4, V5}
  if x == i.is[1].x
    # i.is[1].v(ω)
    applyintervention(i.is[1], ω, intctx)
  elseif x == i.is[2].x
    # i.is[2].v(ω)
    applyintervention(i.is[2], ω, intctx)
  elseif length(i.is) >= 3 && x == i.is[3].x
    # i.is[3].v(ω)
    applyintervention(i.is[3], ω, intctx)
  elseif length(i.is) >= 4 && x == i.is[4].x
    # i.is[4].v(ω)
    applyintervention(i.is[4], ω, intctx)
  elseif length(i.is) >= 5 && x == i.is[5].x
    # i.is[5].v(ω)
    applyintervention(i.is[5], ω, intctx)
  else
    prepostapply(traits, x, ω)
  end
end

function passintervene(traits,
                       i::MultiIntervention{XS},
                       tags,
                       x,
                       ω) where XS
  index = 1;
  while (index <= length(i.is) && x != i.is[index].x)
    index += 1
  end
  if (index <= length(i.is))
    i.is[index].v(ω)
  else
    prepostapply(traits, x, ω)
  end
end

# We only consider intervention if the intervention types match
@inline passintervene(traits, i::AbstractIntervention, intctx, x, ω) =
  prepostapply(traits, x, ω)

Var.dispatch(traits::trait(Intervene), f, ω::AbstractΩ) =
  passintervene(traits, ω.tags.intervene.intervention, ω.tags.intervene.intctx, f, ω)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using ExcelUtilities
using Test

@testset "ExcelUtilities.jl" begin
    @testset "vectors" begin
        excel_row = "1\t2\t3\t4"
        excel_col = "1\n2\n3\n4"
        d = ExcelUtilities.xlclip_reader(excel_row)
        @test d == [1,2,3,4]
        @test ExcelUtilities.xlclip_writer(d) == excel_col
        @test ExcelUtilities.xlclip_writer(d') == excel_row

        d = ExcelUtilities.xlclip_reader(excel_col)
        @test d == [1,2,3,4]
        @test ExcelUtilities.xlclip_writer(d) == excel_col

    end

    @testset "arrays" begin
        excel_arr = "1\t2\t3\t4\n5\t6\t7\t8\n9\t10\t11\t12"

        d = ExcelUtilities.xlclip_reader(excel_arr)
        @test d == [1 2 3 4;5 6 7 8;9 10 11 12]
        @test ExcelUtilities.xlclip_writer(d) == excel_arr

    end

end

max_issues_repo_name: alecloudenback/ExcelUtilities

max_forks_repo_name: alecloudenback/ExcelUtilities


text:
using SeisIO, SeisNoise, JLD2, PlotlyJS

# This is a script to test the functionality of the MWCS functions of Noise.jl
# This assumes synthetic data generated via genSynth.jl

function testMWCS(finame::String, foname::String, InputDict::Dict, windowlenlist::Array{Float64,1}, winsteplist::Array{Float64,1}, smoothlist::Array{Int64,1}, type::String; plot::Bool=false)
    valData = jldopen(finame)
    dvVlist = valData["info/dvVlist"]
    noiselist = valData["info/noiselist"]

    f_err = jldopen(foname, "a+")
    # hacky way to write metadata only once
    if type == "realData"
        f_err["info/dvV"] = dvVlist
        f_err["info/noise"] = noiselist
        f_err["info/winstep"] = winsteplist
        f_err["info/winlen"] = windowlenlist
        f_err["info/smooth"] = smoothlist
    end
    for dvV in dvVlist
        println("dv/v: $dvV")
        for noiselvl in noiselist
            for win_len in windowlenlist, win_step in windowsteplist, smooth in smoothlist
                signals = valData["$type/$dvV.$noiselvl"]
                u1 = signals["ref"]
                u2 = signals["cur"]

                dist = 0.0
                true_dtt = -dvV
                # compute dt
                time_axis, dt, err, coh = try
                    mwcs(u1,
                         u2,
                         InputDict["freqmin"],
                         InputDict["freqmax"],
                         InputDict["fs"],
                         InputDict["mintimelag"],
                         win_len,
                         win_step,
                         smooth)
                catch y
                    #println("Failed on first step")
                    #println(y)
                    continue
                end

                dt .*= 2π
                err .*= 2π
                # compute WLS regression for dt/t
                m, em, a, ea, m0, em0 = try
                    mwcs_dvv(time_axis,
                             dt,
                             err,
                             coh,
                             InputDict["dtt_lag"],
                             dist/1000,
                             InputDict["dtt_v"],
                             InputDict["dtt_minlag"],
                             InputDict["dtt_width"],
                             InputDict["dtt_sides"])

                catch y
                    #println("Failed on second step")
                    #println(y)
                    continue
                end

                #=
                # Report performance
                println("Noise level: $(noiselvl*100)%")
                println("Smoothing: $smooth")
                println("Window_length: $win_len")
                println("True: $(-true_dtt*100)%")
                println("Estimated: $(-m0*100)%")
                error = abs((m0 - true_dtt)/true_dtt)*100
                println("Error: $error%")
                =#

                if plot==true
                    plot_dtt(time_axis, dt, coh, m0, 0, err, true_dtt)
                end

                f_err["$type/$dvV.$noiselvl/$win_len.$win_step.$smooth"] = [m, em]
            end
        end
    end
    close(f_err)
    close(valData)
end

function plot_dtt(time_axis, dt, coh, m, a, err, true_dtt)
    # Plot dt vs t, error bars, true dt vs t, and the best fit line
    trace1 = scatter(;x=time_axis, y=dt, mode="lines+markers", marker_color=coh, line_color="black", name="Calculated dt/t: $((round(m, sigdigits=3))*100)%")
    true_shift = scatter(;x=time_axis, y=true_dtt*time_axis, mode="lines", line_color="blue", name="True dt/t: $(true_dtt*100)%")
    trace2 = scatter(;x=time_axis, y=(m*time_axis) .+ a, mode="lines", line_color="green", name="Linear Regression")
    plots = [trace1, true_shift, trace2]

    # plot error bars from mwcs()
    err1 = dt .+ err
    err2 = dt .- err
    for i=1:length(dt)
        t = scatter(;x=[time_axis[i], time_axis[i]], y=[err1[i], err2[i]], mode="lines", line_color="black", showlegend=false)
        push!(plots, t)
    end

    layout = Layout(;title="MWCS dt/t", xaxis_title="t (s)", yaxis_title="dt (s)")
    p = PlotlyJS.plot(plots, layout)
    display(p)
    readline()
end

# MWCS input parameters
InputDict_real = Dict( "freqmin"            => 0.1,
                       "freqmax"            => 9.9,
                       "fs"                 => 20.0,
                       "mintimelag"         => -100.0,
                       "dtt_width"          => 200.0,
                       "dtt_lag"            => "static",
                       "dtt_v"              => 1.0,
                       "dtt_minlag"         => 0.0,
                       "dtt_sides"          => "both")

InputDict_synth = Dict( "freqmin"            => 0.1,
                        "freqmax"            => 9.9,
                        "fs"                 => 20.0,
                        "mintimelag"         => 0.0,
                        "dtt_width"          => 200.0,
                        "dtt_lag"            => "static",
                        "dtt_v"              => 1.0,
                        "dtt_minlag"         => 0.0,
                        "dtt_sides"          => "both")

finame = "verificationData.jld2"
foname = "verificationDataError_MWCS_real.jld2"
windowlenlist = collect(20.:10.:60.)
windowsteplist = [10.0]
smoothlist = [3, 9, 15, 21]
type = "realData"

for type in ["rickerConv", "dampedSinusoid", "realData"]
    if type == "realData"
        InputDict = InputDict_real
    else
        InputDict = InputDict_synth
    end
    testMWCS(finame, foname, InputDict, windowlenlist, windowsteplist, smoothlist, type)
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
#-----------------------------------------------------------------------# common
abstract type HistogramStat{T} <: OnlineStat{T} end

# Index of `edges` that `y` belongs in, depending on if bins are `left` closed and if
# the bin on the end is `closed` instead of half-open.
function binindex(edges::AbstractVector, y, left::Bool, closed::Bool)
    a, b = extrema(edges)
    if y < a
        0
    elseif y > b
        length(edges)
    elseif y == a && (left || closed)
        1
    elseif y == b && (!left || closed)
        length(edges) - 1
    elseif left
        if isa(edges, AbstractRange)
            floor(Int, (y - a) / step(edges)) + 1
        else
            searchsortedlast(edges, y)
        end
    else
        if isa(edges, AbstractRange)
            ceil(Int, (y - a) / step(edges))
        else
            searchsortedfirst(edges, y) - 1
        end
    end
end

# requires: edges(o), midpoints(o), counts(o)
split_candidates(o::HistogramStat) = midpoints(o)
Statistics.mean(o::HistogramStat) = mean(midpoints(o), fweights(counts(o)))
Statistics.var(o::HistogramStat) = var(midpoints(o), fweights(counts(o)); corrected=true)
Statistics.std(o::HistogramStat) = sqrt(var(o))
Statistics.median(o::HistogramStat) = quantile(o, .5)
function Base.extrema(o::HistogramStat)
    x, y = midpoints(o), counts(o)
    x[findfirst(x -> x > 0, y)], x[findlast(x -> x > 0, y)]
end
function Statistics.quantile(o::HistogramStat, p = [0, .25, .5, .75, 1])
    x, y = midpoints(o), counts(o)
    inds = findall(x -> x != 0, y)
    quantile(x[inds], fweights(y[inds]), p)
end

#-----------------------------------------------------------------------# Hist
"""
    Hist(edges; left = true, closed = true)

Create a histogram with bin partition defined by `edges`.

- If `left`, the bins will be left-closed.
- If `closed`, the bin on the end will be closed.
    - E.g. for a two bin histogram ``[a, b), [b, c)`` vs. ``[a, b), [b, c]``

# Example

    o = fit!(Hist(-5:.1:5), randn(10^6))

    # approximate statistics
    using Statistics

    mean(o)
    var(o)
    std(o)
    quantile(o)
    median(o)
    extrema(o)
"""
struct Hist{T, R} <: HistogramStat{T}
    edges::R
    counts::Vector{Int}
    out::Vector{Int}
    left::Bool
    closed::Bool

    function Hist(edges::R, T::Type = eltype(edges); left::Bool=true, closed::Bool = true) where {R<:AbstractVector}
        new{T,R}(edges, zeros(Int, length(edges) - 1), [0,0], left, closed)
    end
end
nobs(o::Hist) = sum(o.counts) + sum(o.out)
value(o::Hist) = (x=o.edges, y=o.counts)

midpoints(o::Hist) = midpoints(o.edges)
counts(o::Hist) = o.counts
edges(o::Hist) = o.edges

function area(o::Hist)
    c = o.counts
    e = o.edges
    if isa(e, AbstractRange)
        return step(e) * sum(c)
    else
        return sum((e[i+1] - e[i]) * c[i] for i in 1:length(c))
    end
end

function pdf(o::Hist, y)
    i = binindex(o.edges, y, o.left, o.closed)
    if i < 1 || i > length(o.counts)
        return 0.0
    else
        return o.counts[i] / area(o)
    end
end

function _fit!(o::Hist, y)
    i = binindex(o.edges, y, o.left, o.closed)
    if 1 ≤ i < length(o.edges)
        o.counts[i] += 1
    else
        o.out[1 + (i > 0)] += 1
    end
end

function _merge!(o::Hist, o2::Hist)
    if o.edges == o2.edges
        for j in eachindex(o.counts)
            o.counts[j] += o2.counts[j]
        end
    else
        @warn("Histogram edges do not align.  Merging is approximate.")
        for (yi, wi) in zip(midpoints(o2.edges), o2.counts)
            for k in 1:wi
                _fit!(o, yi)
            end
        end
    end
end

#-----------------------------------------------------------------------# ExpandingHist
const ExpandableRange = Union{StepRange, StepRangeLen, LinRange}

mutable struct ExpandingHist{T, R <: StepRangeLen} <: HistogramStat{T}
    edges::R
    counts::Vector{Int}
    left::Bool
    n::Int
    function ExpandingHist(init::R, T::Type=Number; left::Bool = true) where {R <: ExpandableRange}
        new{T, R}(init, zeros(Int, length(init) - 1), left, 0)
    end
end
function ExpandingHist(b::Int; left::Bool=true)
    ExpandingHist(range(0, stop = 0, length = b + 1), Number; left=left)
end

midpoints(o::ExpandingHist) = midpoints(o.edges)
counts(o::ExpandingHist) = o.counts
edges(o::ExpandingHist) = o.edges

function Base.in(y, o::ExpandingHist)
    a, b = extrema(o.edges)
    o.left ? (a ≤ y < b) : (a < y ≤ b)
end

function _fit!(o::ExpandingHist, y)
    o.n += 1

    # init
    if nobs(o) == 1
        o.edges = range(y, stop=y, length=length(o.edges))
    elseif nobs(o) == 2
        a, b = extrema(o.edges)
        ey = eps(float(y))
        o.edges = range(min(a,y) - ey, stop=max(b, y) + ey, length=length(o.edges))
    end

    expand!(o, y)
    o.counts[binindex(o.edges, y, o.left, false)] += 1
end

function expand!(o::ExpandingHist, y)
    a, b = extrema(o.edges)
    if y > b  # find C such that y <= a + 2^C * (b - a)
        C = ceil(Int, log2((y - a) / (b - a)))
        o.edges = range(a, stop = a + 2^C * (b - a), length = length(o.edges))
    elseif y < a  # find C such that y >= b - 2^C * (b - a)
        C = ceil(Int, log2((b - y) / (b - a)))
        o.edges = range(b - 2^C * (b - a), stop = b, length = length(o.edges))
    end
end

# function adapt!(o::Hist, y)
#     y in o && return nothing

#     a, b = extrema(o.edges)
#     w = b - a

#     # number of widths to extend to the left
#     vl = (a - y) / w
#     n_widths_l = if vl < 0
#         0
#     elseif o.left || o.closed
#         ceil(Int, (a - y) / w)
#     else
#         v = (a - y) / w
#         ceil(Int, v) + isinteger(v)
#     end

#     # number of widths to extend to the right
#     vr = (y - b) / w
#     n_widths_r = if vr < 0
#         0
#     elseif !o.left || o.closed
#         ceil(Int, (y - b) / w)
#     else
#         v = (y - b) / w
#         ceil(Int, v) + isinteger(v)
#     end

#     for i in 1:n_widths_l
#         o.edges = extendleft(o.edges)
#         collapseright!(o.counts)
#     end
#     for i in 1:n_widths_r
#         o.edges = extendright(o.edges)
#         collapseleft!(o.counts)
#     end
#     nothing
# end

# function Base.in(y, o::Hist)
#     a, b = extrema(o.edges)
#     if o.left
#         return o.closed ? (a ≤ y ≤ b) : (a ≤ y < b)
#     else
#         return o.closed ? (a ≤ y ≤ b) : (a < y ≤ b)
#     end
# end

# function collapseleft!(x)
#     for i in eachindex(x)
#         j = 2i - 1
#         if j <= length(x)
#             x[i] = x[j]
#             if j < length(x)
#                 x[i] += x[j + 1]
#             end
#         else
#             x[i] = 0
#         end
#     end
# end
# function collapseright!(x)
#     for i in reverse(eachindex(x))
#         j = 2i - length(x)
#         if j >= 1
#             x[i] = x[j]
#             if j > 1
#                 x[i] += x[j - 1]
#             end
#         else
#             x[i] = 0
#         end
#     end
# end

# function extendleft(rng::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})
#     a, b = extrema(rng)
#     range(2a - b, stop = b, length = length(rng))
# end
# function extendright(rng::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})
#     a, b = extrema(rng)
#     range(a, stop = 2b - a, length = length(rng))
# end

# function extendleft(rng::StepRange)
#     a, b = extrema(rng)
#     (2a - b):step(rng):b
# end
# function extendright(rng::StepRange)
#     a, b = extrema(rng)
#     a:step(rng):(2b - a)
# end

# extendleft(rng) = error("Histogram cannot adapt with bin edges of type ", typeof(rng))
# extendright(rng) = extendleft(rng)

#-----------------------------------------------------------------------# KHist
struct KHistBin{T}
    loc::T
    count::Int
end
Base.isless(a::KHistBin, b::KHistBin) = isless(a.loc, b.loc)
function Base.merge(a::KHistBin, b::KHistBin)
    n = a.count + b.count
    KHistBin(smooth(a.loc, b.loc, b.count / n), n)
end
xy(o::KHistBin) = o.loc, o.count

@deprecate Hist(b::Int) KHist(b::Int)

"""
    KHist(k::Int)

Estimate the probability density of a univariate distribution at `k` approximately
equally-spaced points.

Ref: [http://www.jmlr.org/papers/volume11/ben-haim10a/ben-haim10a.pdf](http://www.jmlr.org/papers/volume11/ben-haim10a/ben-haim10a.pdf)

# Example

    o = fit!(KHist(25), randn(10^6))

    # Approximate statistics
    using Statistics
    mean(o)
    var(o)
    std(o)
    quantile(o)
    median(o)

    using Plots
    plot(o)
"""
struct KHist{T, E<:Extrema{T}} <: HistogramStat{Number}
    bins::Vector{KHistBin{T}}
    b::Int
    ex::E
end
KHist(b::Int, T::Type = Float64) = KHist(KHistBin{T}[], b, Extrema(T))

midpoints(o::KHist) = getfield.(o.bins, :loc)
counts(o::KHist) = getfield.(o.bins, :count)
edges(o::KHist) = vcat(minimum(o.ex), midpoints(getfield.(o.bins, :loc)), maximum(o.ex))

nobs(o::KHist) = isempty(o.bins) ? 0 : sum(x -> x.count, o.bins)

xy(o::KHist) = getfield.(o.bins, :loc), getfield.(o.bins, :count)

Base.extrema(o::KHist) = extrema(o.ex)

function value(o::KHist)
    x, y = xy(o)
    a, b = extrema(o.ex)
    if !isempty(x)
        if x[1] != a && y[1] != 1
            pushfirst!(x, a)
            pushfirst!(y, 1)
        end
        if x[end] != b && y[end] != 1
            push!(x, b)
            push!(y, 1)
        end
    end
    (x=x, y=y)
end

_fit!(o::KHist{T}, y) where {T} = push!(o, KHistBin(T(y), 1))

function Base.push!(o::KHist, y::KHistBin)
    fit!(o.ex, y.loc)
    insert!(o.bins, searchsortedfirst(o.bins, y), y)
    if length(o.bins) > o.b
        mindiff, i = Inf, 0
        for k in Base.OneTo(length(o.bins) - 1)
            diff = o.bins[k + 1].loc - o.bins[k].loc
            if diff < mindiff
                mindiff, i = diff, k
            end
        end
        o.bins[i] = merge(o.bins[i], o.bins[i + 1])
        deleteat!(o.bins, i + 1)
    end
end

function _merge!(a::KHist, b::KHist)
    merge!(a.ex, b.ex)
    for bin in b.bins
        push!(a, bin)
    end
end

function area(o::KHist, ind = length(o.bins) - 1)
    x, y = value(o)
    out = 0.0
    for i in Base.OneTo(ind)
        out += (x[i+1] - x[i]) * middle(y[i+1], y[i])
    end
    out
end

# based on linear interpolation
function pdf(o::KHist, x::Number)
    a, b = extrema(o.ex)
    firstbin = o.bins[1]
    lastbin = o.bins[end]

    # a ≤ firstbin.loc ≤ lastbin.loc ≤ b
    if firstbin.loc < x < lastbin.loc
        i = searchsortedfirst(o.bins, KHistBin(x, 0))
        x1, y1 = xy(o.bins[i - 1])
        x2, y2 = xy(o.bins[i])
        return smooth(y1, y2, (x - x1) / (x2 - x1)) / area(o)
    elseif a < x ≤ firstbin.loc
        x1, y1 = a, 1
        x2, y2 = xy(firstbin)
        return smooth(y1, y2, (x - x1) / (x2 - x1)) / area(o)
    elseif lastbin.loc ≤ x < b
        x1, y1 = xy(lastbin)
        x2, y2 = b, 1
        return smooth(y1, y2, (x - x1) / (x2 - x1)) / area(o)
    elseif x == a || x == b
        return 1 / area(o)
    else
        return 0.0
    end
end

function cdf(o::KHist, x::Number)
    a, b = extrema(o.ex)
    if x < a
        return 0.0
    elseif x == a
        return o.bins[1].count / area(o)
    elseif x ≥ b
        return 1.0
    else
        i = searchsortedfirst(o.bins, KHistBin(x, 0))
        x1, y1 = xy(o.bins[i - 1])
        x2, y2 = xy(o.bins[i])
        h = smooth(y1, y2, (x2 - x) / (x2 - x1))
        return (area(o, i - 2) + (x - x1) * h) / area(o)
    end
end


text: abstract type AbstractSampler end


text:
using Tables

sym(ptr) = ccall(:jl_symbol, Ref{Symbol}, (Ptr{UInt8},), ptr)

struct Query
    stmt::Stmt
    status::Base.RefValue{Cint}
    names::Vector{Symbol}
    types::Vector{Type}
    lookup::Dict{Symbol, Int}
end

struct Row <: Tables.AbstractRow
    q::Query
end

getquery(r::Row) = getfield(r, :q)

Tables.isrowtable(::Type{Query}) = true

Base.IteratorSize(::Type{Query}) = Base.SizeUnknown()
Base.eltype(q::Query) = Row

function reset!(q::Query)
    sqlite3_reset(q.stmt.handle)
    q.status[] = execute(q.stmt)
    return
end

function done(q::Query)
    st = q.status[]
    if st == SQLITE_DONE
        sqlite3_reset(q.stmt.handle)
        return true
    end
    st == SQLITE_ROW || sqliteerror(q.stmt.db)
    return false
end

function getvalue(q::Query, col::Int, ::Type{T}) where {T}
    handle = q.stmt.handle
    t = sqlite3_column_type(handle, col)
    if t == SQLITE_NULL
        return missing
    else
        TT = juliatype(t) # native SQLite Int, Float, and Text types
        return sqlitevalue(ifelse(TT === Any && !isbitstype(T), T, TT), handle, col)
    end
end

Tables.getcolumn(r::Row, ::Type{T}, i::Int, nm::Symbol) where {T} = getvalue(getquery(r), i, T)

Tables.getcolumn(r::Row, i::Int) = Tables.getcolumn(r, getquery(r).types[i], i, getquery(r).names[i])
Tables.getcolumn(r::Row, nm::Symbol) = Tables.getcolumn(r, getquery(r).lookup[nm])
Tables.columnnames(r::Row) = getquery(r).names

function Base.iterate(q::Query)
    done(q) && return nothing
    return Row(q), nothing
end

function Base.iterate(q::Query, ::Nothing)
    q.status[] = sqlite3_step(q.stmt.handle)
    done(q) && return nothing
    return Row(q), nothing
end

"Return the last row insert id from the executed statement"
DBInterface.lastrowid(q::Query) = last_insert_rowid(q.stmt.db)

"""
    DBInterface.prepare(db::SQLite.DB, sql::AbstractString)

Prepare an SQL statement given as a string in the sqlite database; returns an `SQLite.Stmt` compiled object.
See `DBInterface.execute`(@ref) for information on executing a prepared statement and passing parameters to bind.
A `SQLite.Stmt` object can be closed (resources freed) using `DBInterface.close!`(@ref).
"""
DBInterface.prepare(db::DB, sql::AbstractString) = Stmt(db, sql)

"""
    DBInterface.execute(db::SQLite.DB, sql::String, [params])
    DBInterface.execute(stmt::SQLite.Stmt, [params])

Bind any positional (`params` as `Vector` or `Tuple`) or named (`params` as `NamedTuple` or `Dict`) parameters to an SQL statement, given by `db` and `sql` or
as an already prepared statement `stmt`, execute the query and return an iterator of result rows.

Note that the returned result row iterator only supports a single-pass, forward-only iteration of the result rows.
Calling `SQLite.reset!(result)` will re-execute the query and reset the iterator back to the beginning.

The resultset iterator supports the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface, so results can be collected in any Tables.jl-compatible sink,
like `DataFrame(results)`, `CSV.write("results.csv", results)`, etc.
"""
function DBInterface.execute(stmt::Stmt, params=())
    status = execute(stmt, params)
    cols = sqlite3_column_count(stmt.handle)
    header = Vector{Symbol}(undef, cols)
    types = Vector{Type}(undef, cols)
    for i = 1:cols
        header[i] = sym(sqlite3_column_name(stmt.handle, i))
        types[i] = Union{juliatype(stmt.handle, i), Missing}
    end
    return Query(stmt, Ref(status), header, types, Dict(x=>i for (i, x) in enumerate(header)))
end

"""
    SQLite.createtable!(db::SQLite.DB, table_name, schema::Tables.Schema; temp=false, ifnotexists=true)

Create a table in `db` with name `table_name`, according to `schema`, which is a set of column names and types, constructed like `Tables.Schema(names, types)`
where `names` can be a vector or tuple of String/Symbol column names, and `types` is a vector or tuple of sqlite-compatible types (`Int`, `Float64`, `String`, or unions of `Missing`).

If `temp=true`, the table will be created temporarily, which means it will be deleted when the `db` is closed.
If `ifnotexists=true`, no error will be thrown if the table already exists.
"""
function createtable!(db::DB, nm::AbstractString, ::Tables.Schema{names, types}; temp::Bool=false, ifnotexists::Bool=true) where {names, types}
    temp = temp ? "TEMP" : ""
    ifnotexists = ifnotexists ? "IF NOT EXISTS" : ""
    typs = [types === nothing ? "BLOB" : sqlitetype(fieldtype(types, i)) for i = 1:length(names)]
    columns = [string(esc_id(String(names[i])), ' ', typs[i]) for i = 1:length(names)]
    return execute(db, "CREATE $temp TABLE $ifnotexists $nm ($(join(columns, ',')))")
end

"""
    source |> SQLite.load!(db::SQLite.DB, tablename::String; temp::Bool=false, ifnotexists::Bool=false, analyze::Bool=false)
    SQLite.load!(source, db, tablename; temp=false, ifnotexists=false, analyze::Bool=false)

Load a Tables.jl input `source` into an SQLite table that will be named `tablename` (will be auto-generated if not specified).

`temp=true` will create a temporary SQLite table that will be destroyed automatically when the database is closed
`ifnotexists=false` will throw an error if `tablename` already exists in `db`
`analyze=true` will execute `ANALYZE` at the end of the insert
"""
function load! end

load!(db::DB, table::AbstractString="sqlitejl_"*Random.randstring(5); kwargs...) = x->load!(x, db, table; kwargs...)

function load!(itr, db::DB, name::AbstractString="sqlitejl_"*Random.randstring(5); kwargs...)
    # check if table exists
    nm = esc_id(name)
    status = execute(db, "pragma table_info($nm)")
    rows = Tables.rows(itr)
    sch = Tables.schema(rows)
    return load!(sch, rows, db, nm, name, status == SQLITE_DONE; kwargs...)
end

checkdupnames(names) = length(unique(map(x->lowercase(String(x)), names))) == length(names) || error("duplicate case-insensitive column names detected; sqlite doesn't allow duplicate column names and treats them case insensitive")

function load!(sch::Tables.Schema, rows, db::DB, nm::AbstractString, name, shouldcreate; temp::Bool=false, ifnotexists::Bool=false, analyze::Bool=false)
    # check for case-insensitive duplicate column names (sqlite doesn't allow)
    checkdupnames(sch.names)
    # create table if needed
    shouldcreate && createtable!(db, nm, sch; temp=temp, ifnotexists=ifnotexists)
    # build insert statement
    params = chop(repeat("?,", length(sch.names)))
    stmt = Stmt(db, "INSERT INTO $nm VALUES ($params)")
    # start a transaction for inserting rows
    transaction(db) do
        for row in rows
            Tables.eachcolumn(sch, row) do val, col, _
                bind!(stmt, col, val)
            end
            sqlite3_step(stmt.handle)
            sqlite3_reset(stmt.handle)
        end
    end
    analyze && execute(db, "ANALYZE $nm")
    return name
end

# unknown schema case
function load!(::Nothing, rows, db::DB, nm::AbstractString, name, shouldcreate; temp::Bool=false, ifnotexists::Bool=false, analyze::Bool=false)
    state = iterate(rows)
    state === nothing && return nm
    row, st = state
    names = propertynames(row)
    sch = Tables.Schema(names, nothing)
    checkdupnames(sch.names)
    # create table if needed
    shouldcreate && createtable!(db, nm, sch; temp=temp, ifnotexists=ifnotexists)
    # build insert statement
    params = chop(repeat("?,", length(names)))
    stmt = Stmt(db, "INSERT INTO $nm VALUES ($params)")
    # start a transaction for inserting rows
    transaction(db) do
        while true
            Tables.eachcolumn(sch, row) do val, col, _
                bind!(stmt, col, val)
            end
            sqlite3_step(stmt.handle)
            sqlite3_reset(stmt.handle)
            state = iterate(rows, st)
            state === nothing && break
            row, st = state
        end
    end
    analyze && execute(db, "ANALYZE $nm")

    return name
end


text:
abstract AbstractMHRunner <: AbstractRunner #an abstract super type for both Standard and Generalized Metropolis-Hastings

const maxinitialcounter = 100::Int

_runner(p::MHRuntimePolicy,args...;keyargs...) = _runner(traittype(p.runner),p,args...;keyargs...)

_numsamplesinchain(::Type{Val{:burnin}},nburnin::Int,niterations::Int,nindicatorsamples::Int) = nburnin*nindicatorsamples,nburnin*nindicatorsamples
_numsamplesinchain(::Type{Val{:main}},nburnin::Int,niterations::Int,nindicatorsamples::Int) = niterations*nindicatorsamples,0
_numsamplesinchain(::Type{Val{:all}},nburnin::Int,niterations::Int,nindicatorsamples::Int) = (nburnin+niterations)*nindicatorsamples+1,nburnin*nindicatorsamples+1

_storeduring(phase::Symbol,p::MHRuntimePolicy) = in(traitvalue(p.store),[phase,:all])

function createcommon(runner_::AbstractMHRunner,tuner_::AbstractTuner,nparas::Int,nproposals::Int,nindicatorsamples::Int)
    p = runner_.policy
    nchainsamples,nburninsamples = _numsamplesinchain(traittype(p.store),runner_.numburnin,runner_.numiterations,nindicatorsamples)
    i = indicator(traitvalue(p.indicator),nproposals,nindicatorsamples,p.calculationtype)
    a = tunerstate(tuner_,runner_.numburnin,p.calculationtype)
    c = chain(traitvalue(p.chain),nparas,nchainsamples,p.sampletype,p.calculationtype;numburnin=nburninsamples)
    i,a,c
end

function initialize!(runner_::AbstractMHRunner,model_::AbstractModel,state_::AbstractSamplerState,chain_::AbstractChain,store::Bool)
    initialcounter = 0
    #the values is first generated in the proposals sample
    #then in prepare! function below, copied to the from sample
    sample_ = proposals(state_)
    #initialize the proposals field with values generated from the defaults or from the priors
    while ~isfinite((initialize!(runner_.policy.initialize,model_,sample_) ; geometry!(model_,sample_)).logprior[1])
        if (initialcounter+=1) > maxinitialcounter
            error("Problems generating initial proposal, giving up. Please check the parameter priors.")
        end
    end
    #prepare the distribution for the first proposal
    #this should result in the proposals field to be copied to the from field
    #plus other sampler-dependent setup
    prepare!(state_,true)
    store?store!(chain_,sample_,1):nothing
    state_
end


text:
# Basic list

abstract type List{T} end

struct ListNull{T} <: List{T} end

copy(l::ListNull{T}) where {T} = l

mutable struct ListNode{T} <: List{T}
    head::T
    tail::List{T}
end

copy(l::ListNode{T}) where {T} = ListNode{T}(l.head, l.tail)

head(l::ListNode) = l.head
tail(l::ListNode) = l.tail

# iterator
function Base.iterate(l::List, state=l)
    return typeof(state) <: ListNull ? nothing : (state.head, state.tail)
end
cons(value, l::List{T}) where T = ListNode{T}(value, l)

function length(l::List{T}) where T
    n = 0
    while typeof(l) != ListNull{T}
        n += 1
        l = l.tail
    end
    n
end

function cat(a::List{T}, b::List{T}) where T
    if a === nothing
        b
    else
        a = copy(a)
        u = a
        while u.tail !== nothing
            u.tail = copy(u.tail)
            u = u.tail
        end
        u.tail = b
        a
    end
end

function show(io::IO, a::List{T}) where T
    print(io, "List([")
    while typeof(a) != ListNull{T}
        print(io, a.head)
        if typeof(a.tail) != ListNull{T}
            print(io, ", ")
        end
        a = a.tail
    end
    print(io, "])")
end


text:
using Distributed

# to run in parallel set environment variable JULIA_NUM_THREADS to > 1
if haskey(ENV, "JULIA_NUM_THREADS")
    nw = parse(Int, ENV["JULIA_NUM_THREADS"])
    addprocs(nw; exeflags="--project")
end

using VortexCollisions

const nt = 10000

const M = 128
const N = 128

const ℳcut = 0
const 𝒩cut = 0

r(x,y) = sqrt(x^2 + y^2)
Θ(x,y) = atan(y,x)
r₀(x,y,ϵ) = 1 + ϵ * cos(2 * Θ(x,y))

function u_init_exp_r(x,y,ϵ=0.4)
    x̅ = (x-π)
    y̅ = (y-π)
    exp(- (r(x̅,y̅) / r₀(x̅,y̅,ϵ))^4)
end

gr = Grid2d(M,N)
ft = FourierTransform(gr, ℳcut=ℳcut, 𝒩cut=𝒩cut)
op = FokkerPlanckOperator(gr, ft)

# test
run_simulation(op, 10, 1E-6, u_init_exp_r,         "FokkerPlanck_exp_r_dt1E-6_nt10.h5")

# production
# run_simulation(op, nt, 1E-6, u_init_exp_r,         "FokkerPlanck_exp_r_dt1E-6_nt"        * string(nt) * ".h5")


text:
"""
    AffineDistribution(μ, σ, ρ)

A shifted and scaled (affinely transformed) version of `ρ`.

If ``Z`` is a random variable with distribution `ρ`, then `AffineDistribution(μ, σ, ρ)` is
the distribution of the random variable
```math
X = μ + σ Z
```

If `ρ` is a discrete univariate distribution, the probability mass function of the
transformed distribution is given by
```math
P(X = x) = P\\left(Z = \\frac{x-μ}{σ} \\right).
```

If `ρ` is a continuous univariate distribution with probability density function ``f_Z``,
the probability density function of the transformed distribution is given by
```math
f_X(x) = \\frac{1}{|σ|} f_Z\\left( \\frac{x-μ}{σ} \\right).
```

We recommend against using the `AffineDistribution` constructor directly. Instead, use
`+`, `-`, `*`, and `/`. These are optimized for specific distributions and will fall back
on `AffineDistribution` only when they need to.

Affine transformations of discrete variables are easily affected by rounding errors. If you
are getting incorrect results, try using exact `Rational` types instead of floats.

```julia
d = μ + σ * ρ       # Create location-scale transformed distribution
params(d)           # Get the parameters, i.e. (μ, σ, ρ)
```
"""
struct AffineDistribution{T<:Real, S<:ValueSupport, D<:UnivariateDistribution{S}} <: UnivariateDistribution{S}
    μ::T
    σ::T
    ρ::D
    function AffineDistribution{T,S,D}(μ::T, σ::T, ρ::D; check_args=true) where {T<:Real, S<:ValueSupport, D<:UnivariateDistribution{S}}
        check_args && @check_args(AffineDistribution, σ > zero(σ))
        new{T, S, D}(μ, σ, ρ)
    end
end

function AffineDistribution(μ::T, σ::T, ρ::UnivariateDistribution; check_args::Bool=true) where {T<:Real}
    _T = promote_type(eltype(ρ), T)
    D = typeof(ρ)
    S = value_support(D)
    return AffineDistribution{_T,S,D}(_T(μ), _T(σ), ρ; check_args=check_args)
end

function AffineDistribution(μ::Real, σ::Real, ρ::UnivariateDistribution; check_args::Bool=true)
    return AffineDistribution(promote(μ, σ)..., ρ; check_args=check_args)
end

# aliases
const LocationScale{T,S,D} = AffineDistribution{T,S,D}
function LocationScale(μ::Real, σ::Real, ρ::UnivariateDistribution; check_args::Bool=true)
    Base.depwarn("`LocationScale` is deprecated, use `AffineDistribution` instead", :LocationScale)
    if check_args && σ ≤ 0  # preparation for future PR where I remove σ > 0 check
        throw(ArgumentError("σ must be strictly positive."))
    end
    return AffineDistribution(μ, σ, ρ; check_args=false)
end

const ContinuousAffineDistribution{T<:Real,D<:ContinuousUnivariateDistribution} = AffineDistribution{T,Continuous,D}
const DiscreteAffineDistribution{T<:Real,D<:DiscreteUnivariateDistribution} = AffineDistribution{T,Discrete,D}

Base.eltype(::Type{<:AffineDistribution{T}}) where T = T

minimum(d::AffineDistribution) = d.μ + d.σ * minimum(d.ρ)
maximum(d::AffineDistribution) = d.μ + d.σ * maximum(d.ρ)
support(d::AffineDistribution) = affinedistribution_support(d.μ, d.σ, support(d.ρ))
function affinedistribution_support(μ::Real, σ::Real, support::RealInterval)
    return RealInterval(μ + σ * support.lb, μ + σ * support.ub)
end
affinedistribution_support(μ::Real, σ::Real, support) = μ .+ σ .* support

AffineDistribution(μ::Real, σ::Real, d::AffineDistribution) = AffineDistribution(μ + d.μ * σ, σ * d.σ, d.ρ)

#### Conversions

convert(::Type{AffineDistribution{T}}, μ::Real, σ::Real, ρ::D) where {T<:Real, D<:UnivariateDistribution} = AffineDistribution(T(μ),T(σ),ρ)
convert(::Type{AffineDistribution{T}}, d::AffineDistribution{S}) where {T<:Real, S<:Real} = AffineDistribution(T(d.μ),T(d.σ),d.ρ, check_args=false)

#### Parameters

location(d::AffineDistribution) = d.μ
scale(d::AffineDistribution) = d.σ
params(d::AffineDistribution) = (d.μ,d.σ,d.ρ)
partype(::AffineDistribution{T}) where {T} = T

#### Statistics

mean(d::AffineDistribution) = d.μ + d.σ * mean(d.ρ)
median(d::AffineDistribution) = d.μ + d.σ * median(d.ρ)
mode(d::AffineDistribution) = d.μ + d.σ * mode(d.ρ)
modes(d::AffineDistribution) = d.μ .+ d.σ .* modes(d.ρ)

var(d::AffineDistribution) = d.σ^2 * var(d.ρ)
std(d::AffineDistribution) = d.σ * std(d.ρ)
skewness(d::AffineDistribution) = skewness(d.ρ)
kurtosis(d::AffineDistribution) = kurtosis(d.ρ)

isplatykurtic(d::AffineDistribution) = isplatykurtic(d.ρ)
isleptokurtic(d::AffineDistribution) = isleptokurtic(d.ρ)
ismesokurtic(d::AffineDistribution) = ismesokurtic(d.ρ)

entropy(d::ContinuousAffineDistribution) = entropy(d.ρ) + log(d.σ)
entropy(d::DiscreteAffineDistribution) = entropy(d.ρ)

mgf(d::AffineDistribution,t::Real) = exp(d.μ*t) * mgf(d.ρ,d.σ*t)

#### Evaluation & Sampling

pdf(d::ContinuousAffineDistribution, x::Real) = pdf(d.ρ,(x-d.μ)/d.σ) / d.σ
pdf(d::DiscreteAffineDistribution, x::Real) = pdf(d.ρ,(x-d.μ)/d.σ)

logpdf(d::ContinuousAffineDistribution,x::Real) = logpdf(d.ρ,(x-d.μ)/d.σ) - log(d.σ)
logpdf(d::DiscreteAffineDistribution, x::Real) = logpdf(d.ρ,(x-d.μ)/d.σ)

for f in (:cdf, :ccdf, :logcdf, :logccdf)
    @eval $f(d::AffineDistribution, x::Real) = $f(d.ρ, (x - d.μ) / d.σ)
end

quantile(d::AffineDistribution,q::Real) = d.μ + d.σ * quantile(d.ρ,q)

rand(rng::AbstractRNG, d::AffineDistribution) = d.μ + d.σ * rand(rng, d.ρ)
cf(d::AffineDistribution, t::Real) = cf(d.ρ,t*d.σ) * exp(1im*t*d.μ)
gradlogpdf(d::ContinuousAffineDistribution, x::Real) = gradlogpdf(d.ρ,(x-d.μ)/d.σ) / d.σ

#### Syntactic sugar for simple transforms of distributions, e.g., d + x, d - x, and so on

Base.:+(d::UnivariateDistribution, x::Real) = AffineDistribution(x, one(x), d)
Base.:+(x::Real, d::UnivariateDistribution) = d + x
Base.:*(x::Real, d::UnivariateDistribution) = AffineDistribution(zero(x), x, d)
Base.:*(d::UnivariateDistribution, x::Real) = x * d
Base.:-(d::UnivariateDistribution, x::Real) = d + -x
Base.:/(d::UnivariateDistribution, x::Real) = inv(x) * d


text:
_isapproxzero(x::T, ztol) where {T<:Real} = x == zero(T)
_isapproxzero(x::T, ztol) where {T<:AbstractFloat} = -ztol < x < ztol
_isapproxzero(x::Vector{T}, ztol) where {T<:Real} = _isapproxzero(sum(abs, x), ztol)

_isapprox(x::Real, y::Real, ztol) = x == y
# I check with zero because isapprox(0, 1e-100) is false...
# but isapprox(1e-100, 2e-100) should be false
_isapprox(x::AbstractFloat, y::AbstractFloat, ztol) = (_isapproxzero(x, ztol) ? _isapproxzero(y, ztol) : (_isapproxzero(y, ztol) ? false : isapprox(x, y)))
_isapprox(x::Vector{S}, y::Vector{T}, ztol) where {S<:AbstractFloat, T<:AbstractFloat} = (x == zero(x) ? _isapproxzero(y, ztol) : (y == zero(y) ? _isapproxzero(x, ztol) : isapprox(x, y)))

_lt(x::T, y::T, ztol) where {T<:Real} = x < y
_lt(x::S, y::T, ztol) where {S<:Real,T<:Real} = _lt(promote(x, y)..., ztol)
_lt(x::AbstractFloat, y::AbstractFloat, ztol) = x < y && !_isapprox(x, y, ztol)


text:
# ------------------------------------------------------------------------------
# Chebyshev preconditioner
# ------------------------------------------------------------------------------

"""
```julia
Chebyshev(operator)
```

Chebyshev polynomial preconditioner for finite element operators.
    The Chebyshev semi-iterative method is applied to the matrix ``D^{-1} A``,
    where ``D^{-1}`` is the inverse of the operator diagonal.

# Arguments:
- `operator`: finite element operator to precondition

# Returns:
- Chebyshev preconditioner object

# Example:
```jldoctest
# setup
mesh = Mesh2D(1.0, 1.0);
mass = GalleryOperator("mass", 4, 4, mesh);

# preconditioner
chebyshev = Chebyshev(mass);

# verify
println(chebyshev)

# output
chebyshev preconditioner:
eigenvalue estimates:
  estimated minimum 0.2500
  estimated maximum 1.3611
estimate scaling:
  λ_min = a * estimated min + b * estimated max
  λ_max = c * estimated min + d * estimated max
  a = 0.0000
  b = 0.1000
  c = 0.0000
  d = 1.0000
```
"""
mutable struct Chebyshev <: AbstractPreconditioner
    # data never changed
    operator::Operator
    eigenvaluebounds::Array{Float64,1}

    # data empty until assembled
    eigenvalueestimates::Array{Float64,1}
    operatordiagonalinverse::AbstractArray{Float64}

    # inner constructor
    Chebyshev(operator::Operator) = new(operator, [0.0, 0.1, 0.0, 1.0])
end

# printing
# COV_EXCL_START
function Base.show(io::IO, preconditioner::Chebyshev)
    print(io, "chebyshev preconditioner:")

    # eigenvalue estimates
    print(io, "\neigenvalue estimates:")
    @printf(io, "\n  estimated minimum %.4f", preconditioner.eigenvalueestimates[1])
    @printf(io, "\n  estimated maximum %.4f", preconditioner.eigenvalueestimates[2])

    # estimate scaling
    print(io, "\nestimate scaling:")
    print(io, "\n  λ_min = a * estimated min + b * estimated max")
    print(io, "\n  λ_max = c * estimated min + d * estimated max")
    @printf(io, "\n  a = %.4f", preconditioner.eigenvaluebounds[1])
    @printf(io, "\n  b = %.4f", preconditioner.eigenvaluebounds[2])
    @printf(io, "\n  c = %.4f", preconditioner.eigenvaluebounds[3])
    @printf(io, "\n  d = %.4f", preconditioner.eigenvaluebounds[4])
end
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# data for computing symbols
# ------------------------------------------------------------------------------

"""
```julia
getoperatordiagonalinverse(preconditioner)
```

Compute or retrieve the inverse of the symbol matrix diagonal for a Chebyshev
    preconditioner

# Returns:
- Symbol matrix diagonal inverse for the operator

# Example:
```jldoctest
# setup
mesh = Mesh1D(1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);

# preconditioner
chebyshev = Chebyshev(diffusion)

# note: either syntax works
diagonalinverse = LFAToolkit.getoperatordiagonalinverse(chebyshev);
diagonalinverse = chebyshev.operatordiagonalinverse;

# verify
@assert diagonalinverse ≈ [6/7 0; 0 3/4]

# output

```
"""
function getoperatordiagonalinverse(preconditioner::Chebyshev)
    # assemble if needed
    if !isdefined(preconditioner, :operatordiagonalinverse)
        # retrieve diagonal and invert
        diagonalinverse = preconditioner.operator.diagonal^-1

        # store
        preconditioner.operatordiagonalinverse = diagonalinverse
    end

    # return
    return getfield(preconditioner, :operatordiagonalinverse)
end

"""
```julia
geteigenvalueestimates(preconditioner)
```

Compute or retrieve the eigenvalue estimates for a Chebyshev preconditioner

# Returns:
- Eigenvalue estimates for the operator

# Example:
```jldoctest
# setup
mesh = Mesh1D(1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);

# preconditioner
chebyshev = Chebyshev(diffusion)

# estimate eigenvalues
eigenvalueestimates = LFAToolkit.geteigenvalueestimates(chebyshev);

# verify
@assert eigenvalueestimates ≈ [0, 15/7]

# output

```
"""
function geteigenvalueestimates(preconditioner::Chebyshev)
    # assemble if needed
    if !isdefined(preconditioner, :eigenvalueestimates)
        dimension = preconditioner.operator.dimension
        λ_min = 1
        λ_max = 0
        θ_step = 2π / 4
        θ_range = -π/2:θ_step:3π/2-θ_step

        # compute eigenvalues
        if dimension == 1
            for θ_x in θ_range
                A = computesymbols(preconditioner.operator, [θ_x])
                eigenvalues = abs.(eigvals(preconditioner.operatordiagonalinverse * A),)
                λ_min = min(λ_min, eigenvalues...)
                λ_max = max(λ_max, eigenvalues...)
            end
        elseif dimension == 2
            for θ_x in θ_range, θ_y in θ_range
                A = computesymbols(preconditioner.operator, [θ_x, θ_y])
                eigenvalues = abs.(eigvals(preconditioner.operatordiagonalinverse * A),)
                λ_min = min(λ_min, eigenvalues...)
                λ_max = max(λ_max, eigenvalues...)
            end
        elseif dimension == 3
            for θ_x in θ_range, θ_y in θ_range, θ_z in θ_range
                A = computesymbols(preconditioner.operator, [θ_x, θ_y, θ_z])
                eigenvalues = abs.(eigvals(preconditioner.operatordiagonalinverse * A),)
                λ_min = min(λ_min, eigenvalues...)
                λ_max = max(λ_max, eigenvalues...)
            end
        end

        # store
        preconditioner.eigenvalueestimates = [λ_min, λ_max]
    end

    # return
    return getfield(preconditioner, :eigenvalueestimates)
end

"""
```julia
seteigenvalueestimatescaling(preconditioner, eigenvaluebounds)
```

Set the scaling of the eigenvalue estimates for a Chebyshev preconditioner

# Arguments:
- `eigenvaluebounds`: array of 4 scaling factors to use when setting ``\\lambda_{\\text{min}}``
    and ``\\lambda_{\\text{max}}`` based on eigenvalue estimates

``\\lambda_{\\text{min}}`` = a * estimated min + b * estimated max

``\\lambda_{\\text{max}}`` = c * estimated min + d * estimated max

# Example:
```jldoctest
# setup
mesh = Mesh1D(1.0);
diffusion = GalleryOperator("diffusion", 3, 3, mesh);

# preconditioner
chebyshev = Chebyshev(diffusion)

# set eigenvalue estimate scaling
# PETSc default is to use 1.1, 0.1 of max eigenvalue estimate
#   https://www.mcs.anl.gov/petsc/petsc-3.3/docs/manualpages/KSP/KSPChebyshevSetEstimateEigenvalues.html
seteigenvalueestimatescaling(chebyshev, [0.0, 0.1, 0.0, 1.1]);
println(chebyshev)

# output
chebyshev preconditioner:
eigenvalue estimates:
  estimated minimum 0.0000
  estimated maximum 2.1429
estimate scaling:
  λ_min = a * estimated min + b * estimated max
  λ_max = c * estimated min + d * estimated max
  a = 0.0000
  b = 0.1000
  c = 0.0000
  d = 1.1000
```
"""
function seteigenvalueestimatescaling(
    preconditioner::Chebyshev,
    eigenvaluebounds::Array{Float64,1},
)
    if length(eigenvaluebounds) != 4
        Throw(error("exactly four transformation arguments are required")) # COV_EXCL_LINE
    end

    preconditioner.eigenvaluebounds[1] = eigenvaluebounds[1]
    preconditioner.eigenvaluebounds[2] = eigenvaluebounds[2]
    preconditioner.eigenvaluebounds[3] = eigenvaluebounds[3]
    preconditioner.eigenvaluebounds[4] = eigenvaluebounds[4]
    preconditioner
end

# ------------------------------------------------------------------------------
# get/set property
# ------------------------------------------------------------------------------

function Base.getproperty(preconditioner::Chebyshev, f::Symbol)
    if f == :operatordiagonalinverse
        return getoperatordiagonalinverse(preconditioner)
    elseif f == :eigenvalueestimates
        return geteigenvalueestimates(preconditioner)
    else
        return getfield(preconditioner, f)
    end
end

function Base.setproperty!(preconditioner::Chebyshev, f::Symbol, value)
    if f == :operator
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    elseif f == :eigenvaluebounds
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(preconditioner, f, value)
    end
end

# ------------------------------------------------------------------------------
# compute symbols
# ------------------------------------------------------------------------------

"""
```julia
computesymbols(preconditioner, ω, θ)
```

Compute or retrieve the symbol matrix for a Chebyshev preconditioned operator

# Arguments:
- `preconditioner`: Chebyshev preconditioner to compute symbol matrix for
- `ω`:              Smoothing parameter array
                      [degree], [degree, ``\\lambda_{\\text{max}}``], or
                      [degree, ``\\lambda_{\\text{min}}``, ``\\lambda_{\\text{max}}``]
- `θ`:              Fourier mode frequency array (one frequency per dimension)

# Returns:
- Symbol matrix for the Chebyshev preconditioned operator

# Example:
```jldoctest
using LinearAlgebra

for dimension in 1:3
    # setup
    mesh = []
    if dimension == 1
        mesh = Mesh1D(1.0);
    elseif dimension == 2
        mesh = Mesh2D(1.0, 1.0);
    elseif dimension == 3
        mesh = Mesh3D(1.0, 1.0, 1.0);
    end
    diffusion = GalleryOperator("diffusion", 3, 3, mesh);

    # preconditioner
    chebyshev = Chebyshev(diffusion);

    # compute symbols
    A = computesymbols(chebyshev, [1], π*ones(dimension));

    # verify
    using LinearAlgebra;
    eigenvalues = real(eigvals(A));
    if dimension == 1
        @assert min(eigenvalues...) ≈ 0.15151515151515105
        @assert max(eigenvalues...) ≈ 0.27272727272727226
    elseif dimension == 2
        @assert min(eigenvalues...) ≈ -0.25495098334134725
        @assert max(eigenvalues...) ≈ -0.17128758445192374
    elseif dimension == 3
        @assert min(eigenvalues...) ≈ -0.8181818181818181
        @assert max(eigenvalues...) ≈ -0.357575757575757
    end
end

# output

```
"""
function computesymbols(preconditioner::Chebyshev, ω::Array, θ::Array)
    # validate number of parameters
    if length(ω) < 1
        Throw(error("at least one parameter required for Chebyshev smoothing")) # COV_EXCL_LINE
    elseif length(ω) > 3
        Throw(error("no more than three parameters allowed for Chebyshev smoothing")) # COV_EXCL_LINE
    end
    if (ω[1] % 1) > 1E-14 || ω[1] < 1
        Throw(error("first parameter must be degree of Chebyshev smoother")) # COV_EXCL_LINE
    end

    # get operator symbol
    A = computesymbols(preconditioner.operator, θ)

    # set eigenvalue estimates
    λ_min = 0
    λ_max = 0
    if length(ω) == 1
        λ_min = preconditioner.eigenvalueestimates[1]
        λ_max = preconditioner.eigenvalueestimates[2]
    elseif length(ω) == 2
        λ_max = ω[2]
    else
        λ_min = ω[2]
        λ_max = ω[3]
    end
    lower =
        λ_min * preconditioner.eigenvaluebounds[1] +
        λ_max * preconditioner.eigenvaluebounds[2]
    upper =
        λ_min * preconditioner.eigenvaluebounds[3] +
        λ_max * preconditioner.eigenvaluebounds[4]

    # compute Chebyshev smoother of given degree
    D_inv = preconditioner.operatordiagonalinverse
    D_inv_A = D_inv * A
    k = ω[1] # degree of Chebyshev smoother
    α = (upper + lower) / 2
    c = (upper - lower) / 2
    β_0 = -c^2 / (2 * α)
    γ_1 = -(α + β_0)
    E_0 = I
    E_1 = I - (1 / α) * D_inv_A * E_0
    E_n = I
    for _ = 2:k
        E_n = (D_inv_A * E_1 - α * E_1 - β_0 * E_0) / γ_1
        β_0 = (c / 2)^2 / γ_1
        γ_1 = -(α + β_0)
        E_0 = E_1
        E_1 = E_n
    end

    # return
    return E_1
end

# ------------------------------------------------------------------------------


text:
using ImagineInterface, Unitful, IntervalSets, Statistics
using Test

##################################BIDIRECTIONAL STACK########################################3
sample_rate = 50000*Unitful.s^-1
pmin = 0.0*Unitful.μm
pmax = 200.0*Unitful.μm
stack_time = 1.0*Unitful.s
posfwd, posback = gen_bidi_pos(pmin, pmax, stack_time, sample_rate)
exp_time = 0.011*Unitful.s
flash_frac = 0.1 #fraction of time to keep laser on during exposure
z_spacing = 3.1*Unitful.μm
z_pad = 5.0*Unitful.μm
nsamps_stack = ceil(Int, stack_time*sample_rate)

#offset by one sample going forward so that we don't use the end points of the triangle
delay1samp = 1/sample_rate
exp_dist = exp_time * abs(pmax - pmin) / stack_time

exp_intervals_fwd = spaced_intervals(posfwd, z_spacing, exp_time, sample_rate; delay=delay1samp, z_pad = max(z_pad, exp_dist), alignment=:start, rig="ocpi-2")
exp_nsamps = width(exp_intervals_fwd[1]) + 1
flash_nsamps = ImagineInterface.calc_num_samps(flash_frac * exp_time, sample_rate)
offset_nsamps = exp_nsamps - flash_nsamps
exp_intervals_back = map(x-> ClosedInterval(length(posfwd)-(maximum(x)+offset_nsamps-2), length(posfwd)-(maximum(x)-flash_nsamps-1)), exp_intervals_fwd)
las_intervals_fwd = map(x-> ClosedInterval(maximum(x)-flash_nsamps+1, maximum(x)), exp_intervals_fwd)
las_intervals_back = map(x-> ClosedInterval(maximum(x)-flash_nsamps+1, maximum(x)), exp_intervals_back)

samps_las_fwd = gen_pulses(nsamps_stack, las_intervals_fwd)
samps_las_back = gen_pulses(nsamps_stack, las_intervals_back)
samps_cam_fwd = gen_pulses(nsamps_stack, exp_intervals_fwd)
samps_cam_back = gen_pulses(nsamps_stack, exp_intervals_back)

@test all(map(IntervalSets.width, exp_intervals_fwd) .== 549)
@test all(map(IntervalSets.width, exp_intervals_back) .== 549)
#exposure duration
@test isapprox(IntervalSets.duration(exp_intervals_fwd[1]) / sample_rate, exp_time; atol=1/sample_rate)
#laser duration
@test isapprox(IntervalSets.duration(las_intervals_fwd[1]) / sample_rate, flash_frac*exp_time; atol=(1/sample_rate))
@test posfwd[1] == pmin
@test posfwd[end] < pmax
@test posback[1] == pmax
@test posback[end] > pmin
@test length(samps_cam_fwd) == nsamps_stack
#count pulses
nexp = length(exp_intervals_fwd)
@test nexp == length(exp_intervals_back) == 61
@test length(find_pulse_starts(samps_cam_fwd)) == nexp
@test length(find_pulse_stops(samps_cam_fwd)) == nexp
@test length(find_pulse_starts(samps_cam_back)) == nexp
@test length(find_pulse_stops(samps_cam_back)) == nexp
#TODO: test padding

d = gen_bidirectional_stack(pmin, pmax, z_spacing, stack_time, exp_time, sample_rate, flash_frac; z_pad = z_pad)
@test d["positioner"] == vcat(posfwd, posback)
@test d["camera"] == vcat(samps_cam_fwd, samps_cam_back)
@test d["laser"] == vcat(samps_las_fwd, samps_las_back)
@test all(find_pulse_stops(d["camera"]) .== find_pulse_stops(d["laser"]))
@test d["nframes"] == length(exp_intervals_fwd) * 2

cam = d["camera"]
las = d["laser"]
cam1 = cam[1:div(length(cam),2)]
cam2 = cam[div(length(cam),2)+1:end]
las1 = las[1:div(length(las),2)]
las2 = las[div(length(las),2)+1:end]
lexp1 = cam1.&las1 #overlap of laser with exposure
lexp2= reverse(cam2.&las2) #reversed to align with lexp1
@test all(lexp1[2:end].==lexp2[1:end-1]) #exclude endpoints of the triangle

#write it
ocpi2 = rigtemplate("ocpi-2"; sample_rate = sample_rate)
nstacks = 5
pos = first(getpositioners(ocpi2))
cam = first(getcameras(ocpi2))
las = getname(ocpi2, "488nm laser")
append!(pos, "pos", d["positioner"])
append!(cam, "cam", d["camera"])
append!(las, "las", d["laser"])
nframes = d["nframes"]
replicate!(pos, nstacks-1)
replicate!(cam, nstacks-1)
replicate!(las, nstacks-1)
outname = splitext(tempname())[1] *".json"
write_commands(outname, [cam;las;pos], nstacks, nframes, exp_time; isbidi = true)
#should throw an error when passed the wrong number of stacks/frames
@test_throws Exception write_commands(outname, [cam;las;pos], nstacks+1, nframes, exp_time; isbidi = true)
@test_throws Exception write_commands(outname, [cam;las;pos], nstacks, nframes+1, exp_time; isbidi = true)

#read it back in
_ocpi2 = parse_commands(outname)
lasers = getlasers(_ocpi2)
@test length(lasers) == 2
las_all = getname(lasers, "all lasers") #should have been added automatically
@test all(get_samples(las_all) .= true)

#smoothed bidi
smfwd, smbck = ImagineInterface.fit_smoothed_bidi(pmin, pmax, inv(2*stack_time), sample_rate; cutoff = 15.0*Unitful.Hz)
@test abs((maximum(smfwd) - pmax)) <= 0.01*Unitful.μm
@test abs((maximum(smbck) - pmax)) <= 0.01*Unitful.μm
@test abs((minimum(smfwd) - pmin)) <= 0.01*Unitful.μm
@test abs((minimum(smbck) - pmin)) <= 0.01*Unitful.μm
@test std(ustrip.(posfwd.-smfwd)) > 0.5 #make sure it's not just a triangle wave
@test std(ustrip.(posback.-smbck)) > 0.5

#slice_positions
ps = ImagineInterface.slice_positions(pmin, pmax, 5.0*Unitful.μm, 1.0*Unitful.μm)
#with 1um padding we have 198um total to place the slices, can fit 40 slices in a span of 395um
@test length(ps) == 40
#after centering expect an extra 3/2 = 1.5um padding
@test isapprox(first(ps), 2.5*Unitful.μm)
@test isapprox(last(ps), 197.5*Unitful.μm)
@test all(isapprox.(diff(ps), 5.0*Unitful.μm))

##################################UNIDIRECTIONAL STACK########################################
#set reset time equal to stack time, so the piezo waveform should be the same as in the bidi test, with half of the frames
posuni, posreset = gen_sawtooth(pmin, pmax, stack_time, stack_time, sample_rate)
@test posuni == posfwd
@test posreset == posback
flash_frac_ocpi1 = 1.1 #flashing per-exposure doesn't work well on ocpi1
d2 = gen_unidirectional_stack(pmin, pmax, z_spacing, stack_time, stack_time, exp_time, sample_rate, flash_frac_ocpi1; z_pad = z_pad)
@test length(posreset) == length(posback)
@test length(findall(x->x==1, diff(d2["camera"]))) == length(exp_intervals_fwd) #count pulses
@test d2["nframes"] == length(exp_intervals_fwd)

#write it
ocpi1 = rigtemplate("ocpi-1"; sample_rate = sample_rate)
nstacks = 5
pos = first(getpositioners(ocpi1))
cam = first(getcameras(ocpi1))
las = getname(ocpi1, "488nm laser shutter")
append!(pos, "pos", d2["positioner"])
append!(cam, "cam", d2["camera"])
append!(las, "las", d2["laser"])
nframes = d2["nframes"]
replicate!(pos, nstacks-1)
replicate!(cam, nstacks-1)
replicate!(las, nstacks-1)
outname = splitext(tempname())[1] *".json"
write_commands(outname, [cam;las;pos], nstacks, nframes, exp_time; isbidi = false)

#read it back in
_ocpi1 = parse_commands(outname)

##################################STEPPED UNIDIRECTIONAL STACK########################################
#set reset time equal to stack time, so the piezo waveform should be the same as in the bidi test, with half of the frames
pause_time = exp_time + 0.1 * Unitful.s
flash_frac_ocpi1 = 1.1 #flashing per-exposure doesn't work well on ocpi1
d2 = gen_stepped_stack(pmin, pmax, z_spacing, pause_time, stack_time, exp_time, sample_rate, flash_frac_ocpi1)

#write it
ocpi1 = rigtemplate("ocpi-1"; sample_rate = sample_rate)
nstacks = 5
pos = first(getpositioners(ocpi1))
cam = first(getcameras(ocpi1))
las = getname(ocpi1, "488nm laser shutter")
append!(pos, "pos", d2["positioner"])
append!(cam, "cam", d2["camera"])
append!(las, "las", d2["laser"])
nframes = d2["nframes"]
replicate!(pos, nstacks-1)
replicate!(cam, nstacks-1)
replicate!(las, nstacks-1)
outname = splitext(tempname())[1] *".json"
write_commands(outname, [cam;las;pos], nstacks, nframes, exp_time; isbidi = false)
#read it back in
_ocpi1 = parse_commands(outname)

########################### NO PIEZO MOTION ########################
pset = 0.0 * Unitful.μm
inter_exp_time = 0.0001 * Unitful.s #The time between exposure pulses
d3 = gen_2d_timeseries(pset, 10, exp_time, inter_exp_time, sample_rate, flash_frac)
@test all(d3["positioner"] .== pset)
@test length(findall(x->x==1, diff(d3["camera"]))) == 10
@test length(findall(x->x==1, diff(d3["laser"]))) == 10
@test (exp_time + inter_exp_time) * 10 * sample_rate ≈ length(d3["camera"])

max_forks_repo_name: HolyLab/ImagineInterface


text:
"""
    planar push block
        particle with contacts at each corner
"""
struct BimanipulationBlock{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    mass
    J
    μ
    g

    contact_corner_offset
	control_input_offset

    nq
    nu
    nc
    nf
    nb

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# Dimensions
nq = 3 # configuration dimension
nu = 2 * 4 + 2 * 2 # control dimension
nc = 4 # number of contact points
nf = 4 # number of faces for friction cone pyramid
nb = nc * nf
ns = 1

n = 2 * nq
m = nu + nc + nb + nc + nb + 1
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Kinematics
r = 0.1

# contact corner
cc1 = @SVector [r, r]
cc2 = @SVector [-r, r]
cc3 = @SVector [r, -r]
cc4 = @SVector [-r, -r]

# control location
# cu1 = @SVector [r, r]
# cu2 = @SVector [-r, r]
# cu3 = @SVector [-r, -r]
# cu4 = @SVector [r, -r]
cu5 = @SVector [r, 0.0]
cu6 = @SVector [0.0, r]
cu7 = @SVector [-r, 0.0]
cu8 = @SVector [0.0, -r]

contact_corner_offset = @SVector [cc1, cc2, cc3, cc4]
# control_input_offset = @SVector [cu1, cu2, cu3, cu4, cu5, cu6, cu7, cu8]
control_input_offset = @SVector [cu5, cu6, cu7, cu8]

# Parameters
μ = 1.0  # coefficient of friction
g = 9.81
mass = 1.0   # mass
J = 1.0 / 12.0 * mass * ((2.0 * r)^2 + (2.0 * r)^2)

# Methods
M_func(model::BimanipulationBlock, q) = Diagonal(@SVector [model.mass, model.mass, model.J])

function C_func(model::BimanipulationBlock, q, q̇)
	SVector{3}([0.0, model.mass * model.g, 0.0])
end

function rotation_matrix(x)
	SMatrix{2,2}([cos(x) -sin(x); sin(x) cos(x)])
end

function ϕ_func(model::BimanipulationBlock, q)
    p = view(q, 1:2)
	θ = q[3]
    R = rotation_matrix(θ)

	@SVector [(p + R * model.contact_corner_offset[1])[2],
              (p + R * model.contact_corner_offset[2])[2],
              (p + R * model.contact_corner_offset[3])[2],
              (p + R * model.contact_corner_offset[4])[2]]
end

function control_kinematics_func(model::BimanipulationBlock, q)
    p = q[1:2]
	θ = q[3]
    R = rotation_matrix(θ)

    SVector{8}([p + R * model.control_input_offset[1];
              p + R * model.control_input_offset[2];
              p + R * model.control_input_offset[3];
              p + R * model.control_input_offset[4]])
              # p + R * model.control_input_offset[5];
              # p + R * model.control_input_offset[6];
              # p + R * model.control_input_offset[7];
              # p + R * model.control_input_offset[8]])
end

function B_func(model::BimanipulationBlock, q)
	tmp(z) = control_kinematics_func(model, z)
	ForwardDiff.jacobian(tmp, q)
end

function N_func(model::BimanipulationBlock, q)
    tmp(z) = ϕ_func(model, z)
    ForwardDiff.jacobian(tmp, q)
end

function P_func(model::BimanipulationBlock, q)
    map = [1. 0.;
           0. 1.;
           -1. 0.;
           0. -1.]

    function p(x)
        pos = view(x, 1:2)
		θ = x[3]
        R = rotation_matrix(θ)

        [map * (pos + R * model.contact_corner_offset[1])[1:2];
         map * (pos + R * model.contact_corner_offset[2])[1:2];
         map * (pos + R * model.contact_corner_offset[3])[1:2];
         map * (pos + R * model.contact_corner_offset[4])[1:2]]
    end

    ForwardDiff.jacobian(p, q)
end

function friction_cone(model::BimanipulationBlock,u)
    λ = u[model.idx_λ]
    b = u[model.idx_b]

    @SVector [model.μ[1] * λ[1] - sum(b[1:4]),
              model.μ[2] * λ[2] - sum(b[5:8]),
              model.μ[3] * λ[3] - sum(b[9:12]),
              model.μ[4] * λ[4] - sum(b[13:16])]
end

function maximum_dissipation(model::BimanipulationBlock, x⁺, u, h)
    q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]

    ψ = u[model.idx_ψ]
    ψ_stack = [ψ[1] * ones(4);
               ψ[2] * ones(4);
               ψ[3] * ones(4);
               ψ[4] * ones(4)]

    η = u[model.idx_η]

    P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function lagrangian_derivatives(model, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

function fd(model::BimanipulationBlock{Discrete, FixedTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

	qm1 = 0.5 * (q1 + q2⁺)
    vm1 = (q2⁺ - q1) / h
    qm2 = 0.5 * (q2⁺ + q3)
    vm2 = (q3 - q2⁺) / h

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	[q2⁺ - q2⁻;
     (0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2
     + transpose(B_func(model, qm2)) * SVector{8}(u_ctrl[1:8])
     + transpose(N_func(model, q3)) * SVector{4}(λ)
     + transpose(P_func(model, q3)) * SVector{16}(b))]
end

model = BimanipulationBlock{Discrete, FixedTime}(n, m, d,
			mass, J, [[μ for i = 1:nc]..., μ], g,
			contact_corner_offset,
			control_input_offset,
            nq, nu, nc, nf, nb,
            idx_u,
            idx_λ,
            idx_b,
            idx_ψ,
            idx_η,
            idx_s)

function control_input(q, u)
	idx = [(i - 1) * 2 .+ (1:2) for i = 1:4]

    k = control_kinematics_func(model, q)
	k_input1 = u[model.idx_u][9:10]
	k_input2 = u[model.idx_u][11:12]

	d1 = norm(k[1:2] - k_input1)
	d2 = norm(k[3:4] - k_input1)
	d3 = norm(k[5:6] - k_input1)
	d4 = norm(k[7:8] - k_input1)

	e1 = norm(k[1:2] - k_input2)
	e2 = norm(k[3:4] - k_input2)
	e3 = norm(k[5:6] - k_input2)
	e4 = norm(k[7:8] - k_input2)

	_, min_idx1 = findmin([d1, d2, d3, d4])
	_, min_idx2 = findmin([e1, e2, e3, e4])

	return u[idx[min_idx1]], u[9:10], min_idx1, u[idx[min_idx2]], u[11:12], min_idx2
end

function visualize!(vis, model::BimanipulationBlock, q, u;
		r = r,
		u_mag = 1.0,
        Δt = 0.1)

	default_background!(vis)

    setobject!(vis["box"], GeometryBasics.Rect(Vec(-1.0 * r,
		-1.0 * r,
		-1.0 * r),
		Vec(2.0 * r, 2.0 * r, 2.0 * r)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    for i = 1:model.nc
        setobject!(vis["contact$i"], GeometryBasics.Sphere(Point3f0(0),
            convert(Float32, 0.02)),
            MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
    end

	for i = 1:4
        setobject!(vis["control$i"], GeometryBasics.Sphere(Point3f0(0),
            convert(Float32, 0.02)),
            MeshPhongMaterial(color = RGBA(51.0 / 255.0, 1.0, 1.0, 1.0)))
    end

	force_vis1 = ArrowVisualizer(vis[:force1])
	setobject!(force_vis1, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))

	force_vis2 = ArrowVisualizer(vis[:force2])
	setobject!(force_vis2, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))

	uf1, up1, min_idx1, uf2, up2, min_idx2 = control_input(q[2], u[1])
	uf1_norm = uf1 / u_mag
	uf2_norm = uf2 / u_mag

	settransform!(force_vis1,
				Point(up1[1] - uf1_norm[1], 0.0, up1[2] - uf1_norm[2]),
				Vec(uf1_norm[1], 0.0, uf1_norm[2]),
				shaft_radius=0.01,
				max_head_radius=0.025)

	settransform!(force_vis2,
				Point(up2[1] - uf2_norm[1], 0.0, up2[2] - uf2_norm[2]),
				Vec(uf2_norm[1], 0.0, uf2_norm[2]),
				shaft_radius=0.01,
				max_head_radius=0.025)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q)
    for t = 1:T-1
        MeshCat.atframe(anim, t) do
			if t < T-1
				uf1, up1, min_idx1, uf2, up2, min_idx2 = control_input(q[t+1], u[t])
				println("t = $t, control_input 1: $min_idx1, control_input 2: $min_idx2")

				if norm(uf1) < 1.0e-6
					setvisible!(vis[:force1], false)
				else
					setvisible!(vis[:force1], true)
					uf1_norm = uf1 / u_mag

					settransform!(force_vis1,
								Point(up1[1] - uf1_norm[1], 0.0, up1[2] - uf1_norm[2]),
								Vec(uf1_norm[1], 0.0, uf1_norm[2]),
								shaft_radius=0.01,
								max_head_radius=0.025)
				end

				if norm(uf2) < 1.0e-6
					setvisible!(vis[:force2], false)
				else
					setvisible!(vis[:force2], true)
					uf2_norm = uf2 / u_mag

					settransform!(force_vis2,
								Point(up2[1] - uf2_norm[1], 0.0, up2[2] - uf2_norm[2]),
								Vec(uf2_norm[1], 0.0, uf2_norm[2]),
								shaft_radius=0.01,
								max_head_radius=0.025)
				end

			end

            settransform!(vis["box"],
				compose(Translation(q[t+1][1], 0.0, q[t+1][2]), LinearMap(RotY(q[t+1][3]))))

            for i = 1:model.nc
                settransform!(vis["contact$i"],
                    Translation(([q[t+1][1]; 0.0; q[t+1][2]] + RotY(q[t+1][3]) * [contact_corner_offset[i][1]; 0.0; contact_corner_offset[i][2]])...))
            end

			for i = 1:4
				settransform!(vis["control$i"],
					Translation(([q[t+1][1]; 0.0; q[t+1][2]] + RotY(q[t+1][3]) * [control_input_offset[i][1]; 0.0; control_input_offset[i][2]])...))
			end
        end
    end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, -90.0, -1.0),LinearMap(RotZ(pi / 2.0))))
	setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 75)

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using SpecialMatrices

k = 4000
B = Block(k)
for i = 1:k
  B[i] = randn(2,2)
end
#BF = sparse(full(B))

x = randn(k*2,1)

B*x;
Profile.clear()
@time B*x;

#println("Sparse")
#@time BF*x;


text:
# Let's test latex commands to death... especially in light of #444

@testset "lx++1" begin
    # 444
    s = "\\newcommand{\\note}[1]{#1} \\note{A `B` C} D" |> fd2html_td
    @test isapproxstr(s, """

A B C D

        """)
end


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0402 - Implementation of Neural network workflows

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (the guarantee of reproducibility will however be lost)
## Pkg.resolve()
Pkg.instantiate()
using Random
Random.seed!(123)
ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

# We will _not_ run cross validation here to find the optimal hyper-parameters. The process will not be different than those we saw in the lesson on the Perceptron. Instead we focus on creating neural network models, train them based on data and evaluating their predictions.
# For feed-forward neural networks (both for classification and regression) we will use [BetaML](https://github.com/sylvaticus/BetaML.jl), while for Convolutional Neural Networks example we will use the [Flux.jl](https://github.com/FluxML/Flux.jl) package.

# ## Feed-forward neural networks

# ### Binary classification

# Data loading...
using BetaML, DelimitedFiles
data  = readdlm(joinpath(dirname(pathof(BetaML)),"..","test","data","binary2DData.csv"),'\t')
nR   = size(data,1)
idx  = shuffle(1:nR)
data = data[idx,:]
X    = copy(data[:,[2,3]])
y    = max.(0,convert(Array{Int64,1},copy(data[:,1]))) # Converting labels from {-1,1} to {0,1}
((xtrain,xtest),(ytrain,ytest)) = partition([X,y],[0.7,0.3])

# #### Using defaults - hidding complexity

# Model definition...
l1   = DenseLayer(2,5,f=tanh)
l2   = DenseLayer(5,5,f=relu)
l3   = DenseLayer(5,1,f=sigmoid)
mynn = buildNetwork([l1,l2,l3],squaredCost)

# Training...
train!(mynn,xtrain,ytrain)

ŷtrain         = predict(mynn, xtrain) |> makeColVector .|> round .|> Int
ŷtest          = predict(mynn, xtest)  |> makeColVector .|> round .|> Int
trainAccuracy  = accuracy(ŷtrain,ytrain)
testAccuracy   = accuracy(ŷtest,ytest)

# #### Specifying all options

# Model definition...
l1   = DenseLayer(2,5,f=tanh, df= dtanh,rng=copy(FIXEDRNG))
l2   = DenseLayer(5,5,f=relu,df=drelu,rng=copy(FIXEDRNG))
l3   = DenseLayer(5,1,f=sigmoid,df=dsigmoid,rng=copy(FIXEDRNG))
mynn = buildNetwork([l1,l2,l3],squaredCost,dcf=dSquaredCost,name="A classification task") # or crossEntropy / dCrossEntropy

# Training...

function myOwnTrainingInfo(nn,x,y;n,nBatches,epochs,verbosity,nEpoch,nBatch)
    if verbosity == NONE
        return false # doesn't stop the training
    end
    nMsgDict = Dict(LOW => 0, STD => 10,HIGH => 100, FULL => n)
    nMsgs = nMsgDict[verbosity]
    batchSize = size(x,1)
    if verbosity == FULL || ( nBatch == nBatches && ( nEpoch == 1  || nEpoch % ceil(epochs/nMsgs) == 0))

       ϵ = loss(nn,x,y)
       println("MY Training.. \t avg ϵ on (Epoch $nEpoch Batch $nBatch): \t $(ϵ)")
    end
    return false
 end
train!(mynn,xtrain,ytrain,epochs=300,batchSize=6,sequential=true,verbosity=STD,cb=myOwnTrainingInfo,optAlg=ADAM(η=t -> 0.001, λ=1.0, β₁=0.9, β₂=0.999, ϵ=1e-8),rng=copy(FIXEDRNG))

ŷtrain         = predict(mynn, xtrain) |> makeColVector .|> round .|> Int
ŷtest          = predict(mynn, xtest)  |> makeColVector .|> round .|> Int
trainAccuracy  = accuracy(ŷtrain,ytrain)
testAccuracy   = accuracy(ŷtest,ytest)

# ### Multinomial classification

# We want to determine the plant specie given some bothanic measures of the flower

iris     = readdlm(joinpath(dirname(Base.find_package("BetaML")),"..","test","data","iris.csv"),',',skipstart=1)
iris     = iris[shuffle(axes(iris, 1)), :] # Shuffle the records, as they aren't by default
x        = convert(Array{Float64,2}, iris[:,1:4])
y        = map(x->Dict("setosa" => 1, "versicolor" => 2, "virginica" =>3)[x],iris[:, 5]) # Convert the target column to numbers

((xtrain,xtest),(ytrain,ytest)) = partition([x,y],[0.8,0.2],shuffle=false)

ytrain_oh = oneHotEncoder(ytrain) # Convert to One-hot representation (e.g. 2 => [0 1 0], 3 => [0 0 1])

# Define the Artificial Neural Network model
l1   = DenseLayer(4,10,f=relu) # Activation function is ReLU
l2   = DenseLayer(10,3)        # Activation function is identity by default
l3   = VectorFunctionLayer(3,f=softmax) # Add a (parameterless) layer whose activation function (softMax in this case) is defined to all its nodes at once
mynn = buildNetwork([l1,l2,l3],crossEntropy,name="Multinomial logistic regression Model Sepal") # Build the NN and use the squared cost (aka MSE) as error function (crossEntropy could also be used)

# Training it (default to ADAM)
res = train!(mynn,scale(xtrain),ytrain_oh,batchSize=6) # Use optAlg=SGD() to use Stochastic Gradient Descent instead

# Test it
ŷtrain        = predict(mynn,scale(xtrain))   # Note the scaling function
ŷtest         = predict(mynn,scale(xtest))
trainAccuracy = accuracy(ŷtrain,ytrain)
testAccuracy  = accuracy(ŷtest,ytest,tol=1,ignoreLabels=false)

cm = ConfusionMatrix(ŷtrain,ytrain, labels=["setosa", "versicolor", "virginica"])
println(cm)

# ### Regression

# Data Loading and processing..
using Pipe, HTTP, CSV, Plots, DataFrames
urlData = "https://www4.stat.ncsu.edu/~boos/var.select/diabetes.tab.txt"
data = @pipe HTTP.get(urlData).body |> CSV.File(_, delim='\t') |> DataFrame
sex_oh = oneHotEncoder(data.SEX)
X = hcat(data.AGE, Matrix(data[:,3:10]),sex_oh)
y = data.Y
(xtrain,xval),(ytrain,yval) = partition([X,y],[0.8,0.2])

# Model definition...
l1   = DenseLayer(11,20,f=relu)
l2   = DenseLayer(20,20,f=relu)
l3   = DenseLayer(20,1,f=relu) # y is positive
mynn = buildNetwork([l1,l2,l3],squaredCost)

# Training...
trainingLogs = train!(mynn,scale(xtrain),ytrain,batchSize=6,epochs=600)

ŷtrain   = predict(mynn, scale(xtrain)) |> makeColVector
ŷval     = predict(mynn, scale(xval))  |> makeColVector
trainRME = meanRelError(ŷtrain,ytrain,normRec=false)
testRME  = meanRelError(ŷval,yval,normRec=false)
plot(trainingLogs.ϵ_epochs[10:end])
scatter(yval,ŷval,xlabel="obs",ylabel="est",legend=nothing)

# ## Convolutional neural networks

using LinearAlgebra, Statistics,Flux, MLDatasets, Plots

x_train, y_train = MLDatasets.MNIST.traindata(dir = "data/MNIST")
x_train          = permutedims(x_train,(2,1,3)); # For correct img axis
x_train          = convert(Array{Float32,3},x_train);
x_train          = reshape(x_train,(28,28,1,60000));
y_train          = Flux.onehotbatch(y_train, 0:9)
train_data       = Flux.Data.DataLoader((x_train, y_train), batchsize=128)
x_test, y_test   = MLDatasets.MNIST.testdata(dir = "data/MNIST")
x_test           = permutedims(x_test,(2,1,3)); # For correct img axis
x_test           = convert(Array{Float32,3},x_test);
x_test           = reshape(x_test,(28,28,1,10000));
y_test           = Flux.onehotbatch(y_test, 0:9)

model = Chain(
    ## 28x28 => 14x14
    Conv((5, 5), 1=>8, pad=2, stride=2, relu),
    ## 14x14 => 7x7
    Conv((3, 3), 8=>16, pad=1, stride=2, relu),
    ## 7x7 => 4x4
    Conv((3, 3), 16=>32, pad=1, stride=2, relu),
    ## 4x4 => 2x2
    Conv((3, 3), 32=>32, pad=1, stride=2, relu),
    ## Average pooling on each width x height feature map
    GlobalMeanPool(),
    Flux.flatten,
    Dense(32, 10),
    Flux.softmax
)

myaccuracy(ŷ, y) = (mean(Flux.onecold(ŷ) .== Flux.onecold(y)))
myloss(x, y)     = Flux.crossentropy(model(x), y)

opt = Flux.ADAM()
ps  = Flux.params(model)
number_epochs = 4

Flux.@epochs number_epochs Flux.train!(myloss, ps, train_data, opt)

ŷtrain =   model(x_train)
ŷtest  =   model(x_test)
myaccuracy(ŷtrain, y_train)
myaccuracy(ŷtest, y_test)

plot(Gray.(x_train[:,:,1,2]))
cm = ConfusionMatrix(Flux.onecold(ŷtest) .-1 ,Flux.onecold(y_test) .-1)
println(cm)

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


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

@doc raw"""
    ObjsExternalOrgMigrationsCurrent(;
        date_started=nothing,
        team_id=nothing,
    )

    - date_started::Int32
    - team_id::String
"""
mutable struct ObjsExternalOrgMigrationsCurrent <: SwaggerModel
    date_started::Any # spec type: Union{ Nothing, Int32 } # spec name: date_started
    team_id::Any # spec type: Union{ Nothing, String } # spec name: team_id

    function ObjsExternalOrgMigrationsCurrent(;date_started=nothing, team_id=nothing)
        o = new()
        validate_property(ObjsExternalOrgMigrationsCurrent, Symbol("date_started"), date_started)
        setfield!(o, Symbol("date_started"), date_started)
        validate_property(ObjsExternalOrgMigrationsCurrent, Symbol("team_id"), team_id)
        setfield!(o, Symbol("team_id"), team_id)
        o
    end
end # type ObjsExternalOrgMigrationsCurrent

const _property_map_ObjsExternalOrgMigrationsCurrent = Dict{Symbol,Symbol}(Symbol("date_started")=>Symbol("date_started"), Symbol("team_id")=>Symbol("team_id"))
const _property_types_ObjsExternalOrgMigrationsCurrent = Dict{Symbol,String}(Symbol("date_started")=>"Int32", Symbol("team_id")=>"String")
Base.propertynames(::Type{ ObjsExternalOrgMigrationsCurrent }) = collect(keys(_property_map_ObjsExternalOrgMigrationsCurrent))
Swagger.property_type(::Type{ ObjsExternalOrgMigrationsCurrent }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsExternalOrgMigrationsCurrent[name]))}
Swagger.field_name(::Type{ ObjsExternalOrgMigrationsCurrent }, property_name::Symbol) =  _property_map_ObjsExternalOrgMigrationsCurrent[property_name]

function check_required(o::ObjsExternalOrgMigrationsCurrent)
    (getproperty(o, Symbol("date_started")) === nothing) && (return false)
    (getproperty(o, Symbol("team_id")) === nothing) && (return false)
    true
end

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


text:
# Original matlab code from pmtk3 by Daniel Eaton (danieljameseaton@gmail.com)
# Ported by Levi Boyles
include("ars_defs.jl")
include("logging.jl")
include("probability_utils.jl")

#   ARS - Adaptive Rejection Sampling
#         sample perfectly & efficiently from a univariate log-concave
#         function
#
#   logpdf        a function handle to the log of a log-concave function.
#                  evaluated as: logpdf(x)
#
#   domain      the domain of logpdf. may be unbounded.
#                  ex: [1,inf], [-inf,inf], [-5, 20]
#
#   a,b         two points on the domain of logpdf, a<b. if domain is bounded
#                  then use a=domain[1], b=domain[2]. if domain is
#                  unbounded on the left, the derivative of logpdf for x=<a
#                  must be positive. if domain is unbounded on the right,
#                  the derivative for x>=b must be negative.
#
#                  ex: domain = [1,inf], a=1, b=20 (ensuring that
#                  logpdf'(x>=b)<0
#
#   nSamples    number of samples to draw
#
function ars(logpdf::Function, a::Float64, b::Float64, domain::Vector{Float64}, nSamples::Int)
    debug = false

    if domain[1] >= domain[2]
        error("invalid domain")
    end

    if a>=b || isinf(a) || isinf(b) || a<domain[1] || b>domain[2]
        error("invalid a and b")
    end

    numDerivStep = 1e-3 * (b-a)
    S = [a, a+numDerivStep, b-numDerivStep, b]

    # ensure the derivative there is positive
    if domain[1]==-Inf
        f1 = logpdf(S[1])
        f2 = logpdf(S[2])
        if logpdf(S[2]) <= logpdf(S[1])
            error("derivative at a must be positive, since the domain is unbounded to the left");
        end
    end

    # ensure the derivative there is negative
    if domain[2]==Inf
        f3 = logpdf(S[3])
        f4 = logpdf(S[4])
        if f4 >= f3
                error("derivative at b must be negative, since the domain is unbounded to the right");
        end
    end

    # initialize a mesh on which to create upper & lower hulls
    nInitialMeshPoints = 3;
    S = unique([S[1]; S[2]:(S[3]-S[2])/(nInitialMeshPoints+1):S[3]; S[4]]);
    fS = [logpdf(s) for s in S];

    lowerHull, upperHull = arsComputeHulls(S, fS, domain);

    nSamplesNow = 0;
    iterationNo = 1;

    samples = zeros(nSamples)

    while true

#            if debug
#                    figure(1); clf;
#                    arsPlot(upperHull, lowerHull, domain, S, fS, logpdf);
#            end

            # sample x from Hull
            x = arsSampleUpperHull( upperHull );

            lhVal, uhVal = arsEvalHulls( x, lowerHull, upperHull );

            U = rand();

            meshChanged = 0; # flag to indicate if a new point has been added to the mesh

            # three cases for acception/rejection
            if log(U)<=lhVal-uhVal
                    # accept, u is below lower bound
                    nSamplesNow = nSamplesNow + 1;
                    samples[nSamplesNow] = x;

            elseif log(U)<=logpdf(x )-uhVal
                    # accept, u is between lower bound and f
                    nSamplesNow = nSamplesNow + 1;
                    samples[nSamplesNow] = x;

                    meshChanged = 1;
            else
                    # reject, u is between f and upper bound
                    meshChanged = 1;

            end

            if meshChanged == 1
                    S = sort( [S; x] );
                    fS = [logpdf(s) for s in S];

                    lowerHull, upperHull = arsComputeHulls(S, fS, domain);
            end

            if debug
                    print("iteration $iterationNo, samples collected $nSamplesNow");
            end

            iterationNo = iterationNo + 1;

            if nSamplesNow==nSamples
                    break;
            end

    end

    samples
end

function arsComputeHulls(S, fS, domain)

    # compute lower piecewise-linear hull
    # if the domain of logpdf is unbounded to the left or right, then the lower
    # hull takes on -inf to the left or right of the end points of S

    lowerHull = Array(HullNode, length(S)-1)

    for li=1:length(S)-1
            lowerHull[li] = HullNode()
            lowerHull[li].m = (fS[li+1]-fS[li])/(S[li+1]-S[li]);
            lowerHull[li].b = fS[li] - lowerHull[li].m*S[li];
            lowerHull[li].left = S[li];
            lowerHull[li].right = S[li+1];
    end

    # compute upper piecewise-linear hull

    num_upper_segments = 2*(length(S)-2) + isinf(domain[1]) + isinf(domain[2])
    upperHull = Array(HullNode, num_upper_segments)

    i = 0
    if isinf(domain[1])
            # first line (from -infinity)
            m = (fS[2]-fS[1])/(S[2]-S[1]);
            b = fS[1] - m*S[1];
            #pr = 1/m*( exp(m*S[1]+b) - 0 ); # integrating in from -infinity
            pr = computeSegmentLogProb(-Inf, S[1], m, b)

            i += 1
            upperHull[i] = HullNode()
            upperHull[i].m = m;
            upperHull[i].b = b;
            upperHull[i].pr = pr;
            upperHull[i].left = -Inf;
            upperHull[i].right = S[1];

    end

    # second line
    m = (fS[3]-fS[2])/(S[3]-S[2]);
    b = fS[2] - m*S[2];
    #pr = 1/m * ( exp(m*S[2]+b) - exp(m*S[1]+b) );
    pr = computeSegmentLogProb(S[1], S[2], m, b)

    i += 1
    upperHull[i] = HullNode()
    upperHull[i].m = m;
    upperHull[i].b = b;
    upperHull[i].pr = pr;
    upperHull[i].left = S[1];
    upperHull[i].right = S[2];

    # interior lines
    # there are two lines between each abscissa
    for li=2:length(S)-2

            m1 = (fS[li]-fS[li-1])/(S[li]-S[li-1]);
            b1 = fS[li] - m1*S[li];

            m2 = (fS[li+2]-fS[li+1])/(S[li+2]-S[li+1]);
            b2 = fS[li+1] - m2*S[li+1];

            if !isfinite(m1) && !isfinite(m2)
                lprintln("S: $(S[li-1:li+2])")
                lprintln("fS: $(fS[li-1:li+2])")

                error("both hull slopes are infinite")
            end

            #ix = (b1-b2)/(m2-m1); # compute the two lines' intersection

            dx1 = S[li]-S[li-1]
            df1 = fS[li]-fS[li-1]
            dx2 = S[li+2]-S[li+1]
            df2 = fS[li+2]-fS[li+1]

            f1 = fS[li]
            f2 = fS[li+1]
            x1 = S[li]
            x2 = S[li+1]

            # more numerically stable than above
            #ix = (b1-b2)*dx1*dx2 / ( df2*dx1 - df1*dx2)
            ix = ((f1*dx1-df1*x1)*dx2 - (f2*dx2-df2*x2)*dx1) / ( df2*dx1 - df1*dx2)

            if !isfinite(m1) || abs(m1-m2) < 10.0^8 * eps(m1)
                ix = S[li]
                pr1 = -Inf
                #pr2 = 1.0/m2 * ( exp(m2*S[li+1]+b2) - exp(m2*ix+b2) );
                pr2 = computeSegmentLogProb(ix, S[li+1], m2, b2)
            elseif !isfinite(m2)
                ix = S[li+1]
                #pr1 = 1.0/m1 * ( exp(m1*ix+b1) - exp(m1*S[li]+b1) );
                pr1 = computeSegmentLogProb(S[li], ix, m1, b1)
                pr2 = -Inf
            else

                if !isfinite(ix)
                    lprintln("b1: $b1, b2: $b2, m1: $m1, m2: $m2")
                    error("Non finite intersection")
                end

                if abs(ix - S[li]) < 10.0^12 *eps(S[li])
                    ix = S[li]
                elseif abs(ix-S[li+1]) < 10.0^12 *eps(S[li+1])
                    ix = S[li+1]
                end

                if ix < S[li] || ix > S[li+1]
                    lprintln(fS[li-1:li+2])
                    lprintln(S[li-1:li+2])
                    lprintln("b1: $b1, b2: $b2, m1: $m1, m2: $m2")
                    lprintln("ix: $ix")
                    error("Intersection out of bounds -- logpdf is not concave")
                end

#                pr1 = 1.0/m1 * ( exp(m1*ix+b1) - exp(m1*S[li]+b1) );
#                pr2 = 1.0/m2 * ( exp(m2*S[li+1]+b2) - exp(m2*ix+b2) );
                pr1 = computeSegmentLogProb(S[li], ix, m1, b1)
                pr2 = computeSegmentLogProb(ix, S[li+1], m2, b2)
            end

            i += 1;
            upperHull[i] = HullNode()
            upperHull[i].m = m1;
            upperHull[i].b = b1;
            upperHull[i].pr = pr1;
            upperHull[i].left = S[li];
            upperHull[i].right = ix;

            i += 1;
            upperHull[i] = HullNode()
            upperHull[i].m = m2;
            upperHull[i].b = b2;
            upperHull[i].pr = pr2;
            upperHull[i].left = ix;
            upperHull[i].right = S[li+1];

    end

    # second last line
    m = (fS[end-1]-fS[end-2])/(S[end-1]-S[end-2]);
    b = fS[end-1] - m*S[end-1];
    #pr = 1.0/m * ( exp(m*S[end]+b) - exp(m*S[end-1]+b) );
    pr = computeSegmentLogProb(S[end-1], S[end], m ,b)

    i += 1
    upperHull[i] = HullNode()
    upperHull[i].m = m;
    upperHull[i].b = b;
    upperHull[i].pr = pr;
    upperHull[i].left = S[end-1];
    upperHull[i].right = S[end];

    if isinf(domain[2])
            # last line (to infinity)
            m = (fS[end]-fS[end-1])/(S[end]-S[end-1]);
            b = fS[end] - m*S[end];
            #pr = 1.0/m * ( 0 - exp(m*S[end]+b) );
            pr = computeSegmentLogProb(S[end], Inf, m, b)

            i += 1
            upperHull[i] = HullNode()
            upperHull[i].m = m;
            upperHull[i].b = b;
            upperHull[i].pr = pr;
            upperHull[i].left = S[end];
            upperHull[i].right = Inf;
    end

#    Z = sum([upperHull[i].pr for i = 1:length(upperHull)]);
#    #lprintln([upperHull[i].pr for i = 1:length(upperHull)]);
#    for li=1:length(upperHull)
#            upperHull[li].pr = upperHull[li].pr / Z;
#    end

    probs = exp_normalize([upperHull[i].pr for i = 1:length(upperHull)])
    for li=1:length(upperHull)
            upperHull[li].pr = probs[li]
    end

    (lowerHull, upperHull)
end

function computeSegmentLogProb(l, r, m, b)
    #1.0/m * (exp(m*r+b) - exp(m*l+b))
    if l == -Inf
        return -log(m) + m*r+b
    elseif r == Inf
        return -log(-m) + m*l+b
    else
        M = max(m*r+b, m*l+b)
        return -log(abs(m)) + log(abs(exp(m*r+b-M) - exp(m*l+b-M))) + M
    end
end

function arsSampleUpperHull( upperHull )

    cdf = cumsum([upperHull[i].pr for i = 1:length(upperHull)]);

    # randomly choose a line segment
    U = rand();
    li = 0
    for li = 1:length(upperHull)
            if( U<cdf[li] ) break;
            end
    end

    #lprintln(cdf)

    # sample along that line segment
    U = rand();

    m = upperHull[li].m;
    b = upperHull[li].b;
    left = upperHull[li].left;
    right = upperHull[li].right;

    M = max(m*right, m*left)

    #x = (log( U*(exp(m*right+b) - exp(m*left+b)) + exp(m*left+b) ) - b) / m ;
    x = (log( U*(exp(m*right-M) - exp(m*left-M)) + exp(m*left-M) ) + M) / m ;

    @assert x >= left && x <= right
    if isinf(x) || isnan(x)
            lprintln(upperHull[li])
            error("sampled an infinite or NaN x");
    end
    x
end

function arsEvalHulls( x, lowerHull, upperHull )

    lhVal = 0.0
    # lower bound
    if x<minimum([lowerHull[i].left for i = 1:length(lowerHull)])
            lhVal = -Inf;
    elseif x>maximum([lowerHull[i].right for i = 1:length(lowerHull)]);
            lhVal = -Inf;
    else
            for li=1:length(lowerHull)
                    left = lowerHull[li].left;
                    right = lowerHull[li].right;

                    if x>=left && x<=right
                            lhVal = lowerHull[li].m*x + lowerHull[li].b;
                            break;
                    end
            end
    end

    uhVal = 0.0
    # upper bound
    for li = 1:length(upperHull)
            left = upperHull[li].left;
            right = upperHull[li].right;

            if x>=left && x<=right
                    uhVal = upperHull[li].m*x + upperHull[li].b;
                    break;
            end
    end
    (lhVal, uhVal)
end

#function arsPlot(upperHull, lowerHull, domain, S, fS, logpdf)
#
#
#    Swidth = S[end]-S[1];
#    plotStep = Swidth/1000;
#    ext = 0.15*Swidth; # plot this much before a and past b, if the domain is infinite
#
#    left = S[1]; right = S[end];
#    if isinf(domain[1])
#        left = left - ext
#    end
#    if isinf(domain[2])
#        right = right + ext
#    end
#
#    x = left:plotStep:right;
#    fx = logpdf(x);
#
#    plot(x,fx, "k-"); hold on;
#    plot(S, fS, 'ko');
#    title('ARS');
#
#    # plot lower hull
#    for li=1:length(S)-1
#
#            m = lowerHull[li].m;
#            b = lowerHull[li].b;
#
#            x = lowerHull[li].left:plotStep:lowerHull[li].right;
#            plot( x, m*x+b, 'b-' );
#
#    end
#
#    # plot upper bound
#
#    # first line (from -infinity)
#    if isinf(domain[1])
#            x = (upperHull[1].right-ext):plotStep:upperHull[1].right;
#            m = upperHull[1].m;
#            b = upperHull[1].b;
#            plot( x, x*m+b, 'r-');
#    end
#
#    # middle lines
#    for li=2:length(upperHull)-1
#
#            x = upperHull[li].left:plotStep:upperHull[li].right;
#            m = upperHull[li].m;
#            b = upperHull[li].b;
#            plot( x, x*m+b, 'r-');
#
#    end
#
#    # last line (to infinity)
#    if isinf(domain[2])
#            x = upperHull[end].left:plotStep:(upperHull[end].left+ext);
#            m = upperHull[end].m;
#            b = upperHull[end].b;
#            plot( x, x*m+b, 'r-');
#    end
#


text:
struct ApproximateNonlinearSolver <: NonlinearCombinatorialAlgorithm
    subproblem_algorithm::CombinatorialAlgorithm
end

function solve(
    nli::NonlinearCombinatorialInstance,
    algo::ApproximateNonlinearSolver,
)
    # This technique is valid implemented for maximising. Use Holy traits
    # to support minimisation?
    @assert nli.combinatorial_structure.objective == Maximise()

    # TODO: what if the budgeted term in the objective doesn't have integer components?
    # TODO: what about ε?

    # Transform the enumerated nonlinearity as a true function.
    function sum_over_linear(x)
        return sum(
            nli.linear_coefficients[i] for
            i in eachindex(nli.linear_coefficients) if i in x
        )
    end
    function sum_over_nonlinear(x)
        return sum(
            nli.nonlinear_coefficients[i] for
            i in eachindex(nli.nonlinear_coefficients) if i in x
        )
    end
    if nli.nonlinear_function == Square
        nl_func = (x) -> sum_over_linear(x) + (sum_over_nonlinear(x))^2
    elseif nli.nonlinear_function == SquareRoot
        nl_func = (x) -> sum_over_linear(x) + sqrt(sum_over_nonlinear(x))
    else
        error("Nonlinear function not recognised: $(nli.nonlinear_function).")
    end

    # Get an upper bound on the budget.
    max_budget = _upper_bound_budget(nli)::Int

    # Make a linear instance with the linear coefficients
    # (`nli.combinatorial_structure` does not always have this property).
    li = copy(nli.combinatorial_structure, rewards=nli.linear_coefficients)
    bi = MinimumBudget(
        li,
        _round_coefficients(nli.nonlinear_coefficients, nli),
        max_budget,
    )

    # Compute all interesting budgeted solutions, depending on how it is requested to be done.
    solutions = Dict{Int, Vector}[]
    if nli.all_budgets_at_once
        bbi = copy(bi, compute_all_values=true)
        sols = solve(bbi, nli.linear_algo)::MultipleMinBudgetedSolution
        solutions = sols.solutions
    else
        for budget in 0:max_budget
            bbi = copy(bi, min_budget=budget)
            solutions[budget] = solve(bbi, nli.linear_algo).variables
        end
    end

    # Return the best solution.
    best_sol, _ = _pick_best_solution(solutions, nl_func)
    best_sol_dict = Dict(x => 1.0 for x in best_sol)
    return make_solution(nli.combinatorial_structure, best_sol_dict)
end

function _round_coefficients(x::Vector{Int}, ::NonlinearCombinatorialInstance)
    return x
end

function _round_coefficients(
    x::Vector{T},
    nli::NonlinearCombinatorialInstance,
) where {T <: Real}
    return Int[round(Int, v / nli.ε) for v in x]
end

function _round_coefficients(
    x::Dict{K, Int},
    ::NonlinearCombinatorialInstance,
) where {K}
    return x
end

function _round_coefficients(
    x::Dict{K, T},
    nli::NonlinearCombinatorialInstance,
) where {K, T <: Real}
    return Dict{K, Int}(k => round(Int, v / nli.ε) for (k, v) in x)
end

function _float_coefficients(x::Vector{T}) where {T}
    return Float64.(x)
end

function _float_coefficients(x::Dict{K, T}) where {K, T}
    return Dict{K, Float64}(k => v for (k, v) in x)
end

function _upper_bound_budget(nli::NonlinearCombinatorialInstance)
    # Make a linear instance to maximise the total weight.
    nlw = _float_coefficients(_round_coefficients(nli.nonlinear_coefficients, nli))
    li = copy(nli.combinatorial_structure, rewards=nlw)

    x = if fastest_exact(li) !== nothing
        solve(li, fastest_exact(li)).variables
    else
        solve(li, nli.linear_algo).variables
    end

    return Int(
        sum(
            nlw[i] for
            i in eachindex(nli.nonlinear_coefficients) if i in x
        ),
    )
end

function _pick_best_solution(
    solutions::Dict{Int, Vector{T}},
    nl_func::Function,
) where {T}
    best_solution = Int[]
    best_objective = -Inf
    for (budget, sol) in solutions
        # Ignore infeasible cases.
        # TODO: unify representation of infeasible solutions.
        if length(sol) == 0 || sol == [-1]
            continue
        end

        # Compute the maximum.
        f_x = nl_func(sol)
        if f_x > best_objective
            best_solution = sol
            best_objective = f_x
        end
    end

    return best_solution, best_objective
end