text:
#This file contains a list of the currently available mean functions

abstract type Mean end

# Calculate mean for matrix of observations
function mean(m::Mean, X::MatF64)
    nobs = size(X, 2)
    μ = Array{Float64}(undef, nobs)
    @inbounds for i in 1:nobs
        μ[i] = mean(m, view(X, :, i))
    end
    μ
end

# Calculates the stack [dm / dθᵢ] of mean matrix gradients
function grad_stack(m::Mean, X::MatF64)
    nobs = size(X, 2)
    mat = Array{Float64}(undef, nobs, num_params(m))
    @inbounds for i in 1:nobs
        mat[i, :] = grad_mean(m, view(X, :, i))
    end
    mat
end

#————————————————————————————————————————————

include("mZero.jl")          # Zero mean function
include("mConst.jl")         # Constant mean function
include("mLin.jl")           # Linear mean function
include("mPoly.jl")          # Polynomial mean function
include("composite_mean.jl") # Composite means


text:
module PlantHydraulics

using ConstrainedRootSolvers: NewtonBisectionMethod, ReduceStepMethodND,
      ResidualTolerance, SolutionTolerance, SolutionToleranceND, find_peak,
      find_zero
using DocStringExtensions: TYPEDEF, TYPEDFIELDS
using HypergeometricFunctions: _₂F₁
using PkgUtility: GAS_R, T_25, twarn, ρg_MPa
using SpecialFunctions: gamma
using Statistics: mean
using UnPack: @unpack
using WaterPhysics: relative_surface_tension, relative_viscosity

# export public types --- soil vulnerability
export BrooksCorey, VanGenuchten

# export public types --- xylem vulnerability
export LogisticSingle, PowerSingle, WeibullDual, WeibullSingle

# export public types --- pressure volume curve
export PVCurveLinear, PVCurveSegmented

# export public types --- flow mode
export NonSteadyStateMode, SteadyStateMode

# export public types --- hydraulic tissue
export LeafHydraulics, RootHydraulics, StemHydraulics

# export public types --- hydraulic system
export GrassLikeOrganism, PalmLikeOrganism, TreeLikeOrganism, TreeSimple

# export public functions --- initialize plant
export create_grass, create_palm, create_soil_VC, create_tree, fit_soil_VC!

# export public functions --- curves related
export vc_integral, p_from_volume, soil_erwc, soil_k_ratio_erwc,
       soil_k_ratio_p25, soil_k_ratio_rwc, soil_k_ratio_swc, soil_p_25_erwc,
       soil_p_25_rwc, soil_p_25_swc, soil_rwc, soil_swc, xylem_k_ratio,
       xylem_p_crit

# export public functions
export flow_profile!, pressure_profile!, inititialize_legacy!, critical_flow,
       xylem_risk, plant_conductances!, roots_flow!, xylem_flow, update_PVF!,
       temperature_effects!, end_pressure, fit_xylem_VC

include("types/curves.jl")
include("types/flow.jl"  )
include("types/organ.jl" )
include("types/plant.jl" )

include("initialize/legacy.jl")
include("initialize/soil.jl"  )
include("initialize/plant.jl" )

include("curves/capacity.jl")
include("curves/integral.jl")
include("curves/soil.jl"    )
include("curves/xylem.jl"   )

include("hydraulics/conductance.jl")
include("hydraulics/flow.jl"       )
include("hydraulics/pressure.jl"   )
include("hydraulics/temperature.jl")

include("hydraulics/capacitance.jl")

end # module

ext: jl

max_stars_repo_name: climate-machine/Land

max_forks_repo_name: climate-machine/Land


text:
function getPotentialResNet(nTh=2,T=1.0,nY=nTh,R=Float64)
    L = SingleLayer()
    RN = ResNN(L,R.(range(0.0,stop=T,length=nY)))
    N = NN([L;RN])
    return PotentialNN(N)
end

function initializeWeights(d,m,nTh,ar=x->x)
    (K1,b1) = (ar(0.01*randn(m,d+1)),ar(0.1*randn(m)))
    Θ1 = (K1,b1)
    (K2,b2) = (ar(0.01*repeat(randn(m,m),1,1,nTh)),ar(0.1*repeat(randn(m),1,nTh)))
    Θ2 = (K2,b2)
    ΘN = (Θ1, Θ2)
    (w0,A0,c0,z0) = (ones(m),ar(zeros(d+1,d+1)),ar(zeros(d+1)),ar(zeros(1)))
    return (w0,ΘN,A0,c0,z0)
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
@testset "Test sensitivity equation" begin
    # Problem Dimensions:
    # θ - size (6,); x - size (3,); y - size (2,);
    # dy/dx - size (2, 3); dy/dθ - size(2, 6).
    function test_fun(y, dx, dθ, x, θ)
        # Compute y
        y[1] = θ[1]*x[1] + θ[3]*x[2] + θ[5]*x[3]
        y[2] = θ[2]*x[1] + θ[4]*x[2] + θ[6]*x[3]
        # Compute df/dx
        dx[1, 1] = θ[1]
        dx[2, 1] = θ[2]
        dx[1, 2] = θ[3]
        dx[2, 2] = θ[4]
        dx[1, 3] = θ[5]
        dx[2, 3] = θ[6]
        # Compute df/dθ
        dθ[1, 1] = x[1]
        dθ[1, 2] = 0.0
        dθ[1, 3] = x[2]
        dθ[1, 4] = 0.0
        dθ[1, 5] = x[3]
        dθ[1, 6] = 0.0
        dθ[2, 1] = 0.0
        dθ[2, 2] = x[1]
        dθ[2, 3] = 0.0
        dθ[2, 4] = x[2]
        dθ[2, 5] = 0.0
        dθ[2, 6] = x[3]
        return
    end

    jacobian_buffer = zeros(2, 3)

    θ = Float64[1, 2, 3, 4, 5, 6]

    f = ms.sensitivity_equation(test_fun, jacobian_buffer)

    y = zeros(2)
    dydθ = zeros(2, 6)
    dydϕ = zeros(2, 6)
    x = [1, 2, 3]
    dxdθ = Matrix{Float64}(reshape(1:18, 3, 6))
    dxdϕ = Matrix{Float64}(reshape(1:18, 3, 6))
    k = 1
    y_extended = (y, dydθ, dydϕ)
    x_extended = (x, dxdθ, dxdϕ)
    f(y_extended, x_extended, θ)
    @test y == [22, 28]
    @test dydθ == [23.0 49.0 78.0 103.0 133.0 157.0;
                    28.0 65.0 100.0 138.0 172.0 211.0]
    @test dydϕ == [22.0 49.0 76.0 103.0 130.0 157.0;
                   28.0 64.0 100.0 136.0 172.0 208.0]
    θ .= θ + 2
    f(y_extended, x_extended, θ)
    @test y == [34, 40]
    @test dydθ == [35.0 79.0 126.0 169.0 217.0 259.0;
                   40.0 95.0 148.0 204.0 256.0 313.0]
    @test dydϕ == [34.0 79.0 124.0 169.0 214.0 259.0;
                   40.0 94.0 148.0 202.0 256.0 310.0]

    θ .= θ + 2
    f(y_extended, x_extended, θ)
    @test y == [46, 52]
    @test dydθ == [47.0 109.0 174.0 235.0 301.0 361.0;
                   52.0 125.0 196.0 270.0 340.0 415.0]
    @test dydϕ == [46.0 109.0 172.0 235.0 298.0 361.0;
                   52.0 124.0 196.0 268.0 340.0 412.0]
end


