text:
# lame wrappers of required functions to write mp3 files

"""represents the C pointer lame_global_flags*. used by all LAME functions"""
const LAME = Ptr{Void}

"""meaning of LAME return codes, usually only relevant to lame_encode_buffer families"""
LAME_ERRORS = Dict{Int, String}(
    -1 => "mp3buf was too small",
    -2 => "malloc() problem",
    -3 => "lame_init_params() not called",
    -4 => "psycho acoustic problems"
)

LAME_JOINT_STEREO = 1
LAME_MONO = 3

lame_error_string(ret) = get(LAME_ERRORS, ret, "unknown error code $ret")

"""return the version of LAME library"""
function get_lame_version()
    version = ccall((:get_lame_version, libmp3lame), Ptr{UInt8}, ())
    bytestring(version)
end

"""initialize LAME encoder and return the LAME object"""
function lame_init()
    lame = ccall((:lame_init, libmp3lame), LAME, ())
    if lame == C_NULL
        error("Could not initialize LAME encoder")
    end
    lame
end

"""set the number of samples of the input signal"""
function lame_set_num_samples(lame::LAME, nsamples::Integer)
    ret = ccall((:lame_set_num_samples, libmp3lame), Cint,
                (LAME, Culong), lame, nsamples)

    if ret != 0
        error("Unknown error while calling lame_set_num_samples(): $ret")
    end
end

"""set the sampling rate of the input signal in Hz; default = 44100"""
function lame_set_in_samplerate(lame::LAME, samplerate::Integer)
    ret = ccall((:lame_set_in_samplerate, libmp3lame), Cint,
                (LAME, Cint), lame, samplerate)

    if ret != 0
        error("Unknown error while calling lame_set_in_samplerate(): $ret")
    end
end

"""set the number of channels in the input signal; default = 2"""
function lame_set_num_channels(lame::LAME, nchannels::Integer)
    ret = ccall((:lame_set_num_channels, libmp3lame), Cint,
                (LAME, Cint), lame, nchannels)

    if ret != 0
        error("Error while calling lame_set_num_channels(); should be 1 or 2")
    end
end

"""set the sampling rate of the output signal in Hz"""
function lame_set_out_samplerate(lame::LAME, samplerate::Integer)
    ret = ccall((:lame_set_out_samplerate, libmp3lame), Cint,
                (LAME, Cint), lame, samplerate)

    if ret != 0
        error("Unknown error while calling lame_set_out_samplerate(): $ret")
    end
end

"""set the quality for algorithm selection; should be [0..9]"""
function lame_set_quality(lame::LAME, quality::Integer)
    ret = ccall((:lame_set_quality, libmp3lame), Cint,
                (LAME, Cint), lame, quality)

    if ret != 0
        error("Unknown error while calling lame_set_quality(): $ret")
    end
end

"""set the channel encoding mode; usually 1 (LAME_JOINT_STEREO) or 3 (LAME_MONO)"""
function lame_set_mode(lame::LAME, mode::Integer)
    ret = ccall((:lame_set_mode, libmp3lame), Cint,
                (LAME, Cint), lame, mode)

    if ret != 0
        error("Unknown error while calling lame_set_mode(): $ret")
    end
end

"""set the bitrate of the mp3 output, in kbps"""
function lame_set_brate(lame::LAME, bitrate::Integer)
    ret = ccall((:lame_set_brate, libmp3lame), Cint,
                (LAME, Cint), lame, bitrate)

    if ret != 0
        error("Unknown error while calling lame_set_brate(): $ret")
    end
end

"""enable VBR mode"""
function lame_set_VBR(lame::LAME)
    # call the c function with vbr_mtrh = 4
    ret = ccall((:lame_set_VBR, libmp3lame), Cint,
                (LAME, Cint), lame, 4)

    if ret != 0
        error("Unknown error while calling lame_set_VBR(): $ret")
    end
end

"""set VBR quality; 0.000=highest, 9.999=lowest; """
function lame_set_VBR_quality(lame::LAME, quality::Number)
    # call the c function with vbr_mtrh = 4
    ret = ccall((:lame_set_VBR_quality, libmp3lame), Cint,
                (LAME, Cfloat), lame, Float32(quality))

    if ret != 0
        error("Error while calling lame_set_VBR_quality(); is the quality in [0,9]?")
    end
end

"""set internal parameters based on provided settings"""
function lame_init_params(lame::LAME)
    ccall((:lame_init_params, libmp3lame), Cint, (LAME,), lame)
end

"""print the configuration"""
function lame_print_config(lame::LAME)
    ccall((:lame_print_config, libmp3lame), Cint, (LAME,), lame)
end

"""write 16-bit fixed-point PCM samples"""
function lame_encode_buffer!(lame::LAME, left::Ptr{PCM16Sample}, right::Ptr{PCM16Sample},
                             nsamples::Integer, mp3buf::Ptr{UInt8}, mp3buf_size::Integer)
    ret = ccall((:lame_encode_buffer, libmp3lame), Cint,
                (LAME, Ptr{PCM16Sample}, Ptr{PCM16Sample}, Cint, Ptr{UInt8}, Cint),
                lame, left, right, Cint(nsamples), mp3buf, Cint(mp3buf_size))

    if ret < 0
        error("Error while calling lame_encode_buffer(): ", lame_error_string(ret))
    end

    ret
end

"""write interleaved 16-bit fixed-point PCM samples"""
function lame_encode_buffer_interleaved!(lame::LAME, pcm::Ptr{PCM16Sample},
                                         nsamples::Integer, mp3buf::Ptr{UInt8}, mp3buf_size::Integer)
    ret = ccall((:lame_encode_buffer_interleaved, libmp3lame), Cint,
                (LAME, Ptr{PCM16Sample}, Cint, Ptr{UInt8}, Cint),
                lame, pcm, Cint(nsamples), mp3buf, Cint(mp3buf_size))

    if ret < 0
        error("Error while calling lame_encode_buffer_interleaved(): ", lame_error_string(ret))
    end

    ret
end

"""flush any internal PCM buffers to finalize an MP3 file"""
function lame_encode_flush!(lame::LAME, mp3buf::Ptr{UInt8}, mp3buf_size::Integer)
    ret = ccall((:lame_encode_flush, libmp3lame), Cint,
                (LAME, Ptr{UInt8}, Cint), lame, mp3buf, mp3buf_size)

    if ret < 0
        error("Error while calling lame_encode_flush(): ", lame_error_string(ret))
    end

    ret
end

"""flush any internal PCM buffers to finalize an MP3 file, without ID3v1 tags"""
function lame_encode_flush_nogap(lame::LAME, mp3buf::Ptr{UInt8}, mp3buf_size::Integer)
    ret = ccall((:lame_encode_flush_nogap, libmp3lame), Cint,
                (LAME, Ptr{UInt8}, Cint), lame, mp3buf, mp3buf_size)

    if ret < 0
        error("Error while calling lame_encode_flush_nogap(): ", lame_error_string(ret))
    end

    ret
end

"""free all LAME buffers to finalize the encoder"""
function lame_close(lame::LAME)
    ccall((:lame_close, libmp3lame), Cint, (LAME,), lame)
end

"""initialize ID3 tag"""
function id3tag_init(lame::LAME)
    ccall((:id3tag_init, libmp3lame), Void, (LAME,), lame)
end

"""tell LAME to add ID3v2"""
function id3tag_add_v2(lame::LAME)
    ccall((:id3tag_add_v2, libmp3lame), Void, (LAME,), lame)
end

"""force using ID3v2 and not ID3v1"""
function id3tag_v2_only(lame::LAME)
    ccall((:id3tag_v2_only, libmp3lame), Void, (LAME,), lame)
end

"""set ID3v2 unicode tag for title"""
function id3tag_set_title(lame::LAME, title::AbstractString)
    id3tag_set_textinfo_utf16(lame, "TIT2", title)
end

"""set ID3v2 unicode tag for artist"""
function id3tag_set_artist(lame::LAME, artist::AbstractString)
    id3tag_set_textinfo_utf16(lame, "TPE1", artist)
end

"""set ID3v2 unicode tag for album"""
function id3tag_set_album(lame::LAME, album::AbstractString)
    id3tag_set_textinfo_utf16(lame, "TALB", album)
end

"""set ID3v2 unicode tag for year"""
function id3tag_set_year(lame::LAME, year::AbstractString)
    id3tag_set_textinfo_utf16(lame, "TYER", year)
end

"""set ID3v2 unicode tag for comment"""
function id3tag_set_comment(lame::LAME, comment::AbstractString)
    id3tag_set_comment_utf16(lame, comment)
end

"""get the ID3v2 buffer"""
function lame_get_id3v2_tag(lame::LAME, buffer::Array{UInt8}, size::Integer)
    ccall((:lame_get_id3v2_tag, libmp3lame), Csize_t,
          (LAME, Ptr{UInt8}, Csize_t),
          lame, buffer, size)
end

"""prevent LAME from writing ID3 tags"""
function lame_set_write_id3tag_automatic(lame::LAME, v::Integer)
    ccall((:lame_set_write_id3tag_automatic, libmp3lame), Void,
          (LAME, Cint), lame, v)
end

"""set ID3v2 tag of given id"""
function id3tag_set_textinfo_utf16(lame::LAME, id::AbstractString, text::AbstractString)
    # utf-16 representation of the text, along with the BOM
    data = [0xfeff; utf16(text).data]

    ret = ccall((:id3tag_set_textinfo_utf16, libmp3lame), Cint,
                (LAME, Cstring, Ptr{Cushort}),
                lame, id, data)

    if ret != 0
        error("Unexpected error $ret while calling id3tag_set_textinfo_utf16")
    end
end

"""set ID3v2 comment"""
function id3tag_set_comment_utf16(lame::LAME, text::AbstractString)
    # utf-16 representation of the text, along with the BOM
    data = [0xfeff; utf16(text).data]

    ret = ccall((:id3tag_set_comment_utf16, libmp3lame), Cint,
                (LAME, Ptr{Cchar}, Ptr{Cushort}, Ptr{Cushort}),
                lame, C_NULL, C_NULL, data)

    if ret != 0
        error("Unexpected error $ret while calling id3tag_set_comment_utf16")
    end
end


text:
module ODEPrecompileTest
using ModelingToolkit

function system(; kwargs...)
    # Define some variables
    @parameters t σ ρ β
    @variables x(t) y(t) z(t)
    D = Differential(t)

    # Define a differential equation
    eqs = [D(x) ~ σ * (y - x),
        D(y) ~ x * (ρ - z) - y,
        D(z) ~ x * y - β * z]

    @named de = ODESystem(eqs)
    return ODEFunction(de, [x, y, z], [σ, ρ, β]; kwargs...)
end

# Build an ODEFunction as part of the module's precompilation. These cases
# will not work, because the generated RGFs are put into the ModelingToolkit cache.
const f_bad = system()
const f_noeval_bad = system(; eval_expression = false)

# Setting eval_expression=false and eval_module=[this module] will ensure
# the RGFs are put into our own cache, initialised below.
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)
const f_noeval_good = system(; eval_expression = false, eval_module = @__MODULE__)
end

ext: jl


text:
using LatentOrderBookModel

function main(output = "stdout")
    parsed_args = parse_commandline()
    rdpp = ReactionDiffusionPricePaths(
        parsed_args["num_paths"],
        parsed_args["T"],
        parsed_args["p₀"],
        parsed_args["M"],
        parsed_args["β"],
        parsed_args["L"],
        parsed_args["D"],
        parsed_args["σ"],
        parsed_args["nu"],
        SourceTerm(parsed_args["λ"], parsed_args["μ"]),
    )
    if output == "stdout"
        lob_densities,
        raw_price_paths,
        sample_price_paths,
        P⁺s,
        P⁻s = rdpp(parsed_args["SEED"])
        print(sample_price_paths[:, 1])
    else
        return rdpp(parsed_args["SEED"])
    end

end

if "" != PROGRAM_FILE && occursin(PROGRAM_FILE, @__FILE__)
    main()
end


text:
# implement dhc of degree n where n=0 means random search, n=1 means
# use one previous direction etc.

using Knet,JLD
include("cheb.jl")
@useifgpu CUDArt
# using ProgressMeter: @showprogress, Progress, next!
using Knet: Net, params, axpy!, copysync!, BaseArray
typealias BaseVector{T} BaseArray{T,1}

if !isdefined(:MNIST)
    include(Pkg.dir("Knet/examples/mnist.jl"))
    mnist100 = minibatch(MNIST.xtrn, MNIST.ytrn, 100)
    gpu_ytrn = CudaArray(MNIST.ytrn)
    gpu_xtrn = CudaArray(MNIST.xtrn)
    gpu_ytst = CudaArray(MNIST.ytst)
    gpu_xtst = CudaArray(MNIST.xtst)
end

@knet function model02(x; winit=Gaussian(0,.1), hidden=64, f1=:relu, o...)
    h    = wf(x; out=hidden, f=f1, winit=winit)
    return wf(h; out=10, f=:soft, winit=winit)
end

cosine(x::CudaArray,y::CudaArray)=CUBLAS.dot(x,y)/(vecnorm(x)*vecnorm(y))
cosine(x,y)=dot(x,y)/(vecnorm(x)*vecnorm(y))
setseed(42)
gc()
@gpu !isdefined(:w64) && (w64 = scale!(0.01, randn!(CudaArray(Float32, 794*64))))
@gpu !isdefined(:w32) && (w32 = scale!(0.01, randn!(CudaArray(Float32, 794*32))))
!isdefined(:w64cpu) && (w64cpu = scale!(0.01, randn!(Array(Float32, 794*64))))
!isdefined(:w32cpu) && (w32cpu = scale!(0.01, randn!(Array(Float32, 794*32))))

@gpu function initmodel{T}(w::CudaArray{T}, g=nothing; X=784, Y=10, H=div(length(w),X+Y), lr=0.5)
    global gpumodel
    xgpu = Knet.gpu(); Knet.gpu(true)
    if !isdefined(:gpumodel) || !isa(gpumodel, Knet.Net)
        gpumodel = compile(:model02; hidden=H)
        setp(gpumodel, lr=lr)
        forw(gpumodel, zeros(T, X, 1))
        back(gpumodel, zeros(T, Y, 1), softloss)
    end
    gpumodel.reg[2].out = gpumodel.reg[2].out0 = CudaArray(w.ptr, (H,X), w.dev)
    gpumodel.reg[5].out = gpumodel.reg[5].out0 = CudaArray(w.ptr+H*X*sizeof(T), (Y,H), w.dev)
    if g != nothing
        gpumodel.reg[2].dif = gpumodel.reg[2].dif0 = CudaArray(g.ptr, (H,X), g.dev)
        gpumodel.reg[5].dif = gpumodel.reg[5].dif0 = CudaArray(g.ptr+H*X*sizeof(T), (Y,H), g.dev)
    end
    Knet.gpu(xgpu)
    return gpumodel
end

submatrix(a, i, m, n)=pointer_to_array(pointer(a, i), (m, n))

