text:
module LLW2d

# Linear Long Wave (LLW) tsunami in 2D Cartesian Coordinate

struct Matrices{T,M<:AbstractMatrix{T}}
    absorbing_boundary::M
    ocean_depth::M
    x_averaged_depth::M
    y_averaged_depth::M
    land_filter_m::M
    land_filter_n::M
    land_filter_e::M
end

const g_n = 9.80665

# Set station locations.
function set_stations!(ist::AbstractVector{Int},
                       jst::AbstractVector{Int},
                       station_distance_x::T,
                       station_distance_y::T,
                       station_boundary_x::T,
                       station_boundary_y::T,
                       dx::T,
                       dy::T) where T
    # synthetic station locations
    nst = 0

    # let's check ist and jst have a square number of elements before computing n
    @assert mod(sqrt(length(ist)),1.0) ≈ 0
    @assert mod(sqrt(length(jst)),1.0) ≈ 0
    n = floor(Int, sqrt(length(ist)))

    @inbounds for i in 1:n, j in 1:n
        nst += 1
        ist[nst] = floor(Int, (station_boundary_x + (i - 1) * station_distance_x) / dx)
        jst[nst] = floor(Int, (station_boundary_y + (j - 1) * station_distance_y) / dy)
    end
end

function timestep!(dx_buffer::AbstractMatrix{T},
                   dy_buffer::AbstractMatrix{T},
                   eta1::AbstractMatrix{T},
                   mm1::AbstractMatrix{T},
                   nn1::AbstractMatrix{T},
                   eta0::AbstractMatrix{T},
                   mm0::AbstractMatrix{T},
                   nn0::AbstractMatrix{T},
                   absorbing_boundary::AbstractMatrix{T},
                   x_averaged_depth::AbstractMatrix{T},
                   y_averaged_depth::AbstractMatrix{T},
                   land_filter_m::AbstractMatrix{T},
                   land_filter_n::AbstractMatrix{T},
                   land_filter_e::AbstractMatrix{T},
                   dx::Real,dy::Real,dt::Real) where T
    nx, ny = size(eta1)
    @assert (nx, ny) == size(mm1) == size(nn1) == size(eta0) == size(mm0) ==
        size(nn0) == size(x_averaged_depth) == size(y_averaged_depth) ==
        size(land_filter_m) == size(land_filter_n) == size(land_filter_e) ==
        size(absorbing_boundary)

    # diffs
    for j in 1:ny
        for i in 2:nx
            @inbounds dx_buffer[i,j] = (eta0[i,j] - eta0[i - 1,j]) / dx
        end
        @inbounds dx_buffer[1,j] = (eta0[1,j] - 0) / dx
    end
    for i in 1:nx
        for j in 2:ny
            @inbounds dy_buffer[i,j] = (eta0[i,j] - eta0[i, j - 1]) / dy
        end
        @inbounds dy_buffer[i,1] = (eta0[i,1] - 0) / dy
    end

    # Update Velocity
    for j in 1:ny, i in 1:nx
        @inbounds mm1[i,j] = mm0[i, j] - g_n * x_averaged_depth[i, j] * dx_buffer[i, j] * dt
        @inbounds nn1[i,j] = nn0[i, j] - g_n * y_averaged_depth[i, j] * dy_buffer[i, j] * dt
    end

    # boundary condition
    for j in 1:ny, i in 1:nx
        @inbounds mm1[i, j] = mm1[i, j] * land_filter_m[i, j] * absorbing_boundary[i, j]
        @inbounds nn1[i, j] = nn1[i, j] * land_filter_n[i, j] * absorbing_boundary[i, j]
    end

    # diffs
    for j in 1:ny
        @inbounds dx_buffer[nx, j] = (-mm1[nx, j]) / dx
        for i in 1:(nx-1)
            @inbounds dx_buffer[i, j] = (mm1[i + 1, j] - mm1[i, j]) / dx
        end
    end
    for i in 1:nx
        @inbounds dy_buffer[i, ny] = (-nn1[i, ny]) / dy
        for j in 1:(ny-1)
            @inbounds dy_buffer[i, j] = (nn1[i,j + 1] - nn1[i, j]) / dy
        end
    end

    # Update Wave Heigt
    for j in 1:ny, i in 1:nx
        @inbounds eta1[i, j] = eta0[i, j] - (dx_buffer[i, j] + dy_buffer[i, j]) * dt
    end

    # boundary condition
    for j in 1:ny, i in 1:nx
        @inbounds eta1[i, j] = eta1[i, j] * land_filter_e[i, j] * absorbing_boundary[i, j]
    end
    return eta1, mm1, nn1
end

function timestep!(dx_buffer::AbstractMatrix{T},
                   dy_buffer::AbstractMatrix{T},
                   eta1::AbstractMatrix{T},
                   mm1::AbstractMatrix{T},
                   nn1::AbstractMatrix{T},
                   eta0::AbstractMatrix{T},
                   mm0::AbstractMatrix{T},
                   nn0::AbstractMatrix{T},
                   matrices::Matrices{T},
                   dx::Real,dy::Real,dt::Real) where T
    # Unpack the relevant fields of `matrices`
    return timestep!(dx_buffer, dy_buffer, eta1, mm1, nn1, eta0, mm0, nn0,
                     matrices.absorbing_boundary, matrices.x_averaged_depth,
                     matrices.y_averaged_depth, matrices.land_filter_m,
                     matrices.land_filter_n, matrices.land_filter_e, dx, dy, dt)
end

function setup(nx::Int,
               ny::Int,
               bathymetry_val::Real,
               absorber_thickness_fraction::Real,
               apara::Real,
               cutoff_depth::Real,
               T::DataType = Float64)
    # Memory allocation
    ocean_depth = Matrix{T}(undef, nx, ny)
    x_averaged_depth = Matrix{T}(undef, nx, ny)
    y_averaged_depth = Matrix{T}(undef, nx, ny)
    absorbing_boundary = ones(T, nx, ny)
    land_filter_m = ones(T, nx, ny) # land filters
    land_filter_n = ones(T, nx, ny) # "
    land_filter_e = ones(T, nx, ny) # "

    nxa = floor(Int, nx * absorber_thickness_fraction)
    nya = floor(Int, nx * absorber_thickness_fraction)

    # Bathymetry set-up. Users may need to modify it
    fill!(ocean_depth, bathymetry_val)
    @inbounds for j in 1:ny, i in 1:nx
        if ocean_depth[i,j] < 0
            ocean_depth[i,j] = 0
        elseif ocean_depth[i,j] < cutoff_depth
            ocean_depth[i,j] = cutoff_depth
        end
    end

    # average bathymetry for staabsorbing_boundaryered-grid computation
    for j in 1:ny
        for i in 2:nx
            x_averaged_depth[i, j] = (ocean_depth[i, j] + ocean_depth[i - 1, j]) / 2
            if ocean_depth[i, j] <= 0 || ocean_depth[i - 1, j] <= 0
                x_averaged_depth[i, j] = 0
            end
        end
        x_averaged_depth[1, j] = ocean_depth[1, j]
    end
    for i in 1:nx
        for j in 2:ny
            y_averaged_depth[i, j] = (ocean_depth[i, j] + ocean_depth[i, j - 1]) / 2
            if ocean_depth[i, j] <= 0 || ocean_depth[i, j - 1] <= 0
                y_averaged_depth[i, j] = 0
            end
        end
        y_averaged_depth[i, 1] = ocean_depth[i, 1]
    end

    # Land filter
    @inbounds for j in 1:ny,i in 1:nx
        (x_averaged_depth[i, j] < 0) && (land_filter_m[i, j] = 0)
        (y_averaged_depth[i, j] < 0) && (land_filter_n[i, j] = 0)
        (ocean_depth[i, j] < 0) && (land_filter_e[i, j] = 0)
    end

    # Sponge absorbing boundary condition by Cerjan (1985)
    @inbounds for j in 1:ny, i in 1:nx
        if i <= nxa
            absorbing_boundary[i, j] *= exp(-((apara * (nxa - i)) ^ 2))
        end
        if i >= nx - nxa + 1
            absorbing_boundary[i, j] *= exp(-((apara * (i - nx + nxa - 1)) ^ 2))
        end
        if j <= nya
            absorbing_boundary[i, j] *= exp(-((apara * (nya - j)) ^ 2))
        end
        if j >= ny - nya + 1
            absorbing_boundary[i, j] *= exp(-((apara * (j - ny + nya - 1)) ^ 2))
        end
    end
    return Matrices(absorbing_boundary, ocean_depth, x_averaged_depth, y_averaged_depth,
                    land_filter_m, land_filter_n, land_filter_e)
end

# Calculates the initial surface height at point x,y
function get_initial_height(x::T, y::T, peak_x::T, peak_y::T, peak_height::T, cutoff::T) where T

    dx = peak_x - x
    dy = peak_y - y
    distance_to_peak = sqrt(dx^2 + dy^2)

    if distance_to_peak <= cutoff
        return peak_height / 4.0 * ((1 + cospi(dx / cutoff)) * (1 + cospi(dy / cutoff)))
    else
        return 0.0
    end

end

# Initializes the surface height on all points of the grid
function initheight!(height::AbstractMatrix{T},
                     ocean_depth::AbstractMatrix{T},
                     dx::Real,dy::Real,cutoff_distance::Real,
                     peak_height::Real, peak_position::AbstractVector{T}) where T

    @assert size(height) == size(ocean_depth)

    nx, ny = size(height)

    x = (1:nx) .* dx
    y = (1:ny) .* dy

    height .= get_initial_height.(x', y, peak_position[1], peak_position[2], peak_height, cutoff_distance)
    height[ocean_depth .< eps(T)] .= 0.0

end

# Initializes the surface height on all points of the grid
function initheight!(height::AbstractMatrix{T},
                     matrices::Matrices{T},
                     dx::Real,dy::Real,cutoff_distance::Real,
                     peak_height::Real, peak_position::AbstractVector{T}) where T
    # Unpack the relevant field of `matrices`
    initheight!(height, matrices.ocean_depth, dx, dy, cutoff_distance, peak_height, peak_position)
end

end # module


text:
module test_bukdu_plugs_contentparsers_json

using Test
using JSON2
using Bukdu
using .Bukdu.Plug.ContentParsers

buf = IOBuffer(JSON2.write((k=1,)))
@test ContentParsers.parse(ContentParsers.MergedJSON, buf) == Pair{String,Any}["k"=>1]

buf = IOBuffer(JSON2.write((k=1,)))
@test ContentParsers.parse(ContentParsers.JSONDecoder, buf) == Pair{String,Any}["json"=>(k=1,)]

end

module test_bukdu_plugs_contentparsers_urlencoded_multipart

using Test
using Bukdu
using .Bukdu.Deps
using .Bukdu.Plug.ContentParsers
using .Deps.HTTP: Multipart

struct FC <: ApplicationController; conn::Conn; end

function urlencoded(c::FC)
    c.params.q
end

function multipart(c::FC)
    c.params
end

routes() do
    post("/urlencoded", FC, urlencoded)
    post("/multipart", FC, multipart)
end

# https://discourse.julialang.org/t/http-multipart-form-data-processing-by-server/24076/3
req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=---------------------------182023285717490760841965583652"],
    codeunits("""
-----------------------------182023285717490760841965583652
Content-Disposition: form-data; name="image"; filename="file1.jpg"
Content-Type: image/jpeg

......JFIF.............C..........
-----------------------------182023285717490760841965583652
Content-Disposition: form-data; name="num"

2
-----------------------------182023285717490760841965583652--""")
)
got = Router.call(req).got
image = got["image"]
@test image isa Multipart
@test image.filename == "file1.jpg"
@test image.contenttype == "image/jpeg"
@test String(read(image.data)) == "......JFIF.............C.........."
@test got["num"] == "2"

req = Deps.HTTP.Request(
    "POST",
    "/urlencoded",
    ["Content-Type"=>"application/x-www-form-urlencoded", "Content-Length"=>"3"],
    codeunits("q=5")
)
@test Router.call(req).got == "5"

req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=boundary"],
    codeunits("""
--boundary
Content-Disposition: form-data; name="field1"

value1
--boundary
Content-Disposition: form-data; name="field2"; filename="example.txt"

value2
--boundary--""")
)
got = Router.call(req).got
@test got["field1"] == "value1"
field2 = got["field2"]
@test field2 isa Multipart
@test field2.filename == "example.txt"
@test field2.contenttype == ""
@test String(read(field2.data)) == "value2"

req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=boundary"],
    codeunits("""
--boundary
Content-Disposition: form-data; name="field1"

value1
--boundary
Content-Disposition: form-data; name="field2"; filename="example.txt"
Content-Type: text/plain

value2
--boundary--""")
)
got = Router.call(req).got
@test got["field1"] == "value1"
field2 = got["field2"]
@test field2 isa Multipart
@test field2.filename == "example.txt"
@test field2.contenttype == "text/plain"
@test String(read(field2.data)) == "value2"