text:
"""
**Template Likelihood**

Template file for likelihood creation

```julia
TemplateLikelihood()
```
See all functions you need to implement
---

"""
struct TemplateLikelihood{T<:Real} <: Likelihood{T}
    θ::LatentArray{AbstractVector{T}}
    function TemplateLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function TemplateLikelihood{T}(θ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(θ)
    end
end

function TemplateLikelihood()
    TemplateLikelihood{Float64}()
end

function init_likelihood(likelihood::TemplateLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        TemplateLikelihood{T}([zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        TemplateLikelihood{T}()
    end
end

function pdf(l::TemplateLikelihood,y::Real,f::Real)

end

function Base.show(io::IO,model::TemplateLikelihood{T}) where T
    print(io,"Template Likelihood")
end

function compute_proba(l::TemplateLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    for i in 1:N
        pred[i]  = 0.0
    end
    return pred
end

### Local Updates Section ###

function local_updates!(model::VGP{T,<:TemplateLikelihood,<:AnalyticVI}) where {T}
end

function local_updates!(model::SVGP{T,<:TemplateLikelihood,<:AnalyticVI}) where {T}
end

function sample_local!(model::VGP{T,<:TemplateLikelihood,<:GibbsSampling}) where {T}
    return nothing
end

### Natural Gradient Section ###

function expec_μ(model::VGP{T,<:TemplateLikelihood,<:AnalyticVI},index::Integer) where {T}
end

function ∇μ(model::VGP{T,<:TemplateLikelihood}) where {T}
end

function expec_μ(model::SVGP{T,<:TemplateLikelihood,<:AnalyticVI},index::Integer) where {T}
end

function ∇μ(model::SVGP{T,<:TemplateLikelihood}) where {T}
end

function expec_Σ(model::AbstractGP{T,<:TemplateLikelihood,<:AnalyticVI},index::Integer) where {T}
    return model.likelihood.θ[index]
end

function ∇Σ(model::AbstractGP{T,<:TemplateLikelihood}) where {T}
    return model.likelihood.θ
end

### ELBO Section ###

function ELBO(model::AbstractGP{T,<:TemplateLikelihood,<:AnalyticVI}) where {T}
    return expecLogLikelihood(model) - GaussianKL(model)
end

function expecLogLikelihood(model::VGP{T,<:TemplateLikelihood,<:AnalyticVI}) where {T}
    tot = 0.0
    return tot
end

function expecLogLikelihood(model::SVGP{T,<:TemplateLikelihood,<:AnalyticVI}) wher {T}
    tot = 0.0
    return model.inference.ρ*tot
end

### Gradient Section ###

function gradpdf(::TemplateLikelihood,y::Int,f::T) where {T<:Real}
end

function hessiandiagpdf(::TemplateLikelihood,y::Int,f::T) where {T<:Real}
end


text:
using SequentialMonteCarlo
using RNGPool
using SMCExamples.LinearGaussian: LGTheta, Float64Particle, kalmanlogZ,
  defaultLGModel, makeLGModel
import MonteCarloMarkovKernels: simulateChainProgress, simulateChain,
  makeAMKernel, kde, estimateBM
using StaticArrays
using StatsBase
using LinearAlgebra
using Random
using Plots
Plots.gr()

setRNGs(0)
lgModel, theta, ys, ko = defaultLGModel(100)

const truex0 = theta.x0
const truev0 = theta.v0
const trueC = theta.C

@inline function toLGTheta(v::SVector{3, Float64})
  return LGTheta(v[1], v[2], trueC, v[3], truex0, truev0)
end

t0 = SVector{3, Float64}(theta.A, theta.Q, theta.R)
const sigmaProp = SMatrix{3, 3, Float64}(Matrix{Float64}(I, 3, 3))

@inline function lglogprior(theta::LGTheta)
  if theta.A < 0 || theta.A > 10 return -Inf end
  if theta.Q < 0 || theta.Q > 10 return -Inf end
  if theta.R < 0 || theta.R > 10 return -Inf end
  return 0.0
end

function makelgsmcltd(ys::Vector{Float64}, N::Int64, nthreads::Int64)
  smcio = SMCIO{Float64Particle, Nothing}(N, length(ys), nthreads, false, 0.5)
  function ltd(in::SVector{3, Float64})
    theta::LGTheta = toLGTheta(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    model::SMCModel = makeLGModel(theta, ys)
    smc!(model, smcio)
    return smcio.logZhats[length(ys)]
  end
end

function makelgkalmanltd(ys::Vector{Float64})
  function ltd(in::SVector{3, Float64})
    theta::LGTheta = toLGTheta(in)
    lp::Float64 = lglogprior(theta)
    if lp == -Inf return -Inf end
    return kalmanlogZ(theta, ys)
  end
end

logtargetSMC = makelgsmcltd(ys, 1024, Threads.nthreads())
logtargetKalman = makelgkalmanltd(ys)

PSMC = makeAMKernel(logtargetSMC, sigmaProp)
PKalman = makeAMKernel(logtargetKalman, sigmaProp)

Random.seed!(12345)
chainSMC = simulateChainProgress(PSMC, t0, 2^15)
sar = PSMC(:acceptanceRate)

@time chainKalman = simulateChain(PKalman, t0, 2^20)
kar = PKalman(:acceptanceRate)

savefigures = false

vsKalman = (i->(x->x[i]).(chainKalman)).(1:3)
vsSMC =  (i->(x->x[i]).(chainSMC)).(1:3)

plot(kde(vsSMC[1], sar))
plot!(kde(vsKalman[1], kar))
savefigures && savefig("pmmh_kde1.png")

plot(kde(vsSMC[2], sar))
plot!(kde(vsKalman[2], kar))
savefigures && savefig("pmmh_kde2.png")

plot(kde(vsSMC[3], sar))
plot!(kde(vsKalman[3], kar))
savefigures && savefig("pmmh_kde3.png")

contour(kde(vsSMC[1], vsSMC[2], sar))
contour!(kde(vsKalman[1], vsKalman[2], kar))
savefigures && savefig("pmmh_kde12.png")

contour(kde(vsSMC[1], vsSMC[3], sar))
contour!(kde(vsKalman[1], vsKalman[3], kar))
savefigures && savefig("pmmh_kde13.png")

contour(kde(vsSMC[2], vsSMC[3], sar))
contour!(kde(vsKalman[2], vsKalman[3], kar))
savefigures && savefig("pmmh_kde23.png")

plot(autocor(vsSMC[1]))
plot!(autocor(vsKalman[1]))
savefigures && savefig("pmmh_acf1.png")
plot(autocor(vsSMC[2]))
plot!(autocor(vsKalman[2]))
savefigures && savefig("pmmh_acf2.png")
plot(autocor(vsSMC[3]))
plot!(autocor(vsKalman[3]))
savefigures && savefig("pmmh_acf3.png")

estimateBM.(vsSMC)
estimateBM.(vsKalman)


text:
# Use map_func to avoid float precision differences between the same vertex
function _push_cell_vertex!(vert,cell_verts, used_vertices,vertices,nextvert, map_func)
  token = ht_keyindex2!(used_vertices, vert)
  if token > 0 # reuse dofs
      reuse_vert = used_vertices.vals[token]
      push!(cell_verts, reuse_vert)
  else # token <= 0, use new vertex
      Base._setindex!(used_vertices, nextvert, vert, -token)
      push!(cell_verts, nextvert)
      push!(vertices, map_func(vert))
      nextvert += 1
  end
  return nextvert
end

function _generate_2d_hex_centroids!(centroids::Vector{Tensors.Vec{2,T}},LL, n_centroid_rows, n_centroid_cols, hex_width, hex_heigth) where {T}
  x_coord = LL[1]; y_coord = LL[2]
  x_lim = n_centroid_cols*hex_width
  y_lim = n_centroid_rows*hex_heigth
  ss = 1
  for j in 1:(2*n_centroid_rows)
      for i in 1:(n_centroid_cols+1)
          centroid = Tensors.Vec{2,T}((x_coord,y_coord))
          push!(centroids, centroid)
          (x_coord + hex_width) > x_lim ? break : x_coord = x_coord + hex_width
      end
      (ss > 0) ? x_coord = LL[1] + hex_width/2 : x_coord = LL[1]
      ss = ss*(-1)
      (y_coord + hex_heigth/4) > y_lim ? break : y_coord = y_coord + 3/4*hex_heigth;
  end
  return centroids
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
"""
A gate will open while showing a card.
Otherwise, the gate is locked.
"""

abstract type State end

struct Locked <: State
end

isopen(::Locked) = false

struct Unlocked <: State
end

isopen(::Unlocked) = true

mutable struct Gate
    state::State
    Gate() = new(Locked())
end

isopen(g::Gate) = isopen(g.state)
lock!(g::Gate) = (g.state = Locked())
unlock!(g::Gate) = (g.state = Unlocked())

function open(g::Gate)
    if g.state isa Unlocked
        println("The gate is open.")
    end
end


text:
# On the other hand, it might be wise to try a different strategy: let the
# giant squid win.
#
# You aren't sure how many bingo boards a giant squid could play at once, so
# rather than waste time counting its arms, the safe thing to do is to figure
# out which board will win last and choose that one. That way, no matter which
# boards it picks, it will win for sure.
#
# In the above example, the second board is the last to win, which happens
# after 13 is eventually called and its middle column is completely marked. If
# you were to keep playing until this point, the second board would have a sum
# of unmarked numbers equal to 148 for a final score of 148 * 13 = 1924.
#
# Figure out which board will win last. Once it wins, what would its final
# score be?

# read input
open("input.txt", "r") do io
  # read numbers
  global numbers = parse.(Int, split(readline(io), ','))

  # read boards
  global boards = []
  while !eof(io)
    # blank line
    readline(io)

    board = Array{Int, 2}(undef, 5, 5)
    @inbounds for x in 1:5
      board[x, :] = parse.(Int, split(strip(readline(io)), r"\s+"))
    end
    push!(boards, board)
  end
end

# `results` keeps track of how many rows have matched, and how many columns
# have matched, per board. When either of those counts reach 5, bingo! We also
# keep track of the sum of matched numbers per board. Finally, keep track of
# which boards have won, what the last winning board was, and the last winning
# number.
numboards = length(boards)
results = zeros(UInt8, numboards, 10)
sums = zeros(Int, numboards)
done = falses(numboards)
numdone = 0
lastwin = 0
lastwinnum = 0
for num in numbers, b in 1:numboards
  done[b] && continue
  idx = findfirst(isequal(num), boards[b])
  if !isnothing(idx)
    sums[b] += num
    if (results[b, idx[1]] += 1) == 5 || (results[b, idx[2] + 5] += 1) == 5
      done[b] = true
      global numdone += 1
      global lastwin = b
      global lastwinnum = num
      numdone == numboards && break
    end
  end
end

unmarkedsum = sum(boards[lastwin])
println("Bingo!: ", (unmarkedsum - sums[lastwin]) * lastwinnum)


text:
################################################################################
#                                ____ ____ _  _                                #
#                                |  | |__| |\/|                                #
#                                |_\| |  | |  |                                #
################################################################################

type QAM <: Modulation
    M::Int                              # Number of symbols in the constellation
    m::Int                              # Bits per symbol
    α::Real                             # Symbol scaling factor
    β::Real                             # Symbol offset
    constellation::Vector               # Constellation of symbols in original
    grayConstellation::Vector           # Constellation of symbols inn gray-code order
    bitsMap::Vector                     # bitsMap[i] maps to constellation[i]
end

function QAM( M::Integer )
    isinteger(sqrt(M)) || error( "sqrt(M) must be an integer value" )
    m                 = int(log2( M ))
    width             = int(sqrt(M))
    constellation     = Array( Complex64, M )
    grayConstellation = Array( Complex64, M )
    bitsMap           = Array( Int, M )
    α                 = 2/sqrt( (M-1) * 2/3 ) # scaling factor
    β                 = (-width+1)/2

    # http://www.gaussianwaves.com/2012/10/constructing-a-rectangular-constellation-for-16-qam/

    for idx in 0:M-1
        inPhase                    = idx >> div( m, 2 )
        quadrature                 = idx & int(exp2(m/2)-1)
        constellation[idx + 1]     = α * Complex( inPhase+β, quadrature+β )
        bitsMap[idx + 1]           = (encode(Gray, inPhase) << div( m, 2 )) | encode(Gray, quadrature)
        inPhase                    = decode(Gray, inPhase)
        quadrature                 = decode(Gray, quadrature)
        grayConstellation[idx + 1] = α * Complex( inPhase+β, quadrature+β )
    end

    QAM( M, m, α, β, constellation, grayConstellation, bitsMap )
end

function Base.show( io::IO, modem::QAM )
    @printf( io, "QAM{%d}\n", modem.m )
end

################################################################################
#                                _  _ ____ ___                                 #
#                                |\/| |  | |  \                                #
#                                |  | |__| |__/                                #
################################################################################

function modulate( modem::QAM, bits::Integer )
    modem.grayConstellation[ bits+1 ]
end

function modulate( modem::QAM, data::Vector )
    [ modulate( modem, datum ) for datum in data ]
end

################################################################################
#                           ___  ____ _  _ ____ ___                            #
#                           |  \ |___ |\/| |  | |  \                           #
#                           |__/ |___ |  | |__| |__/                           #
################################################################################

function demodulate( modem::QAM, symbol::Complex )
    symbol     = symbol / modem.α
    inPhase    = int( real( symbol ) - modem.β )
    quadrature = int( imag( symbol ) - modem.β )
    bits       = ( inPhase << div( modem.m, 2 )) | quadrature
    bits       = clamp( bits, 0, modem.M-1 )
    bits       = modem.bitsMap[bits+1]
end

function demodulate( modem::QAM, symbols::Vector )
    Int[ demodulate( modem, symbol ) for symbol in symbols ]
end

################################################################################
#                      ____ _  _ ____ _  _ ___  _    ____                      #
#                      |___  \/  |__| |\/| |__] |    |___                      #
#                      |___ _/\_ |  | |  | |    |___ |___                      #
################################################################################

#=

M         = 16
modem     = QAM( M )
modData   = rand(0:M-1, M*4)
symbols   = modulate( modem, modData)
symbols   = symbols .* 1.05
demodData = demodulate( modem, symbols )

Test.@test_approx_eq modData demodData

using PyPlot

function PyPlot.plot( modem::QAM )
    fig = figure()
    scatter( real(modem.constellation), imag(modem.constellation) )

    for i in 0:modem.M-1
        symbol = modem.constellation[ i + 1 ]
        x      = real( symbol )
        y      = imag( symbol )
        bits   = modem.bitsMap[i+1]
        strVal = modem.M > 16 ? string(bits) : bin(bits)
        annotate( strVal, (x, y+0.05), ha="center" )
    end

    setp( axes(), aspect = 1  )
    xlabel( "In Phase" )
    ylabel( "Quadrature" )

    return fig
end

plot( modem )

=#


text:
using ExponentialAction
using Test

@testset "ExponentialAction.jl" begin
    include("util.jl")
    include("expv.jl")
    include("expv_autodiff.jl")
end


text:
using Graphene
using Test
using GeometricalPredicates:Point2D

test_xy = [0;
           17.4]
test_group = Set.([[1]])
test_atom_groups = collect_atom_groups(test_xy)
@test test_group == Set.(test_atom_groups)

test_xy = [0     5;
           17.4  8.7]
test_group = Set.([[1, 2], [2, 1]])
test_atom_groups = collect_atom_groups(test_xy)
@test test_group == Set.(test_atom_groups)

test_xy = [0     35;
           17.4  8.7]
test_group = Set.([[1],[2]])
test_atom_groups = collect_atom_groups(test_xy)
@test test_group == Set.(test_atom_groups)

test_xy = [0     5     5   15    15   20   20   30   30   35;
           17.4  8.7  26.1  8.7  26.1  0   17.4  0   17.4  8.7]
test_group = Set.([[1, 2, 3],
                    [2, 1, 4],
                    [3, 1, 5],
                    [4, 2, 6, 7],
                    [5, 3, 7],
                    [6, 4, 8],
                    [7, 4, 5, 9],
                    [8, 6, 10],
                    [9, 7, 10],
                    [10, 8, 9]])
test_atom_groups = collect_atom_groups(test_xy)
@test test_group == Set.(test_atom_groups)


text:
using BoundaryValueDiffEq
using DiffEqBase, OrdinaryDiffEq, DiffEqDevTools
using Test, LinearAlgebra

println("Shooting method")

function f(du,u,p,t)
  (x, v) = u
  du[1] = v
  du[2] = -x
end

function bc!(resid,sol,p,t)
  resid[1] = sol[1][1]
  resid[2] = sol[end][1] - 1
end

tspan = (0.,100.)
u0 = [0.,1.]
bvp = BVProblem(f, bc!, u0, tspan)
resid_f = Array{Float64}(undef, 2)
sol = solve(bvp, Shooting(Tsit5()))
bc!(resid_f, sol, nothing, sol.t)
@test norm(resid_f) < 1e-7

#Test for complex values
u0 = [0.,1.].+1im
bvp = BVProblem(f, bc!, u0, tspan)
resid_f = Array{ComplexF64}(undef, 2)
sol = solve(bvp, Shooting(Tsit5()))
bc!(resid_f, sol, nothing, sol.t)
@test norm(resid_f) < 1e-7


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_opt_rho2_7(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_opt_rho2_7!(A_copy)
end

@inline function exp_ps_opt_rho2_7!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=2.373110447113857
    coeff2=-11.496160069358545
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=2.373110447113857
    coeff2=-11.496160069358545
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.8558107709451741
    coeff2=6.654420725242113
    coeff3=-0.1464046050697995
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=2.6186718329845604
    coeff2=4.547282965965041
    coeff3=5.695287571427549
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.7086288697769633
    coeff2=1.5171963077831512
    coeff3=-8.805953709057599
    coeff4=25.592199199005485
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.7929434423854941
    coeff2=5.6454298259935385
    coeff3=16.186403987022313
    coeff4=-0.6370814945587178
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=4.521013785174327e-9
    coeff2=-0.016780854076753563
    coeff3=0.12296702237640218
    coeff4=-1.1554804561626417
    coeff5=0.914091240339738
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.01201228900944982
    coeff2=-1.1065314556032247e-12
    coeff3=-2.0779221004950836e-15
    coeff4=-5.249610150102727e-19
    coeff5=4.267589684021346e-27
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.0035871003621533335
    coeff2=-0.199515898070094
    coeff3=3.3717237028855647
    coeff4=-1.0770978480112308
    coeff5=2.1074592605382616
    coeff6=-0.5282957524727887
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.13628986461873388
    coeff2=-0.027077600537154935
    coeff3=0.004774676777784938
    coeff4=0.012069483753455755
    coeff5=-0.009542640674472407
    coeff6=0.8690669836434851
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.09860921274946392
    coeff2=0.2639396617162851
    coeff3=-4.745046748609644
    coeff4=3.145388920252181
    coeff5=3.6815285908979076
    coeff6=-0.00016028279238559438
    coeff7=-9.362127616405266e-5
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=5.692302453618037
    coeff2=-2.4680601957893162
    coeff3=0.7182801012293478
    coeff4=-0.5030745298888237
    coeff5=0.2802466342872092
    coeff6=-7.214590944431316e-5
    coeff7=0.8753927048889629
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=-4.433844812389262
    coeff2=0.14174985095473627
    coeff3=14.35374095270583
    coeff4=-16.423040024910904
    coeff5=1.6798059142303008
    coeff6=3.7664383775767242e-6
    coeff7=0.04873878019200554
    coeff8=1.3559720150054728
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
#!/usr/bin/env julia
# This is the code which be used as the implementation for the TextMate
# bundle command. You create a new command in the bundle and paste this code in.
# The code can be tested without using TextMate by using the mock.sh
# shell script which mock the environment variables and external executables
# used in this file.
# PList, SimpleParser julie modules must be installed however. And the
# gocode executable

using PList
import GoTMSupport

bundle_support = ENV["TM_BUNDLE_SUPPORT"]
dialog = ENV["DIALOG"]
gocode = ENV["TM_GOCODE"]
filepath = ENV["TM_FILEPATH"]
row = int(ENV["TM_LINE_NUMBER"])
col = int(ENV["TM_LINE_INDEX"]) + 1

function register_icons()
    icons = icon_plist()
    run(`$dialog images --register $icons`)
end

function icon_plist()
    writeplist_string(map({"const", "func", "package", "type", "var"}) do v
        {v => string(bundle_support, "/icons/$v.png")}
    end)
end

suggestions = GoTMSupport.completion_choices_plist(filepath, row, col, gocode)

register_icons()

# suggestions = "( { display = law; }, { display = laws; } )"
print(readall(`$dialog popup --suggestions $suggestions`))


text:
export card_db

const card_db =  eval(Meta.parse(read(joinpath(@__DIR__,"..", "data", "card_db.txt"), String)))

## Decklists often don't have the full name of flip cards
## Add some synonyms
const multicards = filter( x -> occursin(r" // ", x["name"]), card_db)

function splitter(x)
    mobj = match(r" // ", x)
    offset = mobj.offset
    String(strip(x[1:(offset-1)]))
end

function three_slashes(x)
    ans = replace(x, r" // " => " /// ")
    ans
end

function one_slash(x)
    ans = replace(x, r" // " => "/")
    ans
end

function one_slash_with_space(x)
    ans = replace(x, r" // " => " / ")
    ans
end

const multicard_synonyms  = Dict( (splitter(x["name"]), x["name"]) for x in multicards)
for x in multicards
    multicard_synonyms[three_slashes(x["name"])] = x["name"]
    multicard_synonyms[one_slash(x["name"])] = x["name"]
    multicard_synonyms[one_slash_with_space(x["name"])] = x["name"]
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
include("../src/RandomField.jl")
using .RandomField

# choose parameters:
## remember, we need that 4|n_i for i=1,...,d
N = [100,100]		# number pixels per dimension
L = [5.0,5.0]	# real space box is [-l,l]^d
f(X::Vector{T} where T<:Real) = 10*gaus_rf(X)

# derived parameters
d = dim(N) # dimension
H = (pi/2) .* N ./L

# allocate memory
X = Cloud(d,N)
G = Cscalar(X)

init!(X,L)

sample!(G,f)

tovtk(surface(G), "../data/surface.vtk")


text:
using TimerOutputs

@timeit begin
    3*3
    @timeit 2+2
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using ForwardDiff
using DiffResults
#using Sobol
using MappedArrays
using LinearAlgebra
using Roots
using BenchmarkTools
using Unitful
#using Distributions
using PositiveFactorizations

include("utils.jl")
include("core.jl")
include("solver_core.jl")

include("IAPWS95.jl")
include("gerg2008.jl")
include("solver5.jl")

minigerg = GERG2008(:H2S,:CH4) #H2S + Ch4
P0 = 4.53e6
T0 = 190
x0 = [0.05,0.95]

x01 = [1.731e-02,0.98269]
x02 = [0.06618,0.93382]
v0 = 5.690353439153631e-5
v0 = 1.00367585015v0
v01 = 2.08e-4
v02 = 6.62e-05
#m = IAPWS95()
#gerg = GERG2008()
#rho0 = 838.025
#T0 = 500
#v0 = 1.0/molar_to_weight(rho0,[m.molecularWeight],[1.0])

#Tx = 450
#rhox = 890.341250
#hx = 1000*weight_to_molar(0.749161585e3  ,[m.molecularWeight],[1.0])
#sx = 1000*weight_to_molar(0.210865845e1 ,[m.molecularWeight],[1.0])
#datx = (m,rhox*1.0u"kg/m^3",Tx*1.0u"K",[1.0])
#dat0 = (m,rho0*1.0u"kg/m^3",T0*1.0u"K",[1.0])

println("loaded")


text:
using Hiccup, Pkg
import Hiccup.div

export files

Base.joinpath() = ""

function validpath(root, path; dirs = true)
  full = normpath(root, path)
  startswith(full, root) &&
    (isfile(full) || (dirs && isdir(full)))
end

extension(f) = last(splitext(f))[2:end]

fileheaders(f) = d("Content-Type" => get(mimetypes, extension(f), "application/octet-stream"))

fileresponse(f) = d(:file => f,
                    :body => read(f),
                    :headers => fileheaders(f))

fresp(f) =
  isfile(f) ? fileresponse(f) :
  isdir(f) ?  dirresponse(f) :
  error("$f doesn't exist")

files(root, dirs = true) =
  branch(req -> validpath(root, joinpath(req[:path]...), dirs=dirs),
         req -> fresp(joinpath(root, req[:path]...)))

# Directories

files_css = """
  table { width:100%; border-radius:5px; }
  td { padding: 5px; }
  tr:nth-child(odd) { background: #f4f4ff; }
  .size { text-align: right; }
  """

function filelink(root, f)
  isdir(joinpath(root, f)) && (f = "$f/")
  a(d(:href=>f), f)
end

dirresponse(f) =
  html(head(style([mux_css, files_css])),
       body(h1("Files"),
            div(".box", table([tr(td(".file", filelink(f, x)),
                                  td(".size", string(filesize(joinpath(f, x)))))
                               for x in ["..", readdir(f)...]]))))

const ASSETS_DIR = "assets"
function packagefiles(dirs=true)
    loadpaths = LOAD_PATH
    function absdir(req)
        pkg = req[:params][:pkg]
        for p in loadpaths
            dir = joinpath(p, pkg, ASSETS_DIR)
            if isdir(dir)
                return dir
            end
        end
        Pkg.dir(String(pkg), ASSETS_DIR)  # Pkg.dir doesn't take SubString
    end

    branch(req -> validpath(absdir(req), joinpath(req[:path]...), dirs=dirs),
           req -> (Base.@warn("""
                        Relying on /pkg/ is now deprecated. Please use the package
                        `AssetRegistry.jl` instead to register assets directory
                        """, maxlog=1);
                   fresp(joinpath(absdir(req), req[:path]...))))
end

const pkgfiles = route("pkg/:pkg", packagefiles(), Mux.notfound())

using AssetRegistry

function assetserve(dirs=true)
    absdir(req) = AssetRegistry.registry["/assetserver/" * req[:params][:key]]
    branch(req -> (isfile(absdir(req)) && isempty(req[:path])) ||
           validpath(absdir(req), joinpath(req[:path]...), dirs=dirs),
           req -> fresp(joinpath(absdir(req), req[:path]...)))
end

const assetserver = route("assetserver/:key", assetserve(), Mux.notfound())


text:
module ImpulseResponse

using Statistics
using Random
using LinearAlgebra
using Distributed
using SharedArrays
using FFTW
using PyPlot

using Libaudio
using Soundcard
using DeviceUnderTest

"""
simulation
"""
function expsinesweep_simulate(fs=48000., f0=22, f1=22000, ts=10, td=3, eq=[([1.0],[1.0])])
    s = Libaudio.expsinesweep(f0, f1, ts, fs)
    m = length(s)
    n = round(Int, td * fs)
    x = zeros(m+n, 1)
    x[1:m,1] = s
    for i in eq
        x = Libaudio.filt(i[1], i[2], x)
    end
    ellip_b = [0.165069005881145, 0.064728220211450, 0.237771476924023, 0.237771476924022, 0.064728220211450, 0.165069005881146]
    ellip_a = [1.0, -1.544070855644442, 2.377072040756431, -1.638501402700271, 0.950992608325718, -0.210354984704200]
    r = Libaudio.filt(ellip_b, ellip_a, x)
    p = median(abs.(r))
    r[r.>p] .= p
    return Libaudio.impulse(s, n, f0, f1, fs, r)
end

"""
soundcard -> soundcard measurement

# Arguments
- 'ms::Matrix': mixing matrix of loudspeakers
- 'mm::Matrix': mixing matrix of microphones
- 'fs': sample rate of the soundcard
- 'f0': starting frequency
- 'f1': stop frequency
- 'ts': ess signal length in seconds
- 'td': ess decaying length in seconds
- 'eq': filter array for equalization
- 'atten': ess signal attenuation
"""
function expsinesweep_asio(ms::Matrix, mm::Matrix, fs=48000, f0=22, f1=22000, ts=10, td=3, eq=[([1.0],[1.0])], atten=-6)
    @assert f0 < f1
    f0 < 1.0 && (f0 = 1.0)
    f1 > fs/2 && (f1 = fs/2)

    s = Libaudio.expsinesweep(f0, f1, ts, fs)
    m = length(s)
    n = round(Int, td * fs)
    x = zeros(m+n, 1)
    x[1:m,1] = s
    for i in eq
        x = Libaudio.filt(i[1], i[2], x)
    end
    r = Soundcard.playrecord(x * 10^(atten/20), ms, mm, fs)
    return Libaudio.impulse(s, n, f0, f1, fs, r)
end

"""
dut -> dut measurement

# Arguments
- 'f': dut registered function dictionary
- 'ms::Matrix': mixing matrix of loudspeakers
- 'mm::Matrix': mixing matrix of microphones
- 'fs': sample rate of the dut
- 'f0': starting frequency of the ess signal
- 'f1': stop frequency of the ess signal
- 'ts': ess signal length in seconds
- 'td': ess decaying length in seconds
- 'eq': filter array for equalization
- 'atten': attenuation of the ess signal in dB
- 'syncatten': attenuation of the sync symbol in dB
- 'f2': starting frequency of the sync symbol
- 'f3': stop frequency of the sync symbol
- 'tsm': time for syncsymbol generation
- 'tsd': time for syncsymbol decaying
- 'tcs': time for context switching

# Examples
    using DeviceUnderTest
    f = DeviceUnderTest.register()
    expsinesweep_fileio(f[some device code], ...)
"""
function expsinesweep_fileio(f, ms::Matrix, mm::Matrix, fs=47999.6, f0=22, f1=22000, ts=10, td=3, eq=[([1.0],[1.0])], atten=-6, syncatten=-10, f2=800, f3=2000, tsm=0.5, tsd=2, tcs=3)
    @assert f0 < f1
    f0 < 1.0 && (f0 = 1.0)
    f1 > fs/2 && (f1 = fs/2)

    s = Libaudio.expsinesweep(f0, f1, ts, fs)
    m = length(s)
    n = round(Int, td * fs)
    x = zeros(m+n, 1)
    x[1:m,1] = s
    for i in eq
        x = Libaudio.filt(i[1], i[2], x)
    end
    sym = convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fs))
    y = Libaudio.encode_syncsymbol(tcs, sym, tsd, x * 10^(atten/20), fs, 1, syncatten)
    out = randstring() * ".wav"
    av = [8000, 16000, 44100, 48000, 96000, 192000]
    Libaudio.wavwrite(out, DeviceUnderTest.mixer(y, ms), av[findmin(abs.(av.-fs))[2]], 32)

    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    try
        f[:init]()
        f[:readyplayrecord](out, ceil(size(y,1)/fs), true)
        f[:playrecord](true)
        r, rate = Libaudio.wavread("./raw_out_mic_all_16bit_48k_8ch_mic_pcm_before_resample_8ch_48000.wav", Float64)
        nx = size(x,1)
        val, p = Libaudio.decode_syncsymbol(r, sym, tsd, nx/fs, fs)
        val || Libaudio.printl(root, :light_red, Libaudio.nows() * " | ImpulseResponse.expsinesweep_fileio: decode symbol failure")
        c = size(r,2)
        rd = zeros(nx, c)
        tol = 2048
        for i = 1:c
            l = p[i] - tol
            rd[:,i] = r[l:l+nx-1,i]
        end
        return Libaudio.impulse(s, n, f0, f1, fs, rd)
    finally
        rm(out, force=true)
        rm("./raw_out_mic_all_16bit_48k_8ch_mic_pcm_before_resample_8ch_48000.wav", force=true)
    end
end

"""
soundcard -> dut measurement
"""
function expsinesweep_asio_fileio(f, ms::Matrix, mm::Matrix, fs=48000, fm=47999.6, f0=22, f1=22000, ts=10, td=3, eq=[([1.0],[1.0])], atten=-6, syncatten=-10, f2=800, f3=2000, tsm=0.5, tsd=2, tcs=3)
    @assert nprocs() > 1
    @assert f0 < f1
    f0 < 1.0 && (f0 = 1.0)
    f1 > fs/2 && (f1 = fs/2)
    wid = workers()

    s = Libaudio.expsinesweep(f0, f1, ts, fs)
    m = length(s)
    n = round(Int, td * fs)
    x = zeros(m+n, 1)
    x[1:m,1] = s
    for i in eq
        x = Libaudio.filt(i[1], i[2], x)
    end
    sym = convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fs))
    y = Libaudio.encode_syncsymbol(tcs, sym, tsd, x * 10^(atten/20), fs, 1, syncatten)

    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    try
        f[:init]()
        f[:readyrecord](ceil(size(y,1)/fs), true)
        phy = Soundcard.mixer(y, ms)
        pcm = SharedArray{Float32,1}(Soundcard.interleave(phy))
        done = remotecall(Soundcard.play, wid[1], size(phy), pcm, fs)  # latency is low
        f[:record](true)
        fetch(done)
        r,rate = Libaudio.wavread("./raw_out_mic_all_16bit_48k_8ch_mic_pcm_before_resample_8ch_48000.wav", Float64)
        nx = size(x,1)
        tx = nx / fs
        syma = convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fm))
        val, p = Libaudio.decode_syncsymbol(r, syma, tsd, tx, fm)
        val || Libaudio.printl(root, :light_red, Libaudio.nows() * " | ImpulseResponse.expsinesweep_asio_fileio: decode symbol failure")
        sa = Libaudio.expsinesweep(f0, f1, ts, fm)
        nxa = round(Int, tx * fm)
        na = nxa-length(sa)

        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.expsinesweep_asio_fileio: decay samples [td x fm]/[tx x fm] $(round(Int,td * fm))/$(na)")
        c = size(r,2)
        tol = 2048
        rd = zeros(nxa, c)
        for i = 1:c
            l = p[i] - tol
            rd[:,i] = r[l:l+nxa-1,i]
        end
        return Libaudio.impulse(sa, na, f0, f1, fm, rd)
    finally
        rm("./raw_out_mic_all_16bit_48k_8ch_mic_pcm_before_resample_8ch_48000.wav", force=true)
    end
end

"""
dut -> soundcard measurement
"""
function expsinesweep_fileio_asio(f, ms::Matrix, mm::Matrix, fs=48000, fm=47999.6, f0=22, f1=22000, ts=10, td=3, eq=[([1.0],[1.0])], atten=-6, syncatten=-10, f2=800, f3=2000, tsm=0.5, tsd=2, tcs=3)
    @assert nprocs() > 1
    @assert f0 < f1
    f0 < 1.0 && (f0 = 1.0)
    f1 > fs/2 && (f1 = fs/2)
    wid = workers()

    s = Libaudio.expsinesweep(f0, f1, ts, fm)
    m = length(s)
    n = round(Int, td * fm)
    x = zeros(m+n, 1)
    x[1:m,1] = s
    for i in eq
        x = Libaudio.filt(i[1], i[2], x)
    end
    sym = convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fm))
    y = Libaudio.encode_syncsymbol(tcs, sym, tsd, x * 10^(atten/20), fm, 1, syncatten)
    out = randstring() * ".wav"
    Libaudio.wavwrite(out, DeviceUnderTest.mixer(y, ms), fs, 32)

    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    try
        f[:init]()
        f[:readyplay](out)
        done = remotecall(f[:play], wid[1])
        r = convert(Matrix{Float64}, Soundcard.record(round(Int,fs*size(y,1)/fm), mm, fs))
        fetch(done)
        nx = size(x,1)
        tx = nx / fm
        syma = convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fs))
        val, p = Libaudio.decode_syncsymbol(r, syma, tsd, tx, fs)
        val || Libaudio.printl(root, :light_red, Libaudio.nows() * " | ImpulseResponse.expsinesweep_fileio_asio: decode symbol failure")
        sa = Libaudio.expsinesweep(f0, f1, ts, fs)
        nxa = round(Int, tx * fs)
        na = nxa-length(sa)

        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.expsinesweep_fileio_asio: decay samples [td x fs]/[tx x fs] $(round(Int,td * fs))/$(na)")
        c = size(r,2)
        tol = 2048
        rd = zeros(nxa, c)
        for i = 1:c
            l = p[i] - tol
            rd[:,i] = r[l:l+nxa-1,i]
        end
        return Libaudio.impulse(sa, na, f0, f1, fs, rd)
    finally
        rm(out, force=true)
    end
