text:
import Libaudio
using HDF5
using FFTW
using Test

function smoothspectrum_test()
    handelpath = joinpath(Libaudio.modulepath(Libaudio), "test/handel.mat")
    y = h5read(handelpath, "y")
    y = y[:,1]
    Fs = h5read(handelpath, "Fs")
    Fs = Fs[1]
    Z = h5read(handelpath, "Z")
    Z = Z[:,1]

    Y = fft(y)
    NFFT = length(y)
    if iseven(NFFT)
        Nout = (NFFT÷2)+1
    else
        Nout = (NFFT+1)÷2
    end
    Y = Y[1:Nout]
    f = ((0:Nout-1)/NFFT)*Fs
    Y = 20log10.(abs.(Y)/NFFT)
    Noct = 3
    Zt = Libaudio.smoothspectrum(Y, collect(f), Noct)
    (Z, Zt)
end
let (Z, Zt) = smoothspectrum_test()
    @test isapprox(maximum(abs.(Z-Zt)), 3e-12, atol=1e-12)
    @info "==== (0.01) smoothspectrum_test ===="
end

function mps_test()
    mpspath = joinpath(Libaudio.modulepath(Libaudio), "test/mps.mat")
    data = h5read(mpspath, "data")
    mpsr = h5read(mpspath, "mps")

    mpsrn = Vector{ComplexF64}(undef,length(mpsr))
    for (i,j) in enumerate(mpsr)
        mpsrn[i] = j.data[1] + (j.data[2] * im)
    end

    mpst = Libaudio.mps(fft(data[:,1]))
    (mpsrn, mpst)
end
let (mpsr, mpst) = mps_test()
    @test  isapprox(maximum(abs.(abs.(mpsr)-abs.(mpst))), 4e-17, atol=1e-17)
    @info "==== (0.02) mps_test ===="
end

function labelnorm_test()
    foo = joinpath(Libaudio.modulepath(Libaudio), "test/alexa_collect.wav")
    lab = joinpath(Libaudio.modulepath(Libaudio), "test/alexa_collect.TextGrid")
    spl1 = Libaudio.labelnorm(foo, lab, foo[1:end-4]*"_norm.wav", 16, "a", 0, 0)
    @info spl1
    spl2 = Libaudio.labellevel(foo[1:end-4]*"_norm.wav", lab)
    @info spl2
    spl3 = Libaudio.labelnorm(foo, lab, foo[1:end-4]*"_norm_rand40.wav", 16, "a", 5, 5, true, 40)
    spl2
end
let spl2 = labelnorm_test()
    @test all(isapprox.(spl2, 106.2, atol=0.5))
    @info "==== (0.1) labelnorm_test ===="
end

function wav2pcm_test()
    foo = joinpath(Libaudio.modulepath(Libaudio), "test/acqua_ieee_male_250ms_10450ms.wav")
    Libaudio.wav2pcm(foo, foo[1:end-4] * "-16.pcm", 16)
    Libaudio.wav2pcm(foo, foo[1:end-4] * "-24.pcm", 24)
    Libaudio.wav2pcm(foo, foo[1:end-4] * "-32.pcm", 32)
    nothing
end
let y = wav2pcm_test()
    # use Adobe Audition to check the results
    @info "==== (0.2) wav2pcm_test ===="
end

function libwav_test()
    filename = joinpath(Libaudio.modulepath(Libaudio), "test/acqua_ieee_male_250ms_10450ms.wav")
    Libaudio.wavmeta(filename)
    r,fs = Libaudio.wavread(filename)
    # m = Libaudio.wavinfo(filename)
    # x = zeros(Float32, m[1] * m[2])
    # Libaudio.wavread!(filename, x)
    # y = permutedims(reshape(x, m[2], m[1]))
    ys,sr,bps = Libaudio.wavread(filename)
    yd,sr,bps = Libaudio.wavread(filename, Float64)

    Libaudio.wavwrite(filename[1:end-4] * "_s32.wav", ys, Int(fs), 32)
    Libaudio.wavwrite(filename[1:end-4] * "_s24.wav", ys, Int(fs), 24)
    Libaudio.wavwrite(filename[1:end-4] * "_s16.wav", ys, Int(fs), 16)
    Libaudio.wavwrite(filename[1:end-4] * "_d32.wav", yd, Int(fs), 32)
    Libaudio.wavwrite(filename[1:end-4] * "_d24.wav", yd, Int(fs), 24)
    Libaudio.wavwrite(filename[1:end-4] * "_d16.wav", yd, Int(fs), 16)
    (r,ys,yd)
end
let (r,ys,yd) = libwav_test()
    @test isapprox(r, ys, atol=1e-7)
    @test isapprox(r, yd, atol=1e-7)
    @info "==== (0.3) libwav_test ===="
end

function wavio_test()
    path = Libaudio.modulepath(Libaudio)
    x,fs,nb = Libaudio.wavread(joinpath(path,"test/acqua_ieee_male_250ms_10450ms.wav"))
    Libaudio.wavwrite(joinpath(path,"test/foo24.wav"), x, Int(fs), 24)
    Libaudio.wavwrite(joinpath(path,"test/foo16.wav"), x, Int(fs), 16)
    Libaudio.wavwrite(joinpath(path,"test/foo32.wav"), x, Int(fs), 32)

    x16,fs,nb = Libaudio.wavread(joinpath(path,"test/foo16.wav"))
    x24,fs,nb = Libaudio.wavread(joinpath(path,"test/foo24.wav"))
    x32,fs,nb = Libaudio.wavread(joinpath(path,"test/foo32.wav"))
    x, x16, x24, x32
end
let (x, x16, x24, x32) = wavio_test()
    @test isapprox(x16, x, atol=1)
    @test isapprox(x24, x, atol=1e-4)
    @test isapprox(x32, x, atol=1e-7)
    @info "==== (1) wavio_test ===="
end

function wavalid_test()
    path = Libaudio.modulepath(Libaudio)
    x = 0.1 * randn(192000,8)
    y = convert(Matrix{Float32},x)
    Libaudio.wavwrite(joinpath(path,"test/valid_d2f.wav"), x, 48000, 32)
    Libaudio.wavwrite(joinpath(path,"test/valid_f2f.wav"), y, 48000, 32)

    xd,fs,nb = Libaudio.wavread(joinpath(path,"test/valid_d2f.wav"), Float64)
    yf1,fs,nb = Libaudio.wavread(joinpath(path,"test/valid_f2f.wav"), Float32)
    yf2,fs,nb = Libaudio.wavread(joinpath(path,"test/valid_f2f.wav"))
    x,xd,y,yf1,yf2
end
let (x,xd,y,yf1,yf2) = wavalid_test()
    @test y == yf1
    @test y == yf2
    @test isapprox(x, xd, atol=1e-5)
    @info "==== (2) wavalid_test ===="
end

function weighting_a_truth_48000()
    AWEIGHT_48kHz_BA = [0.234301792299513 -0.468603584599025 -0.234301792299515 0.937207169198055 -0.234301792299515 -0.468603584599025 0.234301792299512;
                            1.000000000000000 -4.113043408775872 6.553121752655049 -4.990849294163383 1.785737302937575 -0.246190595319488 0.011224250033231]'
end

function weighting_a_truth_16000()
    AWEIGHT_16kHz_BA = [0.531489829823557 -1.062979659647115 -0.531489829823556 2.125959319294230 -0.531489829823558 -1.062979659647116 0.531489829823559;
    1.000000000000000 -2.867832572992163  2.221144410202311 0.455268334788664 -0.983386863616284 0.055929941424134 0.118878103828561]'
end

let (b, a) = Libaudio.weighting_a(48000), truth = weighting_a_truth_48000()
    @test isapprox(b, truth[:,1], atol = 1e-14)
    @test isapprox(a, truth[:,2], atol = 1e-14)
    @info "==== (3.1) A weighting ===="
end
let (b, a) = Libaudio.weighting_a(16000), truth = weighting_a_truth_16000()
    @test isapprox(b, truth[:,1], atol = 1e-14)
    @test isapprox(a, truth[:,2], atol = 1e-14)
    @info "==== (3.2) A weighting ===="
end

function conv_truth()
    convert(Vector{Float64}, Libaudio.conv_naive(collect(1.0:1.0:4.0), collect(4.0:-1.0:-4.0)))
end
let y = Libaudio.conv(collect(1.0:1.0:4.0), collect(4.0:-1.0:-4.0)), truth = conv_truth()
    @test isapprox(y, truth, atol = 1e-14)
    @info "==== (4) conv ===="
end

function xcorr_truth()
    # octave
    y = [-4.0, -11.0, -20.0, -30.0, -20.0, -10.0, -0.0, 10.0, 20.0, 25.0, 24.0, 16.0, 0.0, -0.0, -0.0, 0.0, 0.0]
end
let y = Libaudio.xcorr(collect(1.0:1.0:4.0), collect(4.0:-1.0:-4.0)), truth = xcorr_truth()
    @test isapprox(y, truth, atol = 1e-14)
    @info "==== (5) xcorr ===="
end

function xcorrcoeff_test()
    x = randn(48000)
    s = x[24000:32000]
    r64 = Libaudio.xcorr(s, x, true)
    t64 = Libaudio.xcorrcoeff_threaded(s, x)

    x = randn(Float32,48000)
    s = x[24000:32000]
    r32 = Libaudio.xcorr(s, x, true)
    t32 = Libaudio.xcorrcoeff_threaded(s, x)

    (r64, t64, r32, t32)
end
let (r64, t64, r32, t32) = xcorrcoeff_test()
    (m64r, i64r) = findmax(r64)
    (m64t, i64t) = findmax(t64)
    (m32r, i32r) = findmax(r32)
    (m32t, i32t) = findmax(t32)
    @info m64r, m64t, m32r, m32t
    @info i64r, i64t, i32r, i32t
    @test i64r == 24001
    @test i64t == 24001
    @test i32r == 24001
    @test i32t == 24001
    @info "==== (6) xcorrcoeff_test ===="
end

function filt_truth()
    y = [3.150000000000000e+01, -1.675000000000000e+01, 5.125000000000000e+00, 5.031250000000000e+01, -4.271875000000000e+01, 3.453125000000000e+00, 1.377578125000000e+02, -1.690429687500000e+02]
    [y y y y]
end
let b = [7,6,5], a = [2,3,4], x = Float32[9,1,5,2,7,4,8,3], truth = filt_truth()
    y = Libaudio.filt(b,a,[x x x x])
    @test isapprox(y, truth, atol = 1e-17)
    @info "==== (7) filt ===="
end

function hamming_truth_periodic()
    w = [8.000000000000002e-02, 2.147308806541882e-01, 5.400000000000000e-01, 8.652691193458119e-01, 1.000000000000000e+00, 8.652691193458120e-01, 5.400000000000001e-01, 2.147308806541882e-01]