function initmodel{T}(w::Array{T}, g=nothing; X=784, Y=10, W=length(w), H=div(W,X+Y), lr=0.5)
    global cpumodel
    xgpu=Knet.gpu(); Knet.gpu(false)
    if !isdefined(:cpumodel) || !isa(cpumodel, Knet.Net)
        cpumodel = compile(:model02; hidden=H)
        setp(cpumodel, lr=lr)
        forw(cpumodel, zeros(T, X, 1))
        back(cpumodel, zeros(T, Y, 1), softloss)
    end
    if pointer(w) != pointer(cpumodel.reg[2].out0)
        cpumodel.reg[2].out = cpumodel.reg[2].out0 = submatrix(w, 1, H, X) # CudaArray(w.ptr, (H,X), w.dev)
        cpumodel.reg[5].out = cpumodel.reg[5].out0 = submatrix(w, 1+H*X, Y, H) # CudaArray(w.ptr+H*X*sizeof(T), (Y,H), w.dev)
    end
    if g != nothing && pointer(g) != pointer(cpumodel.reg[2].dif0)
        cpumodel.reg[2].dif = cpumodel.reg[2].dif0 = submatrix(g, 1, H, X) # CudaArray(g.ptr, (H,X), g.dev)
        cpumodel.reg[5].dif = cpumodel.reg[5].dif0 = submatrix(g, 1+H*X, Y, H) # CudaArray(g.ptr+H*X*sizeof(T), (Y,H), g.dev)
    end
    Knet.gpu(xgpu)
    return cpumodel
end

# function wtest2(w, loss=softloss)
#     f = initmodel(w)
#     y = forw(f, MNIST.xtst)
#     return loss(y, MNIST.ytst)
# end

# function wgrad(w, g=similar(w); x=MNIST.xtrn, y=MNIST.ytrn, loss=softloss)
#     f = initmodel(w, g)
#     ypred = forw(f, x)
#     back(f, y, loss)
#     return loss(ypred, y)
# end

function weval(w, g=nothing; x=nothing, y=nothing, loss=softloss)
    f = initmodel(w, g)
    ypred = forw(f, x)
    g != nothing && back(f, y, loss)
    return loss(ypred, y)
end

weval(w::Array,g=nothing)=(Knet.gpu(false);weval(w,g;x=MNIST.xtrn,y=MNIST.ytrn))
weval(w::CudaArray,g=nothing)=(Knet.gpu(true);weval(w,g;x=gpu_xtrn,y=gpu_ytrn))

function weval100(w, g=nothing)
    global _weval1, _weval2, _weval3
    dims = (794*64, length(w))
    if !isdefined(:_weval1) || size(_weval1) != dims
        _weval1 = scale!(1/sqrt(dims[1]), randn!(similar(w, dims)))
        _weval2 = similar(w, (dims[1],))
        _weval3 = similar(w, (dims[1],))
    end
    A_mul_B!(_weval2, _weval1, w)
    if g == nothing
        return weval(_weval2)
    else
        val = weval(_weval2, _weval3)
        At_mul_B!(g, _weval1, _weval3)
        return val
    end
end

# weval(w::Array;o...)=weval(CudaArray(w);o...)
# weval(w::Array,g::Array;o...)=(gg=CudaArray(g);f=weval(CudaArray(w),gg;o...);copy!(g,gg);f)

function wtest(w, loss=softloss; x=gpu_xtrn, y=gpu_ytrn)
    f = initmodel(w)
    ypred = forw(f, x)
    return loss(ypred, y)
end

function gradest0{T}(f::Function, x::BaseVector{T}; rscale=0.005, nkeep=5)
    g0,g1,x0
    x0 = copy(x)
    g0 = similar(x)
    f0 = f(x0, g0)
    xdims = length(x)
    xkeep = zeros(T, xdims, nkeep)
    fkeep = zeros(T, 1, nkeep)
    r = similar(x)
    for i=1:nkeep
        f1 = f(axpy!(rscale, randn!(r), copy!(x,x0)))
        fkeep[i] = f1
        copysync!(xkeep, 1+(i-1)*xdims, x, 1, xdims)
    end
    xdiff = broadcast(-, xkeep, to_host(x0)) # X,N
    fdiff = fkeep - f0                       # 1,N
    g1 = vec(fdiff / xdiff)
    g0 = to_host(g0)
    n0 = vecnorm(g0)
    n1 = vecnorm(g1)
    cosine = dot(g0,g1)/(n0*n1)
    println((:nkeep,nkeep,:rscale,rscale,:cos,cosine,:n0,n0,:n1,n1))
end

function gradest1{T}(f::Function, x::BaseVector{T}; rscale=0.005, nkeep=5)
    g0,g1,x0
    x0 = copy(x)
    g0 = similar(x)
    f0 = f(x0, g0)
    xdims = length(x)
    xkeep = zeros(T, xdims, nkeep)
    fkeep = zeros(T, 1, nkeep)
    r = similar(x)
    for i=1:nkeep
        f1 = f(axpy!(rscale, randn!(r), copy!(x,x0)))
        fkeep[i] = f1
        copysync!(xkeep, 1+(i-1)*xdims, x, 1, xdims)
    end
    g1 = vec(fkeep / xkeep)
    g0 = to_host(g0)
    n0 = vecnorm(g0)
    n1 = vecnorm(g1)
    cosine = dot(g0,g1)/(n0*n1)
    println((:nkeep,nkeep,:rscale,rscale,:cos,cosine,:n0,n0,:n1,n1))
end

function gradest2{T}(f::Function, x::BaseVector{T}; rscale=0.05, nkeep=5)
    g0,g1,x0
    x0 = copy(x)
    g0 = similar(x)
    f0 = f(x0, g0)
    xdims = length(x)
    xkeep = ones(T, 1+xdims, nkeep)
    fkeep = zeros(T, 1, nkeep)
    r = similar(x)
    for i=1:nkeep
        f1 = f(axpy!(rscale, randn!(r), copy!(x,x0)))
        fkeep[i] = f1
        copysync!(xkeep, 2+(i-1)*(1+xdims), x, 1, xdims)
    end
    g1 = vec(fkeep / xkeep)[2:end]
    g0 = to_host(g0)
    n0 = vecnorm(g0)
    n1 = vecnorm(g1)
    cosine = dot(g0,g1)/(n0*n1)
    println((:nkeep,nkeep,:rscale,rscale,:cos,cosine,:n0,n0,:n1,n1))
end

function gradest3{T}(f::Function, x::BaseVector{T}; rscale=0.005, nkeep=5)
    g0,g1,x0,xdiff,fdiff
    x0 = copy(x)
    g0 = similar(x)
    f0 = f(x0, g0)
    xdims = length(x)
    xkeep = zeros(T, xdims, nkeep+1)
    fkeep = zeros(T, 1, nkeep+1)
    copysync!(xkeep, 1, x0, 1, xdims)
    fkeep[1] = f0
    npair = div(nkeep*(nkeep+1),2)
    xdiff = zeros(T, xdims, npair)
    fdiff = zeros(T, 1, npair)
    idiff = 0
    r = similar(x)
    for i=2:nkeep+1
        f1 = f(axpy!(rscale, randn!(r), copy!(x,x0)))
        fkeep[i] = f1
        copysync!(xkeep, 1+(i-1)*xdims, x, 1, xdims)
        for j=1:i-1
            dx1 = xkeep[:,i] - xkeep[:,j]
            copysync!(xdiff, 1+idiff*xdims, dx1, 1, xdims)
            fdiff[idiff+=1] = f1 - fkeep[j]
        end
    end
    g1 = vec(fdiff / xdiff)
    g0 = to_host(g0)
    n0 = vecnorm(g0)
    n1 = vecnorm(g1)
    cosine = dot(g0,g1)/(n0*n1)
    println((:nkeep,nkeep,:rscale,rscale,:cos,cosine,:n0,n0,:n1,n1))
end

# MC estimation of gradient.  cos(alpha)=sqrt(k/n) almost exactly.

function gradest4{T}(f::Function, x::BaseVector{T}; rscale=0.0001, nsample=256, lr=1.0, l2=0.0)
    x0 = copy(x)
    g0 = similar(x)
    f0 = f(x0, g0)
    n0 = vecnorm(g0)
    nx = vecnorm(x0)
    @show CUBLAS.dot(x0,g0)/(n0*nx)
    r = similar(x)
    g = fill!(similar(x),0)
    cs = 0; n1 = nf = f1 = f2 = 0; np = 1
    for i=1:nsample
        scale!(rscale, randn!(r))
        nr = vecnorm(r)
        f1 = f(axpy!(1, r, copy!(x,x0)))
        f2 = f0 + CUBLAS.dot(g, r)
        scale!(1-lr*l2, g)
        axpy!(-lr*(f2-f1)/(nr*nr), r, g) # NLMS algorithm normalizes with nr*nr (wikipedia, which also has convergence proofs)
        n1 = vecnorm(g)
        cs = CUBLAS.dot(g0,g)/(n0*n1)
        if (nf+=1) >= np
            println((nf,:cos,cs,:n1n0,n1/n0,:nfnd,sqrt(nf/length(x)),:f2f1,f2-f1,:f1f0,f1-f0,:f0,f0,:n0,n0,:n1,n1,:nr,nr))
            np *= 2
        end
    end
    println((nf,:cos,cs,:n1n0,n1/n0,:nfnd,sqrt(nf/length(x)),:f2f1,f2-f1,:f1f0,f1-f0,:f0,f0,:n0,n0,:n1,n1,:nr,nr))
end

# Try random steps biased toward the current gradient estimate rather
# than radially symmetric.  Does not improve.

function gradest5{T}(f::Function, x::BaseVector{T}; rscale=0.0001, nsample=128, lr=1.0, l2=0.0, gbias=1.0)
    x0 = copy(x)
    g0 = similar(x)
    f0 = f(x0, g0)
    ng0 = vecnorm(g0)
    nx0 = vecnorm(x0)
    @show CUBLAS.dot(x0,g0)/(ng0*nx0)
    r = similar(x)
    nr0 = sqrt(length(x))
    g = fill!(similar(x),0)
    ng = cs = nf = f1 = f2 = nr = 0; np = 1
    for i=1:nsample
        randn!(r)
        ng > 0 && axpy!(gbias*nr0/ng, g, r)
        scale!(rscale, r)
        nr = vecnorm(r)
        f1 = f(axpy!(1, r, copy!(x,x0)))
        f2 = f0 + CUBLAS.dot(g, r)
        l2 != 0 && scale!(1-lr*l2, g)
        axpy!(-lr*(f2-f1)/(nr*nr), r, g) # NLMS algorithm normalizes with nr*nr (wikipedia, which also has convergence proofs)
        ng = vecnorm(g)
        cs = CUBLAS.dot(g0,g)/(ng0*ng)
        if (nf+=1) >= np
            println((i,:cos,cs,:ngng0,ng/ng0,:f2f1,f2-f1,:f1f0,f1-f0,:f0,f0,:ng0,ng0,:ng,ng,:nr,nr))
            np *= 2
        end
    end
    println((nsample,:cos,cs,:ngng0,ng/ng0,:f2f1,f2-f1,:f1f0,f1-f0,:f0,f0,:ng0,ng0,:ng,ng,:nr,nr))
end

function train(f, data, loss; epochs=100)
    p = Progress(epochs, dt=1, barlen=40)
    for epoch=1:epochs
        for (x,ygold) in data
            ypred = forw(f, x)
            back(f, ygold, loss)
            update!(f)
        end
        zloss = test(f, data, zeroone)
        sloss = test(f, data, softloss)
        p.desc = @sprintf("%d/%d %.6f %.6f ", epoch, epochs, sloss, zloss)
        next!(p)
    end
end

function test(f, data, loss)
    sumloss = numloss = 0
    for (x,ygold) in data
        ypred = forw(f, x)
        sumloss += loss(ypred, ygold)
        numloss += 1
    end
    return sumloss/numloss
end

# inc!(d,k)=(d[k]=1+get(d,k,0))
# avg!(d,k,v;a=0.01)=(d[k]=a*v+(1-a)*get(d,k,v))

function dhc0{T}(f::Function, x::BaseVector{T}; grow=2.0, rmin=1e-5, rscale=0.05, ftol=1e-6, smin=0.05, smax=0.45, salpha=0.1)
    r = similar(x)
    x0 = copy(x)
    f0 = f(x)
    savg = (smin+smax)/2
    nf = np = 1
    while f0 > ftol
        f1 = f(axpy!(rscale, randn!(r), copy!(x,x0)))
        if f1 <= f0
            f0 = f1
            copy!(x0, x)
            savg = (1-salpha)*savg + salpha
        else
            savg = (1-salpha)*savg
        end
        if savg > smax
            # print("+")
            rscale *= grow
            savg = (smin+smax)/2
        elseif savg < smin
            # print("-")
            rscale /= grow
            rscale < rmin && break
            savg = (smin+smax)/2
        end
        if (nf+=1) >= np
            println((nf,:f0,f0,:rscale,rscale,:savg,savg,:x1,x0[1],:err,wtest(x0, zeroone)))
            np *= 2
        end
        if nf % 1000 == 0
            gc()
            # println((:mem, Knet.gpumem(), :gc, (gc();gpusync();Knet.gpumem())))
        end
    end
    println((nf,:f0,f0,:rscale,rscale,:savg,savg,:x1,x0[1],:err,wtest(x0, zeroone)))
    copy!(x,x0)
    return f0
end

Knet.copysync!{T}(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer)=copy!(dest,doffs,src,soffs,n)
# @gpu import CUDArt: to_host
# to_host(x)=x

function try_step{T}(f::Function, x::BaseVector{T}, step::BaseVector{T}, scale::Float64, x0::BaseVector{T}, f0::Float64, xkeep::Matrix{T}, fkeep::Matrix{T}, ikeep::Int)
    xdims = length(x)
    f1 = f(axpy!(scale, step, copy!(x, x0)))
    if f1 > f0
        fkeep[ikeep] = f1
        copysync!(xkeep, 1+(ikeep-1)*xdims, x, 1, xdims)
    else
        fkeep[ikeep] = f0
        copysync!(xkeep, 1+(ikeep-1)*xdims, x0, 1, xdims)
        copysync!(x0, x)
        f0 = f1
    end
    return f0
end