end

# const sndcard_n_out = 8
# const sndcard_n_in = 8
# const mic_n = 1
# mixplay = zeros(Float32, size(essd,2), sndcard_n_out)
# mixplay[1,2] = 1.0f0
# mixrec = zeros(Float32, sndcard_n_in, mic_n)
# mixrec[2,1] = 1.0f0
#
# note: adjust atten for the level of the ess signal
# note: adjust syncatten for the level of the sync symbol, if too high mixer would prompt with mic clipping error!
# note: use long t_ess when possible, for the build-up of the stimulus energy
# note: use long t_decay if room or system dunamics are reverberant
# note: eq = [(b,a),(b,a)...] is prepending transfer function for filter verification, must be designed according to fs or fsd based on mode
# note: mode[1] is the physical device for playback, mode[2] is the physical device for recording

# mixspk = zeros(1,2)
# mixspk[1,1] = 1.0
# mixmic = zeros(9,1)
# mixmic[9,1] = 1.0
"""
dut@fm -> soundcard@fs
"""
function measureclockdrift(f, ms::Matrix{Float64}, mm::Matrix{Float64}, rep=3, fs=48000, f2=800, f3=2000, tsm=0.5, atten=-6)
    #
    #   +--------+--------+--------+--------+--------+ => 5 samples in digital domain played via dut's speaker, whose sample interval is Td.
    #   +-----+-----+-----+-----+-----+-----+-----+--- => 7 samples captured by the standard sampler of the soundcard,
    #
    #   5 x Td ≈ 7 x Tr
    #   or formly, N/Fd ≈ Nm / Fr
    #   Fd ≈ N / Nm x Fr
    #   Fd/Fr ≈ N/Nm = 5/7
    #
    @assert nprocs() > 1
    wpid = workers()

    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: start measuring device clock drift")

    # fileio -> asio
    sync = 10^(atten/20) * convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fs))
    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: measureclockdrift: sync samples $(length(sync))")

    period = [zeros(round(Int,100fs),1); sync]
    signal = [zeros(round(Int,3fs),1); sync; repeat(period,rep,1); zeros(round(Int,3fs),1)]
    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: stimulus formed")

    out = randstring() * ".wav"
    Libaudio.wavwrite(out, DeviceUnderTest.mixer(signal, ms), fs, 32)
    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: filesize $(filesize(out)/1024/1024) MiB")

    measure = Array{Tuple{Float64, Float64, Float64},1}()
    try
        f[:init]()
        f[:readyplay](out)
        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: stimulus pushed to device")

        done = remotecall(f[:play], wpid[1])
        r = convert(Matrix{Float64}, Soundcard.record(size(signal,1), mm, fs))
        fetch(done)
        # Libaudio.wavwrite("clockdrift.wav", r, fs, 16)
        val = true
        for k = 1:size(r,2)
            flag,lbs,pk,pkf,y = Libaudio.extractsymbol(r[:,k], sync, rep+1)
            val = val && flag
            pkfd = diff(pkf)
            chrodrift_100sec = ((pkfd[end] - pkfd[1]) / (rep-1))/fs
            freqdrift_100sec = (size(period,1) - median(pkfd))/fs
            push!(measure, (fs * size(period,1)/median(pkfd), freqdrift_100sec, chrodrift_100sec))
        end
        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: status $val")
    finally
        rm(out)
    end

    for k in measure
        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: time drift every 100 seconds $(k[2]/100)")
        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: temperature drift every 100 seconds $(k[3]/100)")
        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: estimated sample rate $(k[1])")
    end

    measure