end
function hamming_truth_symmetric()
    w = [8.000000000000002e-02, 2.531946911449826e-01, 6.423596296199047e-01, 9.544456792351128e-01, 9.544456792351128e-01, 6.423596296199048e-01, 2.531946911449827e-01, 8.000000000000002e-02]
end
function hann_truth_periodic()
    w = [0.000000000000000e+00, 1.464466094067262e-01, 4.999999999999999e-01, 8.535533905932737e-01, 1.000000000000000e+00, 8.535533905932738e-01, 5.000000000000001e-01, 1.464466094067263e-01]
end
function hann_truth_symmetric()
    w = [0.000000000000000e+00, 1.882550990706332e-01, 6.112604669781572e-01, 9.504844339512095e-01, 9.504844339512095e-01, 6.112604669781573e-01, 1.882550990706333e-01, 0.000000000000000e+00]
end

let w = Libaudio.hamming(8), truth = hamming_truth_symmetric()
    @test isapprox(w, truth, atol = 1e-15)
    @info "==== (8.1) hamming symmetric ===="
end
let w = Libaudio.hamming(8, "periodic"), truth = hamming_truth_periodic()
    @test isapprox(w, truth, atol = 1e-15)
    @info "==== (8.2) hamming periodic ===="
end
let w = Libaudio.hann(8), truth = hann_truth_symmetric()
    @test isapprox(w, truth, atol = 1e-15)
    @info "==== (8.3) hann symmetric ===="
end
let w = Libaudio.hann(8, "periodic"), truth = hann_truth_periodic()
    @test isapprox(w, truth, atol = 1e-15)
    @info "==== (8.4) hann periodic ===="
end

function zeropend_truth()
    ff = ([1, 2, 3, 4, 5, 6, 7, 8, 9], 3)
    tf = (Float64[0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], 4)
    ft = (Float32[1, 2, 3, 4, 5, 6, 7, 8, 9, 0], 4)
    tt = (ComplexF64[0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0], 5)
    (ff, tf, ft, tt)
end
let (ff, tf, ft, tt) = zeropend_truth()
    @test ff == Libaudio.zeropend(collect(1:9), Libaudio.WindowFrame(4,2), false, false)
    @test tf == Libaudio.zeropend(collect(1.:9.), Libaudio.WindowFrame(4,2), true, false)
    @test ft == Libaudio.zeropend(collect(1.f0:9.f0), Libaudio.WindowFrame(4,2), false, true)
    @test tt == Libaudio.zeropend(ComplexF64.(1:9), Libaudio.WindowFrame(4,2), true, true)
    @info "==== (9) zeropend ===="
end

function buffer_truth()
    ff = (Complex{Float64}[1.0+0.0im 3.0+0.0im 5.0+0.0im; 2.0+0.0im 4.0+0.0im 6.0+0.0im; 3.0+0.0im 5.0+0.0im 7.0+0.0im; 4.0+0.0im 6.0+0.0im 8.0+0.0im], 8)
    tf = ([0.0 1.0 3.0 5.0; 0.0 2.0 4.0 6.0; 1.0 3.0 5.0 7.0; 2.0 4.0 6.0 8.0], 10)
    ft = ([1 3 5 7; 2 4 6 8; 3 5 7 9; 4 6 8 0], 10)
    tt = (Float32[0 1 3 5 7; 0 2 4 6 8; 1 3 5 7 9; 2 4 6 8 0], 12)
    (ff, tf, ft, tt)
end
let (ff, tf, ft, tt) = buffer_truth()
    @test ff == Libaudio.buffer(ComplexF64.(1.:9.), Libaudio.WindowFrame(4,2), false, false)
    @test tf == Libaudio.buffer(collect(1.:9.), Libaudio.WindowFrame(4,2), true, false)
    @test ft == Libaudio.buffer(collect(1:9), Libaudio.WindowFrame(4,2), false, true)
    @test tt == Libaudio.buffer(collect(1.f0:9.f0), Libaudio.WindowFrame(4,2), true, true)
    @info "==== (10) buffer ===="
end

function spectrogram_truth()
    # octave:
    # y = buffer(1:100, 16, 4)
    # fft(y)
    y = [
    7.800000000000000e+01 + 0.000000000000000e+00im   2.640000000000000e+02 + 0.000000000000000e+00im  4.560000000000000e+02 + 0.000000000000000e+00im  6.480000000000000e+02 + 0.000000000000000e+00im   8.400000000000000e+02 + 0.000000000000000e+00im   1.032000000000000e+03 + 0.000000000000000e+00im  1.224000000000000e+03 + 0.000000000000000e+00im   1.416000000000000e+03 + 0.000000000000000e+00im   7.720000000000000e+02 + 0.000000000000000e+00im;
    -2.445134153790879e+00 + 3.874624229109006e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im  -8.000000000000000e+00 + 4.021871593700678e+01im   8.386292881545590e+01 - 4.876519307362072e+02im;
    -3.585786437626905e+00 + 1.689949493661167e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -8.000000000000000e+00 + 1.931370849898476e+01im  -4.000000000000000e+00 + 9.656854249492380e+00im;
    -4.941739916456369e+00 + 9.417980255114790e+00im  -8.000000000000000e+00 + 1.197284610132391e+01im  -8.000000000000000e+00 + 1.197284610132391e+01im  -8.000000000000000e+00 + 1.197284610132391e+01im  -8.000000000000000e+00 + 1.197284610132391e+01im  -8.000000000000000e+00 + 1.197284610132391e+01im  -8.000000000000000e+00 + 1.197284610132391e+01im  -8.000000000000000e+00 + 1.197284610132391e+01im   9.538008559557824e+01 - 1.451707589785524e+02im;
    -6.000000000000000e+00 + 6.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -8.000000000000000e+00 + 8.000000000000000e+00im  -4.000000000000000e+00 + 4.000000000000000e+00im;
    -6.472473645916727e+00 + 4.204776819518363e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im  -8.000000000000000e+00 + 5.345429103354389e+00im   9.627676865391413e+01 - 6.481332787817199e+01im;
    -6.414213562373095e+00 + 2.899494936611665e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -8.000000000000000e+00 + 3.313708498984761e+00im  -4.000000000000000e+00 + 1.656854249492381e+00im;
    -6.140652283836026e+00 + 1.533038855493633e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im  -8.000000000000000e+00 + 1.591298939037266e+00im   9.648021693505173e+01 - 1.929449963582681e+01im;
    -6.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -8.000000000000000e+00 + 0.000000000000000e+00im  -4.000000000000000e+00 + 0.000000000000000e+00im
    ]
end
let y = spectrogram_truth()
    x = Libaudio.spectrogram(collect(1.:100.), Libaudio.WindowFrame(16,12), true, true)
    @test isapprox(x[1], y, atol=1e-13)
    # note that we can: isapprox([1.0, 2.0], [1.0+eps(Float32), 2.0+eps(Float32)], atol = 1e-6) -> true
    @info "==== (11) spectrogram ===="
end

function powerspectrum_truth()
    # octave abs(fft(buffer(1:100, 16, 4))).^2 / 16
    y = [
        3.802500000000000e+02   4.356000000000000e+03   1.299600000000000e+04   2.624400000000000e+04   4.410000000000000e+04   6.656400000000000e+04  9.363600000000000e+04   1.253160000000000e+05   3.724900000000000e+04;
        9.420312329436815e+01   1.050965694763527e+02   1.050965694763527e+02   1.050965694763527e+02   1.050965694763527e+02   1.050965694763527e+02  1.050965694763527e+02   1.050965694763527e+02   1.530233727376606e+04;
        1.865317459305202e+01   2.731370849898476e+01   2.731370849898476e+01   2.731370849898476e+01   2.731370849898476e+01   2.731370849898476e+01  2.731370849898476e+01   2.731370849898476e+01   6.828427124746191e+00;
        7.069946592976890e+00   1.295931523537420e+01   1.295931523537420e+01   1.295931523537420e+01   1.295931523537420e+01   1.295931523537420e+01  1.295931523537420e+01   1.295931523537420e+01   1.885744374414300e+03;
        4.499999999999999e+00   8.000000000000002e+00   8.000000000000002e+00   8.000000000000002e+00   8.000000000000002e+00   8.000000000000002e+00  8.000000000000002e+00   8.000000000000002e+00   2.000000000000000e+00;
        3.723316449940346e+00   5.785850768686756e+00   5.785850768686756e+00   5.785850768686756e+00   5.785850768686756e+00   5.785850768686756e+00  5.785850768686756e+00   5.785850768686756e+00   8.418739783176704e+02;
        3.096825406947977e+00   4.686291501015240e+00   4.686291501015240e+00   4.686291501015240e+00   4.686291501015240e+00   4.686291501015240e+00  4.686291501015240e+00   4.686291501015240e+00   1.171572875253810e+00;
        2.503613662714615e+00   4.158264519586321e+00   4.158264519586321e+00   4.158264519586321e+00   4.158264519586321e+00   4.158264519586321e+00  4.158264519586321e+00   4.158264519586321e+00   6.050443735019726e+02;
        2.250000000000000e+00   4.000000000000000e+00   4.000000000000000e+00   4.000000000000000e+00   4.000000000000000e+00   4.000000000000000e+00  4.000000000000000e+00   4.000000000000000e+00   1.000000000000000e+00
    ]
end
let y = powerspectrum_truth()
    x = Libaudio.powerspectrum(collect(1.0:100.0), Libaudio.WindowFrame(16,12), true, true)
    @test isapprox(x[1], y, atol=1e-11)
    @info "==== (12) powerspectrum ===="
end

"""
in this extreme case T=Float32 is insufficient in accuracy but I doubt for normal audio signals
single precision would be enough?
"""
function extractsymbol_test(T)
    m = 10000
    n = 500
    a = randn(T,5m+n)
    t = a[5m+1:5m+n]
    flag1, lbs1, pk1, pkf1, m1 = Libaudio.extractsymbol([a[1:m]; T(0.05)*t; a[m+1:2m]; T(0.001)*t; a[2m+1:3m]; T(0.1)*t; a[3m+1:4m]; T(0.01)*t; a[4m+1:5m]], t, 4, verbose=true, normcoeff=true)
    flag2, lbs2, pk2, pkf2, m2 = Libaudio.extractsymbol([a[1:m]; t; a[m+1:2m]; t; a[2m+1:3m]; t; a[3m+1:4m]; t; a[4m+1:5m]], t, 4, verbose=true, normcoeff=true)
    # @info "extract symbol test: xcorrcoeff enabled" lbs pk pkf
    # @info "extract symbol test: xcorrcoeff disabled" lbsx pkx pkfx
    (flag1 && flag2, lbs1, pk1, pkf1, m1, lbs2, pk2, pkf2, m2)