function dhc1{T}(f::Function, x::BaseVector{T}; nkeep=5, grow=2.0, rmin=1e-5, rscale=0.05, ftol=1e-6, smin=0.05, smax=0.45, salpha=0.1)
    # Keep the last few points and function values:
    xdims = length(x)
    xkeep = zeros(T, xdims, nkeep)
    fkeep = zeros(T, 1, nkeep)

    # Initialize the buffer:
    r = similar(x)
    x0 = copy(x)
    f0 = f(x); nf = np = 1
    for i = 1:nkeep
        f0 = try_step(f, x, randn!(r), rscale, x0, f0, xkeep, fkeep, i); nf+=1
    end

    fdiff = similar(fkeep)
    xdiff = similar(xkeep)
    fgrad = similar(x)
    savg = 0

    while f0 > ftol
        # Estimate the gradient at x0; the / operator solves linear regression:
        broadcast!(-, xdiff, xkeep, to_host(x0))
        fdiff = fkeep - f0
        copy!(fgrad, xdiff / fdiff)

        # It turns out the ideal learning rate is fairly stable around 1.0
        # fscale = rscale * sqrt(xdims) / vecnorm(fgrad)
        fscale = 1.0

        # Try a step opposite the gradient direction
        f0 = try_step(f, x, fgrad, -fscale, x0, f0, xkeep, fkeep, rand(1:nkeep)); nf+=1
        # Try a step in a random direction
        f0 = try_step(f, x, randn!(r), rscale, x0, f0, xkeep, fkeep, rand(1:nkeep)); nf+=1

        # Report
        if (nf+=1) >= np
            println((nf,:f0,f0,:rscale,rscale,:savg,savg,:x1,x0[1],:err,wtest(x0, zeroone)))
            np *= 2
        end

        # Rescale (based on dhc0, TODO: make this adaptive)
        rscale = (nf < 32 ? 0.05 : nf < 512 ? 0.025 : nf < 2048 ? 0.0125 : 0.00625)
    end
    println((nf,:f0,f0,:rscale,rscale,:savg,savg,:x1,x0[1],:err,wtest(x0, zeroone)))
    copy!(x,x0)
    return f0
end

function gd0{T}(f::Function, x::BaseVector{T}, g::BaseVector{T}=similar(x); gscale=1.5, ftol=0.05, alpha=0.1)
    f0 = f(x, g)
    nf = np = 1
    savg = 1.0
    gavg = vecnorm(g)
    while f0 > ftol
        f1 = f(axpy!(-gscale, g, x), g)
        savg = (1-alpha)*savg + alpha * (f1<f0)
        gavg = (1-alpha)*gavg + alpha * vecnorm(g)
        f0 = f1
        if (nf+=1) >= np
            println((nf,:f0,f0,:gscale,gscale,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone)))
            np *= 2
        end
        if nf % 1000 == 0
            gc()
            # println((:mem, Knet.gpumem(), :gc, (gc();gpusync();Knet.gpumem())))
        end
    end
    println((nf,:f0,f0,:gscale,gscale,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone)))
    return f0
end

function gd1{T}(f::Function, fgrad::Function, x::BaseVector{T}, g::BaseVector{T}=similar(x); grow=sqrt(2.0), gmin=1e-5, gscale=1.0, ftol=0.05, lravg=gscale, galpha=0.1)
    x0 = copy(x)
    f0 = fgrad(x, g)
    nf = np = 1
    gnorm = vecnorm(g)
    while f0 > ftol
        f1 = f(axpy!(-gscale, g, copy!(x,x0)))
        while f1 < f0
            gscale *= grow
            f1 = f(axpy!(-gscale, g, copy!(x,x0)))
        end
        while f1 > f0
            gscale /= grow
            f1 = f(axpy!(-gscale, g, copy!(x,x0)))
        end
        copy!(x0, x)
        f0 = fgrad(x0, g)
        @assert f0==f1
        lravg = galpha * gscale + (1-galpha) * lravg
        gnorm = galpha * vecnorm(g) + (1-galpha) * gnorm
        if (nf+=1) >= np
            println((nf,:f0,f0,:lravg,lravg,:gnorm,gnorm,:err,wtest(x0, zeroone)))
            np *= 2
        end
        if nf % 1000 == 0
            gc() # println((:mem, Knet.gpumem(), :gc, (gc();gpusync();Knet.gpumem())))
        end
    end
    println((nf,:f0,f0,:lravg,lravg,:gnorm,gnorm,:err,wtest(x0, zeroone)))
    copy!(x,x0)
    return f0
end

# Test gradient estimation using SGD updates
# using CUBLAS: dot

function gd2{T}(f::Function, x::BaseVector{T}, g::BaseVector{T}=similar(x); gscale=1.5, ftol=0.05, alpha=0.1, lr=0.0001, l2=0)
    f0 = f(x, g)
    nf = np = 1
    savg = 1.0
    gavg = vecnorm(g)

    # sgd update: g := (1-lr*h) g - lr * (g.w-f) * w

    # f(w) = g.w  # assuming loss is linear in w
    # df/dw = g   # g is also the gradient
    # w' = w - lr*g # therefore would be used in updating w in SGD

    # Finding g with online linear regression (sgd)
    # J = (1/2) (g.w-f)^2 + (h/2) |g|^2  # regularized loss from predicting f(w)=g.w
    # dJ/dg = (g.w - f) * w + h * g
    # g' = g - lr * ((g.w - f) * w + h * g)
    # g' = (1-lr*h) g - lr * (g.w-f) * w

    ghat = fill!(similar(g),0)
    axpy!(lr * f0, x, ghat)
    n0_avg = vecnorm(g)
    n1_avg = vecnorm(ghat)
    cs_avg = CUBLAS.dot(g,ghat)/(n0_avg*n1_avg)

    while f0 > ftol
        f1 = f(axpy!(-gscale, g, x), g)
        savg = (1-alpha)*savg + alpha * (f1<f0)
        gavg = (1-alpha)*gavg + alpha * vecnorm(g)
        f0 = f1

        f2 = CUBLAS.dot(ghat, x)
        scale!(1-lr*l2, ghat)
        axpy!(-lr*(f2-f0), x, ghat)

        n0 = vecnorm(g)
        n1 = vecnorm(ghat)
        cs = CUBLAS.dot(g,ghat)/(n0*n1)
        n0_avg = (1-alpha)*n0_avg + alpha * n0
        n1_avg = (1-alpha)*n1_avg + alpha * n1
        cs_avg = (1-alpha)*cs_avg + alpha * cs

        if (nf+=1) >= np
            println((nf,:lr,lr,:l2,l2,:cs,cs_avg,:n0,n0_avg,:n1,n1_avg,:f0,f0,:gscale,gscale,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone)))
            np *= 2
        end
        if nf % 1000 == 0
            gc()
            # println((:mem, Knet.gpumem(), :gc, (gc();gpusync();Knet.gpumem())))
        end
    end
    println((nf,:lr,lr,:l2,l2,:cs,cs_avg,:n0,n0_avg,:n1,n1_avg,:f0,f0,:gscale,gscale,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone)))
    return f0
end

# gd2 does not handle the biases right.  It is using g.x=f.  It should
# be using g.x+b=f with b approximately equal to the value of the
# first point.  Or instead we can use (g.x1+b)-(g.x2+b)=(f1-f2) or
# g.(x1-x2)=(f1-f2) for updates and not worry about biases.  This is gd3.

function gd3{T}(f::Function, x::BaseVector{T}; gscale=1.5, ftol=0.05, alpha=0.01, lr=1.0, l2=0.1)
    g0 = similar(x)
    f0 = f(x, g0)

    nf = np = 1
    savg = 1.0
    gavg = vecnorm(g0)

    g1 = fill!(similar(g0),0)
    n0_avg = vecnorm(g0)
    n1_avg = 0
    cs_avg = 1

    while f0 > ftol
        f1 = f(axpy!(-gscale, g0, x), g0)
        savg = (1-alpha)*savg + alpha * (f1<f0)
        gavg = (1-alpha)*gavg + alpha * vecnorm(g0)

        # Try to estimate the gradient and report the cosine.  Note that we are moving using the real gradient.

        # dx = x - x0 = -gscale*g0
        # f1hat = f0 + g1.(x-x0)
        f1hat = f0 - gscale * CUBLAS.dot(g1, g0)
        # dg = (f1hat-f1)*(x-x0) = (f1hat-f1)*(-gscale)*g0
        # update g = g - lr * dg = g + lr * (f1hat-f1) * gscale * g0
        scale!(1-lr*l2, g1)
        axpy!(lr*(f1hat-f1)*gscale, g0, g1)

        f0 = f1

        n0 = vecnorm(g0)
        n1 = vecnorm(g1)
        cs = CUBLAS.dot(g0,g1)/(n0*n1)
        n0_avg = (1-alpha)*n0_avg + alpha * n0
        n1_avg = (1-alpha)*n1_avg + alpha * n1
        cs_avg = (1-alpha)*cs_avg + alpha * cs

        if (nf+=1) >= np
            println((nf,:lr,lr,:l2,l2,:cs,cs_avg,:n0,n0_avg,:n1,n1_avg,:f0,f0,:gscale,gscale,:savg,savg,:gavg,gavg)) # ,:err,wtest(x, zeroone)))
            np *= 2
        end
        if nf % 1000 == 0
            gc()
            # println((:mem, Knet.gpumem(), :gc, (gc();gpusync();Knet.gpumem())))
        end
    end
    println((nf,:lr,lr,:l2,l2,:cs,cs_avg,:n0,n0_avg,:n1,n1_avg,:f0,f0,:gscale,gscale,:savg,savg,:gavg,gavg)) # ,:err,wtest(x, zeroone)))
    return f0
end

# Run gd, report how much the gradient turns
function gd3turn{T}(f::Function, x::BaseVector{T}; gscale=1.5, ftol=0.05, alpha=0.1)
    g0 = similar(x)
    g1 = similar(x)
    f0 = f(x, g0)
    n0 = gavg = vecnorm(g0)
    cavg = savg = NaN
    nf = np = 1
    while f0 > ftol
        f1 = f(axpy!(-gscale, g0, x), g1)
        n1 = vecnorm(g1)
        gavg = (1-alpha)*gavg + alpha * n1
        savg = isnan(savg) ? 1.0*(f1<f0) : (1-alpha)*savg + alpha * (f1<f0)
        g0g1 = CUBLAS.dot(g0,g1)/(n0*n1)
        cavg = isnan(cavg) ? g0g1 : (1-alpha)*cavg + alpha * g0g1
        (nf+=1) >= np && (np*=2; println((nf,:f0,f0,:cavg,cavg,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone))))
        nf % 100 == 0 && gc()
        copy!(g0,g1); n0,f0 = n1,f1
    end
    println((nf,:f0,f0,:cavg,cavg,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone)))
    return f0
end

using Knet: cpucopy

# Run gd, save points every 2^n iterations
function gd3save{T}(f::Function, x::BaseVector{T}; gscale=0.75, ftol=0.05, alpha=0.1)
    g0 = similar(x)
    f0 = f(x, g0)
    gavg = vecnorm(g0)
    savg = 1.0
    xbuf = Any[]
    nf = np = 1
    while f0 > ftol
        f1 = f(axpy!(-gscale, g0, x), g0)
        gavg = (1-alpha)*gavg + alpha * vecnorm(g0)
        savg = (1-alpha)*savg + alpha * (f1<f0)
        f0 = f1
        (nf+=1) >= np && (np*=2; println((nf,:f0,f0,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone))); push!(xbuf,copy(x)))
        nf % 100 == 0 && gc()
    end
    println((nf,:f0,f0,:savg,savg,:gavg,gavg,:err,wtest(x, zeroone))); push!(xbuf,copy(x))
    return xbuf
end

# gradient descent without all the bells and whistles
function gd3simple{T}(f::Function, x0::BaseVector{T}; gscale=1.0, ftol=0.05, maxnf=typemax(Int), dt=10)
    g0 = similar(x0)
    f0 = f(x0, g0)
    nf = np = nt = 1
    while f0 > ftol
        f0 = f(axpy!(-gscale, g0, x0), g0)
        if (nf+=1) >= np || time() >= nt
            println((nf,:f0,f0,:g0,vecnorm(g0),:x0,vecnorm(x0)))
            nf >= np && (np *= 2)
            time() >= nt && (nt = time()+dt)
            nf >= maxnf && break
            gc()
        end
    end
    println((nf,:f0,f0,:g0,vecnorm(g0),:x0,vecnorm(x0)))
    return f0
end

# gradient descent with adaptive step size
function gd3ada{T}(f::Function, x0::BaseVector{T}; ftol=0.0, maxnf=typemax(Int), dt=10, gmax=1.0)
    g0,x = [ similar(x0) for i=1:2 ]
    @show f0 = f(x0, g0)
    nf = np = nt = nu = 1
    nx,ng = vecnorm(x0),vecnorm(g0)
    gscale = lr = gmax/ng
    while f0 > ftol
        gscale = min(lr, gmax/ng)
        f0 = f(axpy!(-gscale, g0, x0), g0); nf += 1
        nx,ng = vecnorm(x0),vecnorm(g0)
        if nf >= nu
            f1 = f(axpy!(-lr,  g0, copy!(x,x0)))
            lr2 = lr * exp(rand()-1/2)
            f2 = f(axpy!(-lr2, g0, copy!(x,x0)))
            f2 < f1 && (lr = lr2)
            nu = min(nu*1.2,nu+50)
        end
        if nf >= np || time() >= nt
            println((nf,:f0,f0,:g0,ng,:x0,nx,:lr,lr,:gs,gscale))
            gc()
            nf >= maxnf && break
            nt = time()+dt
            nf >= np && (np *= 2)
        end
    end
    println((nf,:f0,f0,:g0,ng,:x0,nx,:lr,lr,:gs,gscale))
    return f0
end

function ego2{T}(f::Function, x::BaseVector{T}; ftol=0.05, alpha=0.01, lr=1.0, l2=0.0, noise=1.0, rmin=1e-6, stepsize=100.0, maxnf=Inf)
    # x0 is the anchor point
    x0 = copy(x)
    g0 = similar(x0)
    f0 = f(x0, g0)
    n0 = vecnorm(g0)
    xdims = length(x0)

    # g1 is the gradient estimate, starts at 0
    g1 = fill!(similar(g0),0)
    n1 = 0

    # arrays for actual gradient (for debugging) and step
    g = similar(x)
    r = similar(x)

    # some statistics
    nf = np = 1
    savg = 0.5
    n0_avg = n0
    n1_avg = n1
    nr_avg = rmin*sqrt(xdims)
    cs_avg = 0.0

    while f0 > ftol
        # rmin is the minimum scale of the random vector components
        # noise is the norm ratio of random vector to gradient vector
        rscale = max(rmin, noise*n1/sqrt(xdims))
        scale!(rscale, randn!(r))

        # add noise to the negative gradient to get the step vector
        axpy!(-1, g1, r)
        # take step scaled with stepsize (default 1), updating f1 and g
        scale!(stepsize, r)
        nr = vecnorm(r)
        isfinite(nr)  || error("!isfinite(r)")
        f1 = f(axpy!(1, r, copy!(x,x0)), g)
        # compute success rate
        savg = (1-alpha)*savg + alpha * (f1<f0)

        # dx = x - x0 = r
        # f2 = f0 + g1.(x-x0) is our estimate
        f2 = f0 + CUBLAS.dot(g1, r)
        # weight decay:
        l2!=0 && scale!(1-lr*l2, g1)
        # dg = (f2-f1)*(x-x0) = (f2-f1)*r
        # nlms update g1 = g1 - (lr/nr^2) * dg = g1 - (lr/nr^2) * (f2-f1) * r
        axpy!(-lr*(f2-f1)/(nr*nr), r, g1)
        n1 = vecnorm(g1)
        isfinite(n1) || error("!isfinite(n1)")

        if f1 < f0
            copy!(x0, x)
            copy!(g0, g)
            f0 = f1
            n0 = vecnorm(g0)
        end

        # Report
        cs = n1==0 ? 0 : CUBLAS.dot(g0,g1)/(n0*n1)
        isfinite(cs) || error("!isfinite(cs)")
        n0_avg = (1-alpha)*n0_avg + alpha * n0
        n1_avg = (1-alpha)*n1_avg + alpha * n1
        nr_avg = (1-alpha)*nr_avg + alpha * nr
        cs_avg = (1-alpha)*cs_avg + alpha * cs
        if (nf+=1) >= np
            # println((i,:cos,cs,:n1n0,n1/n0,:f2f1,f2-f1,:f1f0,f1-f0,:f0,f0,:n0,n0,:n1,n1,:nr,nr))
            println((nf,:f0,f0,:cos,cs_avg,:n0,n0_avg,:n1,n1_avg,:nr,nr_avg,:savg,savg,:err,wtest(x, zeroone),:lr,lr,:l2,l2,:stepsize,stepsize,:noise,noise))
            np *= 2
        end
        if nf % 1000 == 0
            gc()
            # println((:mem, Knet.gpumem(), :gc, (gc();gpusync();Knet.gpumem())))
        end
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:cos,cs_avg,:n0,n0_avg,:n1,n1_avg,:nr,nr_avg,:savg,savg,:err,wtest(x, zeroone),:lr,lr,:l2,l2,:stepsize,stepsize,:noise,noise))
    return f0