end

"""
soundcard@fs -> dut@fm
"""
function measureclockdrift2(f, ms::Matrix{Float64}, mm::Matrix{Float64}, rep=3, fs=48000, f2=800, f3=2000, tsm=0.5, atten=-6)
    #
    #   +--------+--------+--------+--------+--------+ => 5 samples in digital domain played via standard sampler of the soundcard, whose sample interval is Tr.
    #   +-----+-----+-----+-----+-----+-----+-----+--- => 7 samples captured by the imperfect sampler of the device under test,
    #
    #   5 x Tr ≈ 7 x Td
    #   or formly, N/Fr ≈ Nm / Fd
    #   Fd ≈ Nm / N x Fr
    #   Fd/Fr ≈ Nm/N = 5/7
    #
    @assert nprocs() > 1
    wpid = workers()
    root = joinpath(Libaudio.folder(), Libaudio.logfile())

    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift2: start measuring device clock drift")

    sync = 10^(atten/20) * convert(Vector{Float64}, Libaudio.symbol_expsinesweep(f2, f3, tsm, fs))
    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift2: sync samples $(length(sync))")

    period = [zeros(round(Int,100fs),1); sync]
    signal = [zeros(round(Int,3fs),1); sync; repeat(period,rep,1); zeros(round(Int,3fs),1)]
    Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift2: stimulus formed")

    measure = Array{Tuple{Float64, Float64, Float64},1}()
    try
        f[:init]()
        f[:readyrecord](ceil(size(signal,1)/fs), true)
        phy = Soundcard.mixer(signal, ms)
        pcm = SharedArray{Float32,1}(Soundcard.interleave(phy))
        done = remotecall(Soundcard.play, wpid[1], size(phy), pcm, fs)  # latency is low
        f[:record](true)
        fetch(done)
        r, rate = Libaudio.wavread("./raw_out_mic_all_16bit_48k_8ch_mic_pcm_before_resample_8ch_48000.wav", Float64)

        val = true
        for k = 1:size(r,2)
            flag,lbs,pk,pkf,y = Libaudio.extractsymbol(convert(Vector{Float64},r[:,k]), sync, rep+1)
            val = val && flag
            pkfd = diff(pkf)
            chrodrift_100sec = ((pkfd[end] - pkfd[1]) / (rep-1))/fs
            freqdrift_100sec = (size(period,1) - median(pkfd))/fs
            push!(measure, (fs * median(pkfd)/size(period,1), freqdrift_100sec, chrodrift_100sec))
        end
        Libaudio.printl(root, :light_cyan, Libaudio.nows() * " | ImpulseResponse.measureclockdrift: status $val")
    finally
        rm("./raw_out_mic_all_16bit_48k_8ch_mic_pcm_before_resample_8ch_48000.wav")
    end
    measure