req = Deps.HTTP.Request(
    "POST",
    "/multipart",
    ["Content-Type"=>"multipart/form-data; boundary=----WebKitFormBoundaryLeFS8sRvPzeszrWy"],
    codeunits("""
------WebKitFormBoundaryLeFS8sRvPzeszrWy\r\nContent-Disposition: form-data; name=\"user_name\"\r\n\r\n\r\n------WebKitFormBoundaryLeFS8sRvPzeszrWy\r\nContent-Disposition: form-data; name=\"user_famous\"\r\n\r\nfalse\r\n------WebKitFormBoundaryLeFS8sRvPzeszrWy--\r\n""")
)
@test Router.call(req).got == Assoc("user_name"=>"", "user_famous"=>"false")

Routing.empty!()

@test  ContentParsers.rstripcrlf([0x0d, 0x0a]) == []      # CR LF
@test  ContentParsers.rstripcrlf([0x0d])       == [0x0d]  # CR

end # module test_bukdu_plugs_contentparsers_urlencoded_multipart

ext: jl


text:
using egads

width = 5.0
rad   = 1.0

NULL_E = ego(0)

(context, status) = EG_open( )
if (status < 0) error("Can't open, failure code: %i", status) end

# build cylinder 1
datacyl = [0.0, 0.0,    0.0,
           0.0, 0.0,  width,
           rad]
(ebody1, status) = EG_makeSolidBody(context, CYLINDER, datacyl)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# build cylinder 2
datacyl = [0.0, 0.0, -width,
           0.0, 0.0, 0.0,
           rad]
(ebody2, status) = EG_makeSolidBody(context, CYLINDER, datacyl)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# join
toler = Cdouble(0.0)
(emodel, status) = EG_join( ebody1, ebody2, toler )
if (status != EGADS_SUCCESS) cleanup(status, context) end

# save model
status = EG_saveModel(emodel, "test.egads")
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_saveModel(emodel, "test.step")
if (status != EGADS_SUCCESS) cleanup(status, context) end

# delete
status = EG_deleteObject(ebody1)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(ebody2)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(emodel)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# close everything
cleanup( context )


text:
using Requests, ZipFile, DataFrames

datapath = normpath(Pkg.dir(), "Mortchartgen", "data")

dllist = ["documentation"; "availability"; "country_codes"; "notes"; "Pop";
	"morticd07"; "morticd08"; "morticd09"; "Morticd10_part1"; "Morticd10_part2";
	"WPP2017_PopulationByAgeSex_Medium"]

tablelist = ["MortIcd7"; "Morticd8"; "Morticd9"; "Morticd10_part1";
	"Morticd10_part2"; "pop"]

whosuff = ".zip"
wppsuff = ".csv"

function download_who(dlpath = ".", filelist = dllist)
	preurl_def = "http://www.who.int/entity/healthinfo/statistics/"
	preurl_pop = "http://www.who.int/entity/healthinfo/"
	preurl_wpp = "https://esa.un.org/unpd/wpp/DVD/Files/1_Indicators%20(Standard)/CSV_FILES/"
	urlparams = Dict("ua" => "1")
	wppfile = dllist[end]

	for filename in filelist
		if filename == "Pop"
			preurl = preurl_pop
			filesuff = whosuff
		elseif filename == wppfile
			preurl = preurl_wpp
			filesuff = wppsuff
		else
			preurl = preurl_def
			filesuff = whosuff
		end
		filefullname = *(filename, filesuff)
		dlfilepath = normpath(dlpath, filefullname)
		url = *(preurl, filefullname)
		print("$url\n$dlfilepath\n")
		Requests.save(Requests.get(url; query = urlparams), dlfilepath)
		if filesuff == ".zip"
			r = ZipFile.Reader(dlfilepath)
			for f in r.files
				write(normpath(dlpath, f.name), read(f))
			end
			close(r)
		end
	end
end

function table_import(dlpath = ".",
	tablecnf = normpath(datapath, "tableimp.cnf"), filelist = tablelist)
	for filename in tablelist
		dlfilepath = normpath(dlpath, filename)
		if filename == "pop"
			temppath = normpath(dlpath, "Pop")
		else
			temppath = normpath(dlpath, "Deaths")
		end
		mv(dlfilepath, temppath)
		run(`mysqlimport --defaults-extra-file=$tablecnf
			Morticd $temppath`)
		mv(temppath, dlfilepath)
	end
	cp(normpath(dlpath, "country_codes"), normpath(datapath, "country_codes"),
		remove_destination = true)
end

function wpp_convert(dlpath = ".", syr = 1990, eyr = 2016)
	whoages = map((x)->Symbol("Age$x"), [7:25;25])
	wppages = [map((x,y)->"$(x)-$(y)", 5:5:95, 9:5:99); "100+"]
	ages = DataFrame(variable = whoages, AgeGrp = wppages)
	country_codes = readtable(normpath(datapath, "country_codes"))
	rename!(country_codes, [:country; :name], [:Country; :Location])
	wpp = readtable(normpath(dlpath, *(dllist[end], wppsuff)))
	wpp = wpp[((wpp[:Time].>=syr) .& (wpp[:Time].<=eyr)), :]
	wpp = join(wpp, country_codes, on = :Location)
	wpp = join(wpp, ages, on = :AgeGrp)
	wppsexes = Dict()
	for (i, sexcol) in enumerate([:PopMale, :PopFemale])
		wppsexes[i] = DataFrame(variable = wpp[:variable], value = wpp[sexcol] .* 1000,
			Sex = i, Year = wpp[:Time], Country = wpp[:Country])
	end
	wpp = aggregate(vcat(values(wppsexes)...), [:variable; :Year; :Sex; :Country], sum)
	wpp = DataFrame(variable = wpp[:variable], value = wpp[:value_sum], Sex = wpp[:Sex],
		Year = wpp[:Year], Country = wpp[:Country])
	sort!(wpp, cols = [order(:variable, by = (x)->parse("$x"[4:end])); :Country; :Year; :Sex])
	writetable(normpath(datapath, "wppconverted.csv"), wpp)
end


text:
struct NN
    g1
    g2
end

function NN(wordembeds::Matrix{T}, posembeds::Matrix{T}) where T
    w = lookup(Node(zerograd(wordembeds)), Node(name="w"))
    p = lookup(Node(zerograd(posembeds)), Node(name="p"))

    h = concat(1, w, p)
    #h = dropout(h, 0.33)
    d = 300
    batchdims_w = Node(name="batchdims_w")
    h = LSTM(T,150,d,1,0.0,true)(h,batchdims_w)
    # h = dropout(h, 0.33)
    g1 = Graph(h)

    h = Node(name="h")
    hdep = Linear(T, 2d, 300, init_w=OrthoNormal())(h)
    hdep = relu(hdep)
    #hdep = dropout(hdep, 0.33)
    hhead = Linear(T, 2d, 300, init_w=OrthoNormal())(h)
    hhead = relu(hhead)
    #hhead = dropout(hhead, 0.33)

    h = Linear(T,300,300)(hdep)
    h = BLAS.gemm('T', 'N', 1, hhead, h)
    g2 = Graph(h)
    NN(BACKEND(g1), BACKEND(g2))
end

function (nn::NN)(x::Sample, istrain::Bool)
    Merlin.CONFIG.train = istrain

    #=
    if istrain
        p = CPUBackend()(x.p)
        w = CPUBackend()(x.w)
        h = nn.gg1(x.batchdims_w, p, w)
        hs = split(h, x.batchdims_w)
        os = Var[]
        for h in hs
            push!(os, nn.gg2(h))
        end

        heads = map(h -> CPUBackend()(h), x.heads)
        ys = Var[]
        for (h,o) in zip(heads,os)
            y = softmax_crossentropy(h, o)
            push!(ys, y)
        end
        cpu_y = concat(1, ys).data
    end
    =#

    h = nn.g1(x.batchdims_w, x.p, x.w)
    hs = split(h, x.batchdims_w)
    os = Var[]
    for h in hs
        push!(os, nn.g2(h))
    end
    if istrain
        ys = Var[]
        for (h,o) in zip(x.heads,os)
            y = softmax_crossentropy(h, o)
            # y = o
            push!(ys, y)
        end
        # ys[1]
        concat(1, ys)
    else
        cat(1, map(o -> vec(argmax(Array(o.data),1)), os)...)
    end
end


text:
# GPU cross-correlation using files split into Npersplit windows per file
# run splitseis.jl with Npersplit set to 236 before running this
# Note: this will only work if an Nvidia GPU is available
using SeisIO, SeisNoise, Dates, Glob, BenchmarkTools, Printf

Npersplit = 236
cc_len, cc_step = 2^15, 26768
freqmin = 0.003
freqmax = 0.005
maxlag = 12000.
splits = glob("*",joinpath(expanduser("~/SeisNoise_paper/DATA/SPLITSEIS"),string(Npersplit)))
SPLITCORR = joinpath(expanduser("~/SeisNoise_paper/CORR/SPLITSEIS"),string(Npersplit))
XML = expanduser("~/SeisNoise_paper/DATA/XML")

if !isdir(SPLITCORR)
    mkpath(SPLITCORR)
end

# this function reads a file, applies pre-processing, converts to RawData
#  and then returns an FFTData
function seis2fft(file::String,XML::String,cc_len::Int,cc_step::Int,
                  freqmin::Float64,freqmax::Float64)
    S = SeisData()
    read_data!(S,"sac",file)
    detrend!(S)
    taper!(S)
    filtfilt!(S,fl=0.001,fh=0.4,np=3)

    # get xmlfile
    net,sta,loc,chan = split(S[1].id,'.')
    xmlfile = joinpath(XML,net*'.'*sta*".xml")
    Resp = read_meta("sxml",xmlfile)
    translate_resp!(S,Resp)
    remove_resp!(S)
    R = RawData(S[1],cc_len,cc_step) |> gpu
    demean!(R)
    taper!(R)
    bandpass!(R,freqmin,freqmax,corners=3)
    return rfft(R)
end

# this function cross-correlates data for a "split"
function GPUALL(splits,cc_len,cc_step,freqmin,freqmax,maxlag,OUTDIR,XML)
    Nsplits = length(splits)
    for ii = 1:Nsplits
        println("Cross-correlating chunk $ii of $Nsplits ",now())
        files = glob("*",splits[ii])
        Nfiles = length(files)

        # I/O + fft
        Fs = Array{FFTData}(undef,Nfiles)
        for jj = 1:Nfiles
            Fs[jj] = seis2fft(files[jj],XML,cc_len,cc_step,freqmin,freqmax)
        end

        # cross-correlate
        for jj = 1:Nfiles -1
            for kk = jj+1:Nfiles
                C = compute_cc(Fs[jj],Fs[kk],maxlag)
                stack!(C,allstack=true)
                C = cpu(C)
                finalize(C.corr)
                save_corr(C,OUTDIR)
            end
            finalize(Fs[jj].fft)
        end
    end
    return nothing
end

# warm up the GPU
GPUALL(splits[1:1],cc_len,cc_step,freqmin,freqmax,maxlag,SPLITCORR,XML)
toremove = glob("*",SPLITCORR)
rm.(toremove)

tGPU = @elapsed begin
    GPUALL(splits,cc_len,cc_step,freqmin,freqmax,maxlag,SPLITCORR,XML)
end
println("GPU XCORR took $(@sprintf("%.2f",tGPU)) seconds")


text:
using .AxisKeys

"""

# convert_a3d

# Convert the output file(s) created by cmdstan to a KeyedArray.

$(SIGNATURES)

"""
function convert_a3d(a3d_array, cnames, ::Val{:keyedarray})
    psymbols= Symbol.(cnames)
    pa = permutedims(a3d_array, [1, 3, 2])
    wrapdims(pa,
        iteration=1:size(pa, 1),
        chain=1:size(pa, 2),
        param=psymbols
    )
end

function matrix(ka::KeyedArray, sym::Union{Symbol, String})
    n = string.(axiskeys(ka, :param))
    syms = string(sym)
    sel = String[]
    for (i, s) in enumerate(n)
        if length(s) > length(syms) && syms == n[i][1:length(syms)] &&
            n[i][length(syms)+1] in ['[', '.', '_']
            append!(sel, [n[i]])
        end
    end
    length(sel) == 0 && error("$syms not in $n")
    ka(param=Symbol.(sel))
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
export jld22sac
#
using Dates, SeisIO, JLD2, SeisNoise, Sockets, Logging
"""
    Converting JLD2 format to sac files. There are options to read in JLD2 data
    in different types: SeisData (defined in SeisIO) and CorrData (defined in SeisNoise).
"""