end
let (sane, lbs, pk, pkf, m, lbsx, pkx, pkfx, mx) = extractsymbol_test(Float64), (sane32, lbs32, pk32, pkf32, m32, lbsx32, pkx32, pkfx32, mx32) = extractsymbol_test(Float32)
    @info sane, sane32
    p = 10000
    q = 500
    @test lbs[1] == 1+p
    @test lbs[2] == lbs[1]-1+(q+p)+1
    @test lbs[3] == lbs[2]-1+(q+p)+1
    @test lbs[4] == lbs[3]-1+(q+p)+1
    @test lbsx[1] == 1+p
    @test lbsx[2] == lbsx[1]-1+(q+p)+1
    @test lbsx[3] == lbsx[2]-1+(q+p)+1
    @test lbsx[4] == lbsx[3]-1+(q+p)+1

    @test lbs32[1] == 1+p
    @test lbs32[2] == lbs32[1]-1+(q+p)+1
    @test lbs32[3] == lbs32[2]-1+(q+p)+1
    @test lbs32[4] == lbs32[3]-1+(q+p)+1
    @test lbsx32[1] == 1+p
    @test lbsx32[2] == lbsx32[1]-1+(q+p)+1
    @test lbsx32[3] == lbsx32[2]-1+(q+p)+1
    @test lbsx32[4] == lbsx32[3]-1+(q+p)+1
    @info "==== (13.1) extractsymbol_test ===="
end

function spl_test()
    calib250, fs = Libaudio.wavread(joinpath(Libaudio.modulepath(Libaudio), "test\\2018-06-22T15-37-20-913+42AA_114.0_105.4_26XX_12AA_0dB_UFX.wav"))
    calib1000, fs = Libaudio.wavread(joinpath(Libaudio.modulepath(Libaudio), "test\\2018-06-22T15-38-53-76+42AB_114.0__26XX_12AA_0dB_UFX.wav"))
    sp, fs = Libaudio.wavread(joinpath(Libaudio.modulepath(Libaudio), "test\\acqua_ieee_male_250ms_10450ms.wav"))

    s1, aanw = Libaudio.spl(calib250[:,1], sp[:,1:1], sp[:,1], 1, Libaudio.WindowFrame(fs,16384,16384÷4), 0.3, 10.3, 100, 12000, 114.0)
    s2, abnw = Libaudio.spl(calib1000[:,1], sp[:,1:1], sp[:,1], 1, Libaudio.WindowFrame(fs,16384,16384÷4), 0.3, 10.3, 100, 12000, 114.0)
    s3, aaaw = Libaudio.spl(calib250[:,1], sp[:,2:2], sp[:,2], 1, Libaudio.WindowFrame(fs,16384,16384÷4), 0.3, 10.3, 100, 12000, 105.4, weighting="A")
    s4, abaw = Libaudio.spl(calib1000[:,1], sp[:,2:2], sp[:,2], 1, Libaudio.WindowFrame(fs,16384,16384÷4), 0.3, 10.3, 100, 12000, 105.4, weighting="A")
    return (s1 && s2 && s3 && s4, aanw, abnw, aaaw, abaw)
end
let (sane, aanw, abnw, aaaw, abaw) = spl_test()
    # @info "spl test results" aanw abnw aaaw abaw
    @info sane
    @test abs(aanw[1] - abnw[1]) < 0.5
    @test isapprox(aanw[1], 128.166, atol=1e-3)
    @test isapprox(abnw[1], 127.993, atol=1e-3)
    @test isapprox(aaaw[1], 125.182, atol=1e-3)
    @test isapprox(abaw[1], 116.33, atol=1e-3)
    @info "==== (14) spl_test ===="
end

function impulse_test()
    ndecay = 1000
    a = Libaudio.expsinesweep(100, 8000, 1, 16000)
    r = [zeros(BigFloat,17); a; zeros(BigFloat, ndecay-17)]
    r = [r r]
    f,h,d,m = Libaudio.impulse(a, ndecay, 100, 8000, 16000, r)

    u,q = findmax(view(d,:,1))
    v,p = findmax(view(m,:,1))
    (u,v,p,q)
end
let (u,v,p,q) = impulse_test()
    @test p-q == 17
    @test isapprox(u, v, atol=1e-15)
    @info "==== (15) impulse_test ===="
end

function symbol_test()
    fs = 48000
    x = randn(8192, 4)
    s = Libaudio.symbol_expsinesweep(800, 2000, 1, fs)
    se = Libaudio.encode_syncsymbol(0.5, s, 0.1, x, fs, 2, -6)
    Libaudio.decode_syncsymbol(se[:,2:2], s, 0.1, size(x,1)/fs, fs)
end
let (sane, loc) = symbol_test()
    @info sane
    @test loc[1] == convert(Int, 0.5 * 48000 + length(Libaudio.symbol_expsinesweep(800, 2000, 1, 48000)) + 0.1 * 48000) + 1
    @info "==== (16) symbol_test ===="
end

function list_test()
    mp = Libaudio.modulepath(Libaudio)
    shallowfolders = Libaudio.list(mp)
    deepfolders = Libaudio.list(mp, recursive=true)
    shallowfiles = Libaudio.list(mp, ".jl")
    deepfiles = Libaudio.list(mp, ".jl", recursive=true)

    # @info shallowfolders deepfolders shallowfiles deepfiles
    (shallowfolders, deepfolders, shallowfiles, deepfiles)
end
let (shallowfolders, deepfolders, shallowfiles, deepfiles) = list_test()
    @test length(shallowfolders) ≥ 3
    @test length(deepfolders) ≥ 5
    @test length(shallowfiles) == 0
    @test length(deepfiles) == 2
    @info "==== (17) list_test ===="
end

function checksum_test()
    mp = Libaudio.modulepath(Libaudio)
    Libaudio.checksum_randinit(mp, "jl.sha")
    before = Libaudio.checksum_validate(mp, "jl.sha", ".jl", recursive=true)
    Libaudio.checksum_update(mp, "jl.sha", ".jl", recursive=true)
    after = Libaudio.checksum_validate(mp, "jl.sha", ".jl", recursive=true)
    (before, after)
end
let (before, after) = checksum_test()
    @test before == false
    @test after == true
    rm(joinpath(Libaudio.modulepath(Libaudio),"jl.sha"))
    @info "==== (18) checksum_test ===="
end

function sigdistratio_test()
    x = randn(8192)
    sdr = Libaudio.sigdistratio([zeros(1024); x; zeros(1024)], view(x,:))
    # @info "signal to distortion ratio (dB)" sdr
    return sdr
end
let sdr = sigdistratio_test()
    @test sdr > 160
    @info "==== (19) sigdistratio_test ===="
end

function resample_test()
    mp = Libaudio.modulepath(Libaudio)
    x, fs = Libaudio.wavread(joinpath(mp, "test/acqua_ieee_male_250ms_10450ms.wav"), Float64)
    y = Libaudio.resample_vhq(x, fs, 16000)
    # wavwrite(y, joinpath(mp, "test/acqua_ieee_male_250ms_10450ms_16khz.wav"), Fs=16000, nbits=32)
    z = Libaudio.resample_vhq(y, 16000, 48000)
    # wavwrite(z, joinpath(mp, "test/acqua_ieee_male_250ms_10450ms_48khz.wav"), Fs=48000, nbits=32)
    # set_zero_subnormals(true)
    sdr = Libaudio.sigdistratio([zeros(1024);view(z,:,1);zeros(1024)], view(x,:,1))
    # set_zero_subnormals(false)
    # @info "resample sdr value" sdr
    return sdr
end
let sdr = resample_test()
    @test sdr > 26.0
    @info "==== (20) resample_test ===="
end

function idealmask_test()
    path = Libaudio.modulepath(Libaudio)
    x1,fs = Libaudio.wavread(joinpath(path, "test/female_with_1s_silence_48khz.wav"), Float64)
    x2,fs = Libaudio.wavread(joinpath(path, "test/acqua_ieee_male_250ms_10450ms.wav"), Float64)
    binarymask, y1, y2 = Libaudio.idealsoftmask(x1[:,1], x2[:,1], fs)
    Libaudio.wavwrite(joinpath(path, "test/idealmask.wav"), [y1 y2], fs, 32)
    sdr1 = Libaudio.sigdistratio([zeros(8000);y1[:,1];zeros(8000)], x1[:,1])
    sdr2 = Libaudio.sigdistratio([zeros(8000);y2[:,1];zeros(8000)], x2[:,1])
    (sdr1, sdr2)
end
let (sdr1, sdr2) = idealmask_test()
    # @info "ideal mask test" sdr1 sdr2
    @test sdr1 > 13
    @test sdr2 > 11
    @info "==== (21) idealmask_test ===="
end