end

# Consider the points x0, x0+r, x0+g, x0+g+r
# where g is the grad estimate and r is a similar sized random vector
# Evaluate using them for steps and for gradient estimation

function ego2b{T}(f::Function, x::BaseVector{T}; ftol=0.05, alpha=0.1, lr=1.0, l2=0.0, noise=1.0, rmin=1e-6, step=10.0, maxnf=Inf, update=1)
    dg = similar(x)             # -step * g
    dr = similar(x)             # noise * step * (ng/nr) * randn
    dh = similar(x)             # dg + dr
    dz = similar(x)
    xg = similar(x)             # x + dg
    xr = similar(x)             # x + dr
    xh = similar(x)             # x + dh
    xz = similar(x)
    g0 = similar(x)             # real gradient from x
    fx = f(x,g0)
    ng0 = vecnorm(g0)
    g = scale!(rmin, randn!(similar(x))) # gradient estimate
    ng = vecnorm(g)
    nr = sqrt(length(x))        # norm of randn
    dfg = dfr = dfh = 0.0
    nf = nt = cs = 0; np = 1

    while fx > ftol
        scale!(-step, copy!(dg,g))
        axpy!(1, dg, copy!(xg,x))
        fg = f(xg)

        scale!(noise * step * ng / nr, randn!(dr))
        axpy!(1, dr, copy!(xr,x))
        fr = f(xr)

        scale!(0.0001, randn!(dz))
        axpy!(1, dz, copy!(xz,x))
        fz = f(xz)

        axpy!(1, dr, copy!(dh,dg))
        axpy!(1, dh, copy!(xh,x))
        fh = f(xh)

        # Which point is the best to move to?
        dfg = (1-alpha)*dfg + alpha*(fg<fx)*(fg-fx)
        dfr = (1-alpha)*dfr + alpha*(fr<fx)*(fr-fx)
        dfh = (1-alpha)*dfh + alpha*(fh<fx)*(fh-fx)

        if update==1
            fh2 = fx + CUBLAS.dot(g, dh)
            axpy!(-lr*(fh2-fh)/(vecnorm(dh)^2), dh, g)
        elseif update==2
            fr2 = fx + CUBLAS.dot(g, dr)
            axpy!(-lr*(fr2-fr)/(vecnorm(dr)^2), dr, g)
        elseif update==3
            fz2 = fx + CUBLAS.dot(g, dz)
            axpy!(-lr*(fz2-fz)/(vecnorm(dz)^2), dz, g)
        end
        ng = vecnorm(g)

        # Report
        nf += 1
        if (nf>=np) # (time()>=nt)
            cs = CUBLAS.dot(g,g0)/(ng*ng0)
            println((nf,:fx,fx,:ng,ng,:ng0,ng0,"ng/ng0",ng/ng0,:cos,cs,"sqrt(nf/nd)",sqrt(nf/length(x)),:dfg,dfg,:dfr,dfr,:dfh,dfh))
            nt=10+time(); np*=2
        end

        # Update anchor point based on xg
        if fg < fx
            copy!(x, xg)
            fx = f(x, g0)
            ng0 = vecnorm(g0)
        end
        nf >= maxnf && break
        nf % 100 == 0 && gc()
    end
    println((nf,:fx,fx,:ng,ng,:ng0,ng0,:cos,cs,:dfg,dfg,:dfr,dfr,:dfh,dfh))
end

# Make the step size (maybe later also the noise size) adaptive.

function ego3{T}(f::Function, x::BaseVector{T}; ftol=0.05, alpha=0.01, lr=0.5, l2=0.0, noise=2.0, rmin=1e-6, stepsize=128.0, smin=0.2, smax=0.8, salpha=0.1, grow=sqrt(2.0), maxnf=Inf)
    # x0 is the anchor point
    x0 = copy(x)
    g0 = similar(x0)
    f0 = f(x0, g0)
    n0 = vecnorm(g0)
    xdims = length(x0)

    # g1 is the gradient estimate, starts at 0
    g1 = fill!(similar(g0),0)
    n1 = 0

    # arrays for actual gradient (for debugging) and step
    # TODO: we should measure the angle to target rather than the (ever changing) local gradient.
    g = similar(x)
    r = similar(x)

    # some statistics
    nf = 1; t0 = 0
    savg = 0.5
    n0_avg = n0
    n1_avg = n1
    nr_avg = rmin*sqrt(xdims)
    cs_avg = 0.0

    while f0 > ftol
        # rmin is the minimum scale of the noise vector components
        # noise is the norm ratio of noise vector to gradient vector
        rscale = max(rmin, noise*n1/sqrt(xdims))
        scale!(rscale, randn!(r)) # this makes vecnorm(r) ~ noise*vecnorm(g1)
        # add noise to the negative gradient to get the step vector and scale with stepsize
        axpy!(-1, g1, r)
        scale!(stepsize, r)
        nr = vecnorm(r)

        # take step scaled with stepsize, updating f1 and g
        f1 = f(axpy!(1, r, copy!(x,x0)), g)
        nf += 1

        # compute success rate and update stepsize accordingly
        savg = (1-salpha)*savg + salpha * (f1<f0)
        savg < smin && (stepsize /= grow; savg = 0.5)
        savg > smax && (stepsize *= grow; savg = 0.5)

        # update gradient estimate
        # f2 = f0 + g1.(x-x0) is our fval estimate
        # note that r = x - x0
        f2 = f0 + CUBLAS.dot(g1, r)
        # weight decay:
        l2!=0 && scale!(1-lr*l2, g1)
        # J = (f2-f1)^2/2
        # dg = dJ/dg1 = (f2-f1)*(x-x0) = (f2-f1)*r
        # nlms update g1 = g1 - lr * dg * (1/nr^2) = g1 - lr * (f2-f1) * r * (1/nr^2)
        axpy!(-lr*(f2-f1)/(nr*nr), r, g1)
        n1 = vecnorm(g1)

        if f1 < f0
            copy!(x0, x)
            copy!(g0, g)
            f0 = f1
            n0 = vecnorm(g0)
        end

        # Report
        cs = (n1==0) ? 0 : CUBLAS.dot(g0,g1)/(n0*n1)
        n0_avg = (1-alpha)*n0_avg + alpha * n0
        n1_avg = (1-alpha)*n1_avg + alpha * n1
        nr_avg = (1-alpha)*nr_avg + alpha * nr
        cs_avg = (1-alpha)*cs_avg + alpha * cs
        if time() >= t0
            println((nf,:f0,f0,:cos,cs_avg,:n0,n0_avg,:n1,n1_avg,:nr,nr_avg,:savg,savg,:err,wtest(x, zeroone),:lr,lr,:l2,l2,:stepsize,stepsize,:noise,noise))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:cos,cs_avg,:n0,n0_avg,:n1,n1_avg,:nr,nr_avg,:savg,savg,:err,wtest(x, zeroone),:lr,lr,:l2,l2,:stepsize,stepsize,:noise,noise))
    return f0
end

urand(min,max)=rand()*(max-min)+min
lrand(min,max)=exp(urand(log(min),log(max)))
function test_ego3()
    copy!(w64, scale!(0.01, convert(Array{Float32}, randn(MersenneTwister(42), size(w64)))))
    lr=lrand(1.2/2, 1.2*2)
    l2=(rand()<0.5 ? 0 : lrand(0.001/lr,1.0/lr))
    noise=lrand(1.6/2, 1.6*2)
    stepsize=lrand(100/2, 100*2)
    @show (lr,l2,noise,stepsize)
    ego3(weval,w64; maxnf=1000, lr=lr, l2=l2, noise=noise, stepsize=stepsize)
end

# Separate gradient steps and noise steps

function ego4{T}(f::Function, x::BaseVector{T}; ftol=0.05, alpha=0.01, lr=0.5, l2=0.0, noise=2.0, rmin=1e-3, stepsize=128.0, smin=0.2, smax=0.8, salpha=0.1, grow=sqrt(2.0), maxnf=Inf, rstep=2)
    # x0 is the anchor point
    x0 = copy(x)
    g0 = similar(x0)
    f0 = f(x0, g0)
    n0 = vecnorm(g0)
    xdims = length(x0)

    # g1 is the gradient estimate, starts at 0
    g1 = fill!(similar(g0),0)
    n1 = 0

    # arrays for actual gradient (for debugging) and step
    # TODO: we should measure the angle to target rather than the (ever changing) local gradient.
    g = similar(x)
    r = similar(x)

    # some statistics
    nf = 1; t0 = 0
    savg = 0.5
    n0_avg = n0
    n1_avg = n1
    nr_avg = rmin*sqrt(xdims)
    cs_avg = 0.0

    while f0 > ftol
        if n1==0 || nf%rstep==0
            scale!(rmin, randn!(r))
        else
            axpy!(-stepsize, g1, r)
        end
        nr = vecnorm(r)

        # take step scaled with stepsize, updating f1 and g
        f1 = f(axpy!(1, r, copy!(x,x0)), g)
        nf += 1

        # compute success rate and update stepsize accordingly
        savg = (1-salpha)*savg + salpha * (f1<f0)
        savg < smin && (stepsize /= grow; savg = 0.5)
        savg > smax && (stepsize *= grow; savg = 0.5)

        # update gradient estimate
        # f2 = f0 + g1.(x-x0) is our fval estimate
        # note that r = x - x0
        f2 = f0 + CUBLAS.dot(g1, r)
        # weight decay:
        l2!=0 && scale!(1-lr*l2, g1)
        # J = (f2-f1)^2/2
        # dg = dJ/dg1 = (f2-f1)*(x-x0) = (f2-f1)*r
        # nlms update g1 = g1 - lr * dg * (1/nr^2) = g1 - lr * (f2-f1) * r * (1/nr^2)
        axpy!(-lr*(f2-f1)/(nr*nr), r, g1)
        n1 = vecnorm(g1)

        if f1 < f0
            copy!(x0, x)
            copy!(g0, g)
            f0 = f1
            n0 = vecnorm(g0)
        end

        # Report
        cs = (n1==0) ? 0 : CUBLAS.dot(g0,g1)/(n0*n1)
        n0_avg = (1-alpha)*n0_avg + alpha * n0
        n1_avg = (1-alpha)*n1_avg + alpha * n1
        nr_avg = (1-alpha)*nr_avg + alpha * nr
        cs_avg = (1-alpha)*cs_avg + alpha * cs
        if time() >= t0
            println((nf,:f0,f0,:cos,cs_avg,:n0,n0_avg,:n1,n1_avg,:nr,nr_avg,:savg,savg,:err,wtest(x, zeroone),:lr,lr,:l2,l2,:stepsize,stepsize,:noise,noise))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:cos,cs_avg,:n0,n0_avg,:n1,n1_avg,:nr,nr_avg,:savg,savg,:err,wtest(x, zeroone),:lr,lr,:l2,l2,:stepsize,stepsize,:noise,noise))
    return f0
end

# Make stepsize and noise size really adaptive by experimenting and
# comparing their impact on f decrease. We try two moves from a given
# point.  The first with the best step/noise.  The second with a
# mutated step/noise.  If the second one wins we update best
# step/noise.  We use both data points to update grad estimate.

function ego5{T}(f::Function, x::BaseVector{T}; ftol=0.05, rmin=1e-6, lr=0.5, noise=2.0, step=128.0, maxnf=Inf)
    # gradient estimate
    g = fill!(similar(x), 0)

    # x0 is the best point known so far
    x0 = copy(x)
    f0 = f(x); nf = 1

    # if x2 is better than x1 we update step and noise
    # if x1 or x2 is better than x0 we update x0
    # in all cases we update the gradient estimate
    x1 = similar(x)
    x2 = similar(x)
    dx1 = similar(x)
    dx2 = similar(x)
    t0 = 0

    while f0 > ftol

        # x1 = x0 + dx1; dx1 = -step *  (g + noise  * |g| * r/|r|)  where r is a random vector
        rscale = max(rmin, vecnorm(g)/sqrt(length(x)))
        axpy!(-step, g, scale!(step * noise * rscale, randn!(dx1)))
        f1 = f(axpy!(1, dx1, copy!(x1,x0))); nf += 1

        # x2 = x0 - dx2; dx2 = step2 * (g + noise2 * |g| * r/|r|)
        step2 = step * exp(0.2*(rand()-0.5))
        noise2 = noise * exp(0.2*(rand()-0.5))
        axpy!(-step2, g, scale!(step2 * noise2 * rscale, randn!(dx2)))
        f2 = f(axpy!(1, dx2, copy!(x2,x0))); nf += 1

        # update gradient estimate:
        # h = f0 + g.dx
        # J = (h - f)^2/2
        # dg = (h - f) * dx
        # g = g - (lr/|dx|^2) * (h - f) * dx

        h1 = f0 + CUBLAS.dot(g, dx1)
        h2 = f0 + CUBLAS.dot(g, dx2)
        axpy!(-lr*(h1-f1)/vecnorm(dx1)^2, dx1, g)
        axpy!(-lr*(h2-f2)/vecnorm(dx2)^2, dx2, g)

        if f2 < f0 && f2 < f1
            copy!(x0, x2)
            f0 = f2
        end
        if f1 < f0 && f1 < f2
            copy!(x0, x1)
            f0 = f1
        end
        if f2 < f1
            step = step2
            # noise getting too low gets the algorithm stuck.
            noise2 > 1 && (noise = noise2)
        end
        if time() >= t0
            println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
    return f0
end

# ego5b: Update f2 (thus step,noise) every 10 steps instead of every step.