#base function to output usage information and format descriptions.
function jld22sac()
    print("
        Converting JLD2 format to SAC files.
        jld22sac(jldfile,sacrootdir; informat, outformat,subset, verbose)
        Requred  arguments:
            jldfile
            sacrootdir
        Optional:
            informat = infile (jld file) structure code (string).
                TD: JLD2 FILE -> Timestamp group -> SeisData/CorrData; Type of data is detected automatically.
                CTD: JLD2 FILE -> Component -> Timestamp group -> SeisData/CorrData;
            outformat = outfile (sac file) structre code (string)
                SRC: Source -> Receiver -> Component -> SAC (timestamps and component information are
                    embeded in file names) #this is for CorrData only. [DEFAULT]
                    Specifically, this is good/designed for poststack CorrData, e.g.,
                    1 stack per day.
                TG: Timestamp group -> SAC   #good for SeisData. if for CorrData, the
                    station pair and component will be embeded in file names.
            subset = a Dictionary containing subset condtions (lists). DEFAULT is empty
                Current support subset keys: timestamplist, stationlist (in net.sta format), componentlist.
            verbose = true/false. [DEFAULT] is false.

        NOTES:
        1. For CorrData, the freqmin, freqmax, cc_len, cc_step are stored in SAC headers
            USER0, USER1, USER2, USER3, respectively
        2. For CorrData, the maxlag is used to determine the start (B) and end (E) times

        \n")

    return
end

#
function jld22sac(jldfile::String,sacrootdir::String; informat::String="TD", outformat::String="SRC",
    subset::Dict=Dict(), verbose::Bool=false)

    #first check JLD formatcode and sac output format/structrue code
    if informat != "TD" && informat != "CTD"
        error("The input JLDFORMATCODE must be one of: TD (default) or CTD. Run jld22sac() for detailed explanations.\n")
    end
    if outformat != "SRC" && outformat != "TG"
        error("The input JLDFORMATCODE must be one of: SRC (default), or TG. Run jld22sac() for detailed explanations.\n")
    end

    #checking subset options
    subsetflag_ts = false
    subsetflag_sta = false
    subsetflag_comp = false
    subsetkey_ts = "timestamplist"
    subsetkey_sta = "stationlist" #consider both source and receiver in CorrData situation.
    subsetkey_comp = "componentlist" #component list
    if !isempty(subset)
        subsetkeys = collect(keys(subset))
        if verbose
            println("Applying subset. Will only save subset data for: ")
            println(subsetkeys)
        end

        if in(subsetkey_ts,subsetkeys)
            subsetflag_ts = true
        end
        if in(subsetkey_sta,subsetkeys)
            subsetflag_sta = true
        end
        if in(subsetkey_comp,subsetkeys)
            subsetflag_comp = true
        end
    end

    #create SAC data directory
    mkpath(sacrootdir)

    jfile = jldopen(jldfile,"r")

    #get first level group tags
    group1 = keys(jfile)
    if subsetflag_ts
        group1 = subset[subsetkey_ts] #use subset timestamplist
    end

    for g1 = group1
        if g1 != "info" && g1 != "errors"
            if informat == "TD" #only one level group
                println("Working on group: ",g1)
                dlist = keys(jfile[g1]) #data list
                for dfile = dlist
                    d = jfile[joinpath(g1,dfile)]
                    # println("CorrData to SAC")
                    # println(typeof(d))
                    if typeof(d) == CorrData # data file has CorrData type
                        stemp = split(d.name,".")
                        srname = join([stemp[1],stemp[2]],".")  #use the first entry in location as the source name
                        rcname = join([stemp[5],stemp[6]],".")
                        comp = d.comp
                        S = corr2seis(d)
                    elseif typeof(d) == SeisData || typeof(d) == SeisChannel # data file has SeisData type
                        # println(d.id)
                        stemp = split(d.id[1],".")
                        srname = ""
                        rcname = join([stemp[1],stemp[2]],".")
                        comp = stemp[4]
                        S = SeisData(d)
                        if S.n > 1
                            @warn "$dfile has more than 1 channels, only headers from the 1st channel are used in convering!"
                        end
                    else
                        error(typeof(d), " datatype not supported.")
                    end
                    if outformat == "SRC" #Source -> Receiver -> SAC
                        if subsetflag_sta
                            if in(srname,subset[subsetkey_sta]) || in(rcname,subset[subsetkey_sta])
                                if subsetflag_comp
                                    if in(comp,subset[subsetkey_comp])
                                        if verbose
                                            println("  --> converting: ",dfile)
                                        end
                                        writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,srname,rcname,comp)))
                                    end
                                else
                                    if verbose
                                        println("  --> converting: ",dfile)
                                    end
                                    writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,srname,rcname,comp)))
                                end
                            end
                        else
                            if subsetflag_comp
                                if in(comp,subset[subsetkey_comp])
                                    if verbose
                                        println("  --> converting: ",dfile)
                                    end
                                    writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,srname,rcname,comp)))
                                end
                            else
                                if verbose
                                    println("  --> converting: ",dfile)
                                end
                                writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,srname,rcname,comp)))
                            end
                        end
                    elseif outformat == "TG" #Timestamp -> SAC
                        if subsetflag_sta
                            if in(srname,subset[subsetkey_sta]) || in(rcname,subset[subsetkey_sta])
                                if subsetflag_comp
                                    if in(comp,subset[subsetkey_comp])
                                        if verbose
                                            println("  --> converting: ",dfile)
                                        end
                                        writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,g1)))
                                    end
                                else
                                    if verbose
                                        println("  --> converting: ",dfile)
                                    end
                                    writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,g1)))
                                end
                            end
                        else
                            if subsetflag_comp
                                if in(comp,subset[subsetkey_comp])
                                    if verbose
                                        println("  --> converting: ",dfile)
                                    end
                                    writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,g1)))
                                end
                            else
                                if verbose
                                    println("  --> converting: ",dfile)
                                end
                                writesac_rich(S,datatype=typeof(d),outdir=string(joinpath(sacrootdir,g1)))
                            end
                        end
                    end
                end
            else #need to do anotehr level loop for format CTSD or CTCD
                error("CTSD and CTCD are not ready yet. Coming soon.")
            end
        end
    end #end of loop for the first level group tags

    JLD2.close(jfile)
end


text:
TensorIndex(s::TensorIndex) = s
TensorIndex(name::AbstractString,tensor_index_type::TensorIndexType,is_up::Bool=true)::TensorIndex = tensor.TensorIndex(name,tensor_index_type,is_up)
TensorIndex(name::Symbol,tensor_index_type::TensorIndexType,is_up::Bool=true)::TensorIndex = tensor.TensorIndex(name,tensor_index_type,is_up)
-( x::TensorIndex ) = TensorIndex(x.__neg__())

TensorIndexType(name::AbstractString,dummy_name::AbstractString,dim::Int,metric_name::AbstractString)::TensorIndexType = tensor.TensorIndexType(name,dummy_name=dummy_name,dim=dim,metric_name=metric_name)
TensorIndexType(name::AbstractString,dummy_name::AbstractString,dim::Int)::TensorIndexType = tensor.TensorIndexType(name,dummy_name=dummy_name,dim=dim)
TensorIndexType(name::AbstractString,dummy_name::AbstractString)::TensorIndexType = tensor.TensorIndexType(name,dummy_name=dummy_name)
TensorIndexType(name::Symbol,dummy_name::Symbol,dim::Int,metric_name::Symbol)::TensorIndexType = tensor.TensorIndexType(string(name),string(dummy_name),dim,string(metric_name))

function Base.getproperty(A::IndexedTensor, k::Symbol)
    if k == :head
        m = getproperty(PyCall.PyObject(A),k)
        return convert(TensorHead,m)
    #elseif k == :indices
    #    m = getproperty(PyCall.PyObject(A),:indices)
    #    println(m)
    #    return convert(Array{TensorIndex},m)
    elseif k == :free
        m = getproperty(PyCall.PyObject(A),k)
        return convert(Array{Tuple{TensorIndex,Int}},m)
    elseif k == :free_indices
        m = pycall(getproperty(PyCall.PyObject(A),:get_free_indices),Array{TensorIndex})
        return m
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    else
        M1 = getproperty(PyCall.PyObject(A), k)
        M1
    end
end

function Base.getproperty(A::TensorIndex, k::Symbol)
    if k == :is_up
        m = getproperty(PyCall.PyObject(A),k)
        return convert(Bool,m)
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    else
        M1 = getproperty(PyCall.PyObject(A), k)
        M1
    end
end
function Base.getproperty(A::T, k::Symbol) where {T <: Tensor}
    if k == :free
        m = getproperty(PyCall.PyObject(A),k)
        return convert(Array{Tuple{TensorIndex,Int}},m)
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    else
        M1 = getproperty(PyCall.PyObject(A), k)
        M1
    # else error?
    end
end

function Base.getproperty(A::TensorIndexType, k::Symbol)
    if k == :delta || k == :metric
        m = getproperty(PyCall.PyObject(A),k)
        return convert(TensorHead,m)
    elseif k == :set_metric
        m = getproperty(PyCall.PyObject(A),k)
        return m
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    else
        M1 = getproperty(PyCall.PyObject(A), k)
        M1
    end
end

macro indices(xs...)
    ret = Expr(:block)
    TIT = xs[1]
    for x in xs[2:end]
        push!(ret.args,:($(esc(x)) = TensorIndex($(string(x)),$(esc(TIT)))))
    end
    return ret
end


text:
# ---------   Utility Functions ----------- #
getvarid(s::Symbol) = string(":($s)")

getvarid(e::Expr)   = string("$(getvarid(e.args[1])), $(e.args[2])")

macro getvarid(e)
  # usage: @getvarid x[2][1+5], will return a tuple like (:x, 2, 6)
  return parse( getvarid(e) )
end

dot(x) = dot(x, x)

invlogit(x) = 1.0 ./ (exp(-x) + 1.0)

logit(x) = log(x ./ (1.0 - x))

function randcat(p::Vector{Float64}) # More stable, faster version of rand(Categorical)
  # if(any(p .< 0)) error("Negative probabilities not allowed"); end
  r, s = rand(), 1.0
  for j = 1:length(p)
    r -= p[j]
    if(r <= 0.0) s = j; break; end
  end
  return s
end

type NotImplementedException <: Exception end

# Numerically stable sum of values represented in log domain.
function logsum(xs :: Vector{Float64})
  largest = maximum(xs)
  ys = map(x -> exp(x - largest), xs)
  result = log(sum(ys)) + largest
  return result
end

using Distributions
# KL-divergence
function kl(p::Normal, q::Normal)
  return (log(q.σ / p.σ) + (p.σ^2 + (p.μ - q.μ)^2) / (2 * q.σ^2) - 0.5)
end


text:
const CSV_FIELDS = [
        :secao
        :item
        :id_posicao
        :fator_risco
        :local_registro
        :carteira_negoc
        :v01
        :v02
        :v03
        :v04
        :v05
        :v06
        :v07
        :v08
        :v09
        :v10
        :v11
        :v12
        :v12_MaM
]

const CSV_DELIMITER = ";"
const NA_STRING = ""
const DECIMAL_DEMILITER = ","

function write_csv(io::IO, doc::Documento)

    # struct Documento
    #     id_docto::String
    #     id_docto_versao::String
    #     data_base::String
    #     id_inst_financ::Int64
    #     tipo_arq::Symbol
    #     nome_contato::String
    #     fone_contato::String
    #     ativo::Vector{ItemCarteira}
    #     passivo::Vector{ItemCarteira}
    #     derivativo::Vector{ItemCarteira}
    #     ativo_fundo::Vector{ItemCarteira}
    #     atividade_financeira::Vector{ItemCarteira}

    _write_csv_header(io)

    for item_carteira in sorted_keys(doc.ativo)
        _write_csv_item_carteira(io, item_carteira, doc.ativo[item_carteira], "ativo")
    end
    for item_carteira in sorted_keys(doc.passivo)
        _write_csv_item_carteira(io, item_carteira, doc.passivo[item_carteira], "passivo")
    end
    for item_carteira in sorted_keys(doc.derivativo)
        _write_csv_item_carteira(io, item_carteira, doc.derivativo[item_carteira], "derivativo")
    end
    for item_carteira in sorted_keys(doc.ativo_fundo)
        _write_csv_item_carteira(io, item_carteira, doc.ativo_fundo[item_carteira], "ativo_fundo")
    end
    for item_carteira in sorted_keys(doc.passivo_fundo)
        _write_csv_item_carteira(io, item_carteira, doc.passivo_fundo[item_carteira], "passivo_fundo")
    end
    for item_carteira in sorted_keys(doc.derivativo_fundo)
        _write_csv_item_carteira(io, item_carteira, doc.derivativo_fundo[item_carteira], "derivativo_fundo")
    end
    for item_carteira in sorted_keys(doc.atividade_financeira)
        _write_csv_item_carteira(io, item_carteira, doc.atividade_financeira[item_carteira], "atividade_financeira")
    end

end

function _write_csv_header(io::IO)
    println(io, join(CSV_FIELDS, CSV_DELIMITER))
end