text:
# Generated by running
# julia --color=yes generate_buildjl.jl O/OpenBLAS/OpenBLAS32@0.3.9/build_tarballs.jl JuliaBinaryWrappers/OpenBLAS32_jll.jl OpenBLAS32-v0.3.9+2
# in the root the the source tree of https://github.com/JuliaPackaging/Yggdrasil/
# by first replacing the `include` by its content, see https://github.com/JuliaPackaging/Yggdrasil/issues/858
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libopenblas"], :libopenblas),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/OpenBLAS32_jll.jl/releases/download/OpenBLAS32-v0.3.9+2"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.aarch64-linux-gnu-libgfortran3.tar.gz", "70c72d420d3358c885f00ab57b48acbd1af76b3d575fd77071e36fbbaa00892e"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.aarch64-linux-gnu-libgfortran4.tar.gz", "5046fbf88664004298b5033782365ec2481b012c08e501d81aa88f341b00d6d7"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.aarch64-linux-gnu-libgfortran5.tar.gz", "ee396e4681f2582bb1af846399946bc877771e7703b51ded58d943eddac5efc6"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.aarch64-linux-musl-libgfortran3.tar.gz", "b79b8c5b473885c7e4988d1a62b52e90d9f5034c8000b39e1b831a10ba753a5c"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.aarch64-linux-musl-libgfortran4.tar.gz", "9931954d062b0dd603c6b4bb049672404726bdbd4969b22068400accbbd233a4"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.aarch64-linux-musl-libgfortran5.tar.gz", "376f594def5d788f2d235b02f6f54fe7e9956e893ff3a0a42e12abb0ac1cff15"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.armv7l-linux-gnueabihf-libgfortran3.tar.gz", "edd34511fe5876f47883eaec62bc2e7b220d9396c32c924e5496a6cbfb86a0ff"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.armv7l-linux-gnueabihf-libgfortran4.tar.gz", "ad10a79218b8bc8ccdb442d640bc392d8c1d41a89d49de116819f022525d6cc5"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.armv7l-linux-gnueabihf-libgfortran5.tar.gz", "04c5b704ae04c209d7c0f9506c29bf870760272096a7f106b4161a803fbbfcf6"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.armv7l-linux-musleabihf-libgfortran3.tar.gz", "668f70e1b4b76feee0ee9d51f2b256dfd2d9a05f3c7936d6ed5e3bb16c6f91f6"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.armv7l-linux-musleabihf-libgfortran4.tar.gz", "03b70ad1cb44c8530c66fbab8adfd6d1ca0b7c66f34a6b4cfe8c93c264c5f881"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.armv7l-linux-musleabihf-libgfortran5.tar.gz", "24d3be8016ca46467d8a207334609c761bc6b7d3a07634ca81ddd2be52625a1b"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-linux-gnu-libgfortran3.tar.gz", "41ce6b10e963824d2cbbbbf31e02ba5edc4784d9ce8adf4edfd611042753b472"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-linux-gnu-libgfortran4.tar.gz", "4ba05d858cde46d35120295e2b0f549d666327b25f15892520dda94af3ce29b9"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-linux-gnu-libgfortran5.tar.gz", "1bd5dd9c01cb418d3c1b30b0d02826c7f193883e5a493b38728b81db52f318cb"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-linux-musl-libgfortran3.tar.gz", "fa5ee765cea55afb5de007c4c9fe28db10c75f426a895f39f57363da3e248a71"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-linux-musl-libgfortran4.tar.gz", "50db036850af08df7b241ab25f5e7f668fc0dbb7ac779b899693f386af21fac8"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-linux-musl-libgfortran5.tar.gz", "e38be09ffa8a48010e395a8cf57d0e2348d8147a4724054962baf2750fabe6c9"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-w64-mingw32-libgfortran3.tar.gz", "7673c3c0eed1e48d0bb3771d73a9b192e95af78b134b68cb40cd69e533849a53"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-w64-mingw32-libgfortran4.tar.gz", "eb21686bbc2e880c5c0fe59da595816b427a023422810c2866ed826626763fa7"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.i686-w64-mingw32-libgfortran5.tar.gz", "433ec40a8496efce9680bc6297efa510374f7e9038ef0eaacd314fd7133fa278"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.powerpc64le-linux-gnu-libgfortran3.tar.gz", "772480f50a235b311a2d8e767339cca46665d0ee4aa39abf24f3ec4d4c261f5c"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.powerpc64le-linux-gnu-libgfortran4.tar.gz", "ae3596cb668b02e7c82302c8e7b3e445260a32252b5d85fdc6bd7cddfc7783a0"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.powerpc64le-linux-gnu-libgfortran5.tar.gz", "c94e277724c5ef483775cceac9f16b778c8c4a9a823353dbd4bfaffd137d7a58"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-apple-darwin14-libgfortran3.tar.gz", "bb41f81ec7aae04a773c188fe9683bc6ad47af211e87f3ede1195ab83961f951"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-apple-darwin14-libgfortran4.tar.gz", "99d88d0a7dc55bafe071c61a231fe36279f89317fef196954dbea80c566be96d"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-apple-darwin14-libgfortran5.tar.gz", "56ee6d750897c91f29d701ead5b2f44a2a12a72d5c1eb9b09dc2c0eb3b65d1d1"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-linux-gnu-libgfortran3.tar.gz", "13861d28a1832ce0d5f1b8ecb84de2f19d7970ab9188d07da6da209ac4c4e5d0"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-linux-gnu-libgfortran4.tar.gz", "e2786ab3c111e3cb33f4e4eda50de5ea536aea99edcc952e7d6a87c82b802c36"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-linux-gnu-libgfortran5.tar.gz", "639b8f40d0c8f9177574b413b9e475eb9045abe2e82078c5a4dde4efbfe39a07"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-linux-musl-libgfortran3.tar.gz", "0bba4279e2e524a6504c16b7cc4266621ba09f62880b6847e291955c90ec5c6c"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-linux-musl-libgfortran4.tar.gz", "cfcc99aa3a9124e54225e01299e86e8b7599b0de75835db12d88e643cd1ef83a"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-linux-musl-libgfortran5.tar.gz", "a9dd680a5bcefad112d03993a35d6de434baafc055a71e2f4961d5f8177b983b"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-unknown-freebsd11.1-libgfortran3.tar.gz", "9c877be5d99c2fb5d479e72c5945d076994b720431e11bfc1426ff69c9147183"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-unknown-freebsd11.1-libgfortran4.tar.gz", "90bb169c6fc582243eabf6f2849d3fb53b0b7e17a3e39f4693c004a4f3249f4c"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-unknown-freebsd11.1-libgfortran5.tar.gz", "0e7c2f8125cda30d1e6cae556caea06a7db3b98479119b361826f5973316632f"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-w64-mingw32-libgfortran3.tar.gz", "21949a80fde51fda975a423c893843c7b94557f9286bb73ba079792e33be3c8d"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-w64-mingw32-libgfortran4.tar.gz", "417fe312177ed42a6a95cf4b5bfd39aa202199ca1e210d9a5bd25b622f511944"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/OpenBLAS32.v0.3.9.x86_64-w64-mingw32-libgfortran5.tar.gz", "b945a95923389da1ec7efd76bbfa1337ac525a37581053dc6aed0cf113c5a1d7"),
)

# Install unsatisfied or updated dependencies:
# We added `, isolate=true` as otherwise, it would segfault when closing `OpenBLAS32`,
# probably because it is conflicting with Julia openblas.
unsatisfied = any(!satisfied(p; verbose=verbose, isolate=true) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

ext: jl


text:
function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didOpen")}}, params)
    return DidOpenTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/didOpen")},DidOpenTextDocumentParams}, server)
    server.isrunning = true
    uri = r.params.textDocument.uri
    if URI2(uri) in keys(server.documents)
        doc = server.documents[URI2(uri)]
        doc._content = r.params.textDocument.text
        doc._version = r.params.textDocument.version
        get_line_offsets(doc, true)
    else
        server.documents[URI2(uri)] = Document(uri, r.params.textDocument.text, false, server)
        doc = server.documents[URI2(uri)]
        doc._version = r.params.textDocument.version
        if any(i->startswith(uri, filepath2uri(i)), server.workspaceFolders)
            doc._workspace_file = true
        end
        set_open_in_editor(doc, true)
        if is_ignored(uri, server)
            doc._runlinter = false
        end
    end
    parse_all(doc, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("julia/reloadText")}}, params)
    return DidOpenTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("julia/reloadText")},DidOpenTextDocumentParams}, server)
    server.isrunning = true
    uri = r.params.textDocument.uri
    if URI2(uri) in keys(server.documents)
        doc = server.documents[URI2(uri)]
        doc._content = r.params.textDocument.text
        doc._version = r.params.textDocument.version
    else
        doc = server.documents[URI2(uri)] = Document(uri, r.params.textDocument.text, false, server)
        doc = server.documents[URI2(uri)]
        doc._version = r.params.textDocument.version
        if any(i->startswith(uri, filepath2uri(i)), server.workspaceFolders)
            doc._workspace_file = true
        end
        set_open_in_editor(doc, true)
        if is_ignored(uri, server)
            doc._runlinter = false
        end
    end
    get_line_offsets(doc)
    parse_all(doc, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didClose")}}, params)
    return DidCloseTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/didClose")},DidCloseTextDocumentParams}, server)
    uri = r.params.textDocument.uri
    !haskey(server.documents, URI2(uri)) && return
    doc = server.documents[URI2(uri)]
    empty!(doc.diagnostics)
    publish_diagnostics(doc, server)
    if !is_workspace_file(doc)
        delete!(server.documents, URI2(uri))
    else
        set_open_in_editor(doc, false)
    end
end

function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didSave")}}, params)
    return DidSaveTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/didSave")},DidSaveTextDocumentParams}, server)
    uri = r.params.textDocument.uri
    doc = server.documents[URI2(uri)]
    parse_all(doc, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/willSave")}}, params)
    return WillSaveTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/willSave")},WillSaveTextDocumentParams}, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/willSaveWaitUntil")}}, params)
    return WillSaveTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/willSaveWaitUntil")},WillSaveTextDocumentParams}, server)
    response = JSONRPC.Response(r.id, TextEdit[])
    send(response, server)
end

function JSONRPC.parse_params(::Type{Val{Symbol("textDocument/didChange")}}, params)
    return DidChangeTextDocumentParams(params)
end

function process(r::JSONRPC.Request{Val{Symbol("textDocument/didChange")},DidChangeTextDocumentParams}, server)
    if !haskey(server.documents, URI2(r.params.textDocument.uri))
        server.documents[URI2(r.params.textDocument.uri)] = Document(r.params.textDocument.uri, "", true, server)
    end
    doc = server.documents[URI2(r.params.textDocument.uri)]
    if r.params.textDocument.version < doc._version
        write_transport_layer(server.pipe_out, JSON.json(Dict("jsonrpc" => "2.0", "method" => "julia/getFullText", "params" => doc._uri)), server.debug_mode)
        return
    end
    doc._version = r.params.textDocument.version

    if length(r.params.contentChanges) == 1 && !endswith(doc._uri, ".jmd") && first(r.params.contentChanges).range !== nothing
        tdcce = first(r.params.contentChanges)
        new_cst = _partial_update(doc, tdcce)
        ls_diags = Diagnostic[]
        if server.runlinter && doc._runlinter
            scopepass(getroot(doc))
            mark_errors(doc, ls_diags)
        end
        send(JSONRPC.Request{Val{Symbol("textDocument/publishDiagnostics")},PublishDiagnosticsParams}(nothing, PublishDiagnosticsParams(doc._uri, ls_diags)), server)
    else
        for tdcce in r.params.contentChanges
            applytextdocumentchanges(doc, tdcce)
        end
        parse_all(doc, server)
    end
end