function ego5b{T}(f::Function, x::BaseVector{T}; ftol=0.05, rmin=1e-6, lr=0.5, noise=2.0, step=128.0, maxnf=Inf)
    # gradient estimate
    g = fill!(similar(x), 0)

    # x0 is the best point known so far
    x0 = copy(x)
    f0 = f(x); nf = 1

    # if x2 is better than x1 we update step and noise
    # if x1 or x2 is better than x0 we update x0
    # in all cases we update the gradient estimate
    x1 = similar(x)
    x2 = similar(x)
    dx1 = similar(x)
    dx2 = similar(x)
    t0 = 0

    while f0 > ftol

        # x1 = x0 + dx1; dx1 = -step *  (g + noise  * |g| * r/|r|)  where r is a random vector
        rscale = max(rmin, vecnorm(g)/sqrt(length(x)))
        axpy!(-step, g, scale!(step * noise * rscale, randn!(dx1)))
        f1 = f(axpy!(1, dx1, copy!(x1,x0))); nf += 1

        if nf % 10 == 0
            # x2 = x0 - dx2; dx2 = step2 * (g + noise2 * |g| * r/|r|)
            step2 = step * exp(0.2*(rand()-0.5))
            noise2 = noise * exp(0.2*(rand()-0.5))
            axpy!(-step2, g, scale!(step2 * noise2 * rscale, randn!(dx2)))
            f2 = f(axpy!(1, dx2, copy!(x2,x0))); nf += 1
        else
            f2 = Inf
        end

        # update gradient estimate:
        # h = f0 + g.dx
        # J = (h - f)^2/2
        # dg = (h - f) * dx
        # g = g - (lr/|dx|^2) * (h - f) * dx

        h1 = f0 + CUBLAS.dot(g, dx1)
        nf%10==0 && (h2 = f0 + CUBLAS.dot(g, dx2))
        axpy!(-lr*(h1-f1)/vecnorm(dx1)^2, dx1, g)
        nf%10==0 && (axpy!(-lr*(h2-f2)/vecnorm(dx2)^2, dx2, g))

        if f2 < f0 && f2 < f1
            copy!(x0, x2)
            f0 = f2
        end
        if f1 < f0 && f1 < f2
            copy!(x0, x1)
            f0 = f1
        end
        if f2 < f1
            step = step2
            # noise getting too low gets the algorithm stuck.
            noise2 > 1 && (noise = noise2)
        end
        if time() >= t0
            println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
    return f0
end

# Try separating the grad move and noise move again.

function ego6{T}(f::Function, x::BaseVector{T}; ftol=0.05, rmin=1e-6, lr=0.5, step=128.0, maxnf=Inf, salpha=0.1, smin=0.2, smax=0.8, grow=sqrt(2.0))
    # gradient estimate
    g = fill!(similar(x), 0)

    # x0 is the best point known so far
    x0 = copy(x)
    f0 = f(x); nf = 1

    # if x2 is better than x1 we update step and noise
    # if x1 or x2 is better than x0 we update x0
    # in all cases we update the gradient estimate
    x1 = similar(x)
    x2 = similar(x)
    dx = similar(x)
    r = similar(x)
    savg = 0.5
    t0 = 0

    while f0 > ftol

        # x1 = x0 - step * g
        f1 = f(axpy!(-step, g, copy!(x1,x0))); nf += 1

        # x2 = x1 + noise
        rscale = max(rmin, step*vecnorm(g)/sqrt(length(x)))
        f2 = f(axpy!(rscale, randn!(r), copy!(x2,x1))); nf += 1

        # update gradient estimate:
        # h = f0 + g.dx
        # J = (h - f)^2/2
        # dg = (h - f) * dx
        # g = g - (lr/|dx|^2) * (h - f) * dx

        axpy!(-1, x0, copy!(dx,x2))
        h2 = f0 + CUBLAS.dot(g, dx)
        axpy!(-lr*(h2-f2)/vecnorm(dx)^2, dx, g)

        # compute success rate and update stepsize accordingly
        savg = (1-salpha)*savg + salpha * (f1<f0)
        savg < smin && (step /= grow; savg = 0.5)
        savg > smax && (step *= grow; savg = 0.5)

        if f1 < f0 && f1 < f2
            f0 = f1
            copy!(x0, x1)
        elseif f2 < f0 && f2 < f1
            f0 = f2
            copy!(x0, x2)
        end

        if time() >= t0
            println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:lr,lr,:err,wtest(x0, zeroone)))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:lr,lr,:err,wtest(x0, zeroone)))
    return f0
end

using Knet: cslice!, csize

function wmini(ws...; x=gpu_xtrn, y=gpu_ytrn, loss=softloss, batchsize=100)
    global xmini, ymini
    if !isdefined(:xmini) || size(xmini)!=csize(x,batchsize)
        xmini = similar(x, csize(x, batchsize))
        ymini = similar(y, csize(y, batchsize))
    end
    nd = size(x, ndims(x))
    id = rand(1:nd,batchsize)
    cslice!(xmini, x, id)
    cslice!(ymini, y, id)
    fval = map(ws) do w
        f = initmodel(w)
        ypred = forw(f, xmini)
        loss(ypred, ymini)
    end
    length(fval)==1 && (fval=fval[1])
    return fval
end

# ego7: Trying minibatching.  Each minibatch will change f.  In order to do
# meaningful comparisons we need to compare multiple points on the
# same minibatch.  So we change the calling convention of f to take an
# arbitrary number of w's and return a corresponding number of values.

function ego7{T}(f::Function, x0::BaseVector{T}; ftol=0.05, rmin=1e-6, lr=1.0, noise=1.2, step=50.0, maxnf=Inf, maxdx=0.1)
    # gradient estimate
    g = fill!(similar(x0), 0)

    # if x1 or x2 is better than x0 we update x0
    # if x2 is better than x1 we update step and noise
    # in all cases we update the gradient estimate
    x1 = similar(x0)
    x2 = similar(x0)
    dx1 = similar(x0)
    dx2 = similar(x0)
    nf = t0 = 0; f0 = Inf

    while f0 > ftol

        # x1 = x0 + dx1; dx1 = -step *  (g + noise  * |g| * r/|r|)  where r is a random vector
        rscale = max(rmin, vecnorm(g)/sqrt(length(x0)))
        axpy!(-step, g, scale!(step * noise * rscale, randn!(dx1)))
        n1 = vecnorm(dx1); n1 > maxdx && scale!(maxdx/n1, dx1)
        axpy!(1, dx1, copy!(x1,x0))

        # x2 = x0 - dx2; dx2 = step2 * (g + noise2 * |g| * r/|r|)
        step2 = step * exp(0.2*(rand()-0.5))
        noise2 = noise * exp(0.2*(rand()-0.5))
        axpy!(-step2, g, scale!(step2 * noise2 * rscale, randn!(dx2)))
        n2 = vecnorm(dx2); n2 > maxdx && scale!(maxdx/n2, dx2)
        axpy!(1, dx2, copy!(x2,x0))

        (f0,f1,f2) = f(x0,x1,x2); nf+=3

        # update gradient estimate:
        # h = f0 + g.dx
        # J = (h - f)^2/2
        # dg = (h - f) * dx
        # g = g - (lr/|dx|^2) * (h - f) * dx

        h1 = f0 + CUBLAS.dot(g, dx1)
        h2 = f0 + CUBLAS.dot(g, dx2)
        axpy!(-lr*(h1-f1)/vecnorm(dx1)^2, dx1, g)
        axpy!(-lr*(h2-f2)/vecnorm(dx2)^2, dx2, g)

        if f2 < f0 && f2 < f1
            copy!(x0, x2)
            f0 = f2
        end
        if f1 < f0 && f1 < f2
            copy!(x0, x1)
            f0 = f1
        end
        if f2 < f1
            # step = step2
            # # noise getting too low gets the algorithm stuck.
            # noise2 > 1 && (noise = noise2)
        end
        if time() >= t0
            println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
    return f0
end

# Let us simplify ego7, separating the gradient step