function _write_csv_item_carteira(io::IO, item::ItemCarteira, fluxos::Fluxos, secao::String)
    csv_values = Vector{String}([
        secao                                  # :secao
        to_string(item.item)                   # :item
        to_string(item.id_posicao)             # :id_posicao
        to_string(item.fator_risco)            # :fator_risco
        to_string(item.local_registro)         # :local_registro
        to_string(item.carteira_negoc)         # :carteira_negoc
        to_string(get_valor_alocado(fluxos,  1)) # :v01
        to_string(get_valor_alocado(fluxos,  2)) # :v02
        to_string(get_valor_alocado(fluxos,  3)) # :v03
        to_string(get_valor_alocado(fluxos,  4)) # :v04
        to_string(get_valor_alocado(fluxos,  5)) # :v05
        to_string(get_valor_alocado(fluxos,  6)) # :v06
        to_string(get_valor_alocado(fluxos,  7)) # :v07
        to_string(get_valor_alocado(fluxos,  8)) # :v08
        to_string(get_valor_alocado(fluxos,  9)) # :v09
        to_string(get_valor_alocado(fluxos, 10)) # :v10
        to_string(get_valor_alocado(fluxos, 11)) # :v11
        to_string(get_valor_alocado(fluxos, 12)) # :v12
        to_string(get_valor_mam(fluxos, 12))     # :v12_MaM
    ])
    println(io, join(csv_values, CSV_DELIMITER))
end

to_string(x::Symbol) = "$x"
to_string(x::Nothing) = NA_STRING
function to_string(x::Float64)
    out = "$x"
    if DECIMAL_DEMILITER != "."
        out = replace(out, "." => DECIMAL_DEMILITER)
    end
    return out
end

"Retorna o valor alocado no vertice de um ItemCarteira. Retorna 0.0 caso não haja valor alocado."
function get_valor_alocado(fluxos::Fluxos, codigo_vertice::Symbol)
    if haskey(fluxos.vertices, codigo_vertice)
        return fluxos.vertices[codigo_vertice].valor_alocado
    else
        return 0.0
    end
end
get_valor_alocado(fluxos::Fluxos, i::Int64) = get_valor_alocado(fluxos, get_codigo_vertice(i))

"Retorna o valor MaM no vertice de um ItemCarteira. Retorna 0.0 caso não haja valor MaM."
function get_valor_mam(fluxos::Fluxos, codigo_vertice::Symbol)
    if haskey(fluxos.vertices, codigo_vertice)
        return fluxos.vertices[codigo_vertice].valor_mam
    else
        return 0.0
    end
end
get_valor_mam(fluxos::Fluxos, i::Int64) = get_valor_mam(fluxos, get_codigo_vertice(i))

write_csv(path::String, doc::Documento) = open(path, "w+") do io write_csv(io, doc) end


text:
struct LL{V}
    v::V
    i::Int
end

islist(::Any) = false
islist(l::Union{LL, Term, AbstractArray, Tuple}) = true

Base.empty(l::LL) = empty(l.v)
Base.isempty(l::LL) = l.i > length(l.v)

Base.length(l::LL) = length(l.v)-l.i+1
@inline car(l::LL) = l.v[l.i]
@inline cdr(l::LL) = isempty(l) ? empty(l) : LL(l.v, l.i+1)

Base.length(t::Term) = length(arguments(t)) + 1 # PIRACY
Base.isempty(t::Term) = false
@inline car(t::Term) = operation(t)
@inline cdr(t::Term) = arguments(t)

@inline car(v) = first(v)
@inline cdr(v) = isempty(v) ? empty(l) : LL(v, 2)

@inline take_n(ll::LL, n) = isempty(ll) || n == 0 ? empty(ll) : @views ll.v[ll.i:n+ll.i-1] # @views handles Tuple
@inline take_n(ll, n) = @views ll[1:n]

drop_n(ll, n) = n === 0 ? ll : drop_n(cdr(ll), n-1)

@inline drop_n(ll::Term, n) = drop_n(arguments(ll), n-1)
@inline drop_n(ll::AbstractArray, n) = drop_n(LL(ll, 1), n)
@inline drop_n(ll::LL, n) = LL(ll.v, ll.i+n)

# Fastest possible dict
struct MatchDict{K,V}
    k::K
    v::V
end
function matchdict(names)
    MatchDict((; map(=>, names, 1:length(names))...),
              ((Ref{Any}(nothing) for n in names)...,))
end
@inline _copy(r) = Ref{Any}(r[])
assoc(d::MatchDict, k, v) = MatchDict(d.k, (d1 = map(_copy,d.v); d1[d.k[k]][] = v; d1))
#Base.haskey(d::MatchDict, k) = d.v[d.k[k]] !== nothing
Base.getindex(d::MatchDict, k) = d.v[d.k[k]][]


text:
using MPIMeasurements

scanner = MPIScanner("MPS.toml")
daq = getDAQ(scanner)

params = toDict(daq.params)
params["studyName"]="TestTobi"
params["studyDescription"]="A very cool measurement"
params["scannerOperator"]="Tobi"
params["dfStrength"]=[20e-3]
params["acqNumAverages"]=100

@info "PULL OUT THE SAMPLE!"
readline(stdin)
uBG = measurement(daq, params, controlPhase=true)

@info "PUT THE SAMPLE IN!"
readline(stdin)
# This version does not store the data
u = measurement(daq, params, controlPhase=true)

# This version does store the data in a custom location
#filename = measurement(daq,"/home/labuser/test.mdf", params, controlPhase=true)

# This version does store the data in the MDFStore
#filename = measurement(daq, MDFStore, params,  controlPhase=true)

#u = loadBGCorrData(filename)

showDAQData(daq,u.-uBG)

disconnect(daq)


text:
module LibCImGui

using CImGuiPack_jll
export CImGuiPack_jll

const IS_LIBC_MUSL = occursin("musl", Base.BUILD_TRIPLET)
if Sys.isapple() && Sys.ARCH === :aarch64
    include("../lib/aarch64-apple-darwin20.jl")
elseif Sys.islinux() && Sys.ARCH === :aarch64 && !IS_LIBC_MUSL
    include("../lib/aarch64-linux-gnu.jl")
elseif Sys.islinux() && Sys.ARCH === :aarch64 && IS_LIBC_MUSL
    include("../lib/aarch64-linux-musl.jl")
elseif Sys.islinux() && startswith(string(Sys.ARCH), "arm") && !IS_LIBC_MUSL
    include("../lib/armv7l-linux-gnueabihf.jl")
elseif Sys.islinux() && startswith(string(Sys.ARCH), "arm") && IS_LIBC_MUSL
    include("../lib/armv7l-linux-musleabihf.jl")
elseif Sys.islinux() && Sys.ARCH === :i686 && !IS_LIBC_MUSL
    include("../lib/i686-linux-gnu.jl")
elseif Sys.islinux() && Sys.ARCH === :i686 && IS_LIBC_MUSL
    include("../lib/i686-linux-musl.jl")
elseif Sys.iswindows() && Sys.ARCH === :i686
    include("../lib/i686-w64-mingw32.jl")
elseif Sys.islinux() && Sys.ARCH === :powerpc64le
    include("../lib/powerpc64le-linux-gnu.jl")
elseif Sys.isapple() && Sys.ARCH === :x86_64
    include("../lib/x86_64-apple-darwin14.jl")
elseif Sys.islinux() && Sys.ARCH === :x86_64 && !IS_LIBC_MUSL
    include("../lib/x86_64-linux-gnu.jl")
elseif Sys.islinux() && Sys.ARCH === :x86_64 && IS_LIBC_MUSL
    include("../lib/x86_64-linux-musl.jl")
elseif Sys.isbsd() && !Sys.isapple()
    include("../lib/x86_64-unknown-freebsd11.1.jl")
elseif Sys.iswindows() && Sys.ARCH === :x86_64
    include("../lib/x86_64-w64-mingw32.jl")
else
    error("Unknown platform: $(Base.BUILD_TRIPLET)")
end

# should be reimplement as macros by using Julia 1.5's @ccall macro
function igText(text)
    ccall((:igText, libcimgui), Cvoid, (Cstring, Cstring), "%s", text)
end

function igTextColored(col, text)
    ccall((:igTextColored, libcimgui), Cvoid, (ImVec4, Cstring), col, text)
end

function igTextDisabled(text)
    ccall((:igTextDisabled, libcimgui), Cvoid, (Cstring,), text)
end

function igTextWrapped(text)
    ccall((:igTextWrapped, libcimgui), Cvoid, (Cstring,), text)
end

function igLabelText(label, text)
    ccall((:igLabelText, libcimgui), Cvoid, (Cstring, Cstring), label, text)
end

function igBulletText(text)
    ccall((:igBulletText, libcimgui), Cvoid, (Cstring,), text)
end

function igTreeNodeStrStr(str_id, text)
    ccall((:igTreeNodeStrStr, libcimgui), Bool, (Cstring, Cstring), str_id, text)
end

function igTreeNodePtr(ptr_id, text)
    ccall((:igTreeNodePtr, libcimgui), Bool, (Ptr{Cvoid}, Cstring), ptr_id, text)
end

function igTreeNodeExStrStr(str_id, flags, text)
    ccall((:igTreeNodeExStrStr, libcimgui), Bool, (Cstring, ImGuiTreeNodeFlags, Cstring), str_id, flags, text)
end

function igTreeNodeExPtr(ptr_id, flags, text)
    ccall((:igTreeNodeExPtr, libcimgui), Bool, (Ptr{Cvoid}, ImGuiTreeNodeFlags, Cstring), ptr_id, flags, text)
end

function igSetTooltip(text)
    ccall((:igSetTooltip, libcimgui), Cvoid, (Cstring,), text)
end

function igOpenPopup(str_id)
    ccall((:igOpenPopup, libcimgui), Cvoid, (Cstring,), str_id)
end

function igLogText(text)
    ccall((:igLogText, libcimgui), Cvoid, (Cstring,), text)
end

# exports
const PREFIXES = ["ig", "Im", "IMGUI_", "imnodes_", "ImPlot_", "ImVector_"]
for name in names(@__MODULE__; all=true), prefix in PREFIXES
    if startswith(string(name), prefix)
        @eval export $name
    end
end

end


text:
module InsertState

using Mixtape
import Mixtape: CompilationContext, transform, allow
using MacroTools
using CodeInfoTools
using BenchmarkTools

# This shows an example where we recursively modify method calls to record state. Very monadic, dare I say :)

module Target

function baz(y::Float64)
    return y + 40.0
end

function foo(x::Float64)
    return baz(x + 20.0)
end

end

struct MyMix <: CompilationContext end

allow(ctx::MyMix, m::Module, fn, args...) = m == Target
debug(ctx::MyMix) = false

mutable struct Recorder
    d::Dict
    ret::Any
    Recorder() = new(Dict(), nothing)
    Recorder(d, ret) = new(d, ret)
end

swap(r, e) = e
function swap(r, e::Expr)
    e.head == :call || return e
    return Expr(:call, r, e.args[1:end]...)
end

function transform(::MyMix, src)
    b = CodeInfoTools.Builder(src)
    q = push!(b, Expr(:call, Recorder))
    for (v, st) in b
        b[v] = swap(q, st)
    end
    return CodeInfoTools.finish(b)
end

function (r::Recorder)(f::Function, args...)
    args = map(a -> a isa Recorder ? a.ret : a, args)
    rec = call(f, args...; ctx = MyMix())
    if rec isa Recorder
        merge!(r.d, rec.d)
        r.d[(f, args...)] = rec.ret
        r.ret = rec.ret
    else
        r.d[(f, args...)] = rec
        r.ret = rec
    end
    return r
end

Mixtape.@load_abi()
display(call(Target.foo, 5.0; ctx = MyMix()))
@btime call(Target.foo, 5.0; ctx = MyMix())

end # module


text:
mutable struct OptimizerResult{XT, YT, VT <: AbstractArray{XT}, OST <: OptimizerStatus{XT,YT}}
    status::OST   # iteration number, residuals and convergence info

    x::VT    # current solution
    f::YT    # current function
    g::VT    # current gradient
end

function OptimizerResult(x::VT, y::YT) where {XT, YT, VT <: AbstractVector{XT}}
    status = OptimizerStatus{XT,YT}()
    result = OptimizerResult{XT,YT,VT,typeof(status)}(status, zero(x), zero(y), zero(x))
    clear!(result)
end

status(result::OptimizerResult) = result.status

solution(result::OptimizerResult) = result.x
minimizer(result::OptimizerResult) = result.x
Base.minimum(result::OptimizerResult) = result.f

function clear!(result::OptimizerResult{XT,YT}) where {XT,YT}
    clear!(result.status)

    result.x .= XT(NaN)
    result.f  = YT(NaN)
    result.g .= XT(NaN)

    return result
end

function residual!(result::OptimizerResult, x, f, g)
    status = result.status

    status.rxₐ = sqeuclidean(x, result.x)
    status.rxᵣ = status.rxₐ / norm(x)

    status.rfₐ = norm(f - result.f)
    status.rfᵣ = status.rfₐ / norm(f)

    status.rgₐ = norm(g - result.g)
    status.rg  = norm(g)

    status.Δf  = f - result.f
    status.Δf̃ = result.g ⋅ x - result.g ⋅ result.x

    status.f_increased = abs(f) > abs(result.f)

    status.x_isnan = any(isnan, x)
    status.f_isnan = any(isnan, f)
    status.g_isnan = any(isnan, g)

    return status
end

function update!(result::OptimizerResult, x, f, g)
    residual!(result, x, f, g)

    result.x .= x
    result.f  = f
    result.g .= g

    return result
end

function initialize!(result::OptimizerResult, x, f, g)
    clear!(result)
    update!(result, x, f, g)
end

function next_iteration!(result::OptimizerResult)
    increase_iteration_number!(status(result))