function _partial_update(doc::Document, tdcce::TextDocumentContentChangeEvent)
    cst = getcst(doc)
    # cst = CSTParser.parse(doc._content, true)
    insert_range = get_offset(doc, tdcce.range)
    doc._content = updated_text = edit_string(doc._content, insert_range, tdcce.text)
    doc._line_offsets = nothing

    i1, i2, loc1, loc2 = get_update_area(cst, insert_range)
    is = insert_size(tdcce.text, insert_range)
    if isempty(updated_text)
        empty!(cst.args)
        StaticLint.clear_meta(cst)
        cst.span = cst.fullspan = 0
    elseif 0 < i1 <= i2
        old_span = cst_len(cst, i1, i2)
        ps = ParseState(updated_text, loc1)
        args = EXPR[]
        if i1 == 1 && (ps.nt.kind == CSTParser.Tokens.WHITESPACE || ps.nt.kind == CSTParser.Tokens.COMMENT)
            CSTParser.next(ps)
            push!(args, CSTParser.mLITERAL(ps.nt.startbyte, ps.nt.startbyte, "", Tokens.NOTHING))
        else
            push!(args, CSTParser.parse(ps)[1])
        end
        while ps.nt.startbyte < old_span + loc1 + is && !ps.done
            push!(args, CSTParser.parse(ps)[1])
        end
        new_span = 0
        for i = 1:length(args)
            new_span += args[i].fullspan
        end
        # remove old blocks
        while old_span + is < new_span && i2 < length(cst.args)
            i2+=1
            old_span += cst.args[i2].fullspan
        end
        for i = i1:i2
            StaticLint.clear_meta(cst.args[i])
        end
        deleteat!(cst.args, i1:i2)

        #insert new blocks
        for a in args
            insert!(cst.args, i1, a)
            CSTParser.setparent!(cst.args[i1], cst)
            i1 += 1
        end
    else
        StaticLint.clear_meta(cst)
        cst = CSTParser.parse(updated_text, true)
    end
    CSTParser.update_span!(cst)
    doc.cst = cst
    if doc.cst.typ === CSTParser.FileH
        doc.cst.val = doc.path
        doc.cst.ref = doc
    end
end

insert_size(inserttext, insertrange) = sizeof(inserttext) - max(last(insertrange) - first(insertrange), 0)

function cst_len(x, i1 = 1, i2 = length(x.args))
    n = 0
    @inbounds for i = i1:i2
        n += x.args[i].fullspan
    end
    n
end

function get_update_area(cst, insert_range)
    loc1 = loc2 = 0
    i1 = i2 = 0

    while i1 < length(cst.args)
        i1 += 1
        a = cst.args[i1]
        if loc1 <= first(insert_range) <= loc1 + a.fullspan
            loc2 = loc1
            i2 = i1
            if !(loc1 <= last(insert_range) <= loc1 + a.fullspan)
                while i2 < length(cst.args)
                    i2 += 1
                    a = cst.args[i2]
                    if loc2 <= last(insert_range) <= loc2 + a.fullspan
                        if i2 < length(cst.args) && last(insert_range) <= loc2 + a.fullspan
                            i2+=1
                        end
                        break
                    end
                    loc2 += a.fullspan
                end
            elseif i2 < length(cst.args) && last(insert_range) == loc1 + a.fullspan
                i2 += 1
            end
            break
        end
        loc1 += a.fullspan
    end
    return i1, i2, loc1, loc2
end

function applytextdocumentchanges(doc::Document, tdcce::TextDocumentContentChangeEvent)
    if tdcce.range == tdcce.rangeLength == nothing
        # No range given, replace all text
        doc._content = tdcce.text
    else
        editrange = get_offset(doc, tdcce.range)
        doc._content = edit_string(doc._content, editrange, tdcce.text)
    end
    doc._line_offsets = nothing
end

function edit_string(text, editrange, edit)
    if first(editrange) == last(editrange) == 0
        text = string(edit, text)
    elseif first(editrange) == 0 && last(editrange) == sizeof(text)
        text = edit
    elseif first(editrange) == 0
        text = string(edit, text[nextind(text, last(editrange)):end])
    elseif first(editrange) == last(editrange) == sizeof(text)
        text = string(text, edit)
    elseif last(editrange) == sizeof(text)
        text = string(text[1:first(editrange)], edit)
    else
        text = string(text[1:first(editrange)], edit, text[min(lastindex(text), nextind(text, last(editrange))):end])
    end
end

function check_refs(doc::Document)
    binds, refs = [], []
    for doc1 in doc.server.documents
        if doc1[2].root == doc.root
            StaticLint.collect_bindings_refs(doc1[2].cst, binds, refs)
        end
    end
    urefs = []
    for r in refs
        # if (r.ref isa CSTParser.Binding && r.ref.val isa CSTParser.EXPR) && !(r.ref in binds || r.ref.val in binds)
        if !(any(r.ref == b.binding for b in binds) || r.ref isa SymbolServer.SymStore || r.ref isa String)
            push!(urefs, r)
        end
    end
    binds, refs, urefs
end


text:
using DataFrames

## copied from NanoSoldier.jl
# from submisson.jl
# `x` can only be Expr, Symbol, QuoteNode, T<:Number, or T<:AbstractString
phrase_argument(x::Union{Expr, Symbol, QuoteNode}) = string(x)
phrase_argument(x::Union{AbstractString, Number})  = repr(x)

function parse_submission_string(submission_string)
    fncall = match(r"`.*?`", submission_string).match[2:end-1]
    argind = findfirst(isequal('('), fncall)
    name = fncall[1:(argind - 1)]
    parsed_args = Meta.parse(replace(fncall[argind:end], ";" => ","))
    args, kwargs = Vector{String}(), Dict{Symbol,String}()
    if isa(parsed_args, Expr) && parsed_args.head == :tuple
        started_kwargs = false
        for x in parsed_args.args
            if isa(x, Expr) && (x.head == :kw || x.head == :(=)) && isa(x.args[1], Symbol)
                @assert !haskey(kwargs, x.args[1]) "kwargs must all be unique"
                kwargs[x.args[1]] = phrase_argument(x.args[2])
                started_kwargs = true
            else
                @assert !started_kwargs "kwargs must come after other args"
                push!(args, phrase_argument(x))
            end
        end
    else
        push!(args, phrase_argument(parsed_args))
    end
    return name, args, kwargs
end

## format benchmark results

function Base.convert(::Type{Dict}, t::BenchmarkTools.Trial)
    data = Dict{String, Float64}()
    data["times"] = sum(t.times) / length(t.times)
    data["gctimes"] = sum(t.gctimes) / length(t.gctimes)
    data["memory"] = t.memory
    data["allocs"] = t.allocs
    return data
end

function flatten_benchmark_results(res, data::Dict{String, Dict}, prefix = "")
    for key in keys(res)
        if res[key] isa BenchmarkTools.Trial
            data["$prefix.$key"] = convert(Dict, res[key])
        elseif res[key] isa BenchmarkGroup
            flatten_benchmark_results(res[key], data, "$prefix.$key")
        end
    end
end

function bmresults_to_dataframe(res)
    flatten_res = Dict{String, Dict}()
    flatten_benchmark_results(res, flatten_res)

    names = collect(keys(flatten_res))
    df= DataFrame(Name=names,
                  Times=map(x->flatten_res[x]["times"], names),
                  GCtimes=map(x->flatten_res[x]["gctimes"], names),
                  Memory=map(x->flatten_res[x]["memory"], names),
                  Allocs=map(x->flatten_res[x]["allocs"], names),
                  )
    return df
end


text:
function codegen_residual(model, mass_matrix, dynamics_bias, contact_kinematics, contact_kinematics_jacobians)
    # dimensions
    nz = num_var(model)
    nθ = num_data(model, nf=length(friction_coefficients(model)))

    # residual
    @variables z[1:nz] θ[1:nθ] μ[1:1]
    r = residual(model, mass_matrix, dynamics_bias, contact_kinematics, contact_kinematics_jacobians, z, θ, μ)
    rz = Symbolics.jacobian(r, z)
    rθ = Symbolics.jacobian(r, θ)

    # methods
    r_model = Symbolics.build_function(r, z, θ, μ)[2]
    rz_model = Symbolics.build_function(rz, z, θ)[2]
    rθ_model = Symbolics.build_function(rθ, z, θ)[2]

    return r_model, rz_model, rθ_model
end


text:
# Input Struct

struct Weight
    node::Int64
    weight::Int64
end

struct NodeLabel
    node::Int64
    label::String
    sequence::String
    info::String
end

struct EdgeLabel
    src::String
    dst::String
    overlap::String
    info::String
end

struct Path
    name::String
    path::Array{Int64}
end

struct GFAResult
    g::SimpleDiGraph
    w::Array{Weight}
    l::Array{NodeLabel}
    e::Array{EdgeLabel}
    p::Array{Path}
    weight_inherited::Bool
end

# Graph Component

struct SubGraphMap
    node::Int64
    master_node::Int64
end

mutable struct GraphResult
    g::SimpleDiGraph
    sg_map::Array{SubGraphMap}
    node_label::Array{NodeLabel}
    weight::Array{Weight}
    edge_label::Array{EdgeLabel}
    source_node::Array{NodeLabel}
    sink_node::Array{NodeLabel}
    path::Array{Path}
    weight_inherited::Bool
    GraphResult() = new()
end

mutable struct GraphComponent
    graph::Array{GraphResult}
    lone_node::Array{GraphResult}
    lone_cycle::Array{GraphResult}
    GraphComponent() = new()
end

# Longest Path

mutable struct LongestPath
    graph::GraphResult
    g_opt::SimpleDiGraph
    path::Array{Int64}
    label_path::Array{String}
    obj::Int64
end


text:
using GZip
using Test

##########################
# test_context("GZip tests")
##########################

#for epoch in 1:10

tmp = mktempdir()

test_infile = @__FILE__
test_compressed = joinpath(tmp, "runtests.jl.gz")
test_empty = joinpath(tmp, "empty.jl.gz")

@static if Sys.iswindows()
    gunzip = "gunzip.exe"
elseif Sys.isunix()
    gunzip = "gunzip"
end

test_gunzip = true
try
    run(pipeline(`which $gunzip`, devnull))
catch
    global test_gunzip
    test_gunzip = false
end