function ego8{T}(f::Function, x0::BaseVector{T}; ftol=0.05, rmin=1e-6, lr=1.0, noise=1.2, step=50.0, maxnf=Inf, smin=0.5, smax=0.95, alpha=0.01, grow=sqrt(2))
    g = scale!(rmin, randn!(similar(x0)))
    r = similar(g)
    x1 = similar(x0)
    nf = t0 = 0; favg = f0 = f(x0)
    savg = 0.5
    while f0 > ftol
        axpy!(-step, g, copy!(x1, x0)) # x1 = x0 - step * g
        (f0,f1) = f(x0,x1); nf+=2
        favg = (1-alpha)*favg + alpha * f0
        # compute success rate and update stepsize accordingly
        savg = (1-alpha)*savg + alpha * (f1<f0)
        #savg < smin && (step /= grow; savg = 0.5)
        #savg > smax && (step *= grow; savg = 0.5)
        if f1 < f0
            copy!(x0, x1)
            f0 = f1
        else
            # axpy!(noise*step*vecnorm(g)/sqrt(length(g)), randn!(r), x1)
            axpy!(0.1/sqrt(length(g)), randn!(r), x1)
            (f0,f1) = f(x0,x1); nf+=2
            dx = axpy!(-1, x0, x1)
            f2 = f0 + CUBLAS.dot(g, dx)
            axpy!(-lr*(f2-f1)/vecnorm(dx)^2, dx, g)
        end
        if time() >= t0
            println((nf,:favg,favg,:savg,savg,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
            t0 = 10+time()
        end
        nf % 100 == 0 && gc()
        nf >= maxnf && break
    end
    println((nf,:favg,favg,:savg,savg,:gnorm,vecnorm(g),:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
    return wtest(x0, softloss)
end

# This is back to batch version with separate grad update and point update steps:
# Failed yet another time.  Cannot beat ego2 using twice as many evals.

function ego9{T}(f::Function, x0::BaseVector{T}; ftol=0.05, rstep=0.0001, gstep=100.0, lr=0.1, maxnf=typemax(Int))
    x = similar(x0)
    d = similar(x0)
    g = similar(x0)
    fill!(g, 0)
    f0 = f(x0)
    dims = length(x0)
    nf = np = 1
    while f0 > ftol
        # Update the gradient
        axpy!(rstep, randn!(d), copy!(x,x0))            # dx = rstep * d
        err = f0 + rstep * CUBLAS.dot(d,g) - f(x)       # fhat = f0 + g.dx
        axpy!(-(lr*err)/(rstep*dims), d, g)             # g -= lr*err*dx/|dx|^2; |dx|=rstep*|d|=rstep*sqrt(dims)
        # Update the point
        axpy!(-gstep, g, copy!(x,x0))
        fx = f(x)
        if fx < f0
            copy!(x0, x)
            f0 = fx
        end
        (nf+=2)>=np && (np*=2;println((nf,:f0,f0,:gnorm,vecnorm(g),:gstep,gstep,:rstep,rstep,:lr,lr,:err,wtest(x0, zeroone))))
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:gnorm,vecnorm(g),:gstep,gstep,:rstep,rstep,:lr,lr,:err,wtest(x0, zeroone)))
end

# ego9 did not work.  Here is the simplified version of ego2 for comparison:

function ego2c{T}(f::Function, x0::BaseVector{T}; ftol=0.05, noise=1.0, step=100.0, lr=1.0, gmin=0.1, maxnf=typemax(Int))
    x = similar(x0)
    d = similar(x0)
    g = fill!(similar(x0),0)
    gnorm = 0.0
    f0 = f(x0)
    rnorm = sqrt(length(x0))
    dnorm = 0; nf = np = 1
    while f0 > ftol
        axpy!(-step, g, scale!(noise*step*max(gmin,gnorm)/rnorm, randn!(d)))
        dnorm = vecnorm(d)
        fx = f(axpy!(1, d, copy!(x,x0)))                # dx = -step*g + step*noise*(gnorm/rnorm)*r
        err = f0 + CUBLAS.dot(d,g) - fx                 # fhat = f0 + g.dx
        axpy!(-(lr*err)/(dnorm^2), d, g)                # g -= lr*err*dx/|dx|^2
        gnorm = vecnorm(g)
        if fx < f0
            copy!(x0, x)
            f0 = fx
        end
        (nf+=1)>=np && (np*=2;println((nf,:f0,f0,:gnorm,gnorm,:dnorm,dnorm,:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone))))
        nf >= maxnf && break
    end
    println((nf,:f0,f0,:gnorm,gnorm,:dnorm,dnorm,:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
end

# And this is ego2 adapted to minibatches

function ego2mini{T}(f::Function, x0::BaseVector{T}; ftol=0.0, noise=1.0, step=100.0, lr=1.0, maxnf=typemax(Int))
    x = similar(x0)
    d = similar(x0)
    g = fill!(similar(x0),0)
    gnorm = 0
    f0 = Inf
    rnorm = sqrt(length(x0))
    dnorm = nf = 0; np = 1
    while f0 > ftol
        axpy!(-step, g, scale!(noise*step*max(gnorm,.0001)/rnorm, randn!(d)))
        dnorm = vecnorm(d)
        axpy!(1, d, copy!(x,x0))                        # dx = -step*g + step*noise*(gnorm/rnorm)*r
        (f0,fx) = f(x0,x)
        err = f0 + CUBLAS.dot(d,g) - fx                 # fhat = f0 + g.dx
        axpy!(-(lr*err)/(dnorm^2), d, g)                # g -= lr*err*dx/|dx|^2
        gnorm = vecnorm(g)
        if fx < f0
            copy!(x0, x)
            f0 = fx
        end
        (nf+=2)>=np && (np*=2;println((nf,:f0,wtest(x0,softloss),:gnorm,gnorm,:dnorm,dnorm,:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone))))
        nf >= maxnf && break
    end
    println((nf,:f0,wtest(x0,softloss),:gnorm,gnorm,:dnorm,dnorm,:step,step,:noise,noise,:lr,lr,:err,wtest(x0, zeroone)))
end

# minibatch version of ego9

function ego9mini{T}(f::Function, x0::BaseVector{T}; ftol=0.0, rstep=0.0001, gstep=100.0, lr=1.0, l2=0.1, maxnf=typemax(Int))
    x1 = similar(x0)
    x2 = similar(x0)
    d = similar(x0)
    g = similar(x0)
    fill!(g, 0)
    f0 = Inf
    dims = length(x0)
    nf = np = 1
    while f0 > ftol
        # Update the gradient
        axpy!(rstep, randn!(d), copy!(x1,x0))           # dx = rstep * d
        (f0,f1) = f(x0,x1)
        err = f0 + rstep * CUBLAS.dot(d,g) - f1         # fhat = f0 + g.dx
        l2!=0 && scale!(1-lr*l2, g)
        axpy!(-(lr*err)/(rstep*dims), d, g)             # g -= lr*err*dx/|dx|^2; |dx|=rstep*|d|=rstep*sqrt(dims)
        # Update the point
        axpy!(-gstep, g, copy!(x2,x0))
        (f0,f2) = f(x0,x2)
        if f2 < f0
            copy!(x0, x2)
            f0 = f2
        end
        (nf+=1)>=np && (np*=2;println((nf,:f0,wtest(x0,softloss),:gnorm,vecnorm(g),:gstep,gstep,:rstep,rstep,:lr,lr,:l2,l2,:err,wtest(x0, zeroone))))
        nf >= maxnf && break
    end
    println((nf,:f0,wtest(x0,softloss),:gnorm,vecnorm(g),:gstep,gstep,:rstep,rstep,:lr,lr,:l2,l2,:err,wtest(x0, zeroone)))
end

# Adaptive step size version of ego9mini

function ego9ada{T}(f::Function, x0::BaseVector{T}; ftol=0.0, rstep=0.0001, gstep=10.0, lr=1.0, l2=0.1, maxnf=typemax(Int))
    x1 = similar(x0)
    x2 = similar(x0)
    x3 = similar(x0)
    d = similar(x0)
    g = similar(x0)
    fill!(g, 0)
    f0 = Inf
    dims = length(x0)
    nf = nt = 0; np = 1
    for nf = 1:maxnf
        # Update the gradient
        axpy!(rstep, randn!(d), copy!(x1,x0))           # dx = rstep * d
        (f0,f1) = f(x0,x1)
        err = f0 + rstep * CUBLAS.dot(d,g) - f1         # fhat = f0 + g.dx
        l2!=0 && scale!(1-lr*l2, g)
        axpy!(-(lr*err)/(rstep*dims), d, g)             # g -= lr*err*dx/|dx|^2; |dx|=rstep*|d|=rstep*sqrt(dims)
        # Update the step size
        axpy!(-gstep, g, copy!(x2,x0))
        if nf%100 == 0 && false #DBG
            gstep2 = exp(rand()-0.5)*gstep
            axpy!(-gstep2, g, copy!(x3,x0))
            (f0,f2,f3) = f(x0,x2,x3)
            (f3 < f2) && (gstep = gstep2)
        else
            (f0,f2) = f(x0,x2)
        end
        (f2 < f0) && (f0=f2; copy!(x0,x2))
        nf >= maxnf && break
        if (nf+=1)>=np || time()>=nt
            nf >= np && (np*=2)
            time() >= nt && (nt=60+time())
            println((nf,:f0,wtest(x0,softloss),:gnorm,vecnorm(g),:gstep,gstep,:rstep,rstep,:lr,lr,:l2,l2,:err,wtest(x0, zeroone)))
        end
    end
    println((nf,:f0,wtest(x0,softloss),:gnorm,vecnorm(g),:gstep,gstep,:rstep,rstep,:lr,lr,:l2,l2,:err,wtest(x0, zeroone)))
end

function wsample(N=1; xdim=794*64, f=weval, gscale=0.8, ftol=0.05, xscale=1/sqrt(xdim), maxnf=Inf)
    xbuf = Array(Float32, 0)
    x = CudaArray(Float32, xdim)
    g = CudaArray(Float32, xdim)
    for n=1:N
        scale!(xscale, randn!(x))
        f0 = f(x, g)
        nf = 1; np = 16
        append!(xbuf, to_host(x)); println((nf,:f,f0,:x,vecnorm(x),:g,vecnorm(g)))
        while f0 > ftol
            f0 = f(axpy!(-gscale, g, x), g)
            (nf+=1) >= np && (np*=2; append!(xbuf, to_host(x)); println((nf,:f,f0,:x,vecnorm(x),:g,vecnorm(g))))
            nf % 100 == 0 && gc()
            nf >= maxnf && break
        end
        append!(xbuf, to_host(x)); println((nf,:f,f0,:x,vecnorm(x),:g,vecnorm(g)))
    end
    return reshape(xbuf, (xdim, div(length(xbuf),xdim)))
end

function dik(w0,ws,wi)
    ds = ws - w0; ns = vecnorm(ds)
    di = wi - w0; ni = vecnorm(di)
    cs = ni==0 ? 0 : dot(ds,di)/(ns*ni)
    sn = 1-cs*cs
    return (cs*ni, sn*ni)
end

function derivative(f,x,v; dx=0.01)
    nv = vecnorm(v)
    x1 = similar(x)
    f0 = f(x)
    println((0,f0,0))
    for i=1:10
        f1 = f(axpy!(dx*i/nv, v, copy!(x1,x)))
        println((i,f1,f1-f0)); f0=f1
    end
end

# Approximate the second derivative in direction d from x for function
# f.

if !isdefined(:wsample_jld) || !isa(wsample_jld,Dict)
    wsample_jld = load("wsample.jld")
    x1cpu = wsample_jld["x"][:,1]; x1gpu = CudaArray(x1cpu)
    x5cpu = wsample_jld["x"][:,5]; x5gpu = CudaArray(x5cpu)
    x9cpu = wsample_jld["x"][:,9]; x9gpu = CudaArray(x9cpu)
    d1cpu = axpy!(-1, x1cpu, copy(x9cpu)); d1gpu = CudaArray(d1cpu)
    d5cpu = axpy!(-1, x5cpu, copy(x9cpu)); d5gpu = CudaArray(d5cpu)
    t5cpu = similar(x5cpu); t5gpu = CudaArray(t5cpu)
end

function deriv1(d::CudaArray; x=x5gpu, t=t5gpu, f=weval, n=5, a=-0.1, b=0.1)
    f0 = chebft(s->f(axpy!(s/vecnorm(d),d,copy!(t,x))); n=n, a=a, b=b)
    f1 = chder(f0; n=n, a=a, b=b)
    chebev(f1,0.0; n=n, a=a, b=b)
end

function deriv2(d::CudaArray; x=x5gpu, t=t5gpu, f=weval, n=5, a=-0.1, b=0.1)
    f0 = chebft(s->f(axpy!(s/vecnorm(d),d,copy!(t,x))); n=n, a=a, b=b)
    f1 = chder(f0; n=n, a=a, b=b)
    f2 = chder(f1; n=n, a=a, b=b)
    chebev(f2,0.0; n=n, a=a, b=b)
end

function deriv2(d::Array; x=x5cpu, t=t5cpu, f=weval, n=5, a=-0.1, b=0.1)
    f0 = chebft(s->f(axpy!(s/vecnorm(d),d,copy!(t,x))); n=n, a=a, b=b)
    f1 = chder(f0; n=n, a=a, b=b)
    f2 = chder(f1; n=n, a=a, b=b)
    chebev(f2,0.0; n=n, a=a, b=b)
end

function deltaf(d::CudaArray; x=x5gpu, t=t5gpu, f=weval, n=5, a=-0.1, b=0.1)
    f0 = chebft(s->f(axpy!(s/vecnorm(d),d,copy!(t,x))); n=n, a=a, b=b)
    f1 = chder(f0; n=n, a=a, b=b)
    f2 = chder(f1; n=n, a=a, b=b)
    d1 = chebev(f1,0.0; n=n, a=a, b=b)
    d2 = chebev(f2,0.0; n=n, a=a, b=b)
    -d1^2/d2
end

normalize(x)=scale!(1/vecnorm(x),x)

function deriv12{T}(f::Function, x::CudaArray{T}, d::CudaArray{T}; n=5, a=-0.1, b=0.1, dnorm=vecnorm(d))
    global deriv12tmp
    if !isdefined(:_deriv12tmp) || !issimilar(deriv12tmp,x)
        deriv12tmp = similar(x)
    end
    f0 = chebft(s->f(axpy!(s/dnorm,d,copy!(deriv12tmp,x))); n=n, a=a, b=b)
    f1 = chder(f0; n=n, a=a, b=b)
    f2 = chder(f1; n=n, a=a, b=b)
    map(T, (chebev(f1,0.0; n=n, a=a, b=b),
            chebev(f2,0.0; n=n, a=a, b=b)))
end

# x0 is our best point, u our best direction (unit vector) so far,
# u1 and u2 are the first and second derivatives in the u direction.
# Try direction v = u + noise

# If u2 is underestimated, u1^2/u2 will be larger than f0.  We know
# minf=0 so we can correct u1^2/u2 < f0; u1^2/f0 < u2; set
# u2=max(u2,u1^2/f0).  If the limit is hit, step size will be u1/u2 =
# u1/(u1^2/f0) = f0/u1.

function enewton1{T}(f::Function, x0::BaseVector{T}; ftol=0.05, maxnf=typemax(Int), dt=3600, noise=1.0)
    f0 = f1 = f(x0)
    n0 = vecnorm(x0)
    x = similar(x0)
    v = similar(x0)
    u = normalize(randn!(similar(x0)))
    u1 = u2 = 0
    nf = np = nt = 1
    nv = 1/sqrt(length(x0))
    while true
        normalize(axpy!(1, u, scale!(noise/nv, randn!(v))))
        (u1,u2) = deriv12(f, x0, u, dnorm=1.0)
        (v1,v2) = deriv12(f, x0, v, dnorm=1.0)
        if ((v2 < 0 && v2 < u2) ||
            (u2 > 0 && v2 > 0 && v1*v1/v2 > u1*u1/u2))
            copy!(u,v); (u1,u2) = (v1,v2)
        end
        if u2 < 0 || u1*u1/u2 > f0
            dx = -f0/u1; df = -f0
        else
            dx = -u1/u2; df = -u1*u1/u2
        end
        abs(dx) > 0.1*n0 && (dx = T(sign(dx)*0.1*n0))
        f1 = f(axpy!(dx, u, copy!(x,x0))); nf+=1
        nf % 100 == 0 && gc()
        if nf >= np || time() >= nt || nf >= maxnf || f0 <= ftol
            println((nf,:f0,f0,:x0,n0,:u1,u1,:u2,u2,:dx,dx,:df,df,f1-f0,:err,T(wtest(x0, zeroone))))
            nf >= np && (np *= 2); nt = dt+time()
            (nf >= maxnf || f0 <= ftol) && break
        end
        if f1 < f0
            copy!(x0,x); f0 = f1; n0 = vecnorm(x0)
        end
    end
end


text:
using GFlops
using Test

function my_axpy!(a, x, y)
    @inbounds @simd for i in eachindex(x)
        y[i] += a*x[i]
    end
end

function my_prod(m, v)
    res = similar(v, size(m, 1))

    for i in 1:size(m, 1)
        acc = zero(eltype(m))
        for j in 1:length(v)
            acc += m[i,j]*v[j]
        end
        res[i] = acc
    end
    res
end

import BenchmarkTools: @benchmark
struct FakeResults
    times
end
macro benchmark(e)
    quote
        FakeResults([2.0])
    end
end

@testset "GFlops" begin
    @testset "Counter" begin
        let
            cnt = GFlops.Counter()
            iob = IOBuffer()
            show(iob, cnt)
            @test String(take!(iob)) == """
Flop Counter:
 add32: 0
 sub32: 0
 mul32: 0
 div32: 0
 add64: 0
 sub64: 0
 mul64: 0
 div64: 0
 sqrt32: 0
 sqrt64: 0
"""
        end
    end

    @testset "@count_ops" begin
        @testset "mul+add 64" begin
            let N = 100
                a = 2.5
                x = rand(N)
                y = similar(x)

                cnt = @count_ops my_axpy!(a, x, y)
                @test cnt.add64 == N
                @test cnt.mul64 == N
                @test GFlops.flop(cnt) == 2*N

                m = rand(N, N)
                v = rand(N)

                cnt = @count_ops(my_prod(m, v))
                @test cnt.add64 == N*N
                @test cnt.mul64 == N*N
                @test GFlops.flop(cnt) == 2*N*N
            end
        end

        @testset "mul+add 32" begin
            let N = 100
                a = 2.5f0
                x = rand(Float32, N)
                y = similar(x)

                cnt = @count_ops my_axpy!(a, x, y)
                @test cnt.add32 == N
                @test cnt.mul32 == N
                @test GFlops.flop(cnt) == 2*N

                m = rand(Float32, N, N)
                v = rand(Float32, N)

                cnt = @count_ops(my_prod(m, v))
                @test cnt.add32 == N*N
                @test cnt.mul32 == N*N
                @test GFlops.flop(cnt) == 2*N*N
            end
        end

        @testset "sqrt" begin
            let cnt = @count_ops sqrt(4.2)
                @test cnt.sqrt64 == 1
                @test GFlops.flop(cnt) == 1
            end

            let cnt = @count_ops sqrt(4.2f0)
                @test cnt.sqrt32 == 1
                @test GFlops.flop(cnt) == 1
            end
        end

        @testset "interpolated arguments" begin
            let N = 100

                T = Float64
                cnt = @count_ops my_axpy!(pi, $(rand(T, N)), $(rand(T, N)))
                @test cnt.add64 == N
                @test cnt.mul64 == N
                @test GFlops.flop(cnt) == 2*N

                T = Float32
                cnt = @count_ops my_axpy!(pi, $(rand(T, N)), $(rand(T, N)))
                @test cnt.add32 == N
                @test cnt.mul32 == N
                @test GFlops.flop(cnt) == 2*N
            end
        end

        @testset "broadcast" begin
            let N = 100

                x = 42.0
                cnt1 = @count_ops sin(x)
                cnt2 = @count_ops sin.($(fill(x, N)))
                @test GFlops.flop(cnt1) != 0
                @test cnt2 == N*cnt1

                x = 42.0f0
                cnt1 = @count_ops sin(x)
                cnt2 = @count_ops sin.($(fill(x, N)))
                @test GFlops.flop(cnt1) != 0
                @test cnt2 == N*cnt1
            end
        end
    end

    @testset "@gflops" begin
        let
            N = 100
            a = 2.5
            x = rand(N)
            y = Vector{Float64}(undef, N)

            @test @gflops(my_axpy!(a, x, y))          == N
            @test @gflops(my_axpy!(π, $(rand(N)), y)) == N
        end

        let
            N = 100
            m = rand(N, N)
            v = rand(N)
            @test @gflops(my_prod(m, v)) == N*N
        end
    end
end


text:
mutable struct QueueRecord{E}
    frames::Vector{Frame{E}}
    timestep::Float64
    nframes::Int # number of active Frames
end
function QueueRecord(::Type{E}, capacity::Int, timestep::Float64, frame_capacity::Int=100) where {E}
    frames = Array{Frame{E}}(undef, capacity)
    for i in 1 : length(frames)
        frames[i] = Frame(E, frame_capacity)
    end
    QueueRecord{E}(frames, timestep, 0)
end

Base.show(io::IO, rec::QueueRecord) = print(io, "QueueRecord(nframes=", rec.nframes, ")")

capacity(rec::QueueRecord) = length(rec.frames)
nframes(rec::QueueRecord) = rec.nframes
function nstates(rec::QueueRecord)
    retval = 0
    for frame_index in 1 : nframes(rec)
        retval += length(rec.frames[frame_index])
    end
    return retval
end

function Base.deepcopy(rec::QueueRecord)
    retval = QueueRecord(capacity(rec), rec.timestep, capacity(rec.frames[1]))
    for i in 1 : rec.nframes
        copyto!(retval.frames[i], rec.frames[i])
    end
    retval
end

pastframe_inbounds(rec::QueueRecord, pastframe::Int) = 1 ≤ 1-pastframe ≤ rec.nframes

"""
Indexed by pastframe, so pastframe == 0 is the current scene, -1 is the previous frame, etc.
"""
Base.getindex(rec::QueueRecord, pastframe::Int) = rec.frames[1 - pastframe]

get_time(rec::QueueRecord, pastframe::Int) = -get_elapsed_time(rec, pastframe)
get_timestep(rec::QueueRecord) = rec.timestep
get_elapsed_time(rec::QueueRecord, pastframe::Int) = -pastframe*rec.timestep
function get_elapsed_time(
    rec::QueueRecord,
    pastframe_farthest_back::Int,
    pastframe_most_recent::Int,
    )

    (pastframe_most_recent - pastframe_farthest_back)*rec.timestep
end

function Base.empty!(rec::QueueRecord)
    rec.nframes = 0
    return rec
end

function push_back_records!(rec::QueueRecord)
    for i in min(rec.nframes+1, capacity(rec)) : -1 : 2
        copyto!(rec.frames[i], rec.frames[i-1])
    end
    return rec
end
function Base.insert!(rec::QueueRecord{E}, frame::Frame{E}, pastframe::Int=0) where {E}
    copyto!(rec[pastframe], frame)
    return rec
end
function Base.get!(frame::Frame{E}, rec::QueueRecord{E}, pastframe::Int=0) where {E}
    copyto!(frame, rec[pastframe])
    frame
end
function update!(rec::QueueRecord{E}, frame::Frame{E}) where {E}
    push_back_records!(rec)
    insert!(rec, frame, 0)
    rec.nframes = min(rec.nframes+1, capacity(rec))
    return rec
end

function allocate_frame(rec::QueueRecord{E}) where {E}
    max_n_objects = maximum(length(rec[j]) for j in 0 : 1-length(rec))
    return Frame(E, max_n_objects)
end

const EntityQueueRecord{S,D,I} = QueueRecord{Entity{S,D,I}}

Base.length(record::QueueRecord) = nframes(record)
function Base.iterate(record::QueueRecord, state::Int64=(-nframes(record)+1))
    return state<=0 ? (record[state], state+1) : nothing
end

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


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

# ToDo: Add literature references to AdaptiveMHTuning docstring.

"""
    struct AdaptiveMHTuning <: MHProposalDistTuning

Adaptive MCMC tuning strategy for Metropolis-Hastings samplers.

Adapts the proposal function based on the acceptance ratio and covariance
of the previous samples.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct AdaptiveMHTuning <: MHProposalDistTuning
    "Controls the weight given to new covariance information in adapting the
    proposal distribution."
    λ::Float64 = 0.5

    "Metropolis-Hastings acceptance ratio target, tuning will try to adapt
    the proposal distribution to bring the acceptance ratio inside this interval."
    α::IntervalSets.ClosedInterval{Float64} = ClosedInterval(0.15, 0.35)

    "Controls how much the spread of the proposal distribution is
    widened/narrowed depending on the current MH acceptance ratio."
    β::Float64 = 1.5

    "Interval for allowed scale/spread of the proposal distribution."
    c::IntervalSets.ClosedInterval{Float64} = ClosedInterval(1e-4, 1e2)

    "Reweighting factor. Take accumulated sample statistics of previous
    tuning cycles into account with a relative weight of `r`. Set to
    `0` to completely reset sample statistics between each tuning cycle."
    r::Real = 0.5
end

export AdaptiveMHTuning

mutable struct ProposalCovTuner{
    S<:MCMCBasicStats
} <: AbstractMCMCTunerInstance
    config::AdaptiveMHTuning
    stats::S
    iteration::Int
    scale::Float64
end

(tuning::AdaptiveMHTuning)(chain::MHIterator) = ProposalCovTuner(tuning, chain)

function ProposalCovTuner(tuning::AdaptiveMHTuning, chain::MHIterator)
    m = totalndof(getdensity(chain))
    scale = 2.38^2 / m
    ProposalCovTuner(tuning, MCMCBasicStats(chain), 1, scale)
end

function _cov_with_fallback(d)
    rng = bat_determ_rng()
    smplr = bat_sampler(d)
    T = float(eltype(rand(rng, smplr)))
    n = totalndof(varshape(d))
    C = fill(T(NaN), n, n)
    try
        C[:] = cov(d)
    catch err
        if err isa MethodError
            C[:] = cov(nestedview(rand(rng, smplr, 10^5)))
        else
            throw(err)
        end
    end
    return C
end

_approx_cov(target::Distribution) = _cov_with_fallback(target)
_approx_cov(target::DistLikeDensity) = _cov_with_fallback(target)
_approx_cov(target::AbstractPosteriorDensity) = _approx_cov(getprior(target))
_approx_cov(target::BAT.TransformedDensity{<:Any,<:BAT.DistributionTransform}) =
    BAT._approx_cov(target.trafo.target_dist)
_approx_cov(target::RenormalizedDensity) = _approx_cov(parent(target))
_approx_cov(target::DensityWithDiff) = _approx_cov(parent(target))

function tuning_init!(tuner::ProposalCovTuner, chain::MHIterator, max_nsteps::Integer)
    Σ_unscaled = _approx_cov(getdensity(chain))
    Σ = Σ_unscaled * tuner.scale

    chain.proposaldist = set_cov(chain.proposaldist, Σ)

    nothing
end

tuning_reinit!(tuner::ProposalCovTuner, chain::MCMCIterator, max_nsteps::Integer) = nothing

function tuning_postinit!(tuner::ProposalCovTuner, chain::MHIterator, samples::DensitySampleVector)
    # The very first samples of a chain can be very valuable to init tuner
    # stats, especially if the chain gets stuck early after:
    stats = tuner.stats
    append!(stats, samples)
end

function tuning_update!(tuner::ProposalCovTuner, chain::MHIterator, samples::DensitySampleVector)
    stats = tuner.stats
    stats_reweight_factor = tuner.config.r
    reweight_relative!(stats, stats_reweight_factor)
    # empty!.(stats)
    append!(stats, samples)

    config = tuner.config

    α_min = minimum(config.α)
    α_max = maximum(config.α)

    c_min = minimum(config.c)
    c_max = maximum(config.c)

    β = config.β

    t = tuner.iteration
    λ = config.λ
    c = tuner.scale
    Σ_old = Matrix(get_cov(chain.proposaldist))

    S = convert(Array, stats.param_stats.cov)
    a_t = 1 / t^λ
    new_Σ_unscal = (1 - a_t) * (Σ_old/c) + a_t * S

    α = eff_acceptance_ratio(chain)

    max_log_posterior = stats.logtf_stats.maximum

    if α_min <= α <= α_max
        chain.info = MCMCIteratorInfo(chain.info, tuned = true)
        @debug "MCMC chain $(chain.info.id) tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"
    else
        chain.info = MCMCIteratorInfo(chain.info, tuned = false)
        @debug "MCMC chain $(chain.info.id) *not* tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"

        if α > α_max && c < c_max
            tuner.scale = c * β
        elseif α < α_min && c > c_min
            tuner.scale = c / β
        end
    end

    Σ_new = new_Σ_unscal * tuner.scale

    chain.proposaldist = set_cov(chain.proposaldist, Σ_new)
    tuner.iteration += 1

    nothing
end

tuning_finalize!(tuner::ProposalCovTuner, chain::MCMCIterator) = nothing

tuning_callback(::ProposalCovTuner) = nop_func


text:
include("src/Blake3Hash.jl")

using BenchmarkTools
using .Blake3Hash

@benchmark begin
    Blake3Hash.Simd128.update!(h, d)
    Blake3Hash.Simd128.digest(h, o)
end setup=begin
    h = Blake3Hash.Simd128.Blake3Ctx()
    d = rand(UInt8, 2*8192)
    o = Vector{UInt8}(undef, 32)
end


text:
import YAML

# convenience extensions of YAML library
import Base.getindex
import Base.keys
import Base.values

Base.length(d::YAML.SequenceNode) = length(d.value)
Base.getindex(d::YAML.SequenceNode, k::Integer) = d.value[k]
Base.iterate(d::YAML.SequenceNode, state=1)  = state>length(d) ? nothing : (d.value[state], state+1)
Base.keys(s::YAML.MappingNode) =  [e[1].value for e=s.value]
Base.values(s::YAML.MappingNode) =  [e[2].value for e=s.value]
Base.getindex(d::YAML.MappingNode, s::AbstractString) = d.value[findfirst(isequal(s),keys(d))][2]
Base.iterate(d::YAML.MappingNode, state=1)  = state>length(d) ? nothing : ((d.value[state][1].value=>d.value[state][2]), state+1)
Base.getindex(d::YAML.MappingNode, s::Symbol) = d[string(s)]
Base.length(d::YAML.MappingNode) = length(d.value)

const symbols_from_greek = Dict(
    :μ=>:mu,
    :σ=>:sigma,
    :Σ=>:Sigma,
    :ρ=>:rho,
    :λ=>:lambda,
    :β=>:beta,
    :δ=>:delta
)
const symbols_to_greek = Dict(v=>k for (k,v) in symbols_from_greek)

#
#  Language object: essentially a table from tag to a constructor
#

struct Language
    objects::Dict{String, Any}
end

Language() = Language(Dict{String, Any}())

function add_language_element!(lang, fun, tag::AbstractString)
    lang.objects[tag] = fun
end
function add_language_element!(lang, fun, tags::Vector{AbstractString})
    for tag in tags
        add_language_element!(lang, fun, tag)
    end
end
function add_language_elements!(lang, elts::AbstractDict{<:AbstractString, <:Any})
    for (k,v) in elts
        add_language_element!(lang, v, k)
    end
end

####
#### Functions to import yaml tree (with locations and tags)
####

const yaml_basis_tags = [
    "tag:yaml.org,2002:null",
    "tag:yaml.org,2002:bool",
    "tag:yaml.org,2002:int",
    "tag:yaml.org,2002:float",
    "tag:yaml.org,2002:binary",
    "tag:yaml.org,2002:timestamp",
    "tag:yaml.org,2002:omap",
    "tag:yaml.org,2002:pairs",
    "tag:yaml.org,2002:set",
    "tag:yaml.org,2002:str",
    "tag:yaml.org,2002:seq",
    "tag:yaml.org,2002:map",
]

# with upcoming version of YAML.jl
function yaml_node_from_string(txt::AbstractString)
    cons = YAML.Constructor()
    YAML.add_multi_constructor!((c,s,m)->m, cons, "tag:yaml.org")
    YAML.add_multi_constructor!((c,s,m)->m, cons, "!")
    data = YAML.load(txt, cons)
    return data
end

# function yaml_node_from_string(txt::AbstractString, minilanguage::Language=Language())
#     # Didn't find how to access the top node more easily.
#     #
#     yml_types = Dict{String,Function}()
#     for tag in yaml_basis_tags
#         yml_types[tag] = (c,n)->n
#     end
#     for tag in keys(minilanguage.objects)
#         yml_types[tag] = (c,n)->n
#     end
#     return YAML.load(txt, yml_types)
# end

function yaml_node_from_file(fn::AbstractString, minilanguage::Language=Language())
    txt = open(f->read(f,String), fn)
    txt = replace(txt, "\r"=>"")
    return yaml_node_from_string(txt) # , minilanguage)
end

##
# Function to evaluate string from calibration
#
# TODO: not safe and hazardous error policy
##

function eval_string(s::AbstractVector{Expr}, d:: AbstractDict{Symbol, <:Number})

        # extract expressions and variable names in proper order
        nms = collect(keys(d))
        exprs = collect(values(d))

        # build expression to Core.evaluate system in correct order
        to_eval = Expr(:block)
        to_eval.args = [:($(i[1])=$(i[2])) for i in zip(nms, exprs)]

        # add one line to return a tuple of all data
        ret = Expr(:tuple); ret.args = s

        # now Core.evaluate and get data
        expr = :(
            let
                $to_eval
                $ret
            end
        )
        data = eval(expr)
end

eval_string(s::Symbol, d::AbstractDict{Symbol, <:Number}) = get(d,s,s)

eval_string(s::Expr, d::AbstractDict{Symbol, <:Number}) = eval_string([s], d)[1]

function eval_string(s::AbstractString, d::AbstractDict{Symbol, <:Number}=Dict{Symbol, Float64}())
        expr = Meta.parse(s)
        eval_string(expr, d)
end

###
# Some errors potentially occuring when parsing yaml structure
##

struct Position
    l1::Int
    c1::Int
    l2::Int
    c2::Int
end

abstract type LanguageError <: Exception end

struct ParsingError{T<:LanguageError} <: Exception
    pos::Position
    error::T
end

function ParsingError(node::YAML.Node, error)
    ms, em = node.start_mark, node.end_mark
    l1, c1 = Int(ms.line), Int(ms.column)
    l2, c2 = Int(em.line), Int(em.column)
    pos = Position(l1, c1, l2, c2)
    return ParsingError(pos, error)
end

struct LangInvalidCall <: LanguageError
    tag::String
    err::MethodError
end

struct LangUnknownTag <: LanguageError
    tag::String
end

function message(e::LangInvalidCall)
    s = IOBuffer()
    showerror(s, e.err)
    msg = String(take!(s))
    return msg
end
message(e::LangUnknownTag) = "Unknown tag: \"$(e.tag)\""

Base.showerror(io::IO, e::LanguageError) = print(io, message(e))
Base.showerror(io::IO, e::ParsingError) = print(io, "Language Error: $(e.pos.l1), $(e.pos.l2): $(message(e.error)))")

#
# construct object from tag and arguments
#

function construct_object(lang, tag, args::AbstractDict{<:Any,<:Any})
    if !( tag in keys(lang.objects))
        throw(LangUnknownTag(tag))
    end
    constructor = lang.objects[tag]
    try
        return constructor(;args...)
    catch err
        throw( LangInvalidCall(tag, err) )
    end
end

function construct_object(lang, tag, args)
    if !(tag in keys(lang.objects))
        throw(LangUnknownTag(tag))
    end
    constructor = lang.objects[tag]
    try
        return constructor(args...)
    catch err
        throw( LangInvalidCall(tag, err) )
    end
end

#
# option types
# maybe this is not needed anymore thanks to constant propagation

abstract type GreekTolerance end
struct FromGreek<:GreekTolerance
    # greek characters are automatically converted to alphanumerical
end
struct ToGreek<:GreekTolerance
    # alphanumerical characters are automatically converted to greek
end
struct NoGreek<:GreekTolerance
    # no conversion
end

###
# evaluate yaml structures
#

# scalar: text, float, int, ...
function eval_node(s::YAML.ScalarNode, calibration::AbstractDict{Symbol, <:Number}, minilang::Language=Language(), greek_tol::GreekTolerance=NoGreek())
    v = s.value
    tag = s.tag
    if tag == "tag:yaml.org,2002:str"
        return try
            eval_string(v, calibration)
        catch err
            # TODO: fix behaviour when type is not recognized
            v
        end
    elseif tag == "tag:yaml.org,2002:int"
        return parse(Int, v)
    elseif tag == "tag:yaml.org,2002:float"
        return parse(Float64, v)
    end
end

# sequences
function eval_node(node::YAML.SequenceNode, calibration::AbstractDict{Symbol, <:Number}, minilang::Language=Language(), greek_tol::GreekTolerance=NoGreek())
    children = [eval_node(ch, calibration, minilang, greek_tol) for ch in node]
    tag = node.tag
    if tag == "tag:yaml.org,2002:seq"
        if typeof(children) <: Vector{Vector{Float64}}
            res = vcat([x' for x in children]...)
            return res
        else
            return children
        end
    else # custom object with positional arguments
        childs = tuple(children...)
        try
            return construct_object(minilang, tag, childs)
        catch err
            throw(ParsingError(node, err))
        end
    end
end

# map
function eval_node(node::YAML.MappingNode, calibration::AbstractDict{Symbol, <:Number}, minilang::Language=Language(), greek_tol::GreekTolerance=NoGreek())
    d = Dict()
    for i=1:length(node)
        k = Symbol(node.value[i][1].value)
        if typeof(greek_tol)<:FromGreek
            if haskey(symbols_from_greek, k)
                k = symbols_from_greek[k]
            end
        elseif typeof(greek_tol)<:ToGreek
            if haskey(symbols_to_greek, k)
                k = symbols_to_greek[k]
            end
        end
        v = eval_node(node.value[i][2], calibration, minilang, greek_tol)
        d[k] = v
    end
    if node.tag == "tag:yaml.org,2002:map"
        return d
    else # custom object with positional arguments
        try
            return construct_object(minilang, node.tag, d)
        catch err
            throw(ParsingError(node, err))
        end
    end
end

max_issues_repo_name: EconForge/Dolang


text:
"""
    xbrs(RR::Vector{<:Real}, SBP::Vector{<:Real};...)

Calculate the xBRS index for the assesment of baroreflex sensitivity. Based on Westerhof, B. E. et al. (2004). Time-domain cross-correlation baroreflex sensitivity: performance on the EUROBAVAR data set. Journal of hypertension, 22(7), 1371-1380

# Keyword Arguments

* minCor: Minimal significant correlation, dafaults to 0.632 (p = 0.05, two-sided for 10 s window)
* tExcerpt: Length of the sliding window in seconds, defaults to 10 s
* delays: Which delays, shift of the RR window shall be considered, defaults to 0:5

# Return

Returns an xBRS structure, that can be be plotted for visual inspection.
"""
function xbrs(RR::Vector{<:Real}, SBP::Vector{<:Real}; minCor::Real = 0.632, tExcerpt::Int = 10, delays::Union{Vector{Int}, UnitRange} = 0:5)
    1 > minCor > 0 || throw(DomainError("The minimal correlation hast to be postive and ≤1"))

    # Resample signals at 1 Hz through cubic splines
    time = cumsum(RR./1000) # ms to s
    time1Hz = 1:round(Int, time[end])
    S1Hz = Spline1D(time, SBP; k = 3)(time1Hz)
    RR1Hz = Spline1D(time, RR; k = 3)(time1Hz)

    # Correlation and xBRS
    len = time1Hz[end] - tExcerpt - maximum(delays)
    maxCors = zeros(Float64, len)
    delayMaxCors = zeros(Int, len)
    xBRSs = zeros(Float64, len)
    xBRSs2 = zeros(Float64, len)
    for i in 1:len
        maxCor, delayMaxCor = findmax(map(delay -> cor(S1Hz[i:i+tExcerpt-1], RR1Hz[(i:i+tExcerpt-1) .+ delay]), delays))
        maxCors[i] = maxCor
        delayMaxCors[i] = delayMaxCor - 1
        xBRSs[i] = std(RR1Hz[(i:i+tExcerpt-1) .+ delayMaxCors[i]]) / std(S1Hz[i:i+tExcerpt-1]) # 2017
        xBRSs2[i] = linReg(S1Hz[i:i+tExcerpt-1], RR1Hz[(i:i+tExcerpt-1) .+ delayMaxCors[i]])[1] / maxCors[i] # 2004
    end

    valid = findall(x -> x >= minCor, maxCors)
    xBRSg = gmean(xBRSs[valid])
    xBRSg2 = gmean(xBRSs2[valid])

    return xBRS(value = xBRSg,
    xBRSg2 = xBRSg2,
    xBRSs = xBRSs,
    xBRSs2 = xBRSs2,
    valid = valid,
    nValid = length(valid),
    maxCors = maxCors,
    delayMaxCors = delayMaxCors,
    RR = RR,
    SBP = SBP)
end

"""
Struct that stores all information regarting the xBRS etimation. The final result is stored in 'xBRSg'. It can be plotted for visual inspection.
"""
@with_kw mutable struct xBRS
    value::Real = 0
    xBRSg2::Real = 0
    nValid::Int = 0
    xBRSs::Vector{<:Real} = Float64[]
    xBRSs2::Vector{<:Real} = Float64[]
    maxCors::Vector{<:Real} = Float64[]
    delayMaxCors::Vector{Int} = Int[]
    valid::Vector{Int} = Int[]
    RR::Vector{<:Real} = Float64[]
    SBP::Vector{<:Real} = Float64[]
end

"""
linReg(x, y)

Solves the linear regression problem y = mx + b using least squares method through '\\'. Returns (m, b)
"""
linReg(x::Vector{<:Real}, y::Vector{<:Real}) = ([x ones(length(x))] \ y)

"""
    gmean(x::Vector{<:Real})

Compute the geometric mean for the positive Vector x.

"""
gmean(x::Vector{<:Real}) = exp(1/length(x) * sum(log.(x)))

@recipe function f(res::xBRS)
    layout := (4,1)

    @series begin
        subplot := 1
        yguide := "SBP [mmHg]"
        label := ""
        seriescolor := :grey

        res.SBP
    end

    @series begin
        subplot := 2
        yguide := "RR [ms]"
        label := ""
        seriescolor := :grey

        res.RR
    end

    @series begin
        subplot := 3
        marker := :d
        yguide := "correlation"
        label := ""
        seriescolor := :grey

        res.valid, res.maxCors[res.valid]
    end

    @series begin
        subplot := 3
        seriestype := :hline
        yguide := "correlation"
        label := "mean valid cor: $(round(mean(res.maxCors[res.valid]), digits = 2))"
        legend := :topright
        seriescolor := :red

        [mean(res.maxCors[res.valid])]
    end

    @series begin
        subplot := 4
        # seriestype := :scatter
        marker := :d
        yguide := "xBRS [ms/mmHg]"
        label := ""
        seriescolor := :grey

        res.valid, res.xBRSs[res.valid]
    end

    @series begin
        subplot := 4
        seriestype := :hline
        yguide := "xBRS [ms/mmHg]"
        xguide := "beat #"
        label := "xBRS: $(round(res.value, digits = 2)) ms/mmHg"
        legend := :topright
        seriescolor := :steelblue

        [res.xBRSg]
    end
end


text:
function max_value_changed(obj::ParameterObserver, arg0::jdouble, arg1::ParameterDriver)
    return jcall(obj, "maxValueChanged", void, (jdouble, ParameterDriver), arg0, arg1)
end

function min_value_changed(obj::ParameterObserver, arg0::jdouble, arg1::ParameterDriver)
    return jcall(obj, "minValueChanged", void, (jdouble, ParameterDriver), arg0, arg1)
end

function name_changed(obj::ParameterObserver, arg0::JString, arg1::ParameterDriver)
    return jcall(obj, "nameChanged", void, (JString, ParameterDriver), arg0, arg1)
end

function reference_date_changed(obj::ParameterObserver, arg0::AbsoluteDate, arg1::ParameterDriver)
    return jcall(obj, "referenceDateChanged", void, (AbsoluteDate, ParameterDriver), arg0, arg1)
end

function reference_value_changed(obj::ParameterObserver, arg0::jdouble, arg1::ParameterDriver)
    return jcall(obj, "referenceValueChanged", void, (jdouble, ParameterDriver), arg0, arg1)
end

function scale_changed(obj::ParameterObserver, arg0::jdouble, arg1::ParameterDriver)
    return jcall(obj, "scaleChanged", void, (jdouble, ParameterDriver), arg0, arg1)
end

function selection_changed(obj::ParameterObserver, arg0::jboolean, arg1::ParameterDriver)
    return jcall(obj, "selectionChanged", void, (jboolean, ParameterDriver), arg0, arg1)
end

function value_changed(obj::ParameterObserver, arg0::jdouble, arg1::ParameterDriver)
    return jcall(obj, "valueChanged", void, (jdouble, ParameterDriver), arg0, arg1)
end


text:
println("Testing...")

using myJuliaUtils
using Test
using SafeTestsets

@time @safetestset "Testing vech" begin include("TestingVech.jl") end
@time @safetestset "Testing Wish" begin include("TestWish/TestingWish.jl") end
@time @safetestset "Testing cholx and randx" begin include("TestingCholXRandX.jl") end
@time @safetestset "Testing some functions in extfunc" begin include("TestingExtfunc.jl") end
@time @safetestset "Testing ACF function" begin include("testAcf/testAcf.jl") end
@time @safetestset "Testing decVCV function" begin include("testdecVCV.jl") end
@time @safetestset "Testing the fitting of distributions" begin include("testFitDist.jl") end
@time @safetestset "Testing the data manipulation function" begin include("test_data_manipulation.jl") end
@time @safetestset "Testing generate data functions" begin include("test_generate.jl") end
@time @safetestset "Testing PIT estimation" begin include("testPIT.jl") end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
lines = readlines("input.txt")
octomap = parse.(Int, reduce(hcat, collect.(lines)))
idx = CartesianIndices(octomap)
minval = typemin(Int)

o = CartesianIndex(1,1)
V = CartesianIndices((-1:1, -1:1))
padded = fill(minval, size(octomap) .+ 2)
padded[idx .+ o] = octomap[idx]

for step = 1:1000
    padded .+= 1
    flashmap = fill(false, size(padded))
    while (yxs = findall(padded .>= 10)) |> length > 0
        flashmap[yxs] .= true
        padded[yxs] .= minval # Prevent any further flashing for this step
        [padded[yx .+ V] .+= 1 for yx = yxs]
    end
    padded[findall(flashmap)] .= 0

    if sum(flashmap) == 100
        println(step)
        break
    end
end


text:
@testset "Hello Function Tests" begin

  @test isdefined(Julz, :hello) == true

  originalSTDOUT = STDOUT

  (outRead, outWrite) = redirect_stdout()

  Julz.hello()

  close(outWrite)

  data = readavailable(outRead)

  close(outRead)

  redirect_stdout(originalSTDOUT)

  @test String(data) == "Hello, World!\n"

end


text:
using Test
using DynamicalBilliards
using DynamicalBilliards.Testing

function test_no_escape(p, bd, N = 1e4)
    xmin, ymin, xmax, ymax = cellsize(bd)

    xt, yt, vyt, vxt, tvec = timeseries(p, bd, Int(N); dt = Inf)

    @test length(xt) == N + 1
    @test typeof(xt) == typeof(tvec) == Vector{eltype(p)}

    @testset "$(tag(p, bd)) correct tvector" begin
        for i in 2:length(tvec)
            @test tvec[i] > tvec[i-1]
        end
    end

    dx1 = minimum(xt) - xmin
    dx2 = xmax - maximum(xt)
    dy1 = minimum(yt) - ymin
    dy2 = ymax - maximum(yt)

    println(tag(p, bd))
    println("dx1: ", dx1, " dx2: ", dx2)
    println("dy1: ", dy1, " dy2: ", dy2)
    println()
    @testset "$(tag(p, bd)) no escape" begin
        for d in (dx1, dx2, dy1, dy2)
            @test d ≥ -acclevel(p, bd)
        end
    end
end

billiards_testset("No escape", test_no_escape; caller = ergodic_tests)

function test_visited_obstacles(p, bd, N = 50)

    ts, obst = visited_obstacles(p, bd, N)
    cts, = evolve(p, bd, N)

    @test length(ts) == length(cts)
    cc = cumsum(cts)
    for i in eachindex(cc)
        @test cc[i] ≈ ts[i]  atol = 1e-12
    end
end
billiards_testset("Visited obstacles", test_visited_obstacles; caller = ergodic_tests)

function test_movin_periodic(p, bd, N = 1e3)
    xmin, ymin, xmax, ymax = cellsize(bd)

    xt, yt = timeseries(p, bd, Int(N), dt = Inf)
    @test length(xt) == N + 1
    @test typeof(xt) == typeof(xt) == Vector{eltype(p)}

    dx1 = minimum(xt) - xmin
    dx2 = xmax - maximum(xt)
    dy1 = minimum(yt) - ymin
    dy2 = ymax - maximum(yt)

    @test maximum(abs.((dx1, dx2))) > xmax - xmin
    @test maximum(abs.((dy1, dy2))) > ymax - ymin

end

billiards_testset("movin periodic", test_movin_periodic; caller = periodic_tests)

function test_raysplit_ts_inside(p, bd, ray, N = 1e4)

    xmin, ymin, xmax, ymax = cellsize(bd)

    xt, yt = timeseries(p, bd, Int(N), ray, dt = Inf)
    @test length(xt) == N + 1
    @test typeof(xt) == typeof(xt) == Vector{eltype(p)}

    dx1 = minimum(xt) - xmin
    dx2 = xmax - maximum(xt)
    dy1 = minimum(yt) - ymin
    dy2 = ymax - maximum(yt)

    for d in (dx1, dx2, dy1, dy2)
        @test d ≥ -acclevel(p, bd)
    end

end

billiards_testset("timeseries raysplit", test_raysplit_ts_inside; caller = ray_tests)


text:
# ---
# title: 419. Battleships in a Board
# id: problem419
# author: Indigo
# date: 2021-06-21
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/battleships-in-a-board/description/>
# hidden: true
# ---
#
# Given an 2D board, count how many battleships are in it. The battleships are
# represented with `'X'`s, empty slots are represented with `'.'`s. You may
# assume the following rules:
#
#   * You receive a valid board, made of only battleships or empty slots.
#   * Battleships can only be placed horizontally or vertically. In other words, they can only be made of the shape `1xN` (1 row, N columns) or `Nx1` (N rows, 1 column), where N can be of any size.
#   * At least one horizontal or vertical cell separates between two battleships - there are no adjacent battleships.
#
# **Example:**
#
#
#
#     X..X
#     ...X
#     ...X
#
#
# In the above board there are 2 battleships.
#
# **Invalid Example:**
#
#
#
#     ...X
#     XXXX
#     ...X
#
#
# This is an invalid board that you will not receive - as battleships will
# always have a cell separating between them.
#
# **Follow up:**
# Could you do it in **one-pass** , using only **O(1) extra memory** and
# **without modifying** the value of the board?
#
#
## @lc code=start
using LeetCode

function count_battleships(board::Matrix{Char})
    res = 0
    for I in CartesianIndices(board)
        board[I] == '.' && continue
        i, j = I.I
        (i == 1 || board[i - 1, j] == '.') &&
            (j == 1 || board[i, j - 1] == '.') &&
            (res += 1)
    end
    return res
end
## @lc code=end