end


text:
struct PhoneTrackingParams
  isusingapp::BitVector
  detection_prob::Float64
  detection_delay::Float64
end

PhoneTrackingParams(N::Integer, detection_prob::Real=1.0, detection_delay::Real=0.5) =
  PhoneTrackingParams(BitVector(undef, N), detection_prob, detection_delay)

function resample!(rng::AbstractRNG, params::PhoneTrackingParams, prob::Real)
  for i in 1:length(params.isusingapp)
    params.isusingapp[i] = rand(rng) < prob
  end
  params
end

function resamplebyhouseholds!(rng::AbstractRNG, params::PhoneTrackingParams, prob::Real, household_ptrs::Vector{Tuple{Ti,Ti}} where Ti<:Integer)

  num_individuals = length(household_ptrs)
  num_users = rand(rng, Binomial(num_individuals, prob)) # makes the number same as in the independent case

  randomized_households = shuffle!(rng, unique(household_ptrs)) # not optimal but gets the job done

  fill!(params.isusingapp, false)
  num_assigned = 0
  for (household_start, household_end) in randomized_households
    if num_assigned >= num_users
      break
    end

    household = UnitRange(household_start, household_end)

    assignements_in_household = min(length(household), num_users-num_assigned)
    for member in 1:assignements_in_household
      params.isusingapp[household[member]] = true
    end
    num_assigned += assignements_in_household
  end
end

function PhoneTrackingParams(rng::AbstractRNG,
                             N::Integer,
                             usage::Real,
                             detection_delay::Real=0.25,
                             detection_prob::Real=1.0,
                             household_ptrs::Union{Nothing,Vector{Tuple{Ti,Ti}} where Ti <: Integer} = nothing)
  params = PhoneTrackingParams(N, detection_prob, detection_delay)
  if household_ptrs === nothing
    resample!(rng, params, usage)
  else
    resamplebyhouseholds!(rng, params, usage, household_ptrs)
  end
  params
end

uses_phone_tracking(params::PhoneTrackingParams, person_id::Integer) = params.isusingapp[person_id]

function saveparams(dict, p::PhoneTrackingParams, prefix::AbstractString="")
  dict[prefix*"isusingapp"] = p.isusingapp
  dict[prefix*"detection_prob"] = p.detection_prob
  dict[prefix*"detection_delay"] = p.detection_delay
end


text:
input = (split.(readlines("day3.txt"),""))
input_mat = hcat(input...) .|> x -> parse(Int64,x)

function get_rate(rate)
    keep = zeros(Int64,size(input_mat)[1])

    for i in 1:size(input_mat)[1]
        # println(i)
        dd = [[j,count(x->x==j,input_mat[i,:])] for j in u]
        # println(dd)
        ddd = hcat(dd...)
        if rate == "gamma"
            keep[i] = ddd[1,argmax(ddd[2,:])]
        elseif rate == "epsilon"
            keep[i] = ddd[1,argmin(ddd[2,:])]
        end
    end
    return "0b"*join(keep) |> x -> parse(Int64,x)
end

get_rate("gamma")*get_rate("epsilon")


text:
@testset "Neuron value tests" begin

    ml(l, lfun=LazyMutable; insize=nin(l)[]) = mutable(l, inputvertex("in", insize, layertype(l)), layerfun = lfun)

    function tr(l, data)
        outshape = collect(size(data))
        outshape[[actdim(ndims(data))]] .= nout(l)
        example = [(param(data), ones(Float32, outshape...))];
        Flux.train!((x,y) -> Flux.mse(l(x), y), params(l), example, Descent(0.1))
    end

    function tr(l, output, inputs...)
        example = [(param.(inputs), output)];
        Flux.train!((x,y) -> Flux.mse(l(x...), y), params(l), example, Descent(0.1))
    end

    @testset "Ewma" begin
        import NaiveNASflux:agg
        m = Ewma(0.3)
        @test agg(m, missing, [1,2,3,4]) == [1,2,3,4]
        @test agg(m, [1,2,3,4], [5,6,7,8]) ≈ [3.8, 4.8, 5.8, 6.8]
    end

    @testset "Neuron value Dense default" begin
        l = ml(Dense(3,5))
        @test size(neuron_value(l)) == (5,)
    end

    @testset "Neuron value RNN default" begin
        l = ml(RNN(3,5))
        @test size(neuron_value(l)) == (5,)
    end

    @testset "Neuron value Conv default" begin
        l = ml(Conv((2,3), 4=>5))
        @test size(neuron_value(l)) == (5,)
    end

    @testset "Neuron value Dense act contrib" begin
        l = ml(Dense(3,5), ActivationContribution)
        @test neuron_value(l) == zeros(5)
        tr(l, ones(Float32, 3, 4))
        @test size(neuron_value(l)) == (5,)
    end

    @testset "Neuron value RNN act contrib" begin
        l = ml(RNN(3,5), ActivationContribution)
        @test neuron_value(l) == zeros(5)
        tr(l, ones(Float32, 3, 8))
        @test size(neuron_value(l)) == (5,)
    end

    @testset "Neuron value Conv act contrib" begin
        l = ml(Conv((3,3), 2=>5, pad=(1,1)), ActivationContribution)
        @test neuron_value(l) == zeros(5)
        tr(l, ones(Float32, 4,4,2,5))
        @test size(neuron_value(l)) == (5,)
    end

    @testset "Neuron value MaxPool act contrib" begin
        l = ml(MaxPool((3,3)), ActivationContribution, insize=2)
        @test ismissing(neuron_value(l))
        tr(l, ones(Float32, 4,4,2,5))

        @test ismissing(minΔninfactor(l))
        @test ismissing(minΔnoutfactor(l))
        @test size(neuron_value(l)) == (2,)
    end

    @testset "Elem add ActivationContribution" begin
        ac(l) = ActivationContribution(l)
        v = ac >> ml(Dense(2,3)) + ml(Dense(4,3))
        tr(v, [1 1 1]', [1 2 3]', [4 5 6]')
        @test size(neuron_value(v)) == (3,)

        @test minΔninfactor(v) == 1
        @test minΔnoutfactor(v) == 1

        g = CompGraph(vcat(inputs.(inputs(v))...), v)
        @test size(g(ones(Float32, 2,2), ones(Float32, 4, 2))) == (nout(v), 2)
    end

    @testset "Concat ActivationContribution" begin
        v = concat(ml(Dense(2,3)), ml(Dense(4,5)), layerfun=ActivationContribution)
        tr(v,ones(nout(v), 1), [1 2 3]', [4 5 6 7 8]')
        @test size(neuron_value(v)) == (nout(v),)

        @test minΔninfactor(v) == 1
        @test minΔnoutfactor(v) == 1

        g = CompGraph(vcat(inputs.(inputs(v))...), v)
        @test size(g(ones(Float32, 2,2), ones(Float32, 4, 2))) == (nout(v), 2)
    end

    @testset "Function ActivationContribution" begin
        # Not really an absorbvertex, but ActivationContribution should work on stuff which is not FLux layers.
        v = absorbvertex(ActivationContribution(x -> 2x), 3, ml(Dense(2,3)))
        tr(v, ones(nout(v), 1))
        @test size(neuron_value(v)) == (nout(v),)

        @test minΔninfactor(v) == 1
        @test minΔnoutfactor(v) == 1

        g = CompGraph(vcat(inputs.(inputs(v))...), v)
        @test size(g(ones(Float32, 2,2))) == (nout(v), 2)
    end

    @testset "Mutate ActivationContribution" begin
        l = ml(Dense(3,5), ActivationContribution ∘ LazyMutable)
        Δnout(l, -1)
        Δoutputs(l, v -> 1:nout_org(v))
        apply_mutation(l)
        @test size(l(ones(Float32, 3,2))) == (4, 2)
        @test size(neuron_value(l)) == (4,)
    end

    @testset "Mutate ActivationContribution MaxPool" begin
        l1 = ml(Conv((3,3), 2=>5, pad=(1,1)), ActivationContribution ∘ LazyMutable)
        l2 = mutable(MaxPool((3,3), pad=(1,1)), l1, layerfun = ActivationContribution ∘ LazyMutable)
        g = CompGraph(inputs(l1), l2)

        # Mutate before activation contribution for l2 has been initialized
        Δnout(l1, 1)
        Δoutputs(l1, v -> 1:nout_org(v))
        apply_mutation(g)
        @test size(neuron_value(l1)) == (6,)
        @test ismissing(neuron_value(l2))

        # This will initialize it
        @test size(g(ones(Float32, 4,4,2,3))) == (2,2,6,3)

        Δnout(l1, -2)
        Δoutputs(l1, v -> 1:nout_org(v))
        apply_mutation(g)
        tr(g, ones(Float32, 2,2,4,3), ones(Float32, 4,4,2,3))
        @test size(neuron_value(l1)) == size(neuron_value(l2)) == (4,)

        Δnin(l2, -1)
        Δoutputs(l1, v -> 1:nout_org(v))
        apply_mutation(g)
        tr(g, ones(Float32, 2,2,3,3), ones(Float32, 4,4,2,3))
        @test size(neuron_value(l1)) == size(neuron_value(l2)) == (3,)
    end
end


text:
using YaoExperiment
using Yao, Yao.ConstGate, YaoExtensions
using Test

# x, y, z, h, s, sdag, t, tdag
@testset "basic gates" begin
    for G in [:X, :Y, :Z, :H, :S, :Sdag, :T, :Tdag]
        g = string(G)
        @eval @test quil($G) == "$($(g))"
        @eval @test quil(put(3,1=>$G)) == "$($g) 0"
    end

    for G in [:SWAP, :CNOT]
        g = string(G)
        @eval @test quil($G) == "$($(g))"
        @eval @test quil(put(3,(2,1)=>$G)) == "$($g) 1 0"
    end
end
# id
@test quil(put(4, 2=>I2)) == "I 1"

# cx, cy, cz, ch, ct
# ccx, ccy, ccz, cch, cct
@testset "basic controled gates" begin
    for G in [:X, :Y, :Z, :H, :T]
        g = G == :X ? "NOT" : string(G)
        @eval @test quil(control(3, 1, 2=>$G)) == "C$($g) 0 1"
        @eval @test quil(control(3, -1, 2=>$G)) == "¬C$($g) 0 1"
        @eval @test quil(control(7, (1,3), 2=>$G)) == "CC$($g) 0 2 1"
    end
end

# rx, ry, rz
# crz, crx, cry
@testset "(controled) rotation gates" begin
    for G in [:X, :Y, :Z]
        g = string(G)
        @eval @test quil(control(3, 1, 2=>rot($G, 0.3))) == "CR$($g)(0.3) 0 1"
        @eval @test quil(rot($G, 0.3)) == "R$($g)(0.3)"
        @eval @test quil(put(3, 2=>rot($G, 0.3))) == "R$($g)(0.3) 1"
    end
    @test quil(rot(SWAP, 0.3)) == "PSWAP(0.3)"
    @test quil(put(3, (3,1)=>rot(SWAP, 0.3))) == "PSWAP(0.3) 2 0"
end

# measure
@testset "measure" begin
    @test quil(Measure(3), YaoExperiment.QuilAddress([1,4,3])) == "MEASURE 0 3 2"
end

@testset "composite" begin
    # kron
    @test quil(kron(5, 2=>X, 4=>T)) == "X 1\nT 3"
    @test quil(PauliString(I2, I2, X, X)) == "X 2\nX 3"
    # chain
    @test quil(chain(put(3,1=>Z), put(3, 2=>X))) == "Z 0\nX 1"
    # Concentrator
    @test quil(subroutine(5, put(3,1=>Z), [3,4,1])) == "Z 2"
    # cache and diff
    @test quil(kron(5, 2=>X, 4=>T) |> cache) == "X 1\nT 3"
    @test quil(rot(X, 0.3)) == "RX(0.3)"
end

@testset "show" begin
    # show
    @test repr(MIME("quil"), control(3, 2, 1=>X)) |> String == "CNOT 1 0\n"
end

@testset "system test" begin
    circuit = chain(qft_circuit(4), Measure())
    println(quil(circuit))
end


text:
using DashHtmlComponents

function Header(title::String)
    return html_div(
        html_h2(
            title,
            style = (borderBottom = "thin lightgrey solid", marginRight = 20,)
        )
    )
end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
include("../src/global/ImportFolders.jl")

import Config
import Server
import Connection
import Event
import Helper

function onBeforeResponse(sock::Base.TCPSocket, data::Array{Any, 1}, ans::Connection.Answer)
  println("before-response\n", data)
  ans.id   = UInt8(1)
  ans.data = [UInt64(rand(1:300))]
end

function onAfterResponse(data::Array{Any, 1})
  println("after-response\n", data)
end

con = Server.create(ip"127.0.0.1", Config.val(:conServerPort), true)
Event.on(con.observer, Connection.EVENT_BEFORE_RESPONSE, onBeforeResponse)
Event.on(con.observer, Connection.EVENT_AFTER_REQUEST, onAfterResponse)
Server.run(con)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    VERSION

const VERSION = begin
    project = joinpath(dirname(@__FILE__), "..", "..", "Project.toml")
    v = []
    for line in eachline(project)
        if startswith(line, "version")
            push!(v, strip(strip(split(line, " = ")[2]), ['"']))
            break
        end
    end
    v[1]
end


text:
dirname = @__DIR__
dirname[end] == '/' || (dirname = dirname * "/")
dir_netlib = dirname * "../dat/netlib/"

include("./benchmark.jl")

solvers = [SOLVERS_OS; SOLVERS_COMM]

# First compilation round
@info "Compilation rounds... This may take a few minutes"
results = benchmark(
    dir_netlib,
    ["AFIRO.SIF"],
    solvers,
    verbose=false
);

# Compilation round for Tulip
benchmark(
    dir_netlib,
    readdir(dir_netlib),
    [sTLP()],
    verbose=false
);

# Real deal
@info "Launching Netlib benchmark"
flush(stdout)
results = benchmark(
    dir_netlib,
    readdir(dir_netlib),
    solvers,
    verbose=true
);

# extract results
@info "Parsing results"
df = extract_results(results)

# Export individual results
@info "Exporting to csv file"
CSV.write(dirname*"res_netlib.csv", df)


text:
"""
```julia
randDiagonal(d, n)

randDiagonal(n)
```
- `d` : default `Normal(0，1)`, entry distribution
- `n` : dimension

# Examples

Generates a 3 by 3 diagonal matrix, with non-zero elements from `Normal(0,1)`
```julia
randDiagonal(3)

3×3 Diagonal{Float64, Vector{Float64}}:
 0.440359   ⋅         ⋅
  ⋅        1.94832    ⋅
  ⋅         ⋅       -0.52536
```

Generates a 5 by 5 diagonal matrix, with non-zero elements from `Poisson(2)`
```julia
randDiagonal(Poisson(2),5)

5×5 Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅  ⋅  ⋅
 ⋅  0  ⋅  ⋅  ⋅
 ⋅  ⋅  0  ⋅  ⋅
 ⋅  ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  ⋅  3
```
"""
function randDiagonal(d::D, n::Int)  where D<:S
    Diagonal(rand(d,n))
end

randDiagonal(n::Int) =  randDiagonal(Normal(), n::Int)

"""
```julia
randTriangular(d , n ;  diag , Diag, upper )

randTriangular(n;diag, upper)
```
- `d` : entry distribution
- `n` : dimension
- `diag` : default `diag = d`, diagonal entry distribution
- `Diag` : default `Diag = true`, `true` includes diagonal, `false` with diagonal entries 0
- `upper` : default `upper = true`, `true` gives upper triangular, `false` gives lower triangular

# Examples
Generate an upper triangular matrix with entries Standard Normal
```julia
randTriangular(3)

3×3 UpperTriangular{Float64, Matrix{Float64}}:
 -0.572757  -0.459518   -1.60622
   ⋅         0.0216834  -0.416529
   ⋅          ⋅         -1.00807
```

Generate a 3 by 3 strictly lower triangular matrix, with nonzero entries uniform from ``\\{1,2,3\\}``
```julia
randTriangular(1:3,3,upper=false,Diag=false)

3×3 LowerTriangular{Int64, Transpose{Int64, Matrix{Int64}}}:
 0  ⋅  ⋅
 3  0  ⋅
 3  2  0
```
"""
function randTriangular(d::D, n::Int; diag=d::D, Diag=true, upper = true)  where D<:S
    M = UpperTriangular(zeros(eltype(d),n,n))
    for i in 1:n, j in i+1:n
        M[i,j] = rand(d)
    end
    if Diag
        M += randDiagonal(diag,n)
    end

    if !upper
        M = transpose(M)
    end
    M
end
randTriangular(n::Int; Diag=true, upper = true) = randTriangular(Normal(),n,Diag=Diag,upper=upper)

"""
```julia
randMatrix(d::D, n::Int, m = n::Int; norm = false::Bool) where D<:S

randMatrix(n::Int, m = n::Int; norm = false::Bool)
```
- `d` : entry distribution
- `n`,`m` : default `m = n` , dimensions
- `norm` : default `false`, if `norm` set to `true`, then the matrix will be normlaized with ``\\operatorname{min}(n,m)^{-1/2}``.
# Examples
Generates a 2 by 2 random  matrix with entries from the Standard  Gaussian.
```julia
randMatrix(2)

2×2 Matrix{Float64}:
 1.74043  -1.30317
 0.72765   0.639943
```

Generates a 3 by 2 random  matrix with entries uniformly from {1,2,3,...,10}.
```julia
randMatrix(1:10,3,2)

3×2 Matrix{Int64}:
  1  3
  6  4
 10  1
```
Generate a normalized random 2 by 2  Matrix with entries  `Poisson(2)` rvs.
Need to import the `Distributions` package for `Poisson(2)`
```julia
using Distributions
randMatrix(Poisson(2),2,norm = true)

2×2 Matrix{Float64}:
 1.41421   0.0
 0.707107  1.41421
```
"""
function randMatrix(d::D, n::Int, m = n::Int ; norm = false::Bool)  where D<:S

    M =  rand(d,n,m)

    if norm
         M/=sqrt(min(n,m))
    end

    return M
end

function randMatrix(n::Int, m = n::Int; norm = false::Bool)
    return randMatrix(Normal(), n,m ; norm = norm)
end

max_stars_repo_path: src/MatrixModels/IID Matrices.jl

max_issues_repo_path: src/MatrixModels/IID Matrices.jl

max_forks_repo_path: src/MatrixModels/IID Matrices.jl


text:
#coefficients for two dimensional example

####### dimension

d = 2

######## starting point

u = [1.,0.]

####### drift

bB = [ -0.2   -1;   0.5  -0.4]
B = copy(bB)
#beta = [0.5,-0.5]
bbeta = [0., 0.]
beta = copy(bbeta)

b(t,x) = exp(-0.2*t)*bB*x + bbeta
#integral of exp(-0.2*t) to be used
ph(t,s) = 5.*exp(-0.2*s)-5.*exp(-0.2*t)

####### diffusion coefficient

# uniformly elliptic
function sigma0(s,x)
 	rho = 2atan(2*x[2])/pi
 	return( 0.7* si* [(2-rho)*[1, 1]  (2+rho)*[1,-1]])

end
a0(s,x) = sigma0(s,x)*sigma0(s,x)'

# [0.5,1]-banded elliptic
function sigma1(s,x)
 	m = norm(x-[-0.8, 0.25])
	[1 1; -1  (1 + 1atan(2m)/pi)]
end
a1(s,x) = sigma1(s,x)*sigma1(s,x)'

# only depending on time
sigma2(t,x) = sigma0(0., [sin(t), cos(t)])
a2(s,x) = sigma2(s,x)*sigma2(s,x)'

# choose one
sigma(t,x) = sigma0(t,x)
a(t,x) = a0(t,x)

# fucntion to ignorantly test, if a is uniformly elliptic
function testa(a, N)
	Ra(a,x) = scalar(x'a(0,x)x/(x'x))
	xmax = [NaN, NaN]
	xmin = [NaN, NaN]
	rmax = -Inf
	rmin = Inf
	for i in 1:N
		x = 5*randn(2)
		r = Ra(a, x)
		if (r > rmax)
			rmax = r
			xmax = x
		end
		if (r < rmin)
			rmin = r
			xmin = x
		end

	end
	return rmin, xmin,rmax, xmax