try
    #########################
    # test_group("Compress Test1: gzip.jl")
    ##########################

    data = open(x->read(x, String), test_infile);

    first_char = data[1]

    gzfile = gzopen(test_compressed, "wb")
    @test write(gzfile, data) == sizeof(data)
    @test close(gzfile) == Z_OK
    @test close(gzfile) != Z_OK

    @test_throws EOFError write(gzfile, data)

    if test_gunzip
        data2 = read(`$gunzip -c $test_compressed`, String)
        @test data == data2
    end

    data3 = gzopen(x->read(x, String), test_compressed)
    @test data == data3

    # Test gzfdio
    raw_file = open(test_compressed, "r")
    gzfile = gzdopen(fd(raw_file), "r")
    data4 = read(gzfile, String)
    close(gzfile)
    close(raw_file)
    @test data == data4

    # Test peek
    gzfile = gzopen(test_compressed, "r")
    @test peek(gzfile) == UInt(first_char)
    read(gzfile, String)
    @test peek(gzfile) == -1
    close(gzfile)

    # Screw up the file
    raw_file = open(test_compressed, "r+")
    seek(raw_file, 3) # leave the gzip magic 2-byte header
    write(raw_file, zeros(UInt8, 10))
    close(raw_file)

    try
        gzopen(x->read(x, String), test_compressed)
        throw(ErrorException("Expecting ArgumentError or similar"))
    catch ex
        @test typeof(ex) <: Union{ArgumentError,ZError,GZError} ||
              contains(ex.msg, "too many arguments")
    end

    ##########################
    # test_group("gzip file function tests (writing)")
    ##########################
    gzfile = gzopen(test_compressed, "wb")
    write(gzfile, data) == sizeof(data)
    @test flush(gzfile) == Z_OK

    NEW = GZip.GZLIB_VERSION > "1.2.3.9"
    pos = position(gzfile)
    NEW && (pos2 = position(gzfile,true))
    @test_throws ErrorException seek(gzfile, 100)   # can't seek backwards on write
    @test position(gzfile) == pos
    NEW && (@test position(gzfile,true) == pos2)
    @test skip(gzfile, 100)
    @test position(gzfile) == pos + 100
    NEW && (@test position(gzfile,true) == pos2)

    @test_throws MethodError truncate(gzfile, 100)
    @test_throws MethodError seekend(gzfile)

    @test close(gzfile) == Z_OK

    gzopen(test_empty, "w") do io
        a = UInt8[]
        @test gzwrite(io, pointer(a), length(a)*sizeof(eltype(a))) == Int32(0)
    end

    ##########################
    # test_group("gzip file function tests (strategy read/write)")
    ##########################

    # rewrite the test file
    modes = "fhR "
    if GZip.ZLIB_VERSION >= (1,2,5,2)
        modes = "fhRFT "
    end
    for ch in modes
        if ch == ' '
            ch = ""
        end
        for level = 0:9
            gzfile = gzopen(test_compressed, "wb$level$ch")
            @test write(gzfile, data) == sizeof(data)
            @test close(gzfile) == Z_OK

            file_size = filesize(test_compressed)

            #println("wb$level$ch: ", file_size)

            if ch == 'T'
                @test(file_size == sizeof(data))
            elseif level == 0
                @test(file_size > sizeof(data))
            else
                @test(file_size < sizeof(data))
            end

            # readline test
            gzf = gzopen(test_compressed)
            s = IOBuffer()
            while !eof(gzf)
                write(s, readline(gzf))
            end
            data2 = String(take!(s));

            # readuntil test
            seek(gzf, 0)
            while !eof(gzf)
                write(s, readuntil(gzf, 'a'; keep=true))
            end
            data3 = String(take!(s));
            close(gzf)

            @test(data == data2)
            @test(data == data3)

        end
    end

    # Test to create and read an empty file
    gzfile = gzopen(test_compressed, "wb")
    @test write(gzfile, "") == 0
    @test close(gzfile) == Z_OK
    gzfile = gzopen(test_compressed, "r")
    @test eof(gzfile) == true
    @test close(gzfile) == Z_OK

    ##########################
    # test_group("gzip array/matrix tests (write/read)")
    ##########################

    let BUFSIZE = 65536
        for level = 0:3:6
            for T in [Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128,
                      Float32,Float64,ComplexF32,ComplexF64]

                minval = 34567
                try
                    minval = min(typemax(T), 34567)
                catch
                    # do nothing
                end

                # Ordered array
                b = zeros(T, BUFSIZE)
                if !isa(T, Complex)
                    for i = 1:length(b)
                        b[i] = (i-1)%minval;
                    end
                else
                    for i = 1:length(b)
                        b[i] = (i-1)%minval - (minval-(i-1))%minval * im
                    end
                end

                # Random array
                if isa(T, AbstractFloat)
                    r = (T)[rand(BUFSIZE)...];
                elseif isa(T, ComplexF32)
                    r = Int32[rand(BUFSIZE)...] + Int32[rand(BUFSIZE)...] * im
                elseif isa(T, ComplexF64)
                    r = Int64[rand(BUFSIZE)...] + Int64[rand(BUFSIZE)...] * im
                else
                    r = b[rand(1:BUFSIZE, BUFSIZE)];
                end

                # Array file
                b_array_fn = joinpath(tmp, "b_array.raw.gz")
                r_array_fn = joinpath(tmp, "r_array.raw.gz")

                gzaf_b = gzopen(b_array_fn, "w$level")
                write(gzaf_b, b)
                close(gzaf_b)

                #println("$T ($level) ordered: $(filesize(b_array_fn))")

                gzaf_r = gzopen(r_array_fn, "w$level")
                write(gzaf_r, r)
                close(gzaf_r)

                #println("$T ($level) random: $(filesize(r_array_fn))")

                b2 = zeros(T, BUFSIZE)
                r2 = zeros(T, BUFSIZE)

                b2_infile = gzopen(b_array_fn)
                read(b2_infile, b2);
                close(b2_infile)

                r2_infile = gzopen(r_array_fn)
                read(r2_infile, r2);
                close(r2_infile)

                @test b == b2
                @test r == r2
            end
        end
    end
finally
    rm(tmp, recursive=true)
end

#end  # for epoch


text:
include(joinpath("src", "ED-model", "glu_ed_model.jl"))
import .Glu_ED_Model

res_ed = []
for var in [0.0001:0.0005:0.005;]
    push!(res_ed, Glu_ED_Model.glu_ed_model(var))
end


text:
require("ArgParse")
require("Iterators")
require("Requests")

using ArgParse, Iterators, Requests
import JSON
import GZip
import Requests

push!(LOAD_PATH, dirname(@__FILE__()))

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input-file", "-i"
            help = "file name with list of pdb ids to process"
            default = "list_short.txt"
            arg_type = String
            #required = true
        "--pdb-directory", "-d"
            help = "directory name to collect pdb files from rcsb server"
            default = "files/"
            arg_type = String
        "--mesh-size", "-m"
            help = "mesh size for saving structures and group data"
            arg_type = Float64
            default = 0.3
            #required = false
        "--lattice-size", "-l"
            help = "lattice size for saving structures and group data"
            arg_type = Float64
            default = 1.22
            #required = false
        "--threshold", "-t"
            help = "threshold for rotamer's clustering in sidechains' database"
            arg_type = Float64
            default = 1.7
            #required = false
        "--verbose"
            help = "show debug messages while processing data"
            action = :store_true
        "--output-backbone", "-b"
            help = "file in .json format for saving backbone database"
            default = "backbone.json.gz"
            #required = true
        "--output-sidechain", "-s"
            help = "file in .json format for saving sidechain database"
            default = "sidechains.json.gz"
            #required = true
    end

    return parse_args(s)
end

#every aminoacid should have specific number of heavy atoms (there should be more specific check but this would be fine by now)
atomsCounter = {
  "ALA" => 5,
  "ARG" => 11,
  "ASN" => 8,
  "ASP" => 8,
  "CYS" => 6,
  "GLN" => 9,
  "GLU" => 9,
  "GLY" => 4,
  "HIS" => 10,
  "ILE" => 8,
  "LEU" => 8,
  "LYS" => 9,
  "MET" => 8,
  "PHE" => 11,
  "PRO" => 7,
  "SER" => 6,
  "THR" => 7,
  "TRP" => 14,
  "TYR" => 12,
  "VAL" => 7
}

function loadFromRCSB(code :: String, destination_path :: String)
  if length(code) < 4
    return
  end
  pdb_file_name = string(destination_path, code, ".pdb")
  if isfile(pdb_file_name)
    #file already loaded, do nothing
    return
  end
  temp_file_name = string(destination_path, code, ".pdb.gz")
  res = get(string("http://www.rcsb.org/pdb/cgi/export.cgi/", code,
      ".pdb.gz?format=PDB&pdbId=", code, "&compression=gz"))
  if (statuscode(res) != 200)
    println(string("couldn't load PDB file, ID=", code))
    return
  end

  open(temp_file_name, "w") do gzfile
    write(gzfile, Requests.bytes(res))
  end

  gzstream = GZip.gzopen(temp_file_name)
    open(pdb_file_name, "w") do destFile
      write(destFile, readall(gzstream))
    end
  close(gzstream)
  rm(temp_file_name)
end

type PDBAtomInfo
  serial :: Int
  atom :: String
  altLoc :: Char
  resName :: String
  chainID :: Char
  resSeq :: Int
  iCode :: Char
  x :: Float64
  y :: Float64
  z :: Float64
  occupancy :: Float64
  tempFactor :: Float64
  segmentID :: String
  element :: String
  charge :: String
end

function parseAtomInfoFromString(line :: String)
  #TODO: check correctness
  #println(line)
  result = PDBAtomInfo(
      int(line[7:12]),
      strip(line[13:16]),
      line[17],
      strip(line[18:20]),
      line[22],
      int(line[23:26]),
      line[27],
      float(strip(line[31:38])),
      float(strip(line[39:46])),
      float(strip(line[47:54])),
      float(strip(line[55:60])),
      float(strip(line[61:66])),
      strip(line[73:76]),
      strip(line[77:78]),
      strip(line[79:80]))
  #println(result)
  result
end

immutable GeometryVector
  coordinates :: Array{Number, 1}
end
GeometryVector() = GeometryVector([0, 0, 0]) #by default is 3-dimensional

GeometryVectorOp2(a :: GeometryVector, b :: GeometryVector, op) = GeometryVector(map(op, zip(a.coordinates, b.coordinates)))

+(a :: GeometryVector, b :: GeometryVector) = GeometryVectorOp2(a, b, x-> x[1] + x[2])
-(a :: GeometryVector, b :: GeometryVector) = GeometryVectorOp2(a, b, x-> x[1] - x[2])
/(a :: GeometryVector, b :: Number) = GeometryVector(map(x -> x / b, a.coordinates))

round2(a :: GeometryVector) = GeometryVector(map(x::Float64->round(x), a.coordinates))
round2(a :: GeometryVector, meshSize :: Float64) = round2(a / meshSize)

*(a :: GeometryVector, b :: Number) = GeometryVector(map(x -> x*b, a.coordinates))
*(a :: GeometryVector, b :: GeometryVector) = sum(map(x -> x[1] * x[2], zip(a.coordinates, b.coordinates)))

len(a :: GeometryVector) = sqrt(a*a)
#println(len(GeometryVector([1,2,3])))

normalize(a :: GeometryVector) = a / len(a)

projection(projected :: GeometryVector, whereToProject :: GeometryVector) = whereToProject * (projected * whereToProject)

function projectToAxes(v :: GeometryVector, x :: GeometryVector, y :: GeometryVector, z :: GeometryVector)
  v_x = v*x
  v_y = (v - x*v_x)*y
  v_z = (v - x*v_x - y*v_y)*z
  GeometryVector([v_x, v_y, v_z])
end

function cross3d(a :: GeometryVector, b :: GeometryVector)
  GeometryVector([
    a.coordinates[2]*b.coordinates[3] - a.coordinates[3]*b.coordinates[2],
    -(a.coordinates[1]*b.coordinates[3] - a.coordinates[3]*b.coordinates[1]),
    a.coordinates[1]*b.coordinates[2] - a.coordinates[2]*b.coordinates[1]
    ])