end

"""
function hFIR = eq_calibration(mix_spk, mix_mic, f_anchor, f_start, f_stop, atten)

    fs = 48000;
    ess_f0 = 10;
    ess_f1 = fs/2;
    ess_time = 30;
    ess_decay = 5;
    %atten = -20;
    nfft = 65536;
    nnyq = nfft/2+1;

    %@ measure the un-eq'd impulse response in the anechoic
    [fundamental, harmonics, response_t] = impulse_response_exponential_sine_sweep(mix_spk, mix_mic, ess_f0, ess_f1, ess_time, ess_decay, 'asio', [1], [1], atten);

    %@ cut the window with the impulse response out for analysis
    p = 4.3e4;
    x = fundamental(p+1:p+nfft);
    x_spec = fft(x)./nfft;
    x_phase = angle(x_spec);
    f = ((0:nnyq-1)'./nfft).*fs;
    x_spec_db = 20*log10(abs(x_spec(1:nnyq)));

    % smooth and plots
    Noct = 3;
    x_spec_db_sm = smoothSpectrum(x_spec_db,f,Noct);
    figure(1);
    subplot(2,1,1); plot(20*log10(abs(x))); title('impulse response in time domain'); xlabel('samples'); hold on;
    %subplot(2,1,2); semilogx(f,x_spec_db,f,x_spec_db_sm); title('impulse response in frequency domain, with 1/3 octave smoothing'); xlabel('Hz'); hold on;
    subplot(2,1,2); semilogx(f,x_spec_db_sm); title('impulse response in frequency domain, with 1/3 octave smoothing'); xlabel('Hz'); hold on;

    % construct a minimal phase filter of the smoothed impulse response
    x_spec_sm = 10.^(x_spec_db_sm/20);
    x_spec_sm = [x_spec_sm; flipud(x_spec_sm(2:end-1))];
    x_spec_sm = x_spec_sm .* exp(x_phase*1i);
    x_sm = real(ifft(x_spec_sm));

    zs = fft(x_sm);
    figure(2); semilogx(f, 20*log10(abs(zs(1:nnyq))), 'k'); hold on;

    zms = mps(fft(x_sm));
    zm = real(ifft(zms)); % it is not symmetrical
    zms = fft(zm);
    semilogx(f, 20*log10(abs(zms(1:nnyq))), 'c--');
    xlabel('Hz'); title('minimal phase filter constructed based on smoothed frequency response'); grid on;

    % do the actual flatten work
    %f1 = 22;
    f1 = f_start;
    f2 = f_stop;
    f1 = ceil(f1 * nfft / fs);
    f2 = floor(f2 * nfft / fs);
    target = x_spec_db_sm(round(f_anchor * nfft / fs));
    %target = x_spec_db(round(f_anchor * nfft / fs));
    H = zeros(nnyq,1);
    H(f1:f2) = target - x_spec_db_sm(f1:f2);
    %H(f1:f2) = target - x_spec_db(f1:f2);
    H(1:f1-1) = H(f1);
    figure(3); semilogx(f, H, 'r'); hold on;

    % compensation filter in time and frequency domain
    H = 10.^(H/20) * nfft;
    H = [H; flipud(H(2:end-1))];
    H = H .* exp(x_phase*1i);
    h = real(ifft(H));
    hs = fft(h)./nfft;
    semilogx(f, 20*log10(abs(hs(1:nnyq))), 'b--');

    hms = mps(fft(h));
    hm = real(ifft(hms)); % it is not symmetrical
    hms = fft(hm)./nfft;
    semilogx(f, 20*log10(abs(hms(1:nnyq))), 'c--');
    xlabel('Hz'); title('Compensation filter and its minimal-phase realization'); grid on;
    figure(2); semilogx(f, 20*log10(abs(hms(1:nnyq)))+x_spec_db_sm, 'b'); grid on;

    [fundamental, harmonics, response_t] = impulse_response_exponential_sine_sweep(mix_spk, mix_mic, ess_f0, ess_f1, ess_time, ess_decay, 'asio', hm/nfft, [1], atten);

    p = 4.3e4;
    x = fundamental(p+1:p+nfft);
    x_spec = fft(x)./nfft;
    x_phase = angle(x_spec);
    f = ((0:nnyq-1)'./nfft).*fs;
    x_spec_db = 20*log10(abs(x_spec(1:nnyq)));

    % smooth and plots
    Noct = 3;
    x_spec_db_sm = smoothSpectrum(x_spec_db,f,Noct);
    figure(1);
    subplot(2,1,1); plot(20*log10(abs(x)), 'r'); grid on;
    %subplot(2,1,2); semilogx(f,x_spec_db,f,x_spec_db_sm); grid on;
    subplot(2,1,2); semilogx(f,x_spec_db_sm); grid on;

    hFIR = hm / nfft;
end
"""
function minimalphase_fir_original(ms::Matrix, mm::Matrix, fs, f0, f1=fs/2, attenuate=-50, unitcircle=1024, tess=30, tdecay=3, noct=3)
    PyPlot.ion()

    fundamental, harmonic, dirac, measure = expsinesweep_asio(ms, mm, fs, f0, f1, tess, tdecay, [([1.0],[1.0])], attenuate)

    f, x = Libaudio.freqz(fundamental[:,1], fs, unitcircle)
    m = length(f)
    ϕ = [angle.(x); angle.(conj(reverse(x[2:end-1])))]
    xd = 20log10.(abs.(x))
    xds = Libaudio.smoothspectrum(xd, f, noct)

    PyPlot.figure("impulse response and 1/3 oct smoothing")
    PyPlot.semilogx(f, xd)
    PyPlot.semilogx(f, xds)
    PyPlot.xlabel("Hz")
    PyPlot.grid(true)

    ω, y = Libaudio.freqz(harmonic[:,1], fs, unitcircle)
    yd = 20log10.(abs.(y))

    PyPlot.figure("harmonic response")
    PyPlot.semilogx(ω, yd)
    PyPlot.xlabel("Hz")
    PyPlot.grid(true)

    return (f, m, ϕ, xds)
end

function minimalphase_fir_design(f, m, ϕ, xds, fs, f0, f1, fx, unitcircle=1024)
        # construct a minimal phase filter of the smoothed impulse response
        xs = 10 .^ (xds/20)
        xs = [xs; reverse(xs[2:end-1])]
        fundamental_sm = real(ifft(xs .* exp.(ϕ*im)))

        y = fft(fundamental_sm)
        z = Libaudio.mps(y)
        fundamental_mp = real(ifft(z))    # not symmetrical

        PyPlot.figure("minimal-phase design")
        PyPlot.semilogx(f, 20log10.(abs.(y[1:m])))
        PyPlot.semilogx(f, 20log10.(abs.(fft(fundamental_mp)[1:m])))
        PyPlot.xlabel("Hz")
        PyPlot.grid(true)

        # do the actual flatten work
        nfft = 2unitcircle
        ω0 = max(round(Int, f0*nfft/fs), 1)
        ω1 = min(round(Int, f1*nfft/fs), m)
        ωx = round(Int, fx*nfft/fs)
        anchor = xds[ωx]

        H = zeros(m)
        H[ω0:ω1] .= anchor .- xds[ω0:ω1]
        H[1:ω0-1] .= H[ω0]
        H[ω1+1:end] .= H[ω1]
        PyPlot.figure("minimal-phase compensation filter")
        PyPlot.semilogx(f, H)

        H = 10 .^ (H/20) * nfft
        H = [H; reverse(H[2:end-1])]
        H = H .* exp.(ϕ*im)
        h = real(ifft(H))
        hs = fft(h)/nfft
        PyPlot.semilogx(f, 20log10.(abs.(hs[1:m])))

        hms = Libaudio.mps(fft(h))
        hm = real(ifft(hms))    # not symmetrical
        hms = fft(hm)/nfft
        PyPlot.semilogx(f, 20log10.(abs.(hms[1:m])))
        PyPlot.xlabel("Hz")
        PyPlot.grid(true)

        PyPlot.figure("minimal-phase design")
        PyPlot.semilogx(f, 20log10.(abs.(hms[1:m]))+xds)
        hfir = hm / nfft
end

function minimalphase_fir_verification(h, ms::Matrix, mm::Matrix, fs, f0, f1=fs/2, attenuate=-50, unitcircle=1024, tess=30, tdecay=3, noct=3)

    fundamental, harmonic, dirac, measure = expsinesweep_asio(ms, mm, fs, f0, f1, tess, tdecay, [(h,[1.0])], attenuate)

    f, x = Libaudio.freqz(fundamental[:,1], fs, unitcircle)
    xd = 20log10.(abs.(x))
    xds = Libaudio.smoothspectrum(xd, f, noct)

    PyPlot.figure("impulse response and 1/3 oct smoothing")
    PyPlot.semilogx(f, xd)
    PyPlot.semilogx(f, xds)

    ω, y = Libaudio.freqz(harmonic[:,1], fs, unitcircle)
    yd = 20log10.(abs.(y))

    PyPlot.figure("harmonic response")
    PyPlot.semilogx(ω, yd)
    PyPlot.xlabel("Hz")
    PyPlot.grid(true)

    return nothing
end

end # module


text:
@testset "discretediag.jl" begin
    nparams = 4
    nchains = 2
    samples = rand(-100:100, 100, nparams, nchains)

    @testset "results" begin
        for method in
            (:weiss, :hangartner, :DARBOOT, :MCBOOT, :billingsley, :billingsleyBOOT)
            between_chain, within_chain = @inferred(discretediag(samples; method=method))

            @test between_chain isa NamedTuple{(:stat, :df, :pvalue)}
            for name in (:stat, :df, :pvalue)
                x = getfield(between_chain, name)
                @test x isa Vector{Float64}
                @test length(x) == nparams
            end

            @test within_chain isa NamedTuple{(:stat, :df, :pvalue)}
            for name in (:stat, :df, :pvalue)
                x = getfield(within_chain, name)
                @test x isa Matrix{Float64}
                @test size(x) == (nparams, nchains)
            end
        end
    end

    @testset "exceptions" begin
        @test_throws ArgumentError discretediag(samples; method=:somemethod)
        for x in (-0.3, 0, 1, 1.2)
            @test_throws ArgumentError discretediag(samples; frac=x)
        end
    end
end


text:
# BigWig Statistics
# =================

"""
    coverage(reader, chrom, chromstart, chromend; usezoom=false)::Int

Compute the coverage of values in `[chromstart, chromend]` of `chrom`.

If `usezoom` is `true`, this function tries to use precomputed statistics (zoom) in the file.
This is often faster but not exact in most cases.
"""
function coverage(reader::Reader, chrom::AbstractString, chromstart::Integer, chromend::Integer; usezoom=false)::Int
    chromid = reader.chroms[chrom][1]
    chromstart -= 1
    if usezoom
        zoom = BBI.find_best_zoom(reader.zooms, UInt32(chromend - chromstart))
        if zoom !== nothing
            return BBI.coverage(zoom, chromid, UInt32(chromstart), UInt32(chromend))
        end
    end
    return exact_coverage(reader, chromid, UInt32(chromstart), UInt32(chromend))
end