end


text:
#Julia Control Flow

function larger(x, y)
    if (x>y)
        return x
    end
    return y
end
println(larger(7,8))

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
using KrylovMethods
using Base.Test

println("=== Testing blockBiCGSTB for real matrix === ")
nrhs = 4;
n = 100
A  = sprandn(n,n,.1) + 10*speye(n)
D  = diag(A)
Af(x) = A*x
M(x)  = D.\x
rhs = randn(n,nrhs)

# test flag for early stopping
xt = blockBiCGSTB(A,rhs,tol=1e-6,maxIter=3,out=2)
@test xt[2]==-1

# test handling of zero rhs
xt = blockBiCGSTB(A,zeros(n,nrhs),tol=1e-6,maxIter=3,out=2)
@test  all(xt[1].==0)
@test xt[2]==-9

x0 = blockBiCGSTB(full(A),rhs,tol=1e-6)
x1 = blockBiCGSTB(A,rhs,tol=1e-6,out=1)
x2 = blockBiCGSTB(Af,rhs,tol=1e-6)
x3 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=100,x=randn(size(rhs)))
x5 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=100,M1=M)

@test vecnorm(A*x0[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x1[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x2[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x3[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x5[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(x2[1]-x1[1])/vecnorm(x1[1]) < 1e-12
@test vecnorm(x3[1]-x1[1])/vecnorm(x1[1]) < 1e-5
@test vecnorm(x0[1]-x1[1])/vecnorm(x1[1]) < 1e-5

println("Testing BlockBiCGSTAB for complex matrix")
A  = sprandn(n,n,.1) + 10*speye(n) + im*(sprandn(n,n,.1) + 10*speye(n) )
D  = diag(A)
Af(x) = A*x
M(x)  = D.\x
rhs = randn(n,nrhs) + 1im * randn(n,nrhs)

x1 = blockBiCGSTB(A,rhs,tol=1e-6)
x2 = blockBiCGSTB(Af,rhs,tol=1e-6)
x3 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=200,x=randn(size(rhs))+im*randn(size(rhs)))
x4 = blockBiCGSTB(Af,rhs,tol=1e-6,maxIter=200,M1=M)

@test vecnorm(A*x1[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x2[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x3[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(A*x4[1]-rhs)/vecnorm(rhs) < 1e-6
@test vecnorm(x2[1]-x1[1])/vecnorm(x1[1]) < 1e-12
@test vecnorm(x3[1]-x1[1])/vecnorm(x1[1]) < 1e-5

println("=== BlockBiCGSTAB: All tests passed ===")


text:
# This file holds the definition of the functions pertaining to the
#   CompoundState type.

# The functions of the CompoundState type require the State type.
requiredTypes = [ "state", "compoundState" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export setName!,
       setStateTarget!,
       addStateToCompound!,
       removeStateFromCompound!,
       clearStatesFromCompound!,
       generateHierarchy

"""
```
setName!( compState::CompoundState,
          name::String )
```
This function sets the name of the compound state `compState` to `name`.

This function returns `nothing`.
"""
function setName!( compState::CompoundState, name::String )::Void

    compState.name = name
    return

end  # setName!( compState::CompoundState, name::String )

"""
```
setStateTarget!( compState::CompoundState,
                 target::T )
    where T <: Integer
```
This function sets the target number of personnel members in compound state
`compState` to `target`. If the number is less than zero, it means there's no
target.

This function returns `nothing`.
"""
function setStateTarget!( compState::CompoundState, target::T )::Void where T <: Integer

    compState.stateTarget = target < 0 ? -1 : target
    return

end  # setStateTarget!( compState, target )

"""
```
addStateToCompound!( compState::CompoundState,
                     stateList::String... )
```
This function adds the states in `stateList` to the list of states composing the
compound state `compState`.

This function returns `nothing`.
"""
function addStateToCompound!( compState::CompoundState,
    stateList::String... )::Void

    for stateName in stateList
        if stateName ∉ compState.stateList
            push!( compState.stateList, stateName )
        end  # if stateName ∉ compState.stateList
    end  # for stateName in stateList

    return

end  # addStateToCompound!( compState, stateList )

"""
```
removeStateFromCompound!( compState::CompoundState,
                          stateList::String... )
```
This function removes the states in `stateList` from the list of states
composing the compound state `compState`.

This function returns `nothing`.
"""
function removeStateFromCompound!( compState::CompoundState,
    stateList::String... )::Void

    stateFlags = map( stateName -> stateName ∈ stateList, compState.stateList )
    deleteat!( compState.stateList, stateFlags )
    return

end  # removeStateFromCompound!( compState, stateList )

"""
```
clearStatesFromCompound!( compState::CompoundState )
```
This function clears the list of states composing the compound state
`compState`.

This function returns `nothing`.
"""
function clearStatesFromCompound!( compState::CompoundState )::Void

    empty!( compState.stateList )
    return

end  # clearStatesFromCompound!( compState )

"""
```
generateHierarchy( mpSim::ManpowerSimulation,
                   attrList::String... )
```
This is a shorthand functin for `generateHierarchy( mpSim, true, attrList )`.
"""
function generateHierarchy( mpSim::ManpowerSimulation,
    attrList::String... )::Void

    return generateHierarchy( mpSim, true, attrList... )

end  # generateHierarchy( mpSim, attrList... )

"""
```
generateHierarchy( mpSim::ManpowerSimulation,
                   useCat::Bool,
                   attrList::String... )
```
This function generates a hierarchy of compound states in the manpower
simulation `mpSim`, splitting them up by the values of the attributes in
`attrList`, in the order given. If `useCat` is `true`, the function will attempt
to get the names for the compound states where available, and will generate
descriptive names otherwise.

This function returns `nothing`.
"""
function generateHierarchy( mpSim::ManpowerSimulation, useCat::Bool,
    attrList::String... )::Void

    # Retain only the attributes that actually exist in the simulation.
    attrsInSim = vcat( mpSim.initAttrList, mpSim.otherAttrList )
    attrNames = map( attr -> attr.name, attrsInSim )
    tmpAttrList = unique( collect( attrList ) )
    filter!( attrName -> attrName ∈ attrNames, tmpAttrList )

    if isempty( tmpAttrList )
        return
    end  # if isempty( tmpAttrList )

    attrs = similar( tmpAttrList, PersonnelAttribute )

    # Get the attributes with the retained names.
    for ii in eachindex( tmpAttrList )
        attrInd = findfirst( attr -> attr.name == tmpAttrList[ ii ],
            attrsInSim )
        attrs[ ii ] = attrsInSim[ attrInd ]
    end  # for ii in eachindex( tmpAttrList )

    # Create a list of all the states in the catalogue.
    catStateList = Dict{String, Array{String}}()

    if useCat
        XLSX.openxlsx( mpSim.catFileName ) do catXF
            stateCat = catXF[ "States" ]
            sLine = 2

            while isa( stateCat[ "A$sLine" ], String )
                stateName = stateCat[ "A$sLine" ]
                nAttrs = stateCat[ "G$sLine" ]
                catStateList[ stateName ] = Array{String}( nAttrs, 2 )

                for ii in 1:nAttrs
                    jj = 6 + ii * 2
                    catStateList[ stateName ][ ii, : ] = [
                        stateCat[ XLSX.CellRef( sLine, jj ) ],
                        stateCat[ XLSX.CellRef( sLine, jj + 1 ) ] ]
                end  # for ii in 1:nAttrs

                sLine += 1
            end  # while isa( catXF[ "A$sLine" ], String )
        end  # XLSX.openxlsx( mpSim.catFileName ) do catXF
    end  # if useCat

    return partitionByAttribute( collect( keys( mpSim.stateList ) ), attrs, 1,
        mpSim, catStateList, collect( keys( catStateList ) ) )

end  # generateHierarchy( mpSim, useCat, attrList )

function Base.show( io::IO, compState::CompoundState )

    print( io, "    Compound state: " * compState.name )

    if isempty( compState.stateList )
        print( io, "\n      No component states in compound state." )
    else
        print( io, "\n      Component states: " *
            join( compState.stateList, ", " ) )
    end  # if isempty( compState.stateList )

    if compState.stateTarget >= 0
        print( io, "\n      Personnel target: $(compState.stateTarget) personnel members" )
    else
        print( io, "\n      No personnel target for compound state." )
    end  # if compState.stateTarget >= 0

end  # show( io, compoundState )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

"""
```
isPersonnelOfState( persAttrs::Dict{String, Any},
                    compState::CompoundState,
                    mpSim::ManpowerSimulation )
```
This function tests if the person with attributes in `persAttrs` belongs to any
of the states making up the compound state `compStates`, where the base states
are defined in the manpowre simulation `mpSim`.

This function returns a `Bool`, the result of the test.
"""
function isPersonnelOfState( persAttrs::Dict{String, Any},
    compState::CompoundState, mpSim::ManpowerSimulation )::Bool

    return any( stateName -> isPersonnelOfState( persAttrs,
        mpSim.stateList[ stateName ] ), compState.stateList )

end  # isPersonnelOfState( persAttrs, compState, mpSim )

"""
```
readCompoundState( sheet::XLSX.Worksheet,
                   sLine::Int )
```
This function reads and processes the custom compound state defined in the Excel
worksheet `sheet` on line `sLine`.

This function returns a `CompoundState` object, the processed compound state.
"""
function readCompoundState( sheet::XLSX.Worksheet, sLine::Int )::CompoundState

    newCompState = CompoundState( sheet[ "H$sLine" ] )
    nStates = sheet[ "J$sLine" ]
    stateList = sheet[ XLSX.CellRange( sLine, 11, sLine, 10 + nStates ) ]
    addStateToCompound!( newCompState, stateList... )

    return newCompState

end  # readCompoundState( sheet::XLSX.worksheet, sLine::Int )::CompoundState

"""
```
processCompoundStates( mpSim::ManpowerSimulation )
```
This function processes the compound states that have been configured based on a
definition in the state catalogue.

This function returns `nothing`.
"""
function processCompoundStates( mpSim::ManpowerSimulation )::Void

    stateList = collect( keys( mpSim.stateList ) )

    for compStateName in keys( mpSim.compoundStates )
        compState = mpSim.compoundStates[ compStateName ]
        compStateReqs = compState.requirements

        isCompState = map( stateList ) do stateName
            stateReqs = mpSim.stateList[ stateName ].requirements
            return all( keys( compStateReqs ) ) do attrName
                return haskey( stateReqs, attrName ) &&
                    ( stateReqs[ attrName ] == compStateReqs[ attrName ] )
            end  # all( keys( compStateReqs ) ) do attrName
        end  # map( stateList ) do stateName

        addCompoundState!( mpSim, compStateName, -1,
            stateList[ isCompState ]... )
    end  # for compState in mpSim.compoundStates

    return

end  # processCompoundStates( mpSim )

"""
```
readHierarchy( sheet::XLSX.Worksheet,
               mpSim::ManpowerSimulation )
```
This function reads the attributes for forming a compound state hierarchy from
the Excel sheet `sheet`, and generates the hierarchy for manpower simulation
`mpSim`.

This function returns `nothing`.
"""
function readHierarchy( sheet::XLSX.Worksheet, mpSim::ManpowerSimulation )::Void

    nAttrs = sheet[ "F4" ]
    attrList = sheet[ XLSX.CellRange( 7, 5, 6+ nAttrs, 5 ) ]
    generateHierarchy( mpSim, attrList... )
    return

end  # readHierarchy( sheet, mpSim )

"""
```
partitionByAttribute( stateList::Vector{String},
                      attrList::Vector{PersonnelAttribute},
                      level::Int,
                      mpSim::ManpowerSimulation,
                      catStateList::Dict{String, Array{String}},
                      catStateNameList::Vector{String},
                      name::String = "",
                      partition::Array{String} = Array{String}( 0, 2 ) ):
```
This function generates a partition of the states in `stateList` and creates
compound states for those parts containing more than one base state. The other
parameters are:
* `stateList`: the ordered attributes for the hierarchy.
* `level`: the level of the hierarchy that needs to be generated.
* `mpSim`: the manpower simulation.
* `catStateList`: a list of all the states in the catalogue with the
  requirements for being in that state.
* `catStateNameList`: a list of all the names of the states in the catalogue.
* `name`: the current auto-generated name for the compound state containing all
  the states in `stateList`. This is not necessarily the name that is used.
* `partition`: the attributes and their values definingthe compound state
  containing all the states in `stateList`.

This function returns `nothing`.
"""
function partitionByAttribute( stateList::Vector{String},
    attrList::Vector{PersonnelAttribute}, level::Int,
    mpSim::ManpowerSimulation, catStateList::Dict{String, Array{String}},
    catStateNameList::Vector{String}, name::String = "",
    partition::Array{String} = Array{String}( 0, 2 ) )::Void

    statePartition = Dict{String, Vector{String}}()
    listOfStates = map( stateName -> mpSim.stateList[ stateName ], stateList )
    attr = attrList[ level ]
    statePartitioned = falses( stateList )

    # Generate the partition.
    for attrVal in attr.possibleValues
        stateInds = map( listOfStates ) do state
            return haskey( state.requirements, attr.name ) &&
                ( attrVal ∈ state.requirements[ attr.name ] )
        end  # map( listOfStates ) do state

        if any( stateInds )
            statePartition[ attrVal ] = stateList[ stateInds ]
        end  # if any( stateInds )

        statePartitioned .|= stateInds
    end  # for attrVal in attrList[ level ].possibleValues

    isPartitionProper = !isempty( statePartition )

    # Generate part for states with attribute not filled in.
    if !all( statePartitioned )
        statePartition[ "undefined" ] = stateList[ .!statePartitioned ]
    end  # if !all( statePartitioned )

    # Create the partition and recurse.
    for attrVal in keys( statePartition )
        newName = name
        newPart = partition

        if isPartitionProper
            newName *= ( name == "" ? "" : "; " ) * "$(attr.name):$attrVal"
            newPart = vcat( partition, [ attr.name attrVal ] )
        end  # if isPartitionProper

        # Check if this compound state already has a name.
        nPart = size( newPart )[ 1 ]
        tmpStateNames = filter( catStateNameList ) do stateName
            nReqs = size( catStateList[ stateName ] )[ 1 ]
            return nPart == nReqs
        end  # filter( catStateNameList ) do stateName

        for ii in 1:nPart
            tmpStateNames = filter( tmpStateNames ) do stateName
                stateReqs = catStateList[ stateName ]
                attrIndex = findfirst( newPart[ ii, 1 ] .== stateReqs[ :, 1 ] )

                if attrIndex == 0
                    return false
                end  # if attrIndex == 0

                return newPart[ ii, 2 ] == stateReqs[ attrIndex, 2 ]
            end  # filter( tmpStateNames ) do stateName
        end  # for ii in 1:nPart

        # Generate compound state.
        if length( statePartition[ attrVal ] ) > 1
            if isPartitionProper
                addCompoundState!( mpSim, isempty( tmpStateNames ) ? newName :
                    tmpStateNames[ 1 ], -1, statePartition[ attrVal ]... )
            end  # if isPartitionProper

            if level < length( attrList )
                partitionByAttribute( statePartition[ attrVal ], attrList,
                    level + 1, mpSim, catStateList, catStateNameList, newName,
                    newPart )
            end  # if level < length( attrList )
        end  # if length( statePartition[ attrVal ] ) > 1
    end  # for attrVal in keys( statePartition )

    return  # for attrVal in attrList[ level ].possibleValues

end  # partitionByAttribute( stateList, attrList, level, mpSim, catStateList,
     #   catStateNameList, name, partition )


text:
using Serialization: serialize, deserialize

"""
    pack(data)

Utility function to serialize `data` into a byte array.
"""
function pack(data)::Vector{UInt8}
    buffer = IOBuffer()
    serialize(buffer, data)
    return take!(buffer)
end

"""
    unpack(bytes)

Utility function to deserialize `bytes` into its original representation.
"""
unpack(bytes::Vector{UInt8}) = IOBuffer(bytes) |> deserialize


text:
# import Singular
import AbstractAlgebra

using Base.Threads

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

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

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

    gb = Groebner.groebner(system)

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

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

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

    ideal_s = Singular.Ideal(R_s, system_s)

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

    bench = @benchmarkable  Singular.std($ideal_s) samples=3
    println(median(run(bench)))
end

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

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

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
    ("siwr", read_SIWR()), # 0.19 vs 0.03
    ("seaijrc", read_SEAIJRC()), # 0.19 vs 0.03
    ("cyclic 7", Groebner.cyclicn(7, ground=ground)), # 0.19 vs 0.03
    ("cyclic 8", Groebner.cyclicn(8, ground=ground)), # 0.60 vs 0.9
    ("katsura 8",Groebner.katsuran(8, ground=ground)), # 1.35 vs 71.121
    ("katsura 9",Groebner.katsuran(9, ground=ground)), # 9.68 vs 775
    ("eco 10",Groebner.eco10(ground=ground)),         # 0.41  vs 31.7
    ("eco 11",Groebner.eco11(ground=ground)),         # 2.5   vs 341
    ("noon 8"    ,Groebner.noonn(8, ground=ground)),  # 0.19  vs 0.36
    ("noon 9"    ,Groebner.noonn(9, ground=ground)),  # 1.8   vs 3.2
    ("henrion 6"    ,Groebner.henrion6(ground=ground)),  # 0.19  vs 0.36
    ("henrion 7"    ,Groebner.henrion7(ground=ground))  # 1.8   vs 3.2
    ]

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

ground = AbstractAlgebra.QQ
run_f4_ff_degrevlex_benchmarks(ground)


text:
"""
# InteractiveCodeSearch.jl --- Interactively search Julia code

Julia has `@edit`, `@less`, etc. which are very handy for reading the
implementation of functions.  However, you need to specify a "good
enough" set of (type) parameters for them to find the location of the
code.

Instead, `InteractiveCodeSearch` provides a few macros to
interactively choose the code you want to read.

## Features

* Interactively choose a method signature before opening the code
  location in your editor.

* Various ways to search methods, such as: by function name `@search show`,
  function call expression `@search show(stdout, "hello")`,
  function call signature `@search show(::IO, ::String)`,
  module name `@search Base`, argument value `@searchmethods 1`,
  argument type `@searchmethods ::Int`, and return type `@searchreturn Int`.

* Interactively search history.  It works in IJulia as well.

## Examples

```julia
using InteractiveCodeSearch
@search show             # search method definitions
@searchmethods 1         # search methods defined for integer
@searchhistory           # search history (Julia ≥ 0.7)
@searchreturn String Pkg # search methods returning a given type (Julia ≥ 0.7)
```

## Requirements

* Interactive matching command.  For example:
  * [peco](https://github.com/peco/peco) (default in terminal)
  * [percol](https://github.com/mooz/percol)
  * [rofi](https://github.com/DaveDavenport/rofi) (GUI; default in IJulia)
"""
module InteractiveCodeSearch
export @search, @searchmethods

import Pkg
using Base
using Base: IOError
using InteractiveUtils: edit, gen_call_with_extracted_types, methodswith

abstract type SearchPolicy end
struct Shallow <: SearchPolicy end
struct Recursive <: SearchPolicy end

mutable struct SearchConfig  # CONFIG
    open
    interactive_matcher::Function
    auto_open::Bool
end

function setmatcher!(cmd::Function, obj::SearchConfig)
    obj.interactive_matcher = cmd
    return
end

function setmatcher!(cmd::Cmd, obj::SearchConfig)
    setmatcher!(convertCmd(cmd), obj)
end

function convertCmd(cmd)
    maybe_warn_matcher(cmd)
    return function (f)
        maybe_warn_matcher(cmd)
        f(cmd)
    end
end

maybe_identifier(s) = !startswith(string(s), "#")

is_locatable(::Any) = false
is_locatable(::Function) = true
is_locatable(t::Type) = !(t <: Vararg)
# https://github.com/JuliaLang/julia/issues/29645

is_defined_in(child, parent) =
    child !== parent && parentmodule(child) === parent

function list_locatables(p::SearchPolicy, m::Module)
    locs = []
    for s in names(m; all=true)
        if maybe_identifier(s)
            x = try
                getfield(m, s)
            catch err
                err isa UndefVarError && continue
                rethrow()
            end
            if is_locatable(x)
                push!(locs, x)
            elseif p isa Recursive && x isa Module && is_defined_in(x, m)
                append!(locs, list_locatables(p, x))
            end
        end
    end
    return locs
end

module_methods(p::SearchPolicy, m::Module) :: Vector{Method} =
    vcat(collect.(methods.(list_locatables(p, m)))...)
# Note: the conversion `:: Vector{Method}` seems to be required only
# for Julia 0.6.

struct _Dummy end

function uninteresting_locs()
    locs = []
    for m in methods(_Dummy)
        path = string(m.file)
        if path != @__FILE__
            push!(locs, (path, m.line))
        end
    end
    return locs
end

"""
    find_source_file(file)

Find source `file` and return its full path.  It just calls
`Base.find_source_file` and return its result for normal Julia
installation.  For nightly Julia build, it tries to guess the right
path when `Base.find_source_file` failed.
"""
function find_source_file(file)
    path = Base.find_source_file(file)
    if path isa AbstractString && ! isfile(path)
        for m in methods(Pkg.add)
            exfile = try
                String(m.file)
            catch err
                continue
            end
            idx = findlast(joinpath(Base.Filesystem.path_separator,
                                    "share", "julia"), exfile)
            if idx isa Nothing
                continue
            end
            prefix = exfile[1:idx[1]]
            if startswith(file, prefix)
                # e.g., relpath = "share/julia/stdlib/v0.7/..."
                relpath = file[length(prefix)+1:end]
                return joinpath(Base.Sys.BINDIR, "..", relpath)
            end
        end
    end
    return path
end

function _readandwrite(matcher)
    proc = matcher() do cmd
        open(`$cmd`, "r+")
    end
    return (proc.out, proc.in, proc)
end

"""
    read_stdout(input::AbstractString, cmd)
    read_stdout(input_provider, cmd)

Julia implementation of "echo {input} | {cmd}".
"""
function read_stdout(input::AbstractString, cmd)
    read_stdout(cmd) do stdin
        write(stdin, input)
    end
end

function read_stdout(input_provider, cmd)
    stdout, stdin, process = _readandwrite(cmd)
    reader = @async read(stdout)
    try
        input_provider(stdin)
    catch err
        if ! (err isa IOError)
            rethrow()
        end
    finally
        close(stdin)
    end
    return fetch(reader)
end

function parse_loc(line)
    rest, lineno = rsplit(line, ":", limit=2)
    _, path = rsplit(rest, " at ", limit=2)
    return String(path), parse(Int, lineno)
end

function run_matcher(input)
    return String(read_stdout(input, CONFIG.interactive_matcher))
end

choose_method(methods::T) where T =
    _choose_method(Base.IteratorSize(T), methods)

function _choose_method(::Base.HasLength, methods)
    if isempty(methods)
        @info "No (interesting) method found"
        return
    end
    if CONFIG.auto_open && length(methods) == 1
        m = first(methods)
        loc = (string(m.file), m.line)
        if loc in uninteresting_locs()
            path, lineno = loc
            @info "Not opening uninteresting location: $path:$lineno"
            return
        end
        return loc
    end
    return _choose_method(Base.SizeUnknown(), methods)
end

function _choose_method(::Base.IteratorSize, methods)
    out = run_matcher() do stdin
        for m in methods
            show(stdin, m)
            println(stdin)
        end
    end
    if isempty(out)
        return
    end
    return parse_loc(out)
end

function run_open(path, lineno)
    @info "Opening $path:$lineno"
    CONFIG.open(find_source_file(path), lineno)
end

maybe_open(::Nothing) = nothing
maybe_open(x::Tuple{String, Integer}) = run_open(x...)

search_methods(methods) = maybe_open(choose_method(methods))

code_search_typed(f, t) = search_methods(methods(f, t))

code_search(::SearchPolicy, f::Base.Callable) = search_methods(methods(f))
code_search(p::SearchPolicy, m::Module) = search_methods(module_methods(p, m))

function code_search(p::SearchPolicy, ::T) where T
    @warn """Cannot search for given value of type $T
             Searching for its type instead..."""
    code_search(p, T)
end

"""
Configuration interface for `InteractiveCodeSearch`.

# Examples

```julia
using InteractiveCodeSearch
InteractiveCodeSearch.CONFIG.interactive_matcher = `peco`  # default in terminal
InteractiveCodeSearch.CONFIG.interactive_matcher = `percol`
InteractiveCodeSearch.CONFIG.interactive_matcher =
    `rofi -dmenu -i -p "🔎"`  # use GUI matcher (default in non-terminal
                              # environment like IJulia)
InteractiveCodeSearch.CONFIG.interactive_matcher =
    `rofi -dmenu -i -p "🔎" -fullscreen`  # bigger screen
InteractiveCodeSearch.CONFIG.open = edit  # default
InteractiveCodeSearch.CONFIG.open = less  # use Base.less to read code
InteractiveCodeSearch.CONFIG.auto_open = true   # default
InteractiveCodeSearch.CONFIG.auto_open = false  # open matcher even when there
                                                # is only one candidate
```

## Using InteractiveCodeSearch.jl by default

Put the following code in your `~/.julia/config/startup.jl`:

```julia
using InteractiveCodeSearch
# InteractiveCodeSearch.CONFIG.interactive_matcher = ...
```
"""
const CONFIG = SearchConfig(
    edit,                       # open
    x->error("uninitialized"),  # interactive_matcher
    true,                       # auto_open
)

should_eval(::Any) = false
should_eval(::Symbol) = true
should_eval(ex::Expr) = ex.head in (:., :ref)
# Given (say) `a.b[c].d[e]` It probably is better to only eval
# `a.b[c].d` and then search for `getindex(a.b[c].d, e)`.  But it's
# (1) a bit harder to implement and (2) evaluating the whole
# expression is still useful.  So let's keep the current
# implementation for a while.

isliteral(::Symbol) = false
isliteral(::Expr) = false
isliteral(::Any) = true

isline(::Any) = false
isline(ex::Expr) = ex.head == :line
isline(::LineNumberNode) = true

single_macrocall(::Any) = nothing
function single_macrocall(x::Expr)
    if x.head == :macrocall && all(isline.(x.args[2:end]))
        return x.args[1]
    elseif x.head == :block
        statements = findall(a -> !isline(a), x.args)
        if length(statements) == 1
            return single_macrocall(x.args[statements[1]])
        end
    end
    return nothing
end

explicitly_typed(::Any) = nothing
function explicitly_typed(ex::Expr)
    if ex.head == :call &&
            all(x isa Expr && x.head == :(::) for x in ex.args[2:end])
        return ex.args[1], [x.args[end] for x in ex.args[2:end]]
    end
    return nothing
end

# Julia >= 0.7:
parse_search_policy(flag::QuoteNode) = parse_search_policy(flag.value)
# Julia 0.6:
function parse_search_policy(flag::Expr)
    @assert flag.head == :quote
    @assert length(flag.args) == 1
    return parse_search_policy(flag.args[1])
end

function parse_search_policy(flag::Symbol)
    if flag in (:shallow, :s)
        return Shallow()
    elseif flag in (:recursive, :r)
        return Recursive()
    end
    error("Invalid flag $flag")
end

"""
    @search x [:shallow | :s | :recursive | :r]

List file locations at which `x` are defined in an interactive matcher
and then open the chosen location in the editor.

When `x` is a module, only the top-level definitions are searched.  To
search all definitions in the submodule, pass `:recursive` or `:r`
flag.

    @search

If no expression is provided, search for the method returned by the
previous execution; i.e., `x` defaults to `ans`.

# Examples
```julia
@search show                      # all method definitions
@search @time                     # all macro definitions
@search Base.Enums                # methods and macros in a module
@search REPL :r                   # search the module recursively
@search *(::Integer, ::Integer)   # methods with specified types
@search dot(π, ℯ)                 # methods with inferred types
```

Note that `@search` evaluates complex expression with `.` and `[]`
such as follows and search the returned value or the type of it:
```julia
@search Base.Multimedia.displays[2].repl
```
"""
macro search(x = :ans, flag = :(:shallow))
    p = parse_search_policy(flag)

    if should_eval(x)
        # Examples:
        #   @search show
        #   @search Base.Enums
        #   @search Base.Multimedia.displays[2].repl
        return :(code_search($p, $(esc(x))))
    end

    macrocall = single_macrocall(x)
    if macrocall !== nothing
        # Examples:
        #   @search @time
        #   @search begin @time end
        return :(code_search($p, $(esc(macrocall))))
    end

    func_type = explicitly_typed(x)
    if func_type !== nothing
        f, ts = func_type
        # Examples:
        #   @search *(::Integer, ::Integer)
        #   @search dot(::AbstractVector, ::SparseVector)
        return :(code_search_typed($(esc(f)), tuple($(esc.(ts)...))))
    end

    # Since `gen_call_with_extracted_types` does not handle literals,
    # let's handle this case here (although there are not much can be
    # done).
    if isliteral(x)
        # Examples:
        #   @search ""
        #   @search 1
        return :(code_search($p, $(esc(x))))
    end

    # Examples:
    #   @search 1 * 2
    #   @search dot([], [])
    gen_call_with_extracted_types(__module__, code_search_typed, x)
end

code_search_methods(T) = search_methods(methodswith(T; supertypes=true))

"""
    @searchmethods x
    @searchmethods ::X

Interactively search through `methodswith(typeof(x))` or
`methodswith(X)`.

# Examples
```julia
@searchmethods 1         # search methods defined for integer
@searchmethods ::Int     # search methods defined for a specified type
```
"""
macro searchmethods(x)
    if x isa Expr && x.head == :(::)
        if length(x.args) > 1
            @info "Ignoring: $(x.args[1:end-1]...) in $x"
        end
        :(code_search_methods($(esc(x.args[end]))))
    else
        :(code_search_methods(typeof($(esc(x)))))
    end
end

################################################################################
#                               matcher binaries                               #
################################################################################

@static if VERSION<v"1.3-"
    const preferred_terminal = Cmd[
        `peco`,
        `percol`,
    ]
else
    using peco_jll
    const preferred_terminal = Union{Function,Cmd}[
        `peco`,
        peco,
        `percol`,
    ]
end

const preferred_gui = Cmd[
    `rofi -dmenu -i -p "🔎"`,
    # what else?
]

function need_gui(stdstreams = [stdout, stdin])
    return !all(isa.(stdstreams, Ref(Base.TTY)))
end

function choose_preferred_command(commands::Vector{<:Union{Function,Cmd}}, default=nothing)
    for cmd in commands
        if !(cmd isa Cmd) || Sys.which(cmd.exec[1]) !== nothing
            return cmd
        end
    end
    return default
end

function choose_interactive_matcher()
    need_gui() && return choose_preferred_command(preferred_gui, preferred_gui[1])
    return choose_preferred_command(vcat(preferred_terminal, preferred_gui), preferred_terminal[1])
end

function matcher_installation_tips(program::AbstractString)
    if program == "peco"
        return """
        See https://github.com/peco/peco for how to install peco.
        """
    elseif program == "rofi"
        msg = """
        See https://github.com/DaveDavenport/rofi for how to install rofi.
        """
    else
        msg = ""
    end
    return """
    $msg
    For terminal usage, `peco` is recommended.
    See https://github.com/peco/peco for how to install peco.
    """
end

function maybe_warn_matcher(cmd::Cmd)
    if Sys.which(cmd.exec[1]) === nothing
        @warn """
        Matcher $(cmd.exec[1]) not installed.
        $(matcher_installation_tips(cmd.exec[1]))
        """
    end
end

function __init__()
    setmatcher!(choose_interactive_matcher(), CONFIG)
end

include("taskmanager.jl")
include("history.jl")
include("return.jl")

end # module