end
#println(normalize(GeometryVector([1,2,3])))

getVector = a :: PDBAtomInfo -> GeometryVector([a.x, a.y, a.z])

type Rotamer
  atoms :: Dict{String, GeometryVector}
  center :: GeometryVector
end
Rotamer() = Rotamer(Dict{String, GeometryVector}(), GeometryVector([0, 0, 0]))

type RotamerInfo
  representatives :: Array{Rotamer, 1}
  total :: Int
  amounts :: Array{Int, 1}
end

RotamerInfo() = RotamerInfo(Rotamer[], 0, Int[])

function getPDBFileNames(input_file_name :: String, directory :: String = "files/")
  input_file = open(input_file_name, "r")
  result = String[]
  if (!isdir(directory))
    mkdir(directory)
  end
  while !eof(input_file)
    code = rstrip(readline(input_file))
    if !isfile(string(directory, code, ".pdb"))
      loadFromRCSB(code, directory)
    end
    if (length(code) > 0)
      push!(result, string(directory, code, ".pdb"))
    end
  end
  close(input_file)
  result
end

# reads and preprocesses data.
# leaves only heavy atoms, and only aminoacids with full representation
#
function readPDB(input_file_name :: String)
  records = Dict{Int, Dict{Int, Dict{String, PDBAtomInfo}}}()
  recordsAA = Dict{Char, Array{Int, 1}}()
  input_file = open(input_file_name, "r")
  lastTimeChainHasEnded = false
  chainNumber = 0
  lastChainId = 0
  while !eof(input_file)
    s = rstrip(readline(input_file), ['\r','\n'])
    if (length(s) < 3)
      continue
    end
    if s[1:3] == "TER"
      lastTimeChainHasEnded = true
    end
    if (length(s) < 5)
      continue
    end
    if s[1:5] == "MODEL"
      lastTimeChainHasEnded = true
    end
    if s[1:4] == "ATOM"
      atom = parseAtomInfoFromString(s)
      if !(atom.resName in ["ALA", "ARG", "ASN", "ASP",
          "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS",
          "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL"])
          continue
      end
      if (atom.element in ["H"])
        continue #ignore hydrogens
      end
      if (atom.chainID != lastChainId)
        lastTimeChainHasEnded = true
      end
      if lastTimeChainHasEnded
        chainNumber = chainNumber + 1
        lastChainId = atom.chainID
        records[chainNumber] = Dict{Int, Dict{String, PDBAtomInfo}}()
        recordsAA[chainNumber] = Int[]
      end
      if !(atom.resSeq in keys(records[chainNumber]))
        records[chainNumber][atom.resSeq] = Dict{String, PDBAtomInfo}()
        push!(recordsAA[chainNumber], atom.resSeq)
      end
      records[chainNumber][atom.resSeq][atom.atom] = atom
      #push!(records, parseAtomInfoFromString(s))
      lastTimeChainHasEnded = false
    end
  end
  close(input_file)
  (records, recordsAA)
end

#helper for aminoacid validation - to skip aminoacids without backbone atoms, or with some info missing
function validateAA(aminoacid :: Dict{String, PDBAtomInfo})
  if (!haskey(aminoacid, "CA"))
    return false
  end
  name = aminoacid["CA"].resName
  if (length(values(aminoacid)) != atomsCounter[name])
    return false
  end
  for v in values(aminoacid)
    if v.resName != name
      #println(string("got invalid aa record named ", name))
      #println(aminoacid)
      return false
    end
  end
  for backbone_key in ["CA", "C", "N", "O"]
    if !(haskey(aminoacid, backbone_key))
      return false
    end
  end
  #TODO: add check for sidechain atoms for specific aminoacid types
  true
end

function splitToFragments(aminoacidIds :: Array{Int, 1}, records :: Dict{Int, Dict{String, PDBAtomInfo}})
  result = Array{Int, 1}[]
  fragment = Int[]
  lastAA = -1
  for aminoacidId in aminoacidIds
    if (lastAA + 1 != aminoacidId)
      if (length(fragment) > 0)
        push!(result, fragment)
        fragment = Int[]
      end
    end
    if validateAA(records[aminoacidId])
      push!(fragment, aminoacidId)
      lastAA = aminoacidId
    end
  end

  if (length(fragment) > 0)
    push!(result, fragment)
  end
  result
end

function processPDB(records :: Dict{Int, Dict{Int, Dict{String, PDBAtomInfo}}},
                recordsAA ::  Dict{Char, Array{Int, 1}})
  recordFragments = Dict{Int, Array{Array{Int, 1}, 1}}()
  for chain in keys(records)
    recordFragments[chain] = Array{Int, 1}[]
    fragments = splitToFragments(recordsAA[chain], records[chain])
    for fragment in fragments
      if (length(fragment) >= 4)
        push!(recordFragments[chain], fragment)
      end
    end
  end
  recordFragments
end

#this function returns d_{i-1, i+1}, d_{i, i+2}, d_{i-1, i+2}
function calculateDistances(aminoacids)
  d1 = len(getVector(aminoacids[3]["CA"]) - getVector(aminoacids[1]["CA"]))
  d2 = len(getVector(aminoacids[4]["CA"]) - getVector(aminoacids[2]["CA"]))
  v1 = getVector(aminoacids[2]["CA"]) - getVector(aminoacids[1]["CA"])
  v2 = getVector(aminoacids[3]["CA"]) - getVector(aminoacids[2]["CA"])
  v3 = getVector(aminoacids[4]["CA"]) - getVector(aminoacids[3]["CA"])
  d3 = sign(cross3d(v1, v2)*v3) * len(v1 + v2 + v3)
  (d1, d2, d3)
end

function calculateDistancesVect(v1, v2, v3)
  d1 = len(v1 + v2)
  d2 = len(v3 + v2)
  d3 = sign(cross3d(v1, v2)*v3) * len(v1 + v2 + v3)
  (d1, d2, d3)
end

typealias AtomPosition GeometryVector
typealias AminoacidInfo Dict{String, AtomPosition}

function getLocalVectors(v1, v2, v3, aminoacid)
  #get local coordinate system
  #i=2

  vp = v2 + v3
  x = normalize(cross3d(v2, vp))
  y = normalize(cross3d(vp, x))
  z = normalize(cross3d(x, y))
  vectors = AminoacidInfo()
  for (s, e) in [("CA", "C"), ("CA", "N"), ("C", "O")]
    if haskey(aminoacid, e)
      vectors[e] = projectToAxes(
        getVector(aminoacid[e]) - getVector(aminoacid[s]),
        x, y, z)
    else
      vectors[e] = GeometryVector([0, 0, 0])
    end
  end
  sidechain = Rotamer()
  for e in keys(aminoacid) #TODO: check if in maincode "CA" position is included or not
    if !(e in ["C", "N", "O"])
      sidechain.atoms[e] = projectToAxes(
        getVector(aminoacid[e]) - getVector(aminoacid["CA"]),
        x, y, z)
    end
  end
  #println(map(x->GeometryVector([x[1], x[2], x[3]]), values(sidechain.atoms)))
  sidechainLength = length(values(sidechain.atoms))
  if (sidechainLength > 0)
    sidechain.center = sum(values(sidechain.atoms)) / sidechainLength
  end
  (x, y, z, aminoacid["CA"].resName, vectors, sidechain)
end

#there is special reason not to take normal average vector, but to take vector with average position and average length for every chain fragment
function getAverage(chainInfo :: Dict{String, Dict{(Int, Int, Int), Array{AminoacidInfo, 1}}})
    result = Dict{String, Dict{(Int, Int, Int), Dict{String, GeometryVector}}}()
    for (aa, aaInfo) in chainInfo
        result[aa] = Dict{(Int, Int, Int), Dict{String, GeometryVector}}()
        for (distances, positions) in aaInfo
            size = length(positions)
            result[aa][distances] = Dict{String, GeometryVector}()
            averageLengths = Dict{String, Float64}()
            for x in positions
                for (k, v) in x
                    if !haskey(result[aa][distances], k)
                        result[aa][distances][k] = GeometryVector([0, 0, 0])
                        averageLengths[k] = 0
                    end
                    result[aa][distances][k] += v
                    averageLengths[k] += len(v)
                end
            end
            for k in keys(result[aa][distances])
                averageLengths[k] = averageLengths[k] / size
                result[aa][distances][k] = normalize(result[aa][distances][k] / size) * averageLengths[k]
            end
            result[aa][distances]["O"] = result[aa][distances]["O"] + result[aa][distances]["C"] #this line to save mean vector from CA to O, with correct C=O average bond length
        end
    end
    #result2 = Dict{String, Dict{(Int, Int, Int), Dict{String, Array{Number, 1}}}}()
    result2 = Dict{String, Dict{Int, Dict{Int, Dict{Int, Dict{String, GeometryVector}}}}}()
    for (aa, aaInfo) in result
        result2[aa] = Dict{Int, Dict{Int, Dict{Int, Dict{String, GeometryVector}}}}()
        for ((d1, d2, d3), positions) in aaInfo
            if !haskey(result2[aa], d1)
              result2[aa][d1] = Dict{Int, Dict{Int, Dict{String, GeometryVector}}}()
            end
            if !haskey(result2[aa][d1], d2)
              result2[aa][d1][d2] = Dict{Int, Dict{String, GeometryVector}}()
            end
            if !haskey(result2[aa][d1][d2], d3)
              result2[aa][d1][d2][d3] = Dict{String, GeometryVector}()
            end
            for (k, v) in positions
                result2[aa][d1][d2][d3][k] = v
            end
        end
    end
    result2
end

#helper method.returns index of nearest rotamer group. or zero if no one is found
function findRotamerGroup(representatives :: Array{Rotamer, 1}, rotamer :: Rotamer, threshold :: Number = 1.7)
  for i in 1:length(representatives)
    if (len(representatives[i].center - rotamer.center) < threshold)
      return i
    end
  end
  0
end

function buildLibraryFragment(positions :: Array{Rotamer, 1}, meshSize :: Float64, threshold :: Float64)
  destination = RotamerInfo()
  #1. group by number of representatives
  # todo: when select greatest. if 1, take unaffected
  # at first - make it simple, but ugly
  if (length(positions)==0)
    println("got zero-length positions")
  end
  for rotamerGroup in sort(collect(groupby(r-> round2(r.center, meshSize), positions)), by=length, rev=true)
    if (length(rotamerGroup[1].atoms) == 0)
      #println("no atoms in rotamer group found")
      continue
    end
    i = findRotamerGroup(destination.representatives, rotamerGroup[1], threshold)
    if (i == 0)
      # add new rotamer to library
      push!(destination.representatives, rotamerGroup[1]) #add only 1 rotamer, they are the same
      push!(destination.amounts, length(rotamerGroup))
    else
      #add rotamer to existing cluster i
      destination.amounts[i] += length(rotamerGroup)
    end
    destination.total += length(rotamerGroup)
  end
  destination