function exact_coverage(reader::Reader, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    cov::Int = 0
    for record in OverlapIterator(reader, chromid, chromstart, chromend)
        cov += BBI.coverage2((record.chromstart, record.chromend), (chromstart, chromend))
    end
    return cov
end

"""
    mean(reader, chrom, chromstart, chromend; usezoom=false)::Float32

Compute the mean of values in `[chromstart, chromend]` of `chrom`.

This function returns `NaN32` if there are no data in that range.

See `coverage` for the `usezoom` keyword argument.
"""
function mean(reader::Reader, chrom::AbstractString, chromstart::Integer, chromend::Integer; usezoom=false)::Float32
    chromid = reader.chroms[chrom][1]
    chromstart -= 1
    if usezoom
        zoom = BBI.find_best_zoom(reader.zooms, UInt32(chromend - chromstart))
        if zoom !== nothing
            return BBI.mean(zoom, chromid, UInt32(chromstart), UInt32(chromend))
        end
    end
    return exact_mean(reader, chromid, UInt32(chromstart), UInt32(chromend))
end

function exact_mean(reader::Reader, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    # compute size-weighted mean
    sum = 0.0f0
    size = 0
    for record in OverlapIterator(reader, chromid, chromstart, chromend)
        cov = BBI.coverage2((record.chromstart, record.chromend), (chromstart, chromend))
        sum += record.value * cov
        size += cov
    end
    return sum / size
end

"""
    minimum(reader, chrom, chromstart, chromend; usezoom=false)::Float32

Compute the minimum of values in `[chromstart, chromend]` of `chrom`.

This function returns `NaN32` if there are no data in that range.

See `coverage` for the `usezoom` keyword argument.
"""
function minimum(reader::Reader, chrom::AbstractString, chromstart::Integer, chromend::Integer; usezoom=false)::Float32
    chromid = reader.chroms[chrom][1]
    chromstart -= 1
    if usezoom
        zoom = BBI.find_best_zoom(reader.zooms, UInt32(chromend - chromstart))
        if zoom !== nothing
            return BBI.minimum(zoom, chromid, UInt32(chromstart), UInt32(chromend))
        end
    end
    return exact_extrema(reader, chromid, UInt32(chromstart), UInt32(chromend))[1]
end

"""
    maximum(reader, chrom, chromstart, chromend; usezoom=false)::Float32

Compute the maximum of values in `[chromstart, chromend]` of `chrom`.

This function returns `NaN32` if there are no data in that range.

See `coverage` for the `usezoom` keyword argument.
"""
function maximum(reader::Reader, chrom::AbstractString, chromstart::Integer, chromend::Integer; usezoom=false)::Float32
    chromid = reader.chroms[chrom][1]
    chromstart -= 1
    if usezoom
        zoom = BBI.find_best_zoom(reader.zooms, UInt32(chromend - chromstart))
        if zoom !== nothing
            return BBI.maximum(zoom, chromid, UInt32(chromstart), UInt32(chromend))
        end
    end
    return exact_extrema(reader, chromid, UInt32(chromstart), UInt32(chromend))[2]
end

function exact_extrema(reader::Reader, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    min = +Inf32
    max = -Inf32
    defined = false
    for record in OverlapIterator(reader, chromid, chromstart, chromend)
        min = Base.min(min, record.value)
        max = Base.max(max, record.value)
        defined = true
    end
    return defined ? (min, max) : (NaN32, NaN32)
end

"""
    std(reader, chrom, chromstart, chromend; usezoom=false)::Float32

Compute the standard deviation of values in `[chromstart, chromend]` of `chrom`.

See `coverage` for the `usezoom` keyword argument.
"""
function std(reader::Reader, chrom::AbstractString, chromstart::Integer, chromend::Integer; usezoom=false)::Float32
    chromid = reader.chroms[chrom][1]
    chromstart -= 1
    if usezoom
        zoom = BBI.find_best_zoom(reader.zooms, UInt32(chromend - chromstart))
        if zoom !== nothing
            return BBI.std(zoom, chromid, UInt32(chromstart), UInt32(chromend))
        end
    end
    return exact_std(reader, chromid, UInt32(chromstart), UInt32(chromend))
end

function exact_std(reader::Reader, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    sum = 0.0f0
    ssq = 0.0f0
    size = 0
    for record in OverlapIterator(reader, chromid, chromstart, chromend)
        cov = BBI.coverage2((record.chromstart, record.chromend), (chromstart, chromend))
        sum += record.value * cov
        ssq += record.value^2 * cov
        size += cov
    end
    return sqrt((ssq - sum^2 / size) / (size - 1))
end


text:
"""
    Optimization solver written by pure Julia
"""
module ArgMin

using LinearAlgebra

include("least_square_solver.jl")
include("convex_optimization_solver.jl")
include("projected_newton_qp_solver.jl")

end


text:
@static if VERSION >= v"1.6"
    const NANSAFE_MODE_ENABLED = @load_preference("nansafe_mode", false)
else
    const NANSAFE_MODE_ENABLED = false
end

const AMBIGUOUS_TYPES = (AbstractFloat, Irrational, Integer, Rational, Real, RoundingMode)

const UNARY_PREDICATES = Symbol[:isinf, :isnan, :isfinite, :iseven, :isodd, :isreal, :isinteger]

const BINARY_PREDICATES = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), :(>=)]

const DEFAULT_CHUNK_THRESHOLD = 12

struct Chunk{N} end

const CHUNKS = [Chunk{i}() for i in 1:DEFAULT_CHUNK_THRESHOLD]

function Chunk(input_length::Integer, threshold::Integer = DEFAULT_CHUNK_THRESHOLD)
    N = pickchunksize(input_length, threshold)
    0 < N <= DEFAULT_CHUNK_THRESHOLD && return CHUNKS[N]
    return Chunk{N}()
end

function Chunk(x::AbstractArray, threshold::Integer = DEFAULT_CHUNK_THRESHOLD)
    return Chunk(length(x), threshold)
end

# Constrained to `N <= threshold`, minimize (in order of priority):
#   1. the number of chunks that need to be computed
#   2. the number of "left over" perturbations in the final chunk
function pickchunksize(input_length, threshold = DEFAULT_CHUNK_THRESHOLD)
    if input_length <= threshold
        return input_length
    else
        nchunks = round(Int, input_length / DEFAULT_CHUNK_THRESHOLD, RoundUp)
        return round(Int, input_length / nchunks, RoundUp)
    end
end

chunksize(::Chunk{N}) where {N} = N

replace_match!(f, ismatch, x) = x

function replace_match!(f, ismatch, lines::AbstractArray)
    for i in eachindex(lines)
        line = lines[i]
        if ismatch(line)
            lines[i] = f(line)
        elseif isa(line, Expr)
            replace_match!(f, ismatch, line.args)
        end
    end
    return lines
end

# This is basically a workaround that allows one to use CommonSubexpressions.cse, but with
# qualified bindings. This is not guaranteed to be correct if the input expression contains
# field accesses.
function qualified_cse!(expr)
    placeholders = Dict{Symbol,Expr}()
    replace_match!(x -> isa(x, Expr) && x.head == :(.), expr.args) do x
        placeholder = Symbol("#$(hash(x))")
        placeholders[placeholder] = x
        placeholder
    end
    cse_expr = CommonSubexpressions.cse(expr, warn=false)
    replace_match!(x -> haskey(placeholders, x), cse_expr.args) do x
        placeholders[x]
    end
    return cse_expr
end


text:
import Singular
import AbstractAlgebra

using Base.Threads

using BenchmarkTools
include("../src/Groebner.jl")

using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

# BenchmarkTools.DEFAULT_PARAMETERS.seconds = 100000
# BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    gb = Groebner.groebner(system)

    # @btime gb = Groebner.groebner($system, reduced=false)

    bench = @benchmarkable Groebner.groebner($system, linalg=:prob) samples=5
    println(median(run(bench)))
end

function benchmark_system_singular(system)
    R = AbstractAlgebra.parent(system[1])
    modulo = AbstractAlgebra.characteristic(R)
    n = AbstractAlgebra.nvars(R)
    ground_s = Singular.N_ZpField(modulo)
    R_s, _ = Singular.PolynomialRing(ground_s, ["x$i" for i in 1:n], ordering=:degrevlex)

    system_s = map(
        f -> AbstractAlgebra.change_base_ring(
                    ground_s,
                    AbstractAlgebra.map_coefficients(c -> ground_s(c.d), f),
                    parent=R_s),
        system)

    ideal_s = Singular.Ideal(R_s, system_s)

    Singular.std(Singular.Ideal(R_s, [system_s[1]]))

    # @btime gb = Singular.std($ideal_s, complete_reduction=false)
    bench = @benchmarkable  Singular.std($ideal_s) samples=5
    println(median(run(bench)))
end

function handle_system_my(name, system)
    println("---------------\n$name")
    println("my")
    benchmark_system_my(system)
    println("---------------")
end

function handle_system_singular(name, system)
    println("---------------\n$name")
    println("singular")
    benchmark_system_singular(system)
    println("---------------")
end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
    ("cyclic 7", Groebner.cyclicn(7, ground=ground)), # 0.06 vs 0.01
    ("cyclic 8", Groebner.cyclicn(8, ground=ground)), # 0.19 vs 0.03
    ("cyclic 9", Groebner.cyclicn(9, ground=ground)), # 0.60 vs 0.9
        ("cyclic 12", Groebner.rootn(12, ground=ground)), # 0.06 vs 0.01
        ("cyclic 13", Groebner.rootn(13, ground=ground)), # 0.19 vs 0.03
        ("cyclic 14", Groebner.rootn(14, ground=ground)), # 0.60 vs 0.9
        ("katsura 11",Groebner.katsura11(ground=ground)), # 1.35 vs 71.121
        ("katsura 12",Groebner.katsura12(ground=ground)), # 9.68 vs 775
        ("katsura 13",Groebner.katsura13(ground=ground)), # 65.22 vs 6802
        ("eco 11",Groebner.eco11(ground=ground)),         # 0.41  vs 31.7
        ("eco 12",Groebner.eco12(ground=ground)),         # 2.5   vs 341
        ("eco 13",Groebner.eco13(ground=ground)),         # 17    vs 4102
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),  # 0.19  vs 0.36
        ("noon 8"    ,Groebner.noonn(8, ground=ground)),  # 1.8   vs 3.2
        ("noon 9"    ,Groebner.noonn(9, ground=ground)),   # 18.1  vs 33.3
        ("henrion 5"    ,Groebner.henrion5(ground=ground)),  # 0.19  vs 0.36
        ("henrion 6"    ,Groebner.henrion6(ground=ground)),  # 1.8   vs 3.2
        ("henrion 7"    ,Groebner.henrion7(ground=ground))   # 18.1  vs 33.3
    ]

    for (name, system) in systems
        task1 = @spawn handle_system_my(name, system)
        task2 = @spawn handle_system_singular(name, system)
        wait(task1)
        wait(task2)
    end
end

ground = AbstractAlgebra.GF(2^31 - 1)
run_f4_ff_degrevlex_benchmarks(ground)

#=

=#


text:
"""
    servedocs_callback!(args...)

Custom callback used in [`servedocs`](@ref) triggered when a file is modified.

If the file is `docs/make.jl`, the callback will check whether any new files
have subsequently been generated in the `docs/src` folder and add them to the
watched files, it will also remove any file that may have been deleted or
renamed.

If the file is either in `docs/src`, a pass of Documenter is triggered to
regenerate the documentation, subsequently the LiveServer will render the
pages produced in `docs/build`.

## Arguments

See the docs of the parent function `servedocs`.
"""
function servedocs_callback!(
            dw::SimpleWatcher,
            fp::AbstractString,
            path2makejl::AbstractString,
            literate::Union{Nothing,String},
            skip_dirs::Vector{String},
            foldername::String,
            buildfoldername::String)
    # ignore things happening in the build folder (generated files etc)
    startswith(fp, joinpath(foldername, buildfoldername)) && return nothing
    # ignore things happening in any skip_dirs
    for dir in skip_dirs
        startswith(fp, dir) && return nothing
    end

    # if the file that was changed is the `make.jl` file, assume that maybe
    # new files have been generated and so refresh the vector of watched files
    if fp == path2makejl
        # it's easier to start from scratch (takes negligible time)
        empty!(dw.watchedfiles)
        scan_docs!(dw, foldername, path2makejl, literate)
    end

    # Run a Documenter pass
    Main.include(abspath(path2makejl))
    file_changed_callback(fp)
    return
end