end

function getRotamerDb(rotamers :: Dict{String, Dict{(Int, Int, Int), Array{Rotamer, 1}}}, meshSize :: Float64, threshold :: Float64)
    result = Dict{String, Dict{(Int, Int, Int), RotamerInfo}}()
    for (aa, aaInfo) in rotamers
        result[aa] = Dict{(Int, Int, Int), RotamerInfo}()
        for (distances, positions) in aaInfo
            size = length(positions)
            rotamerInfo = buildLibraryFragment(positions, meshSize, threshold)
            if (length(rotamerInfo.representatives) > 0)
                result[aa][distances] = rotamerInfo
            end
        end
    end
    result2 = Dict{String, Dict{Int, Dict{Int, Dict{Int, RotamerInfo}}}}()
    for (aa, aaInfo) in result
      result2[aa] = Dict{Int, Dict{Int, Dict{Int, RotamerInfo}}}()
      for (dist, r) in aaInfo
        if !haskey(result2[aa], dist[1])
          result2[aa][dist[1]] = Dict{Int, Dict{Int, RotamerInfo}}()
        end
        if !haskey(result2[aa][dist[1]], dist[2])
          result2[aa][dist[1]][dist[2]] = Dict{Int, RotamerInfo}()
        end
        if !haskey(result2[aa][dist[1]][dist[2]], dist[3])
          result2[aa][dist[1]][dist[2]][dist[3]] = r
        end
      end
    end
    result2
end

function processChainPortion(aminoacids, meshSize = 0.3)
  v1 = getVector(aminoacids[2]["CA"]) - getVector(aminoacids[1]["CA"])
  v2 = getVector(aminoacids[3]["CA"]) - getVector(aminoacids[2]["CA"])
  v3 = getVector(aminoacids[4]["CA"]) - getVector(aminoacids[3]["CA"])
  processChainPortionVec(v1, v2, v3, aminoacid[2], meshSize)
end

function processChainPortionVec(v1, v2, v3, aminoacid, meshSize = 0.3)
  distances = map(x-> convert(Int, round(x/meshSize)), calculateDistancesVect(v1, v2, v3))
  (x, y, z, aa_name, vectors, sidechains) = getLocalVectors(v1, v2, v3, aminoacid)
  (distances, aa_name, vectors, sidechains)
end

function getVectorForSeq(sequence, ks, k, text_file_name, latticeSize = 1.22)
  if (length(ks)<4)
    println("length of keys <=4")
    println(text_file_name)
  end
  for i in [1, 2, 3, k-1, k, k + 1, k + 2, length(ks) - 1, length(ks)]
    if (i < 1 || i > length(ks))
      continue
    end

    if !haskey(sequence[ks[i]], "CA")
      println(sequence[ks[i]])
      println(text_file_name)
    end
  end
  if (k == 1)
    return (round2((getVector(sequence[ks[2]]["CA"]) - getVector(sequence[ks[1]]["CA"]))/latticeSize),
            round2((getVector(sequence[ks[2]]["CA"]) - getVector(sequence[ks[1]]["CA"]))/latticeSize),
            round2((getVector(sequence[ks[3]]["CA"]) - getVector(sequence[ks[2]]["CA"]))/latticeSize))
  end
  if (k == length(ks))
    return (round2((getVector(sequence[ks[length(ks)]]["CA"]) - getVector(sequence[ks[length(ks) - 1]]["CA"]))/latticeSize),
            round2((getVector(sequence[ks[length(ks) - 1]]["CA"]) - getVector(sequence[ks[length(ks) - 2]]["CA"]))/latticeSize),
            round2((getVector(sequence[ks[length(ks)]]["CA"]) - getVector(sequence[ks[length(ks) - 1]]["CA"]))/latticeSize))
  end

  if (k == length(ks) - 1)
    return (round2((getVector(sequence[ks[length(ks) - 1]]["CA"]) - getVector(sequence[ks[length(ks) - 2]]["CA"]))/latticeSize),
            round2((getVector(sequence[ks[length(ks)]]["CA"]) - getVector(sequence[ks[length(ks) - 1]]["CA"]))/latticeSize),
            round2((getVector(sequence[ks[length(ks) - 1]]["CA"]) - getVector(sequence[ks[length(ks) - 2]]["CA"]))/latticeSize))
  end

  return (round2((getVector(sequence[ks[k]]["CA"]) - getVector(sequence[ks[k - 1]]["CA"]))/latticeSize),
          round2((getVector(sequence[ks[k + 1]]["CA"]) - getVector(sequence[ks[k]]["CA"]))/latticeSize),
          round2((getVector(sequence[ks[k + 2]]["CA"]) - getVector(sequence[ks[k + 1]]["CA"]))/latticeSize))
end

function load_atom_info(text_file_name :: String, pdb_dir :: String,
        meshSize :: Float64 = 0.3, latticeSize :: Float64 = 1.22, threshold :: Float64 = 1.7)
  basechainInfo = Dict{String, Dict{(Int, Int, Int), Array{AminoacidInfo, 1}}}()
  sidechainInfo = Dict{String, Dict{(Int, Int, Int), Array{Rotamer, 1}}}()
  pdb_file_names = getPDBFileNames(text_file_name, pdb_dir)
  for pdb_file_name in pdb_file_names
    (atom_infos, atom_info_keys) = readPDB(pdb_file_name)
    chainFragments = processPDB(atom_infos, atom_info_keys)
    for chain in keys(chainFragments)
      for ks in chainFragments[chain]
        if (length(ks) <= 4)
          continue
        end
        for k in 1 : length(ks)
          (v1, v2, v3) = getVectorForSeq(atom_infos[chain], ks, k, pdb_file_name, latticeSize)
          (d, aa, b, s) = processChainPortionVec(v1, v2, v3, atom_infos[chain][ks[k]]) #[atom_infos[chain][i] for i in ks[k - width + 1 : k]])
          if !haskey(basechainInfo, aa)
            basechainInfo[aa] = Dict{(Int, Int, Int), Array{AminoacidInfo, 1}}()
            sidechainInfo[aa] = Dict{(Int, Int, Int), Array{Rotamer, 1}}()
          end
          if !haskey(basechainInfo[aa], d)
            basechainInfo[aa][d] = AminoacidInfo[]
            sidechainInfo[aa][d] = Rotamer[]
          end
          push!(basechainInfo[aa][d], b)
          push!(sidechainInfo[aa][d], s)
        end
      end
    end
  end
  r1 = getAverage(basechainInfo)
  r2 = getRotamerDb(sidechainInfo, meshSize, threshold)
  ( {"data" => r1, "meshSize" => meshSize, "latticeSize" => latticeSize},
    {"data" => r2, "threshold" => threshold, "meshSize"=> meshSize, "latticeSize" => latticeSize})
end

function main()
    parsed_args = parse_commandline()
    (r1, r2) = load_atom_info(parsed_args["input-file"], parsed_args["pdb-directory"], parsed_args["mesh-size"],
        parsed_args["lattice-size"], parsed_args["threshold"])
    output_file = GZip.open(parsed_args["output-backbone"], "w")
    write(output_file, JSON.json(r1))
    close(output_file)
    output_file = GZip.open(parsed_args["output-sidechain"], "w")
    write(output_file, JSON.json(r2))
    close(output_file)
end

main()


text:
using CounterfactualExplanations
using CounterfactualExplanations.Generators

@testset "Construction" begin

    @testset "Generic" begin
        generator = GenericGenerator()
        @test hasfield(GenericGenerator, :loss)
        @test typeof(generator) <: AbstractGradientBasedGenerator
    end

    @testset "Greedy" begin
        generator = GreedyGenerator()
        @test hasfield(GreedyGenerator, :loss)
        @test typeof(generator) <: AbstractGradientBasedGenerator
    end

end

# Data:
using Random
using CounterfactualExplanations.Data
Random.seed!(1234)
xs, ys = Data.toy_data_linear()
X = hcat(xs...)
counterfactual_data = CounterfactualData(X,ys')
x = select_factual(counterfactual_data,rand(1:size(X)[2]))


text:
const BasisSet = Array{CGBF,1}

include("read_basis_g94.jl")
import .g94basis

function build_basis_new( atoms::Atoms )

    data_dict = g94basis.read_basis_set("g94_sto-3g")

    basis = BasisSet()

    Natoms = atoms.Natoms
    atsymbs = atoms.atsymbs

    Nbasis = 0
    for ia = 1:Natoms
        #
        x = atoms.positions[1,ia]
        y = atoms.positions[2,ia]
        z = atoms.positions[3,ia]
        #
        data = data_dict[atsymbs[ia]]
        Nsyms = data.Nsyms
        #
        for isym = 1:Nsyms
            sym = data.syms[isym]
            for (I,J,K) in sym2power[sym]
                cbf = init_CGBF( (x,y,z), (I,J,K) )
                ncontr = data.Ncontr[isym]
                for i = 1:ncontr
                    expn = data.expns[isym][i]
                    coef = data.coeffs[isym][i]
                    push!(cbf,expn,coef)
                end
                push!(basis,cbf)
                Nbasis = Nbasis + 1
            end
        end
    end

    println("Nbasis = ", Nbasis)

    return basis

end

function build_basis( atoms::Atoms, name="sto3g", verbose::Bool=false )

    data = basis_set_data[name]
    basis = BasisSet()

    Natoms = atoms.Natoms
    atsymbs = atoms.atsymbs

    Nbasis = 0
    for ia = 1:Natoms
        #
        atno = ZATOMS[atsymbs[ia]]
        x = atoms.positions[1,ia]
        y = atoms.positions[2,ia]
        z = atoms.positions[3,ia]
        #
        for btuple in data[atno]
            sym, primlist = btuple
            for (I,J,K) in sym2power[sym]
                cbf = init_CGBF( (x,y,z), (I,J,K) )
                for (expn,coef) in primlist
                    push!(cbf,expn,coef)
                end
                push!(basis,cbf)
                Nbasis = Nbasis + 1
            end
        end
    end

    println("Nbasis = ", Nbasis)

    return basis
end

const sym2power = Dict{Any,Any}(
  "S" => [(0,0,0)],
  "P" => [(1,0,0),(0,1,0),(0,0,1)],
  "D" => [(2,0,0),(0,2,0),(0,0,2),(1,1,0),(1,0,1),(0,1,1)],
  "F" => [(3,0,0),(0,3,0),(0,0,3),(2,1,0),(2,0,1),(1,2,0),
          (1,1,1),(1,0,2),(0,2,1),(0,1,2)]
  )