"""
    scan_docs!(dw::SimpleWatcher, args...)

Scans the `docs/` folder and add all relevant files to the watched files.

## Arguments

See the docs of the parent function `servedocs`.
"""
function scan_docs!(dw::SimpleWatcher,
                    foldername::String,
                    path2makejl::String,
                    literate::Union{Nothing,String}
                    )::Nothing
    # Typical expected structure:
    #   docs
    #   ├── make.jl
    #   └── src
    #       ├── *
    #       └── *
    src = joinpath(foldername, "src")

    if !isdir(foldername) || !isdir(src)
        error("I didn't find a $foldername/ or $foldername/src/ folder.")
    end
    # watch the make.jl file as well as

    push!(dw.watchedfiles, WatchedFile(path2makejl))
    if isdir(foldername)
        # add all files in `docs/src` to watched files
        for (root, _, files) ∈ walkdir(joinpath(foldername, "src")), file ∈ files
            push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
        end
    end

    # if the user is not using Literate, return early
    literate === nothing && return

    # If the user gave a specific directory for literate files, add all the
    # files in that directory to the watched files.
    # If the user did not (literate === "") then these files are already watched.
    if !isempty(literate)
        isdir(literate) || error("I didn't find the provided literate folder $literate.")
        for (root, _, files) ∈ walkdir(literate), file ∈ files
            push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
        end
    end

    # When using Literate, each script *.jl is assumed to generate a corresponding *.md
    # file. Only the source (*.jl) file should be watched to avoid an infinite loop where
    # 1. the script is executed, generates the .md file
    # 2. the .md file is seen as modified, triggers the make.jl file
    # 3. the make.jl file executes the script (↩1)
    # So here we remove from the watchlist all files.md that have a files.jl with the same path.
    remove = Int[]
    if isempty(literate)
        # assumption is that the scripts are in `docs/src/...` and that the
        # generated markdown goes in exactly the same spot so for instance:
        #   docs
        #   └── src
        #       ├── index.jl
        #       └── index.md
        for wf ∈ dw.watchedfiles
            spath = splitext(wf.path)
            # ignore non `*.jl` files
            spath[2] == ".jl" || continue
            # if a `.jl` file is found, check if there's any corresponding `.md`
            # file and, if so, remove that file from the watched files to avoid ∞
            k = findfirst(e -> splitext(e.path) == (spath[1], ".md"), dw.watchedfiles)
            k === nothing || push!(remove, k)
        end
    else
        # assumption is that the scripts are in a specific folder and that the
        # generated markdown goes in `docs/src` with the same relative path
        # so for instance:
        #   docs
        #   ├── literate
        #   │   └── index.jl
        #   └── src
        #       └── index.md
        # the logic is otherwise the same as above.
        for (root, _, files) ∈ walkdir(literate), file ∈ files
            spath = splitext(joinpath(root, file))
            spath[2] == ".jl" || continue
            path = replace(spath[1], Regex("^$literate") => joinpath(foldername, "src"))
            k = findfirst(e -> splitext(e.path) == (path, ".md"), dw.watchedfiles)
            k === nothing || push!(remove, k)
        end
    end
    deleteat!(dw.watchedfiles, remove)
    return
end

"""
    servedocs(; kwargs...)

Can be used when developing a package to run the `docs/make.jl` file from
Documenter.jl and then serve the `docs/build` folder with LiveServer.jl.
This function assumes you are in the directory `[MyPackage].jl` with a
subfolder `docs`.

## Keyword Arguments

* `verbose=false`: boolean switch to make the server print information about
                   file changes and connections.
* `literate=nothing`: path to the folder containing the literate scripts; if
                      nothing then it's assumed there's no literate script. If
                      left empty, it's assumed that the literate scripts are in
                      `docs/src`. Any `xxx.jl` file in the literate location is
                      assumed to generate a `xxx.md` file which will not be
                      watched to avoid ∞ loop.
* `doc_env=false`: a boolean switch to make the server start by activating the
                   doc environment or not (i.e. the `Project.toml` in `docs/`).
* `skip_dir=""`: a subpath of `docs/` where modifications should not trigger
                 the generation of the docs, this is useful for instance if
                 you're using Weave and Weave generates some files in
                 `docs/src/examples` in which case you should set
                 `skip_dir=joinpath("docs","src","examples")`.
* `skip_dirs=[]`: same as `skip_dir` but for a list of such dirs. Takes
                  precedence over `skip_dir`.
* `foldername="docs"`: specify a different path for the content.
* `buildfoldername="build"`: specify a different path for the build.
* `makejl="make.jl"`: path of the script generating the documentation relative
                      to `foldername`.
* `host="127.0.0.1"`: where the server will start.
* `port=8000`: port number, an integer between 8000 (default) and 9000.
* `launch_browser=false`: specifies whether to launch a browser at the
                          localhost URL or not.
"""
function servedocs(;
            verbose::Bool=false,
            literate::Union{Nothing,String}=nothing,
            doc_env::Bool=false,
            skip_dir::String="",
            skip_dirs::Vector{String}=String[],
            foldername::String="docs",
            buildfoldername::String="build",
            makejl::String="make.jl",
            host::String="127.0.0.1",
            port::Int=8000,
            launch_browser::Bool=false
            )::Nothing
    # skip_dirs takes precedence over skip_dir
    if isempty(skip_dirs) && !isempty(skip_dir)
        skip_dirs = [skip_dir]
    end

    path2makejl = joinpath(foldername, makejl)

    # The file watcher is a default SimpleWatcher with a custom
    # callback
    docwatcher = SimpleWatcher()
    set_callback!(
        docwatcher,
        fp -> servedocs_callback!(
                docwatcher, fp, path2makejl,
                literate, skip_dirs, foldername, buildfoldername
        )
    )

    # Scan the folder and update the list of files to watch
    scan_docs!(docwatcher, foldername, path2makejl, literate)

    # activate the doc environment if required
    doc_env && Pkg.activate(joinpath(foldername, "Project.toml"))

    # trigger a first pass of Documenter (& possibly Literate)
    Main.include(abspath(path2makejl))

    # note the `docs/build` exists here given that if we're here it means
    # the documenter pass did not error and, therefore that a docs/build
    # has been generated.
    # So we can now serve.
    serve(
        docwatcher,
        host=host,
        port=port,
        dir=joinpath(foldername, buildfoldername),
        verbose=verbose,
        launch_browser=launch_browser
    )

    # when the serve loop is interrupted, de-activate the environment
    doc_env && Pkg.activate()
    return
end

#
# Miscellaneous utils
#

"""
    setverbose(b)

Set the verbosity of LiveServer to either `true` (showing messages upon events) or `false` (default).
"""
setverbose(b::Bool) = (VERBOSE.x = b)

"""
    example()

Simple function to copy an example website folder to the current working directory that can be
watched by the LiveServer to get an idea of how things work.

### Example

```julia
LiveServer.example()
cd("example")
serve()
```
"""
function example(; basedir="")
    isempty(basedir) && (basedir = pwd())
    cp(joinpath(dirname(dirname(pathof(LiveServer))), "example"), joinpath(basedir, "example"))
end

#
# Generate example repo for servedocs
#

const INDEX_MD = raw"""
    # Test
    A link to the [other page](/man/pg1/)
    """
const PG1_JL  = raw"""
    # # Test literate
    # We can include some code like so:
    f(x) = x^5
    f(5)
    """
const MAKE_JL = raw"""
    using Documenter, Literate
    src = joinpath(@__DIR__, "src")
    lit = joinpath(@__DIR__, "literate")
    for (root, _, files) ∈ walkdir(lit), file ∈ files
        splitext(file)[2] == ".jl" || continue
        ipath = joinpath(root, file)
        opath = splitdir(replace(ipath, lit=>src))[1]
        Literate.markdown(ipath, opath)
    end
    makedocs(
        sitename = "testlit",
        pages = ["Home" => "index.md",
                 "Other page" => "man/pg1.md"]
        )
    """

"""
    servedocs_literate_example(dir="servedocs_literate_example")

Generates a folder with the right structure for servedocs+literate example.
Requires Documenter and Literate.
You can then `cd` to that folder and use servedocs:

```
julia> using LiveServer, Documenter, Literate
julia> LiveServer.servedocs_literate_example()
julia> cd("servedocs_literate_example")
julia> servedocs(literate=joinpath("docs", "literate"))
```
"""
function servedocs_literate_example(dirname="servedocs_literate_example")
    isdir(dirname) && rm(dirname, recursive=true)
    mkdir(dirname)
    # folder structure
    src  = joinpath(dirname, "src")
    mkdir(src)
    write(joinpath(src, "$dirname.jl"), "module $dirname\n foo()=1\n end")
    docs = joinpath(dirname, "docs")
    mkdir(docs)
    src = joinpath(docs, "src")
    lit = joinpath(docs, "literate")
    man = joinpath(lit, "man")
    mkdir(src)
    mkdir(lit)
    mkdir(man)
    write(joinpath(src, "index.md"), INDEX_MD)
    write(joinpath(man, "pg1.jl"), PG1_JL)
    write(joinpath(docs, "make.jl"), MAKE_JL)
    return
end


text:
@testset "fallbacks" begin
    x = FallbackType()
    @test getobs(x, 3) == 1234
    @test numobs(x) == 5678
end

@testset "array" begin
    a = rand(2,3)
    @test numobs(a) == 3
    @test @inferred getobs(a, 1) == a[:,1]
    @test @inferred getobs(a, 2) == a[:,2]
    @test @inferred getobs(a, 1:2) == a[:,1:2]
end

@testset "0-dim SubArray" begin
    v = view([3], 1)
    @test @inferred(numobs(v)) === 1
    @test @inferred(getobs(v, 1)) === 3
    @test_throws BoundsError getobs(v, 2)
    @test_throws BoundsError getobs(v, 2:3)
end

@testset "named tuple" begin
    X2, Y2 = rand(2, 3), rand(3)
    dataset = (x=X2, y=Y2)
    @test numobs(dataset) == 3
    o = @inferred getobs(dataset, 2)
    @test o.x == X2[:,2]
    @test o.y == Y2[2]

    o = @inferred getobs(dataset, 1:2)
    @test o.x == X2[:,1:2]
    @test o.y == Y2[1:2]
end

@testset "dict" begin
    dataset = Dict("X" => X, "y" => y)
    @test numobs(dataset) == 15

    @test_broken @inferred getobs(dataset, 2) # not inferred
    o = getobs(dataset, 2)
    @test o["X"] == X[:,2]
    @test o["y"] == y[2]

    o = getobs(dataset, 1:2)
    @test o["X"] == X[:,1:2]
    @test o["y"] == y[1:2]
end

@testset "numobs" begin
    @test_throws MethodError numobs(X,X)
    @test_throws MethodError numobs(X,y)

    @testset "Array, SparseArray, and Tuple" begin
        @test_throws DimensionMismatch numobs((X,XX,rand(100)))
        @test_throws DimensionMismatch numobs((X,X'))
        for var in (Xs, ys, vars...)
            @test @inferred(numobs(var)) === 15
        end
        @test @inferred(numobs(())) === 0
    end

    @testset "SubArray" begin
        @test @inferred(numobs(view(X,:,:))) === 15
        @test @inferred(numobs(view(X,:,:))) === 15
        @test @inferred(numobs(view(XX,:,:,:))) === 15
        @test @inferred(numobs(view(XXX,:,:,:,:))) === 15
        @test @inferred(numobs(view(y,:))) === 15
        @test @inferred(numobs(view(Y,:,:))) === 15
    end
end

@testset "getobs" begin
    @testset "Array and Subarray" begin
        # access outside numobs bounds
        @test_throws BoundsError getobs(X, -1)
        @test_throws BoundsError getobs(X, 0)
        @test_throws BoundsError getobs(X, 16)
        for i in (2, 2:10, [2,1,4])
            @test getobs(XX, i) == XX[:, :, i]
        end
        for i in (2, 1:15, 2:10, [2,5,7], [2,1])
            @test typeof(getobs(Xv, i)) <: Array
            @test typeof(getobs(yv, i)) <: ((i isa Int) ? String : Array)
            @test all(getobs(Xv, i) .== X[:, i])
            @test getobs(Xv,i)  == X[:,i]
            @test getobs(X,i)   == X[:,i]
            @test getobs(XX,i)  == XX[:,:,i]
            @test getobs(XXX,i) == XXX[:,:,:,i]
            @test getobs(y,i)   == y[i]
            @test getobs(yv,i)  == y[i]
            @test getobs(Y,i)   == Y[:,i]
        end
    end

    @testset "SparseArray" begin
        @test typeof(getobs(Xs,2)) <: SparseVector
        @test typeof(getobs(Xs,1:5)) <: SparseMatrixCSC
        @test typeof(getobs(ys,2)) <: Float64
        @test typeof(getobs(ys,1:5)) <: SparseVector
        for i in (2, 1:15, 2:10, [2,5,7], [2,1])
            @test getobs(Xs,i) == Xs[:,i]
            @test getobs(ys,i) == ys[i]
        end
    end

    @testset "Tuple" begin
        # bounds checking correctly
        @test_throws BoundsError getobs((X,y), 16)
        # special case empty tuple
        @test @inferred(getobs((), 10)) === ()
        tx, ty = getobs((Xv, yv), 2:5)
        @test typeof(tx) <: Array
        @test typeof(ty) <: Array
        for i in (1:15, 2:10, [2,5,7], [2,1])
            @test_throws DimensionMismatch getobs((X', y), i)
            @test @inferred(getobs((X,y), i))  == (X[:,i], y[i])
            @test @inferred(getobs((X,yv), i)) == (X[:,i], y[i])
            @test @inferred(getobs((Xv,y), i)) == (X[:,i], y[i])
            @test @inferred(getobs((X,Y), i))  == (X[:,i], Y[:,i])
            @test @inferred(getobs((X,yt), i)) == (X[:,i], yt[:,i])
            @test @inferred(getobs((XX,X,y), i)) == (XX[:,:,i], X[:,i], y[i])
            # compare if obs match in tuple
            x1, y1 = getobs((X1,Y1), i)
            @test all(x1' .== y1)
            x1, y1, z1 = getobs((X1,Y1,sparse(X1)), i)
            @test all(x1' .== y1)
            @test all(x1 .== z1)
        end
        @test getobs((Xv,y), 2) == (X[:,2], y[2])
        @test getobs((X,yv), 2) == (X[:,2], y[2])
        @test getobs((X,Y), 2) == (X[:,2], Y[:,2])
        @test getobs((XX,X,y), 2) == (XX[:,:,2], X[:,2], y[2])
    end
end

@testset "getobs!" begin
    @testset "Array and Subarray" begin
        Xbuf = similar(X)
        getobs!(Xbuf, X)
        @test Xbuf == X
        # access outside numobs bounds
        @test_throws BoundsError getobs!(Xbuf, X, -1)
        @test_throws BoundsError getobs!(Xbuf, X, 0)
        @test_throws BoundsError getobs!(Xbuf, X, 151)

        buff = zeros(4)
        @test @inferred(getobs!(buff, X, 5)) == getobs(X, 5)

        buff = zeros(4, 8)
        @test @inferred(getobs!(buff, X, 3:10)) == getobs(X, 3:10)

        buff = zeros(20,30)
        @test @inferred(getobs!(buff, XX, 5)) == XX[:,:,5]

        buff = zeros(20, 30, 5)
        @test @inferred(getobs!(buff, XX, 6:10)) == XX[:,:,6:10]

        # string vector
        @test getobs!("setosa", y, 1) == "setosa"
        @test getobs!(nothing, y, 1) == "setosa"
    end

    @testset "SparseArray" begin
        # Sparse Arrays opt-out of buffer usage
        @test @inferred(getobs!(nothing, Xs, 1)) == getobs(Xs, 1)
        @test @inferred(getobs!(nothing, Xs, 5:10)) == getobs(Xs, 5:10)
        @test @inferred(getobs!(nothing, ys, 1)) === getobs(ys, 1)
        @test @inferred(getobs!(nothing, ys, 5:10)) == getobs(ys, 5:10)
    end

    @testset "Tuple" begin
        @test_throws MethodError getobs!((nothing,nothing), (X,y))
        @test getobs!((nothing,nothing), (X,y), 1:5) == getobs((X,y), 1:5)
        @test_throws MethodError getobs!((nothing,nothing,nothing), (X,y))
        xbuf = zeros(4,2)
        ybuf = ["foo", "bar"]
        @test_throws MethodError getobs!((xbuf,), (X,y))
        @test_throws MethodError getobs!((xbuf,ybuf,ybuf), (X,y))
        @test_throws DimensionMismatch getobs!((xbuf,), (X,y), 1:5)
        @test_throws DimensionMismatch getobs!((xbuf,ybuf,ybuf), (X,y), 1:5)
        @test @inferred(getobs!((xbuf,ybuf),(X,y), 2:3)) === (xbuf,ybuf)
        @test xbuf == getobs(X, 2:3)
        @test ybuf == getobs(y, 2:3)
        @test @inferred(getobs!((xbuf,ybuf),(X,y), [14,5])) === (xbuf,ybuf)
        @test xbuf == getobs(X, [14,5])
        @test ybuf == getobs(y, [14,5])

        @test getobs!((nothing,xbuf),(Xs,X), 3:4) == (getobs(Xs,3:4),xbuf)
        @test xbuf == getobs(X,3:4)
    end
end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
function call(::Types.Leave, args::Vector, storage, conn::IO)
    ip_addr = string(first(getpeername(conn)))
    client = storage.active_clients[ip_addr]
    channel = storage.active_channels[client.current_channel_id]

    if (client.current_channel_id == "1")
        return (write(conn, "You can't leave the lobby\n"))
    end
    client.current_channel_id = "1"
    write(conn, "You left the channel [$(channel.name)]\n")
end


text:
function check_cutval(A,svec,tvec,S)
    """
    Check the cut value associated with the set of nodes in graph A.
    Source and sink nodes are not explicitly included, but terminal
    edges weights are given in svec and tvec.
    """
    cut = 0
    n = size(A,1)
    Sbar = setdiff(1:n,S)   # complement set

    # terminal edge weight penalties
    cut += sum(tvec[S])
    cut += sum(svec[Sbar])
    AS = A[S,Sbar]
    cut += sum(AS)

    return cut
end

function clique_expansion_penalties(k)
    """
    Symmetric clique splitting penalty for hyperedge of size k.
    Only gives first r = floor(k/2) penalties
    """
    @assert(k > 1)
    r = floor(Int64,(k)/2)
    I = collect(0:k)
    w = I .* (k .- I)
    w = w[1:r+1]
    return w
end

function load_superpixel_cliques(dataset,lambda,numsuperpix=500)
    """
    Load superpixel regions for a certain image dataset from experiements of
    Jegalka et al.

    Save the superpixel regions in hyperedge list format, with clique
    expansion splitting penalties, weighted by lambda
    (use lambda values as given previously by Jegelka et al.)
    """

    # Load superpixel region information
    M = matread("../data/data_reflections/$(dataset)_p$numsuperpix.mat")
    sup_img = M["sup_img"]

    col = size(sup_img,2)                   # number of columns
    row = size(sup_img,1)                   # number of rows
    node = (i,j)->i + (j-1)*row             # go from (i,j) node ID to linear node ID
    region = length(unique(sup_img))        # number of superpixels

    # Save superpixels in alternate format
    Edges = Vector{Vector{Int64}}()
    for e = 1:region
        # each superpixel defines a hyperedge
        push!(Edges,Vector{Int64}())
    end

    EdgeLengths = zeros(Int64,region)
    for i = 1:row
        for j = 1:col
            v = node(i,j)    # linear index of node at gridpoint (i,j)

            # the superpixel, i.e., hyperedge, it's in
            e = round(Int64,sup_img[i,j])
            push!(Edges[e],v)
            EdgeLengths[e] += 1
        end
    end

    # Discard cliques with less than 2 nodes, they do nothing
    keep = findall(x->x>1,EdgeLengths)
    Edges = Edges[keep]
    EdgesW = Vector{Vector{Float64}}()
    for ei = 1:length(Edges)
        # We'll use clique expansion splitting penalties,
        # scaled with previously determined weights of Jegelka et al.
        e = Edges[ei]
        push!(EdgesW,lambda*clique_expansion_penalties(length(e)))
    end

    return Edges, EdgesW
end

function eval_fS_image(Edges,EdgesW,svec,tvec,A,S)
    """
    For a set S, evaluate the objective from a decomposable
    submodular function involving terminal edges as defined
    in svec and tvec, pairwise non-terminal edges
    as defined in A, and higher-order terms from Edges
    and EdgesW.
    """
    n = size(A,1)
    eS = zeros(n)
    eS[S] = ones(length(S))
    Sbar = setdiff(1:n,S)   # complement set

    # terminal edge weight penalties
    cutval = sum(tvec[S])
    cutval += sum(svec[Sbar])

    # interior edge penalties
    AS = A[S,Sbar]
    cutval += sum(AS)

    # higher-order cut penalties
    for i = 1:length(Edges)
        e = Edges[i]
        w = EdgesW[i]
        t = round(Int64,sum(eS[e]))

        # do this because it's a symmetric splitting function,
        # so w is only half the size
        t = min(t,length(e)-t)
        cutval += w[t+1]
    end

    return cutval

end

function sweep_cut_f(Edges,EdgesW,svec,tvec,A,y)
    if length(y) > 10
        println("This code is inefficient. Not recommended for use for vectors this big.")
        return
    end
     Order = sortperm(y,rev = true)
     n = length(y)
     @assert(n == size(A,1))
     Sbest = [Order[1]]
     fbest = eval_fS_image(Edges,EdgesW,svec,tvec,A,Sbest)
     for i = 2:n
         S = Order[1:i]
         f = eval_fS_image(Edges,EdgesW,svec,tvec,A,S)
         if f < fbest
             fbest = f
             Sbest = S
         end
     end

     return Sbest, fbest
 end

function Objective_from_Primals(Xmat,sweepi,Edges,EdgesW,svec,tvec,An)
    """
    Given a set of primal variables from a continuous DSFM method,
    stored as columns of Xmat, find the objective score resulting
    from the best level set S_lam = {i : x_i >= lam}.

    The value of sweepi tells us the best level set as
    determined by the Matlab/C++ code that Xmat comes from.
    """
    sourceside = sum(svec)
    t = size(Xmat,2)
    @assert(t == length(sweepi))
    Fvals = zeros(t)
    for i = 1:t
        x = Xmat[:,i]
        Order = sortperm(x,rev=true)
        wherei = findfirst(x->x==sweepi[i]+1,Order)
        S = Order[1:wherei]
        Fvals[i] = eval_fS_image(Edges,EdgesW,svec,tvec,An,S)
    end
    return Fvals
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
using Blink
using Tachyons
using JSExpr
using CSSUtil

w = Window()

function driver(w)
    dt = get_data_type()
    @show dt

    ui1 = get_mod_mode_network()
    ui2 = get_mod_mode_raster()

    nextui = ui1
    nextui = map(dt["value"]) do v
        if v == "Network"
            ui1
        else
            ui2
        end
    end
    @show nextui

    nextnextui = ""
    map(nextui) do v
        s = v["value"]
        map(s) do x
            if x == "Pairwise"
                nextnextui = "Pairwise UI"
            elseif x == "Advanced"
                nextnextui = "Advanced UI"
            elseif x == "One To All"
                nextnextui = "One To All UI"
            elseif x == "All To One"
                nextnextui = "All To One UI"
            end
        end
    end
    @show nextnextui

    page = vbox(dt,
                 nextui)

    body!(w, page)
end

function get_data_type()

    data_type_prompt = "Step 1: Choose your input data type: "

    # Select between raster and network
    data_type = vbox(Node(:div, data_type_prompt,
                          attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                 Node(:select, "Select Data Type",
                 Node(:option, "Raster"),
                 Node(:option, "Network"), id = "dt",
                 attributes = Dict(:style => "margin-top: 12px; margin-bottom: 12px")))

    s = Scope()
    s.dom = data_type
    onimport(s, JSExpr.@js function ()
                 @var el = this.dom.querySelector("#dt")
                 el.onchange = (function ()
                        $(s["value"])[] = el.value
                    end)
             end)

    s
end

function get_mod_mode_network()

    mod_mode_prompt = "Step 2: Choose a Modelling Mode: "
    mod_mode = vbox(Node(:div, mod_mode_prompt,
                          attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                 Node(:select, "Select Modelling Mode",
                 Node(:option, "Pairwise"),
                 Node(:option, "Advanced"), id = "modelling",
                 attributes = Dict(:style => "margin-top: 12px; margin-bottom: 12px")))

    s = Scope()
    s.dom = mod_mode
    onimport(s, JSExpr.@js function ()
                 @var el = this.dom.querySelector("#modelling")
                 el.onchange = (function ()
                        $(s["value"])[] = el.value
                    end)
             end)

    s
end

function get_mod_mode_raster()

    mod_mode_prompt = "Step 2: Choose a Modelling Mode: "
    mod_mode = vbox(Node(:div, mod_mode_prompt,
                          attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                 Node(:select, "Select Modelling Mode",
                 Node(:option, "Pairwise"),
                 Node(:option, "Advanced"),
                 Node(:option, "One To All "),
                 Node(:option, "All To One"), id = "modelling",
                 attributes = Dict(:style => "margin-top: 12px; margin-bottom: 12px")))

    s = Scope()
    s.dom = mod_mode
    onimport(s, JSExpr.@js function ()
                 @var el = this.dom.querySelector("#modelling")
                 el.onchange = (function ()
                        $(s["value"])[] = el.value
                    end)
             end)
    s
end

driver(w)