text:
# Gross-Nevue ground state using DMRG

include("Schollwock.jl")
using DelimitedFiles
using Plots
using LsqFit
using DataFrames
using CSV

function Hss(m₀::Float64,g₀::Float64,r::Float64,n::Int;ϵ=0.0)
    a = 1.0/(n+1)
    a⁺ = [0.0 0.0;
          1.0 0.0]
    a⁻ = [0.0 1.0;
          0.0 0.0]
    Z = [1.0 0.0;
         0. -1.0]
    # MagneticField = ϵ*randn(Float64)*Diagonal([1.5, 0.5, -0.5, -1.5])
    # return im*(m₀-r/a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))+g₀^2/(2*a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))^2 + m₀*MagneticField
    return im*(-m₀-r/a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))+g₀^2/(2*a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))^2
end

function Hnn(r::Float64,n::Int)
    a = 1.0/(n+1)
    a⁺ = [0.0 0.0;
          1.0 0.0]
    a⁻ = [0.0 1.0;
          0.0 0.0]
    Z = [1.0 0.0;
         0. -1.0]
    II = [1. 0.;
          0. 1.]
    # return im/(2*a)*(-kron(II,a⁻,Z,a⁻)+kron(a⁻,Z,a⁻,II)+kron(II,a⁺,Z,a⁺)-kron(a⁺,Z,a⁺,II)+r*(kron(II,a⁻,a⁺,II)-kron(II,a⁺,a⁻,II)+kron(a⁺,Z,Z,a⁻)-kron(a⁻,Z,Z,a⁺)))
    return im/(2*a)*(kron(II,a⁻,Z,a⁺)-kron(a⁻,Z,a⁺,II)-kron(II,a⁺,Z,a⁻)+kron(a⁺,Z,a⁻,II)+r*(kron(II,a⁻,a⁺,II)-kron(II,a⁺,a⁻,II)+kron(a⁺,Z,Z,a⁻)-kron(a⁻,Z,Z,a⁺)))
end

#Making a Plot#

# function InnerOfAdjoiningMPS(Le::Int,n::Int,BD::Int,r::Float64,m0::Float64,g0::Float64) #Le: length of the smaller system; n: Number of differing sites; BD: Starting Bond dimension
#     mps_long = DMRG(NNHandSSHtoMPO(Hnn(r,Le+n),Hss(m0,g0,r,Le+n),Le+n),Le+n,D=BD,σ=4,ϵ=1.0e-8,maxiteration=30)
#     mps_short_temp = DMRG(NNHandSSHtoMPO(Hnn(r,Le),Hss(m0,g0,r,Le),Le),Le,D=BD,σ=4,ϵ=1.0e-8,maxiteration=30)
#     mps_cut = CutEntanglement(mps_long,Le+1,after=false)
#     tensors = mps_short_temp.tensors
#     for i in Le+1:Le+n
#         tensors[i] = mps_cut.tensors[i]
#     end
#     mps_short = MPS2(tensors, Le+n,mps_short_temp.BD,Le,4)
#     norm_of_MPS = sqrt(abs(inner(mps_short,mps_short)))
#     mps_short.tensors[Le] = mps_short.tensors[Le]/norm_of_MPS
#     return (mps_long.BD,abs(inner(mps_short,mps_long)))
# end

function AbsInnerProductOfMismatchingMPS(M1::MPS2,M2::MPS2)
    M1.length == M2.length && return inner(M1,M2)
    M1.length < M2.length ? (SM,BM) = (M1,M2) : (SM,BM) = (M2,M1)
    l1 = SM.length
    l2 = BM.length
    M_cut = CutEntanglement(BM,l1+1,after=false)
    tensors = SM.tensors
    for l in l1+1:l2
        tensors[l] = M_cut.tensors[l]
    end
    temp_MPS = MPS2(tensors)
    norm_of_MPS = sqrt(abs(inner(temp_MPS,temp_MPS)))
    temp_MPS.tensors[l2] = temp_MPS.tensors[l2]/norm_of_MPS
    return abs(inner(temp_MPS,BM))
end

function UniformInnerProductOfMPS(M1::MPS2,M2::MPS2)
    M1.length == M2.length && return inner(M1,M2)
    M1.length < M2.length ? (SM,BM) = (M1,M2) : (SM,BM) = (M2,M1)
    l1 = SM.length
    l2 = BM.length
    σ = SM.σ
    tensors = SM.tensors
    for l in l1+1:l2
        tensors[l] = ones(ComplexF64,σ,1,1)
    end
    temp_MPS = normalize(MPS2(tensors))
    return abs(inner(temp_MPS,BM))
end

function AddRandomSite(M::MPS2)
    σ = M.σ
    l = M.length
    tensors = M.tensors
    T = typeof(tensors[1][1])
    temp_tensor = rand(T,(σ,1,1))
    tensors[l+1] = temp_tensor/sum(abs2.(temp_tensor))
    Mt = MPS2(tensors,l+1,M.BD,M.CP,M.σ)
    norm_of_MPS = sqrt(abs(inner(Mt,Mt)))
    Mt.tensors[l+1] = Mt.tensors[l+1]/norm_of_MPS
    return Mt
end

function main1(args)
    l0 = parse(Int,args[1+1])
    l = parse(Int,args[2+1])
    r = parse(Float64,args[3+1])
    m0 = parse(Float64,args[4+1])
    g0 = parse(Float64,args[5+1])
    n = parse(Int,args[6+1])
    ϵ0 = parse(Float64,args[7+1])
    D = parse(Int,args[9])
    DMAX = parse(Int,args[10])
    local inners = Array{Float64,1}()
    t = time()
    local SmallGroundStates = Dict{Int,MPS2}()
    local ϵ_Small = zeros(Float64,n)
    local ϵ_Large = zeros(Float64,n)
    local BigGroundStates = Dict{Int,MPS2}()
    local DF_inners = DataFrame(Length = Int[],
                                ϵSmall = Float64[],
                                BD_Small = Int[],
                                ϵLarge = Float64[],
                                BD_Large = Int[],
                                InnerProduct = Float64[],
                                time = Float64[])
    for k in 1:n
        SmallGroundStates[k],ϵ_Small[k] = DMRG(NNHandSSHtoMPO(Hnn(r,l0),Hss(m0,g0,r,l0),l0),D=D,ϵ=ϵ0,maxiteration=3,DMax=DMAX)
    end
    open("Gross-Neveu_r=$(r)_m0=$(m0)_g0=$(g0)_l=$(l).txt","w") do io
        open("MAX_Gross-Neveu_r=$(r)_m0=$(m0)_g0=$(g0)_l=$(l).txt","w") do ioMAX
            for Le in l0+1:1:l
                t1 = time()
                for k in 1:n
                    BigGroundStates[k],ϵ_Large[k] = DMRG(NNHandSSHtoMPO(Hnn(r,Le),Hss(m0,g0,r,Le),Le),D=SmallGroundStates[k].BD,ϵ=ϵ0,maxiteration=3,DMax=DMAX)
                end
                push!(inners,UniformInnerProductOfMPS(SmallGroundStates[1],BigGroundStates[1]))
                for j in 1:n, k in 1:n
                    tmp = UniformInnerProductOfMPS(SmallGroundStates[j],BigGroundStates[k])
                    inners[end] = max(tmp,inners[end])
                    write(io,"$(Le) \t $(round(tmp,digits=5)) \t $(round(time()-t))\n")
                end
                write(ioMAX,"$(Le) \t $(round(inners[end],digits=5)) \t $(round(time()-t))\n")
                push!(DF_inners,(Le,ϵ_Small[1],SmallGroundStates[1].BD,ϵ_Large[1],BigGroundStates[1].BD,inners[end],round(time()-t1,digits=2)))
                SmallGroundStates = deepcopy(BigGroundStates)
                ϵ_Small = deepcopy(ϵ_Large)
            end
            write(ioMAX,"finished computation in $(time()-t).\n")
            CSV.write("Gross-Neveu_r=$(r)_m0=$(m0)_g0=$(g0)_l=$(l).csv",DF_inners)
        end
        write(io,"finished computation in $(time()-t).\n")
    end
    gr()
    p = plot(l0+1:l,inners)
    png("Gross-Neveu_r=$(r)_m0=$(m0)_g0=$(g0)_l=$(l).png")
end

function main2()
    l = parse(Int,ARGS[1+1])
    r = parse(Float64,ARGS[2+1])
    m0 = parse(Float64,ARGS[3+1])
    mf = parse(Float64,ARGS[4+1])
    dm = parse(Float64,ARGS[5+1])
    g0 = parse(Float64,ARGS[6+1])
    D = 25
    local inners = Array{Float64,1}()
    t = time()
    open("Gross-Neveu_Varying-Mass_r=$(r)_l=$(l)_g0=$(g0).txt","w") do io
        for m in m0:dm:mf
            local SmallGroundStates , _ = DMRG(NNHandSSHtoMPO(Hnn(r,l),Hss(m,g0,r,l),l),D=D,ϵ=1.0e-8,maxiteration=5)
            local BigGroundStates , _ = DMRG(NNHandSSHtoMPO(Hnn(r,l+1),Hss(m,g0,r,l+1),l+1),D=SmallGroundStates.BD,ϵ=1.0e-8,maxiteration=5)
            push!(inners,UniformInnerProductOfMPS(SmallGroundStates,BigGroundStates))
            write(io,"$(m) \t $(round(inners[end],digits=5)) \t $(round(time()-t))\n")
        end
        write(io,"finished computation in $(time()-t).\n")
    end
    gr()
    p = plot(m0:dm:mf,inners)
    png("Gross-Neveu_Varying-Mass_r=$(r)_l=$(l)_g0=$(g0).png")
end

function main3(args)
    l = parse(Int,args[1+1])
    r = parse(Float64,args[2+1])
    m0 = parse(Float64,args[3+1])
    rm = parse(Float64,args[4+1])
    n = parse(Int,args[5+1])
    g0 = parse(Float64,args[6+1])
    ϵ0 = parse(Float64,args[7+1])
    D = parse(Int,args[9])
    DMAX = parse(Int,args[10])
    local inners = Array{Float64,1}()
    local DF_inners = DataFrame(m = Float64[],
                                ϵ_Small = Float64[],
                                BD_Small = Int[],
                                ϵ_Large = Float64[],
                                BD_Large = Int[],
                                Inner = Float64[],
                                time = Int[])
    t = time()
    open("Gross-Neveu_Power-Mass2_r=$(r)_l=$(l)_g0=$(g0).txt","w") do io
        local m = m0
        for i in 0:n
            t1 = time()
            local SmallGroundStates,ϵsmall = DMRG(NNHandSSHtoMPO(Hnn(r,l),Hss(m,g0,r,l),l),D=D,ϵ=ϵ0,maxiteration=3,DMax=DMAX)
            local BigGroundStates,ϵbig = DMRG(NNHandSSHtoMPO(Hnn(r,l+1),Hss(m,g0,r,l+1),l+1),D=SmallGroundStates.BD,ϵ=ϵ0,maxiteration=3,DMax=DMAX)
            push!(inners,UniformInnerProductOfMPS(SmallGroundStates,BigGroundStates))
            write(io,"$(m) \t $(round(inners[end],digits=10)) \t $(round(time()-t))\n")
            push!(DF_inners,(m,ϵsmall,SmallGroundStates.BD,ϵbig,BigGroundStates.BD,inners[end],round(Int,time()-t1)))
            m = round(m*rm,digits = 10)
        end
        write(io,"finished computation in $(time()-t).\n")
        CSV.write("Gross-Neveu_Power-Mass2_r=$(r)_l=$(l)_g0=$(g0).csv",DF_inners)
    end
    gr()
    p = plot([m0*rm^i for i in 0:n],inners)
    png("Gross-Neveu_Power-Mass2_r=$(r)_l=$(l)_g0=$(g0).png")
end

function twopoint(j::Int,k::Int) #unnormalized
    Lm = [0.0 0.0;
          1.0 0.0]
    Lp = [0.0 1.0;
          0.0 0.0]
    Z = [1.0 0.0;
        0.0 -1.0]
    LmZ = kron(Lm,Z)
    ZZ = kron(Z,Z)
    ZLp = kron(Z,Lp)
    local O = Dict{Int,Array{Float64,2}}()
    O[j] = LmZ
    O[k] = ZLp
    for s in j+1:k-1
        O[s] = ZZ
    end
    return O
end

function main4(args)
    l = parse(Int,args[1+1])
    r = parse(Float64,args[2+1])
    m0 = parse(Float64,args[3+1])
    rm = parse(Float64,args[4+1]) #multiplicative coefficient for increasing m
    n = parse(Int,args[5+1]) #number of different masses
    g0 = parse(Float64,args[6+1])
    ϵg = parse(Float64,args[7+1]) #precision goal
    D = parse(Int,args[9])
    DMAX = parse(Int,args[10])
    m = m0
    local DF_inners = DataFrame(Distance = collect(0:l-1))
    open("Gross-Neveu_twopoint_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg).txt","w") do io
        t = time()
        for i in 1:n
            exps = Float64[]
            exps2 = ComplexF64[]
            local GS,ϵf = DMRG(NNHandSSHtoMPO(Hnn(r,l),Hss(m,g0,r,l),l),D=D,ϵ=ϵg,maxiteration=3,DMax=DMAX)
            write(io,"precision = $(ϵf)\nBondDimension = $(GS.BD)\nm = $(m)\nTwo point correlators:\n")
            for k in 1:l
                tp = twopoint(1,k)
                ex = (l-1)*1im*expectation(GS,GS,tp)
                write(io,"$(k):\t\t$(ex)\n")
                append!(exps,abs(ex))
                append!(exps2,ex)
            end
            DF_inners[Symbol("ydata_m=",round(m,digits=4))] = exps
            DF_inners[Symbol("ydata_error_m=",round(m,digits=4))] = 2*sqrt(ϵf)*exps
            DF_inners[Symbol("ActualExpectation_m=",round(m,digits=4))] = exps2
            DF_inners[Symbol("BD_m=",round(m,digits=4))] = repeat([GS.BD],l)
            DF_inners[Symbol("mdata_m=",round(m,digits=4))] = repeat([m],l)
            DF_inners[Symbol("ϵ_m=",round(m,digits=4))] = repeat([ϵf],l)
            @. model(x,p) = p[3] + p[2]*exp(-x/p[1])
            write(io,"model(x,p) = p[3] + p[2]*exp(-x/p[1])\n")
            try
                fit = curve_fit(model,0:l-1,exps,[1.0,1.0,0.0])
                write(io,"parameters:\n$(fit.param)\n")
                write(io,"Standard Deviations:\n$(stderror(fit))\n")
            catch
                write(io,"LsqFit Failed\n")
            end
            write(io,"============================================================\n")
            m = m*rm
        end
        CSV.write("Gross-Neveu_twopoint_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg).csv",DF_inners)
        write(io,"finished computation in $(time()-t).")
    end
end

function twopoint2(j::Int,k::Int) #unnormalized
    Lm = [0.0 0.0;
          1.0 0.0]
    I = [1.0 0.0;
         0.0 1.0]
    Lp = [0.0 1.0;
          0.0 0.0]
    Z = [1.0 0.0;
        0.0 -1.0]
    ILp = kron(I,-Lp)
    ZZ = kron(Z,Z)
    LmI = kron(Lm,I)
    local O = Dict{Int,Array{Float64,2}}()
    O[j] = ILp
    O[k] = LmI
    for s in j+1:k-1
        O[s] = ZZ
    end
    return O
end

function main5(args)
    l = parse(Int,args[1+1])
    r = parse(Float64,args[2+1])
    m0 = parse(Float64,args[3+1])
    rm = parse(Float64,args[4+1]) #multiplicative coefficient for increasing m
    n = parse(Int,args[5+1]) #number of different masses
    g0 = parse(Float64,args[6+1])
    ϵg = parse(Float64,args[7+1]) #precision goal
    D = parse(Int,args[9])
    DMAX = parse(Int,args[10])
    m = m0
    local DF_inners = DataFrame(Distance = collect(0:l-1))
    open("Gross-Neveu_twopoint2_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg).txt","w") do io
        t = time()
        for i in 1:n
            exps = Float64[]
            exps2 = ComplexF64[]
            local GS,ϵf = DMRG(NNHandSSHtoMPO(Hnn(r,l),Hss(m,g0,r,l),l),D=D,ϵ=ϵg,maxiteration=5,DMax=DMAX)
            write(io,"precision = $(ϵf)\nBondDimension = $(GS.BD)\nm = $(m)\nTwo point correlators 2:\n")
            for k in 1:l
                tp = twopoint2(1,k)
                ex = (l-1)*1im*expectation(GS,GS,tp)
                write(io,"$(k):\t\t$(ex)\n")
                append!(exps,abs(ex))
                append!(exps2,ex)
            end
            DF_inners[Symbol("ydata_m=",round(m,digits=4))] = exps
            DF_inners[Symbol("ydata_error_m=",round(m,digits=4))] = 2*sqrt(ϵf)*exps
            DF_inners[Symbol("ActualExpectation_m=",round(m,digits=4))] = exps2
            DF_inners[Symbol("BD_m=",round(m,digits=4))] = repeat([GS.BD],l)
            DF_inners[Symbol("mdata_m=",round(m,digits=4))] = repeat([m],l)
            DF_inners[Symbol("ϵ_m=",round(m,digits=4))] = repeat([ϵf],l)
            @. model(x,p) = p[3] + p[2]*exp(-x/p[1])
            write(io,"model(x,p) = p[3] + p[2]*exp(-x/p[1])\n")
            try
                fit = curve_fit(model,0:l-1,exps,[1.0,1.0,0.0])
                write(io,"parameters:\n$(fit.param)\n")
                write(io,"Standard Deviations:\n$(stderror(fit))\n")
            catch
                write(io,"LsqFit Failed\n")
            end
            write(io,"============================================================\n")
            m = m*rm
        end
        CSV.write("Gross-Neveu_twopoint2_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg).csv",DF_inners)
        write(io,"finished computation in $(time()-t).")
    end
end

function main6(args)
    l = parse(Int,args[1+1])
    r = parse(Float64,args[2+1])
    m0 = parse(Float64,args[3+1]) #The m₀ that we know has a good correlation length
    κ = parse(Float64,args[4+1]) #multiplicative span around m₀ we want to study(2 would be a good value)
    n = parse(Int,args[5+1]) #number of different masses
    g0 = parse(Float64,args[6+1])
    ϵg = parse(Float64,args[7+1]) #precision goal
    D = parse(Int,args[9])
    DMAX = parse(Int,args[10])
    rm = (κ)^(1/n)
    m = m0/κ
    local DF_TP = DataFrame(Distance = collect(0:l-1))
    local DF_fit = DataFrame(m0 = Float64[],
        Chi = Float64[],
        b = Float64[],
        error_Chi = Float64[],
        error_b = Float64[],
        inner50 = Float64[],
        error_inner50 = Float64[],
        entanglement1 = Float64[],
        entanglement2 = Float64[],
        entanglement3 = Float64[],
        entanglement4 = Float64[])
    open("Gross-Neveu_TwoPoint_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg)_m0=$(m0).txt","w") do io
        t = time()
        for i in 1:2n+1
            exps = Float64[]
            exps2 = ComplexF64[]
            local GS,ϵf = DMRG(NNHandSSHtoMPO(Hnn(r,l),Hss(m,g0,r,l),l),D=D,ϵ=ϵg,maxiteration=3,DMax=DMAX)
            write(io,"precision = $(ϵf)\nBondDimension = $(GS.BD)\nm = $(m)\nTwo point correlators:\n")
            for k in 1:l
                tp = twopoint(1,k)
                ex = (l-1)*1im*expectation(GS,GS,tp)
                write(io,"$(k):\t\t$(ex)\n")
                append!(exps,abs(ex))
                append!(exps2,ex)
            end
            Ms = GS.tensors
            D = GS.BD
            (σ,j,k) = size(Ms[l])
            F = svd(reshape(permutedims(Ms[l],(2,1,3)),(j,σ*k)))
            local S = F.S
            DF_TP[Symbol("InnerProduct_m=",round(m,digits=4))] = exps
            DF_TP[Symbol("InnerProduct_error_m=",round(m,digits=4))] = 2*sqrt(ϵf)*exps
            DF_TP[Symbol("ActualExpectation_m=",round(m,digits=4))] = exps2
            DF_TP[Symbol("BD_m=",round(m,digits=4))] = repeat([GS.BD],l)
            DF_TP[Symbol("m=",round(m,digits=4))] = repeat([m],l)
            DF_TP[Symbol("ϵ_m=",round(m,digits=4))] = repeat([ϵf],l)
            local small_GS,small_ϵf = DMRG(NNHandSSHtoMPO(Hnn(r,l-1),Hss(m,g0,r,l-1),l-1),D=GS.BD,ϵ=ϵg,maxiteration=3,DMax=DMAX)
            local inner = UniformInnerProductOfMPS(small_GS,GS)
            @. model(x,p) = p[2]*exp(-x/p[1])
            write(io,"model(x,p) =  p[2]*exp(-x/p[1])\n")
            try
                fit = curve_fit(model,2:l-1,exps[3:end],[1.0,1.0])
                write(io,"parameters:\n$(fit.param)\n")
                write(io,"Standard Deviations:\n$(stderror(fit))\n")
                push!(DF_fit,[m,fit.param...,stderror(fit)...,inner,hypot(ϵf,small_ϵf)*inner,S[1:4]...])
            catch e
                write(io,"LsqFit Failed\n")
                @warn e
            end
            write(io,"============================================================\n")
            m = m*rm
        end
        CSV.write("Gross-Neveu_TwoPoint_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg)_m0=$(m0).csv",DF_TP)
        CSV.write("Gross-Neveu_Fitting_l=$(l)_r=$(r)_g0=$(g0)_eps=$(ϵg)_m0=$(m0).csv",DF_fit)
        write(io,"finished computation in $(time()-t).")
    end
end

function Hss2(m₀::Float64,g₀::Float64,r::Float64,a::Float64;ϵ=0.0)
    a⁺ = [0.0 0.0;
          1.0 0.0]
    a⁻ = [0.0 1.0;
          0.0 0.0]
    Z = [1.0 0.0;
         0. -1.0]
    # MagneticField = ϵ*randn(Float64)*Diagonal([1.5, 0.5, -0.5, -1.5])
    # return im*(m₀-r/a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))+g₀^2/(2*a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))^2 + m₀*MagneticField
    return im*(-m₀-r/a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))+g₀^2/(2*a)*(kron(a⁺,a⁻)-kron(a⁻,a⁺))^2
end

function Hnn2(r::Float64,a::Float64)
    a⁺ = [0.0 0.0;
          1.0 0.0]
    a⁻ = [0.0 1.0;
          0.0 0.0]
    Z = [1.0 0.0;
         0. -1.0]
    II = [1. 0.;
          0. 1.]
    # return im/(2*a)*(-kron(II,a⁻,Z,a⁻)+kron(a⁻,Z,a⁻,II)+kron(II,a⁺,Z,a⁺)-kron(a⁺,Z,a⁺,II)+r*(kron(II,a⁻,a⁺,II)-kron(II,a⁺,a⁻,II)+kron(a⁺,Z,Z,a⁻)-kron(a⁻,Z,Z,a⁺)))
    return im/(2*a)*(kron(II,a⁻,Z,a⁺)-kron(a⁻,Z,a⁺,II)-kron(II,a⁺,Z,a⁻)+kron(a⁺,Z,a⁻,II)+r*(kron(II,a⁻,a⁺,II)-kron(II,a⁺,a⁻,II)+kron(a⁺,Z,Z,a⁻)-kron(a⁻,Z,Z,a⁺)))
end

function main7(args)
    l0 = parse(Int,args[2])
    l = parse(Int,args[3])
    r = parse(Float64,args[4])
    m0 = parse(Float64,args[5])
    g0 = parse(Float64,args[6])
    ϵ0 = parse(Float64,args[7])
    D = parse(Int,args[8])
    DMAX = parse(Int,args[9])
    local inners = Array{Float64,1}()
    t = time()
    local a = 1.0/(l+1)
    local DF_inners = DataFrame(Length = Int[],
                                ϵSmall = Float64[],
                                BD_Small = Int[],
                                ϵLarge = Float64[],
                                BD_Large = Int[],
                                InnerProduct = Float64[],
                                time = Float64[])
    local    SmallGroundStates,ϵ_Small = DMRG(NNHandSSHtoMPO(Hnn2(r,a),Hss2(m0,g0,r,a),l0),D=D,ϵ=ϵ0,maxiteration=3,DMax=DMAX)
    open("Fixed-Gross-Neveu_r=$(r)_m0=$(m0)_g0=$(g0)_l=$(l).txt","w") do io
        for Le in l0+1:1:l
            t1 = time()
            local    BigGroundStates,ϵ_Large = DMRG(NNHandSSHtoMPO(Hnn2(r,a),Hss2(m0,g0,r,a),Le),D=SmallGroundStates.BD,ϵ=ϵ0,maxiteration=3,DMax=DMAX)
            push!(inners,UniformInnerProductOfMPS(SmallGroundStates,BigGroundStates))
            write(io,"$(Le) \t $(round(inners[end],digits=5)) \t $(round(time()-t))\n")
            push!(DF_inners,(Le,ϵ_Small,SmallGroundStates.BD,ϵ_Large,BigGroundStates.BD,inners[end],round(time()-t1,digits=2)))
            SmallGroundStates = deepcopy(BigGroundStates)
            ϵ_Small = deepcopy(ϵ_Large)
        end
        CSV.write("Fixed-Gross-Neveu_r=$(r)_m0=$(m0)_g0=$(g0)_l=$(l).csv",DF_inners)
        write(io,"finished computation in $(time()-t).\n")
    end
end

function main()
    fun = parse(Int,ARGS[1])

    if fun == 1 #Inner products of ground states for varying lengths
        #l0,l,r,m0,g0,n=1,ϵ0,D(iniitial BD),DMAX
        main1(ARGS)
    elseif fun == 2
        main2()
    elseif fun == 3 # Inner producs for fixed length but varying mass as a power.
        # l ,r ,m0 ,rm ,n ,g0 ,ϵ0 ,D ,DMAX
        main3(ARGS)
    elseif fun == 4 # Two point function ψ1ψ̄j
        # l, r, m0, rm, n, g0, ϵg, D, DMAX
        main4(ARGS)
    elseif fun == 5 # Two point function ψ̄1ψj
        # l, r, m0, rm, n, g0, ϵg, D, DMAX
        main5(ARGS)
    elseif fun == 6 #calculating both correlation length and asymptotic value of inner products
        # l, r, m0, κ, n, g0, ϵg, D, DMAX
        main6(ARGS)
    elseif fun ==7 #Inner products of ground states for varying lengths with fixed lattice spacing
        #l0,l,r,m0,g0,ϵ0,D(iniitial BD),DMAX
        main7(ARGS)
    end
end

if basename(PROGRAM_FILE) == basename(@__FILE__)
    main()
end

max_stars_repo_name: CNOT/DMRG

max_issues_repo_name: CNOT/DMRG

max_forks_repo_name: CNOT/DMRG


text:
struct iGMRF
    G::GridStructure
    rankDeficiency::Int64
    κ::Float64              # Precision of the field
end

function Base.show(io::IO, obj::iGMRF)

    println(io, "iGMRF")
    println(io, "G :")
    showGridStructure(io, obj.G, prefix = "\t\t\t")
    println(io)
    println(io, "rankDeficiency :\t", obj.rankDeficiency)
    println(io, "κ :\t\t\t", obj.κ)

end

function iGMRF(m₁::Integer, m₂::Integer, order::Integer, κ::Real)::iGMRF

    # Gives the adjacency matrix W for the iGMRF of order 1 or 2 on the regular
    # grid of size (m1 * m2).

    @assert order == 1 || order == 2 "the order should be either 1 or 2."

    if order == 1

        nbs, W = fo_nbs(m₁, m₂)
        condIndSubset = fo_condindsubsets(m₁, m₂)
        rankdef = 1

    else

        nbs, W = so_nbs(m₁, m₂)
        condIndSubset = so_condindsubsets(m₁, m₂)
        rankdef = 3

    end

    W̄ = W - sparse(diagm(length.(nbs)))

    G = GridStructure((m₁, m₂), nbs, condIndSubset, W, W̄)

    return iGMRF(G, rankdef, κ)

end

function fo_nbs(m₁::Integer, m₂::Integer)::Tuple{Vector{Vector{Int64}}, SparseMatrixCSC{Int64,Int64}}

    # 1-off diagonal elements
    v = ones(Int64,m₁)
    v[end] = 0
    V = repeat(v,outer=m₂)
    pop!(V)

    # n-off diagonal elements
    U = ones(Int64,m₁*(m₂-1))

    # get the upper triangular part of the matrix
    m = m₁ * m₂
    D = sparse(1:(m-1), 2:m, V, m, m) + sparse(1:(m-m₁),(m₁+1):m, U, m, m)

    # make D symmetric
    D = D + D'

    # Compute the list of neighbors for each node
    nbs = fill(Int[], m)
    for i = 1:m
        nbs[i] = findall(!iszero, D[:,i])
    end

    # Put the number of neighbors on the diagonal
    W = -D + spdiagm(0 => length.(nbs))

    return (nbs, W)

end

function so_nbs(m₁::Integer, m₂::Integer)::Tuple{Vector{Vector{Int64}}, SparseMatrixCSC{Int64,Int64}}
    # Alternative by adding molecules. There should not be missing values in the grid.
    m = m₁ * m₂
    W = spzeros(Int64,m,m)
    pos = reshape(1:m,m₁,m₂)

    for i=1:m₁
        for j=1:m₂

            S = zeros(Int64,m₁,m₂)

            if (i-2>0)
               S[i-2:i,j] =  S[i-2:i,j] + [1, -2, 1]
            end

            if (i+2<=m₁)
               S[i:i+2,j] =  S[i:i+2,j] + [1, -2, 1]
            end

            if (j-2>0)
               S[i,j-2:j] =  S[i,j-2:j] + [1,-2, 1]
            end

            if (j+2<=m₂)
               S[i,j:j+2] =  S[i,j:j+2] + [1,-2, 1]
            end

            if (i-1>0) && (i+1<=m₁)
                S[i-1:i+1,j] = S[i-1:i+1,j] + [-2, 4, -2]
            end

            if (j-1>0) && (j+1<=m₂)
                S[i,j-1:j+1] = S[i,j-1:j+1] + [-2, 4, -2]
            end

            if (i-1>0) && (j+1<=m₂)
                S[i-1:i,j:j+1] = S[i-1:i,j:j+1] + [-2 2; 2 -2]
            end

            if (i+1<=m₁) && (j+1<=m₂)
                S[i:i+1,j:j+1] = S[i:i+1,j:j+1] + [2 -2; -2 2]
            end

            if (i-1>0) && (j-1>0)
                S[i-1:i,j-1:j] = S[i-1:i,j-1:j] + [2 -2; -2 2]
            end

            if (i+1<=m₁) && (j-1>0)
                S[i:i+1,j-1:j] = S[i:i+1,j-1:j] + [-2 2; 2 -2]
            end

            W[:,pos[i,j]] = S[:]

        end
    end

    # Compute the list of neighbors for each node
    nbs =  Array{Int64,1}[]
    for i = 1:m
        push!(nbs,findall(W[:,i] .< 0))
    end

    return (nbs, W)

end

function fo_condindsubsets(m₁::Integer, m₂::Integer)::Vector{Vector{Integer}}

    condIndSubsetIndex = 2*ones(Int64,m₁,m₂)
    condIndSubsetIndex[1:2:end,1:2:end] .= 1
    condIndSubsetIndex[2:2:end,2:2:end] .= 1

    return Array[findall(vec(condIndSubsetIndex) .==i) for i=1:2]

end

function so_condindsubsets(m₁::Integer, m₂::Integer)::Vector{Vector{Integer}}

    condIndSubsetIndex = zeros(Int64,m₁,m₂)

    condIndSubsetIndex[1:3:end,1:4:end] .= 1
    condIndSubsetIndex[2:3:end,3:4:end] .= 1

    condIndSubsetIndex[1:3:end,2:4:end] .= 2
    condIndSubsetIndex[2:3:end,4:4:end] .= 2

    condIndSubsetIndex[1:3:end,3:4:end] .= 3
    condIndSubsetIndex[2:3:end,1:4:end] .= 3

    condIndSubsetIndex[1:3:end,4:4:end] .= 4
    condIndSubsetIndex[2:3:end,2:4:end] .= 4

    condIndSubsetIndex[3:3:end,1:3:end] .= 5
    condIndSubsetIndex[3:3:end,2:3:end] .= 6
    condIndSubsetIndex[3:3:end,3:3:end] .= 7

    return Array[findall(vec(condIndSubsetIndex) .==i) for i=1:7]

end

function rand(F::iGMRF)::Vector{<:Real}

    @assert F.rankDeficiency == 1 || F.rankDeficiency == 3 "The rank deficiency should be either 1 or 3"

    κ = F.κ
    W = F.G.W
    m₁ = F.G.gridSize[1]
    m₂ = F.G.gridSize[2]
    m = m₁ * m₂

    if F.rankDeficiency == 1

        e₁ = ones(m,1)

        A = e₁

        Q = κ*W + e₁*e₁'

    else

        e₁ = ones(m)
        e₂ = repeat(1:m₁, m₂)
        e₃ = repeat(1:m₂,inner = m₁)

        A = hcat(e₁,e₂,e₃)

        Q = κ*W + e₁*e₁' + e₂*e₂' + e₃*e₃'

    end

    C = cholesky(Q)
    L = C.L

    z = randn(m)

    x = L'\z

#     V = zeros(m,size(A,2))
#     for ii=1:size(A,2)
#        V[:,ii] = C\A[:,ii]
#     end
    V = C\A
    W = A'*V
    U = W\(V')
    c = A' * x
    y = x - U' * c

    return y

end

function logpdf(F::iGMRF, y::Array{<:Real})::Real

    κ = F.κ

    W = F.G.W
    m = F.G.gridSize[1] * F.G.gridSize[2]

    k = F.rankDeficiency

    v = κ*W*y
    q = y'*v

    lpdf =  .5*(m-k)*log(κ) - .5*q

    return lpdf

end

function fullconditionals(F::iGMRF, y::Vector{<:Real})::Vector{NormalCanon}

    κ = F.κ

    W̄ = F.G.W̄
    W = F.G.W

    Q = κ * Array(diag(F.G.W))
    h = -κ*(W̄*y)

    pd = NormalCanon.(h,Q)

    return pd

end

function fullcondlogpdf(F::iGMRF, y::Vector{<:Real})::Vector{<:Real}

    pd = fullconditionals(F::iGMRF,y::Vector{<:Real})

    clpdf = logpdf.(pd,y)

    return clpdf

end

function getconditional(F::GMRF.iGMRF, B::Vector{<:Integer}, x::Vector{<:Real})::MvNormalCanon

    W = F.G.W

    sort!(B)

    A = setdiff(1:(F.G.gridSize[1] * F.G.gridSize[2]), B)

    Waa = W[A,A]
    Wab = W[A,B]

    h = -Wab*x*F.κ

    J = Array(F.κ*Waa)

    pd = MvNormalCanon(h,J)

    return pd

end

ext: jl


text:
using LinQuadOptInterface
using .ClpCInterface

const LQOI = LinQuadOptInterface
const MOI  = LQOI.MOI

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear,
    LQOI.SinVar
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct Optimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    params::Dict{Symbol,Any}
    Optimizer(::Nothing) = new()
end

### Options

# map option name to C function
const optionmap = Dict(
   :PrimalTolerance => set_primal_tolerance,
   :DualTolerance => set_dual_tolerance,
   :DualObjectiveLimit => set_dual_objective_limit,
   :MaximumIterations => set_maximum_iterations,
   :MaximumSeconds => set_maximum_seconds,
   :LogLevel => set_log_level,
   :Scaling => scaling,
   :Perturbation => set_perturbation,
   )
# These options are set by using the ClpSolve object
const solveoptionmap = Dict(
   :PresolveType => set_presolve_type,
   :SolveType => set_solve_type,
   :InfeasibleReturn => set_infeasible_return,
   )

function Optimizer(;kwargs...)
    optimizer = Optimizer(nothing)
    optimizer.params = Dict{String,Any}()
    MOI.empty!(optimizer)
    for (name,value) in kwargs
        optimizer.params[Symbol(name)] = value
    end
    return optimizer
end

LQOI.LinearQuadraticModel(::Type{Optimizer},env) = ClpModel()

LQOI.supported_constraints(optimizer::Optimizer) = SUPPORTED_CONSTRAINTS
LQOI.supported_objectives(optimizer::Optimizer) = SUPPORTED_OBJECTIVES

"""
    replace_inf(x::Real)

Return `Inf` if `x>1e20`, `-Inf` if `x<-1e20`, and `x` otherwise.
"""
function replace_inf(x::Real)
    if x > 1e20
        return Inf
    elseif x < -1e20
        return -Inf
    else
        return x
    end
end

function LQOI.change_variable_bounds!(instance::Optimizer, cols::Vector{Int},
        values::Vector{Float64}, senses::Vector)
    upperbounds = get_col_upper(instance.inner)
    lowerbounds = get_col_lower(instance.inner)
    for (col, value, sense) in zip(cols, values, senses)
        if sense == Cchar('U')
            upperbounds[col] = value
        elseif sense == Cchar('L')
            lowerbounds[col] = value
        else
            error("sense is Cchar('$(Char(sense))'), but only Cchar('U') " *
                  "Cchar('L') are supported.")
        end
    end
    chg_column_upper(instance.inner, upperbounds)
    chg_column_lower(instance.inner, lowerbounds)
end

function LQOI.get_variable_lowerbound(instance::Optimizer, col::Int)
    lower = get_col_lower(instance.inner)
    return replace_inf(lower[col])
end

function LQOI.get_variable_upperbound(instance::Optimizer, col::Int)
    upper = get_col_upper(instance.inner)
    return replace_inf(upper[col])
end

function LQOI.get_number_linear_constraints(instance::Optimizer)
    return get_num_rows(instance.inner)
end

"""
    append_row(instance::Optimizer, row::Int, lower::Float64, upper::Float64,
               rows::Vector{Int}, cols::Vector{Int}, coefs::Vector{Float64})

Given a sparse matrix in the triplet-form `(rows, cols, coefs)`, add row `row`
with upper bound `upper` and lower bound  `lower` to the instance `instance`.
"""
function append_row(instance::Optimizer, row::Int, lower::Float64,
                    upper::Float64, rows::Vector{Int}, cols::Vector{Int},
                    coefs::Vector{Float64})
    indices = if row == length(rows)
        rows[row]:length(cols)
    else
        rows[row]:(rows[row+1]-1)
    end
    add_row(instance.inner, Cint(length(indices)), Cint.(cols[indices] .- 1),
            coefs[indices], lower, upper)
end

function LQOI.add_linear_constraints!(instance::Optimizer, A::LQOI.CSRMatrix{Float64},
        senses::Vector{Cchar}, right_hand_sides::Vector{Float64})
    rows = A.row_pointers
    cols = A.columns
    coefs = A.coefficients
    for (row, (rhs, sense)) in enumerate(zip(right_hand_sides, senses))
        if rhs > 1e20
            error("rhs must always be less than 1e20")
        elseif rhs < -1e20
            error("rhs must always be greater than -1e20")
        end
        lower = -Inf
        upper = Inf
        if sense == Cchar('L')
            upper = rhs
        elseif sense == Cchar('G')
            lower = rhs
        elseif sense == Cchar('E')
            upper = lower = rhs
        else
            error("sense must be Cchar(x) where x is in ['L','G',E']")
        end
        append_row(instance, row, lower, upper, rows, cols, coefs)
    end
end

function LQOI.add_ranged_constraints!(instance::Clp.Optimizer,
        A::LinQuadOptInterface.CSRMatrix{Float64},
        lb::Vector{Float64}, ub::Vector{Float64})
    rows = A.row_pointers
    cols = A.columns
    coefs = A.coefficients
    for row in 1:length(lb)
        append_row(instance, row, lb[row], ub[row], rows, cols, coefs)
    end
end

function LQOI.get_rhs(instance::Optimizer, row::Int)
    lower_bounds = get_row_lower(instance.inner)
    upper_bounds = get_row_upper(instance.inner)
    lower_bound = replace_inf(lower_bounds[row])
    upper_bound = replace_inf(upper_bounds[row])
    if lower_bound > -Inf
        return lower_bound
    elseif upper_bound < Inf
        return upper_bound
    else
        error("Either row_lower or row_upper must be of abs less than 1e20")
    end
end

function LQOI.get_linear_constraint(instance::Optimizer, row::Int)::Tuple{Vector{Int}, Vector{Float64}}
    A = get_constraint_matrix(instance.inner)
    A_row = A[row,:]
    return Array{Int}(A_row.nzind), A_row.nzval
end

function LQOI.change_objective_coefficient!(instance::Optimizer, col, coef)
    objcoefs = get_obj_coefficients(instance.inner)
    objcoefs[col] = coef
    chg_obj_coefficients(instance.inner, objcoefs)
end

function LQOI.change_rhs_coefficient!(instance::Optimizer, row, coef)
    lower_bounds = get_row_lower(instance.inner)
    upper_bounds = get_row_upper(instance.inner)
    lower_bound = replace_inf(lower_bounds[row])
    upper_bound = replace_inf(upper_bounds[row])
    if lower_bound > -Inf
        lower_bounds[row] = coef
        chg_row_lower(instance.inner, lower_bounds)
    elseif upper_bound < Inf
        upper_bounds[row] = coef
        chg_row_upper(instance.inner, upper_bounds)
    else
        error("Either row_lower or row_upper must be of abs less than 1e20")
    end
end

function LQOI.delete_linear_constraints!(instance::Optimizer, start_row::Int, end_row::Int)
    delete_rows(instance.inner, [Cint(i-1) for i in start_row:end_row])
end

function LQOI.change_linear_constraint_sense!(instance::Optimizer, rows::Vector{Int}, senses::Vector{Cchar})
    lower = replace_inf.(get_row_lower(instance.inner))
    upper = replace_inf.(get_row_upper(instance.inner))
    for (sense, row) in zip(senses, rows)
        lb = lower[row]
        ub = upper[row]
        if lb > -Inf
            rhs = lb
        elseif ub < Inf
            rhs = ub
        else
            error("Either row_lower or row_upper must be of abs less than 1e20")
        end
        if sense == Cchar('G')
            lower[row] = rhs
            upper[row] = Inf
        elseif sense == Cchar('L')
            lower[row] = -Inf
            upper[row] = rhs
        elseif sense == Cchar('E')
            lower[row] = rhs
            upper[row] = rhs
        end
    end
    chg_row_upper(instance.inner, upper)
    chg_row_lower(instance.inner, lower)
end

function LQOI.set_linear_objective!(instance::Optimizer, cols::Vector{Int}, coefs::Vector{Float64})
    objective_coefficients = zeros(Float64, get_num_cols(instance.inner))
    for (col, coef) in zip(cols, coefs)
        objective_coefficients[col] += coef
    end
    chg_obj_coefficients(instance.inner, objective_coefficients)
end

function LQOI.change_objective_sense!(instance::Optimizer, sense::Symbol)
    if sense == :min
        set_obj_sense(instance.inner, 1.0)
    elseif sense == :max
        set_obj_sense(instance.inner, -1.0)
    else
        error("sense must be either :min or :max")
    end
end

function LQOI.get_linear_objective!(instance::Optimizer, x::Vector{Float64})
    copyto!(x, get_obj_coefficients(instance.inner))
end

function LQOI.solve_linear_problem!(instance::Optimizer)
    solveroptions = ClpSolve()
    model = instance.inner
    for (name, value) in instance.params
        if haskey(optionmap, name)
            optionmap[name](model, value)
        elseif haskey(solveoptionmap, name)
            solveoptionmap[name](solveroptions,value)
        else
            error("Unrecognized option: $name")
        end
    end
    initial_solve_with_options(instance.inner, solveroptions)
end

function LQOI.get_variable_primal_solution!(instance::Optimizer, x::Vector{Float64})
    copyto!(x, primal_column_solution(instance.inner))
end

function LQOI.get_linear_primal_solution!(instance::Optimizer, x::Vector{Float64})
    copyto!(x, primal_row_solution(instance.inner))
end

function LQOI.get_variable_dual_solution!(instance::Optimizer, x::Vector{Float64})
    copyto!(x, dual_column_solution(instance.inner))
end

function LQOI.get_linear_dual_solution!(instance::Optimizer, x::Vector{Float64})
    copyto!(x, dual_row_solution(instance.inner))
end

function LQOI.get_objective_value(instance::Optimizer)
    return objective_value(instance.inner)
end

function LQOI.get_farkas_dual!(instance::Optimizer, result::Vector{Float64})
    copyto!(result, infeasibility_ray(instance.inner))
    Compat.rmul!(result, -1.0)
end

function LQOI.get_unbounded_ray!(instance::Optimizer, result::Vector{Float64})
    copyto!(result, unbounded_ray(instance.inner))
end

function LQOI.get_termination_status(instance::Optimizer)
    status = ClpCInterface.status(instance.inner)
    if status == 0
        return MOI.Success
    elseif status == 1
        if is_proven_primal_infeasible(instance.inner)
            return MOI.Success
        else
            return MOI.InfeasibleNoResult
        end
    elseif status == 2
        if is_proven_dual_infeasible(instance.inner)
            return MOI.Success
        else
            return MOI.UnboundedNoResult
        end
    elseif status == 3
        return MOI.OtherLimit
    elseif status == 4
        return MOI.OtherError
    else
        error("status returned was $(status), but it must be in [0,1,2,3,4]")
    end
end

function LQOI.get_primal_status(instance::Optimizer)
    if is_proven_dual_infeasible(instance.inner)
        return MOI.InfeasibilityCertificate
    elseif primal_feasible(instance.inner)
        return MOI.FeasiblePoint
    else
        return MOI.UnknownResultStatus
    end
end

function LQOI.get_dual_status(instance::Optimizer)
    if is_proven_primal_infeasible(instance.inner)
        return MOI.InfeasibilityCertificate
    elseif dual_feasible(instance.inner)
        return MOI.FeasiblePoint
    else
        return MOI.UnknownResultStatus
    end
end

function LQOI.get_number_variables(instance::Optimizer)
    return get_num_cols(instance.inner)
end

function LQOI.add_variables!(instance::Optimizer, number_of_variables::Int)
    for i in 1:number_of_variables
        add_column(instance.inner, Cint(0), Int32[], Float64[], -Inf, Inf, 0.0)
    end
end

function LQOI.delete_variables!(instance::Optimizer, start_col::Int, end_col::Int)
    delete_columns(instance.inner, [Cint(i-1) for i in start_col:end_col])
end


text:
# test various training parameters

bigdata = joinpath("..", "data", "big.txt")

vector_size = rand(100:300)
window = rand(1:5)
sample = [0, 1e-5][rand(1:2)]
hs = rand(0:1)
negative = rand(1:5)
threads = rand(1:16)
iter = rand(1:5)
min_count = rand(1:5)
alpha = 0.025
debug = 2
threshold = rand(50:100)
binary = 1
cbow = rand(0:1)
classes = rand(30:100)
save_vocab = "vocab.txt"
read_vocab = "vocab.txt"

word2phrase(bigdata, "bigphrase", min_count=min_count, threshold=threshold,
            debug=debug)

word2vec("bigphrase", "bigvecs.txt", size=vector_size, window=window, sample=sample,
         hs=hs, negative=negative, threads=threads, iter=iter,
         min_count=min_count, alpha=alpha, debug=debug, cbow=cbow,
         save_vocab=save_vocab)
word2vec(bigdata, "bigvecs.bin", size=vector_size, window=window,
         sample=sample, hs=hs, negative=negative, threads=threads,
         iter=iter, binary=binary, read_vocab = read_vocab)

rm("bigphrase")

model = wordvectors("bigvecs.txt")
len_vecs, num_words = size(model)
@test len_vecs == vector_size

word2clusters(bigdata, "bigclus.txt", classes, size=vector_size, window=window,
              sample=sample, hs=hs, negative=negative, threads=threads,
              iter=iter, min_count=min_count, alpha=alpha, debug=debug,
              cbow=cbow, read_vocab=read_vocab, save_vocab=save_vocab)

println("training passed test...")


text:
# Ordinary metrics

###########################################################
#
#   Abstract metric types
#
###########################################################

abstract type UnionPreMetric <: PreMetric end

abstract type UnionSemiMetric <: SemiMetric end

abstract type UnionMetric <: Metric end

###########################################################
#
#   Metric types
#
###########################################################

struct Euclidean <: UnionMetric
    thresh::Float64
end

"""
    Euclidean([thresh])

Create a euclidean metric.

When computing distances among large numbers of points, it can be much
more efficient to exploit the formula

    (x-y)^2 = x^2 - 2xy + y^2

However, this can introduce roundoff error. `thresh` (which defaults
to 0) specifies the relative square-distance tolerance on `2xy`
compared to `x^2 + y^2` to force recalculation of the distance using
the more precise direct (elementwise-subtraction) formula.

# Example:
```julia
julia> x = reshape([0.1, 0.3, -0.1], 3, 1);

julia> pairwise(Euclidean(), x, x)
1×1 Array{Float64,2}:
 7.45058e-9

julia> pairwise(Euclidean(1e-12), x, x)
1×1 Array{Float64,2}:
 0.0
```
"""
Euclidean() = Euclidean(0)

struct WeightedEuclidean{W} <: UnionMetric
    weights::W
end

"""
    PeriodicEuclidean(L)

Create a Euclidean metric on a rectangular periodic domain (i.e., a torus or
a cylinder). Periods per dimension are contained in the vector `L`:
```math
\\sqrt{\\sum_i(\\min\\mod(|x_i - y_i|, p), p - \\mod(|x_i - y_i|, p))^2}.
```
For dimensions without periodicity put `Inf` in the respective component.

# Example
```jldoctest
julia> x, y, L = [0.0, 0.0], [0.75, 0.0], [0.5, Inf];

julia> evaluate(PeriodicEuclidean(L), x, y)
0.25
```
"""
struct PeriodicEuclidean{W} <: UnionMetric
    periods::W
end

struct SqEuclidean <: UnionSemiMetric
    thresh::Float64
end

"""
    SqEuclidean([thresh])

Create a squared-euclidean semi-metric. For the meaning of `thresh`,
see [`Euclidean`](@ref).
"""
SqEuclidean() = SqEuclidean(0)

struct WeightedSqEuclidean{W} <: UnionSemiMetric
    weights::W
end

struct Chebyshev <: UnionMetric end

struct Cityblock <: UnionMetric end
struct WeightedCityblock{W} <: UnionMetric
    weights::W
end

struct TotalVariation <: UnionMetric end
struct Jaccard <: UnionMetric end
struct RogersTanimoto <: UnionMetric end

struct Minkowski{T <: Real} <: UnionMetric
    p::T
end
struct WeightedMinkowski{W,T <: Real} <: UnionMetric
    weights::W
    p::T
end

struct Hamming <: UnionMetric end
struct WeightedHamming{W} <: UnionMetric
    weights::W
end

struct CosineDist <: UnionSemiMetric end
# CorrDist is excluded from `UnionMetrics`
struct CorrDist <: SemiMetric end
struct BrayCurtis <: UnionSemiMetric end

struct ChiSqDist <: UnionSemiMetric end
struct KLDivergence <: UnionPreMetric end
struct GenKLDivergence <: UnionPreMetric end

"""
    RenyiDivergence(α::Real)
    renyi_divergence(P, Q, α::Real)

Create a Rényi premetric of order α.

Rényi defined a spectrum of divergence measures generalising the
Kullback–Leibler divergence (see `KLDivergence`). The divergence is
not a semimetric as it is not symmetric. It is parameterised by a
parameter α, and is equal to Kullback–Leibler divergence at α = 1:

At α = 0, ``R_0(P | Q) = -log(sum_{i: p_i > 0}(q_i))``

At α = 1, ``R_1(P | Q) = sum_{i: p_i > 0}(p_i log(p_i / q_i))``

At α = ∞, ``R_∞(P | Q) = log(sup_{i: p_i > 0}(p_i / q_i))``

Otherwise ``R_α(P | Q) = log(sum_{i: p_i > 0}((p_i ^ α) / (q_i ^ (α - 1))) / (α - 1)``

# Example:
```jldoctest
julia> x = reshape([0.1, 0.3, 0.4, 0.2], 2, 2);

julia> pairwise(RenyiDivergence(0), x, x)
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> pairwise(Euclidean(2), x, x)
2×2 Array{Float64,2}:
 0.0       0.577315
 0.655407  0.0
```
"""
struct RenyiDivergence{T <: Real} <: UnionPreMetric
    p::T # order of power mean (order of divergence - 1)
    is_normal::Bool
    is_zero::Bool
    is_one::Bool
    is_inf::Bool

    function RenyiDivergence{T}(q) where {T}
        # There are four different cases:
        #   simpler to separate them out now, not over and over in eval_op()
        is_zero = q ≈ zero(T)
        is_one = q ≈ one(T)
        is_inf = isinf(q)

        # Only positive Rényi divergences are defined
        !is_zero && q < zero(T) && throw(ArgumentError("Order of Rényi divergence not legal, $(q) < 0."))

        new{T}(q - 1, !(is_zero || is_one || is_inf), is_zero, is_one, is_inf)
    end
end
RenyiDivergence(q::T) where {T} = RenyiDivergence{T}(q)

struct JSDivergence <: UnionSemiMetric end

struct SpanNormDist <: UnionSemiMetric end

# Deviations are handled separately from the other distances/divergences and
# are excluded from `UnionMetrics`
struct MeanAbsDeviation <: Metric end
struct MeanSqDeviation <: SemiMetric end
struct RMSDeviation <: Metric end
struct NormRMSDeviation <: PreMetric end

# Union types
const metrics = (Euclidean,SqEuclidean,PeriodicEuclidean,Chebyshev,Cityblock,TotalVariation,Minkowski,Hamming,Jaccard,RogersTanimoto,CosineDist,ChiSqDist,KLDivergence,RenyiDivergence,BrayCurtis,JSDivergence,SpanNormDist,GenKLDivergence)
const weightedmetrics = (WeightedEuclidean,WeightedSqEuclidean,WeightedCityblock,WeightedMinkowski,WeightedHamming)
const UnionMetrics = Union{UnionPreMetric,UnionSemiMetric,UnionMetric}

###########################################################
#
#  Implementations
#
###########################################################

parameters(::UnionPreMetric) = nothing
parameters(::UnionSemiMetric) = nothing
parameters(::UnionMetric) = nothing
parameters(d::PeriodicEuclidean) = d.periods
for dist in weightedmetrics
    @eval parameters(d::$dist) = d.weights
end

result_type(dist::UnionMetrics, ::Type{Ta}, ::Type{Tb}) where {Ta,Tb} =
    result_type(dist, _eltype(Ta), _eltype(Tb), parameters(dist))
result_type(dist::UnionMetrics, ::Type{Ta}, ::Type{Tb}, ::Nothing) where {Ta,Tb} =
    typeof(_evaluate(dist, oneunit(Ta), oneunit(Tb)))
result_type(dist::UnionMetrics, ::Type{Ta}, ::Type{Tb}, p) where {Ta,Tb} =
    typeof(_evaluate(dist, oneunit(Ta), oneunit(Tb), oneunit(_eltype(p))))

Base.@propagate_inbounds function _evaluate(d::UnionMetrics, a, b)
    _evaluate(d, a, b, parameters(d))
end

# breaks the implementation into eval_start, eval_op, eval_reduce and eval_end

Base.@propagate_inbounds function _evaluate(d::UnionMetrics, a, b, ::Nothing)
    @boundscheck if length(a) != length(b)
        throw(DimensionMismatch("first collection has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    if length(a) == 0
        return zero(result_type(d, a, b))
    end
    s = eval_start(d, a, b)
    for (ai, bi) in zip(a, b)
        s = eval_reduce(d, s, eval_op(d, ai, bi))
    end
    return eval_end(d, s)
end
Base.@propagate_inbounds function _evaluate(d::UnionMetrics, a::AbstractArray, b::AbstractArray, ::Nothing)
    @boundscheck if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    if length(a) == 0
        return zero(result_type(d, a, b))
    end
    begin
        s = eval_start(d, a, b)
        if (IndexStyle(a, b) === IndexLinear() && eachindex(a) == eachindex(b)) || axes(a) == axes(b)
            for I in eachindex(a, b)
                ai = a[I]
                bi = b[I]
                s = eval_reduce(d, s, eval_op(d, ai, bi))
            end
        else
            for (ai, bi) in zip(a, b)
                s = eval_reduce(d, s, eval_op(d, ai, bi))
            end
        end
        return eval_end(d, s)
    end
end

Base.@propagate_inbounds function _evaluate(d::UnionMetrics, a, b, p)
    @boundscheck if length(a) != length(b)
        throw(DimensionMismatch("first collection has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    @boundscheck if length(a) != length(p)
        throw(DimensionMismatch("data collections have length $(length(a)) but parameters have length $(length(p))."))
    end
    if length(a) == 0
        return zero(result_type(d, a, b))
    end
    s = eval_start(d, a, b)
    for (ai, bi, pi) in zip(a, b, p)
        s = eval_reduce(d, s, eval_op(d, ai, bi, pi))
    end
    return eval_end(d, s)
end
Base.@propagate_inbounds function _evaluate(d::UnionMetrics, a::AbstractArray, b::AbstractArray, p::AbstractArray)
    @boundscheck if length(a) != length(b)
        throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b))."))
    end
    @boundscheck if length(a) != length(p)
        throw(DimensionMismatch("arrays have length $(length(a)) but parameters have length $(length(p))."))
    end
    if length(a) == 0
        return zero(result_type(d, a, b))
    end
    begin
        s = eval_start(d, a, b)
        if (IndexStyle(a, b, p) === IndexLinear() && eachindex(a) == eachindex(b) == eachindex(p)) ||
                axes(a) == axes(b) == axes(p)
            for I in eachindex(a, b, p)
                ai = a[I]
                bi = b[I]
                pi = p[I]
                s = eval_reduce(d, s, eval_op(d, ai, bi, pi))
            end
        else
            for (ai, bi, pi) in zip(a, b, p)
                s = eval_reduce(d, s, eval_op(d, ai, bi, pi))
            end
        end
        return eval_end(d, s)
    end
end

_evaluate(dist::UnionMetrics, a::Number, b::Number, ::Nothing) = eval_end(dist, eval_op(dist, a, b))
function _evaluate(dist::UnionMetrics, a::Number, b::Number, p)
    length(p) != 1 && throw(DimensionMismatch("inputs are scalars but parameters have length $(length(p))."))
    eval_end(dist, eval_op(dist, a, b, first(p)))
end

eval_start(d::UnionMetrics, a, b) = _eval_start(d, _eltype(a), _eltype(b))
_eval_start(d::UnionMetrics, ::Type{Ta}, ::Type{Tb}) where {Ta,Tb} =
    _eval_start(d, Ta, Tb, parameters(d))
_eval_start(d::UnionMetrics, ::Type{Ta}, ::Type{Tb}, ::Nothing) where {Ta,Tb} =
    zero(typeof(eval_op(d, oneunit(Ta), oneunit(Tb))))
_eval_start(d::UnionMetrics, ::Type{Ta}, ::Type{Tb}, p) where {Ta,Tb} =
    zero(typeof(eval_op(d, oneunit(Ta), oneunit(Tb), oneunit(_eltype(p)))))
eval_reduce(::UnionMetrics, s1, s2) = s1 + s2
eval_end(::UnionMetrics, s) = s

for M in (metrics..., weightedmetrics...)
    @eval (dist::$M)(a, b) = _evaluate(dist, a, b, parameters(dist))
end

# Euclidean
eval_op(::Euclidean, ai, bi) = abs2(ai - bi)
eval_end(::Euclidean, s) = sqrt(s)
const euclidean = Euclidean()

# Weighted Euclidean
eval_op(::WeightedEuclidean, ai, bi, wi) = abs2(ai - bi) * wi
eval_end(::WeightedEuclidean, s) = sqrt(s)
weuclidean(a, b, w) = WeightedEuclidean(w)(a, b)

# PeriodicEuclidean
function eval_op(::PeriodicEuclidean, ai, bi, p)
    s1 = abs(ai - bi)
    s2 = mod(s1, p)
    s3 = min(s2, p - s2)
    abs2(s3)
end
eval_end(::PeriodicEuclidean, s) = sqrt(s)
peuclidean(a, b, p) = PeriodicEuclidean(p)(a, b)

# SqEuclidean
eval_op(::SqEuclidean, ai, bi) = abs2(ai - bi)
const sqeuclidean = SqEuclidean()

# Weighted Squared Euclidean
eval_op(::WeightedSqEuclidean, ai, bi, wi) = abs2(ai - bi) * wi
wsqeuclidean(a, b, w) = WeightedSqEuclidean(w)(a, b)

# Cityblock
eval_op(::Cityblock, ai, bi) = abs(ai - bi)
const cityblock = Cityblock()

# Weighted City Block
eval_op(::WeightedCityblock, ai, bi, wi) = abs((ai - bi) * wi)
wcityblock(a, b, w) = WeightedCityblock(w)(a, b)

# Total variation
eval_op(::TotalVariation, ai, bi) = abs(ai - bi)
eval_end(::TotalVariation, s) = s / 2
const totalvariation = TotalVariation()

# Chebyshev
eval_op(::Chebyshev, ai, bi) = abs(ai - bi)
eval_reduce(::Chebyshev, s1, s2) = max(s1, s2)
# if only NaN, will output NaN
Base.@propagate_inbounds eval_start(::Chebyshev, a, b) = abs(first(a) - first(b))
const chebyshev = Chebyshev()

# Minkowski
eval_op(dist::Minkowski, ai, bi) = abs(ai - bi)^dist.p
eval_end(dist::Minkowski, s) = s^(1 / dist.p)
minkowski(a, b, p::Real) = Minkowski(p)(a, b)

# Weighted Minkowski
eval_op(dist::WeightedMinkowski, ai, bi, wi) = abs(ai - bi)^dist.p * wi
eval_end(dist::WeightedMinkowski, s) = s^(1 / dist.p)
wminkowski(a, b, w, p::Real) = WeightedMinkowski(w, p)(a, b)

# Hamming
result_type(::Hamming, ::Type, ::Type) = Int # fallback for Hamming
eval_start(d::Hamming, a, b) = 0
eval_op(::Hamming, ai, bi) = ai != bi ? 1 : 0
const hamming = Hamming()

# WeightedHamming
eval_op(::WeightedHamming, ai, bi, wi) = ai != bi ? wi : zero(eltype(wi))
whamming(a, b, w) = WeightedHamming(w)(a, b)

# Cosine dist
eval_start(dist::CosineDist, a, b) =
    zero.(typeof.(eval_op(dist, oneunit(_eltype(a)), oneunit(_eltype(b)))))
eval_op(::CosineDist, ai, bi) = ai * bi, ai * ai, bi * bi
function eval_reduce(::CosineDist, s1, s2)
    a1, b1, c1 = s1
    a2, b2, c2 = s2
    return a1 + a2, b1 + b2, c1 + c2
end
function eval_end(::CosineDist, s)
    ab, a2, b2 = s
    max(1 - ab / (sqrt(a2) * sqrt(b2)), 0)
end
const cosine_dist = CosineDist()

# CorrDist
_centralize(x) = x .- mean(x)
(::CorrDist)(a, b) = CosineDist()(_centralize(a), _centralize(b))
(::CorrDist)(a::Number, b::Number) = CosineDist()(zero(mean(a)), zero(mean(b)))
const corr_dist = CorrDist()

# ChiSqDist
eval_op(::ChiSqDist, ai, bi) = (d = abs2(ai - bi) / (ai + bi); ifelse(ai != bi, d, zero(d)))
const chisq_dist = ChiSqDist()

# KLDivergence
eval_op(dist::KLDivergence, ai, bi) =
    ai > 0 ? ai * log(ai / bi) : zero(eval_op(dist, oneunit(ai), bi))
const kl_divergence = KLDivergence()

# GenKLDivergence
eval_op(dist::GenKLDivergence, ai, bi) =
    ai > 0 ? ai * log(ai / bi) - ai + bi : oftype(eval_op(dist, oneunit(ai), bi), bi)
const gkl_divergence = GenKLDivergence()

# RenyiDivergence
Base.@propagate_inbounds function eval_start(::RenyiDivergence, a, b)
    T = promote_type(_eltype(a), _eltype(b))
    zero(T), zero(T), T(sum(a)), T(sum(b))
end

function eval_op(dist::RenyiDivergence, ai::T, bi::T) where {T <: Real}
    if ai == zero(T)
        return zero(T), zero(T), zero(T), zero(T)
    elseif dist.is_normal
        return ai, ai * ((ai / bi)^dist.p), zero(T), zero(T)
    elseif dist.is_zero
        return ai, bi, zero(T), zero(T)
    elseif dist.is_one
        return ai, ai * log(ai / bi), zero(T), zero(T)
    else # otherwise q = ∞
        return ai, ai / bi, zero(T), zero(T)
    end
end

function eval_reduce(dist::RenyiDivergence,
                                               s1::Tuple{T,T,T,T},
                                               s2::Tuple{T,T,T,T}) where {T <: Real}
    if dist.is_inf
        if s1[1] == zero(T)
            return (s2[1], s2[2], s1[3], s1[4])
        elseif s2[1] == zero(T)
            return s1
        else
            return s1[2] > s2[2] ? s1 : (s2[1], s2[2], s1[3], s1[4])
        end
    else
        return s1[1] + s2[1], s1[2] + s2[2], s1[3], s1[4]
    end
end

function eval_end(dist::RenyiDivergence, s::Tuple{T,T,T,T}) where {T <: Real}
    if dist.is_zero || dist.is_normal
        log(s[2] / s[1]) / dist.p + log(s[4] / s[3])
    elseif dist.is_one
        return s[2] / s[1] + log(s[4] / s[3])
    else # q = ∞
        log(s[2]) + log(s[4] / s[3])
    end
end

renyi_divergence(a, b, q::Real) = RenyiDivergence(q)(a, b)
# Combine docs with RenyiDivergence. Fetching the docstring with @doc causes
# problems during package compilation; see
# https://github.com/JuliaLang/julia/issues/31640
let docstring = Base.Docs.getdoc(RenyiDivergence)
    @doc docstring renyi_divergence
end

# JSDivergence

function eval_op(::JSDivergence, ai::T, bi::T) where {T}
    u = (ai + bi) / 2
    ta = ai > 0 ? ai * log(ai) / 2 : zero(log(one(T)))
    tb = bi > 0 ? bi * log(bi) / 2 : zero(log(one(T)))
    tu = u > 0 ? u * log(u) : zero(log(one(T)))
    ta + tb - tu
end
const js_divergence = JSDivergence()

# SpanNormDist

result_type(dist::SpanNormDist, ::Type{Ta}, ::Type{Tb}) where {Ta,Tb} =
    typeof(eval_op(dist, oneunit(Ta), oneunit(Tb)))
Base.@propagate_inbounds function eval_start(::SpanNormDist, a, b)
    d = first(a) - first(b)
    return d, d
end
eval_op(::SpanNormDist, ai, bi)  = ai - bi
function eval_reduce(::SpanNormDist, s1, s2)
    min_d, max_d = s1
    if s2 > max_d
        max_d = s2
    elseif s2 < min_d
        min_d = s2
    end
    return min_d, max_d
end

eval_end(::SpanNormDist, s) = s[2] - s[1]
(::SpanNormDist)(a::Number, b::Number) = zero(promote_type(typeof(a), typeof(b)))
const spannorm_dist = SpanNormDist()

# Jaccard

eval_start(dist::Jaccard, a, b) =
    zero.(typeof.(eval_op(dist, oneunit(_eltype(a)), oneunit(_eltype(b)))))
function eval_op(::Jaccard, s1, s2)
    abs_m = abs(s1 - s2)
    abs_p = abs(s1 + s2)
    abs_p - abs_m, abs_p + abs_m
end
function eval_reduce(::Jaccard, s1, s2)
    a = s1[1] + s2[1]
    b = s1[2] + s2[2]
    a, b
end
function eval_end(::Jaccard, a)
    v = 1 - (a[1] / a[2])
    return v
end
const jaccard = Jaccard()

# BrayCurtis

eval_start(dist::BrayCurtis, a, b) =
    zero.(typeof.(eval_op(dist, oneunit(_eltype(a)), oneunit(_eltype(b)))))
function eval_op(::BrayCurtis, s1, s2)
    abs_m = abs(s1 - s2)
    abs_p = abs(s1 + s2)
    abs_m, abs_p
end
function eval_reduce(::BrayCurtis, s1, s2)
    a = s1[1] + s2[1]
    b = s1[2] + s2[2]
    a, b
end
function eval_end(::BrayCurtis, a)
    v = a[1] / a[2]
    return v
end
const braycurtis = BrayCurtis()

# Tanimoto

eval_start(::RogersTanimoto, _, _) = 0, 0, 0, 0
function eval_op(::RogersTanimoto, s1, s2)
    tt = s1 && s2
    tf = s1 && !s2
    ft = !s1 && s2
    ff = !s1 && !s2
    tt, tf, ft, ff
end
function eval_reduce(::RogersTanimoto, s1, s2)
    begin
        a = s1[1] + s2[1]
        b = s1[2] + s2[2]
        c = s1[3] + s2[3]
        d = s1[4] + s1[4]
    end
    a, b, c, d
end
function eval_end(::RogersTanimoto, a)
    numerator = 2(a[2] + a[3])
    denominator = a[1] + a[4] + 2(a[2] + a[3])
    numerator / denominator
end
const rogerstanimoto = RogersTanimoto()

# Deviations

(::MeanAbsDeviation)(a, b) = cityblock(a, b) / length(a)
const meanad = MeanAbsDeviation()

(::MeanSqDeviation)(a, b) = sqeuclidean(a, b) / length(a)
const msd = MeanSqDeviation()

(::RMSDeviation)(a, b) = sqrt(MeanSqDeviation()(a, b))
const rmsd = RMSDeviation()

function (::NormRMSDeviation)(a, b)
    amin, amax = extrema(a)
    return RMSDeviation()(a, b) / (amax - amin)
end
const nrmsd = NormRMSDeviation()

###########################################################
#
#  Special method
#
###########################################################

# SqEuclidean/Euclidean
function _pairwise!(r::AbstractMatrix, dist::Union{SqEuclidean,Euclidean},
                    a::AbstractMatrix, b::AbstractMatrix)
    require_one_based_indexing(r, a, b)
    m, na, nb = get_pairwise_dims(r, a, b)
    # the following checks if a'*b can be stored in r directly, it fails for complex eltypes
    inplace = promote_type(eltype(r), typeof(oneunit(eltype(a))'oneunit(eltype(b)))) === eltype(r)
    R = inplace ? mul!(r, a', b) : a'b
    sa2 = sum(abs2, a, dims=1)
    sb2 = sum(abs2, b, dims=1)
    z² = zero(real(eltype(R)))
    if dist.thresh <= 0
        # If there's no chance of triggering the threshold, we can use @simd
        for j = 1:nb
            sb = sb2[j]
            for i = 1:na
                r[i, j] = eval_end(dist, (max(sa2[i] + sb - 2real(R[i, j]), z²)))
            end
        end
    else
        for j = 1:nb
            sb = sb2[j]
            for i = 1:na
                selfterms = sa2[i] + sb
                v = max(selfterms - 2real(R[i, j]), z²)
                if v < dist.thresh * selfterms
                    # The distance is likely to be inaccurate, recalculate directly
                    # This reflects the following:
                    #   while sqrt(x+ϵ) ≈ sqrt(x) + O(ϵ/sqrt(x)) when |x| >> ϵ,
                    #         sqrt(x+ϵ) ≈ O(sqrt(ϵ))             otherwise.
                    v = zero(v)
                    for k = 1:m
                        v += abs2(a[k, i] - b[k, j])
                    end
                end
                r[i, j] = eval_end(dist, v)
            end
        end
    end
    r
end

function _pairwise!(r::AbstractMatrix, dist::Union{SqEuclidean,Euclidean}, a::AbstractMatrix)
    require_one_based_indexing(r, a)
    m, n = get_pairwise_dims(r, a)
    inplace = promote_type(eltype(r), typeof(oneunit(eltype(a))'oneunit(eltype(a)))) === eltype(r)
    # the following checks if a'*b can be stored in r directly, it fails for complex eltypes
    R = inplace ? mul!(r, a', a) : a'a
    sa2 = sum(abs2, a, dims=1)
    safe = dist.thresh <= 0
    z² = zero(real(eltype(R)))
    for j = 1:n
        for i = 1:(j - 1)
            r[i, j] = r[j, i]
        end
        r[j, j] = zero(eltype(r))
        sa2j = sa2[j]
        if safe
            for i = (j + 1):n
                r[i, j] = eval_end(dist, (max(sa2[i] + sa2j - 2real(R[i, j]), z²)))
            end
        else
            for i = (j + 1):n
                selfterms = sa2[i] + sa2j
                v = max(selfterms - 2real(R[i, j]), z²)
                if v < dist.thresh * selfterms
                    v = zero(v)
                    for k = 1:m
                        v += abs2(a[k, i] - a[k, j])
                    end
                end
                r[i, j] = eval_end(dist, v)
            end
        end
    end
    r
end

# Weighted SqEuclidean/Euclidean
function _pairwise!(r::AbstractMatrix, dist::Union{WeightedSqEuclidean,WeightedEuclidean},
                    a::AbstractMatrix, b::AbstractMatrix)
    require_one_based_indexing(r, a, b)
    w = dist.weights
    m, na, nb = get_pairwise_dims(length(w), r, a, b)

    sa2 = wsumsq_percol(w, a)
    sb2 = wsumsq_percol(w, b)
    # the following checks if a'*b can be stored in r directly, it fails for complex eltypes
    inplace = promote_type(eltype(r), typeof(oneunit(eltype(a))'oneunit(eltype(b)))) === eltype(r)
    R = inplace ? mul!(r, a', w .* b) : a'*Diagonal(w)*b
    z² = zero(real(eltype(R)))
    for j = 1:nb
        for i = 1:na
            r[i, j] = eval_end(dist, max(sa2[i] + sb2[j] - 2real(R[i, j]), z²))
        end
    end
    r
end
function _pairwise!(r::AbstractMatrix, dist::Union{WeightedSqEuclidean,WeightedEuclidean},
                    a::AbstractMatrix)
    require_one_based_indexing(r, a)
    w = dist.weights
    m, n = get_pairwise_dims(length(w), r, a)

    sa2 = wsumsq_percol(w, a)
    # the following checks if a'*b can be stored in r directly, it fails for complex eltypes
    inplace = promote_type(eltype(r), typeof(oneunit(eltype(a))'oneunit(eltype(a)))) === eltype(r)
    R = inplace ? mul!(r, a', w .* a) : a'*Diagonal(w)*a
    z² = zero(real(eltype(R)))

    for j = 1:n
        for i = 1:(j - 1)
            r[i, j] = r[j, i]
        end
        r[j, j] = zero(eltype(r))
        for i = (j + 1):n
            r[i, j] = eval_end(dist, max(sa2[i] + sa2[j] - 2real(R[i, j]), z²))
        end
    end
    r
end

# MeanSqDeviation, RMSDeviation, NormRMSDeviation
function _pairwise!(r::AbstractMatrix, dist::MeanSqDeviation, a::AbstractMatrix, b::AbstractMatrix)
    _pairwise!(r, SqEuclidean(), a, b)
    # TODO: Replace by rdiv!(r, size(a, 1)) once julia compat ≥v1.2
    s = size(a, 1)
    for I in eachindex(r)
        r[I] /= s
    end
    return r
end
_pairwise!(r::AbstractMatrix, dist::RMSDeviation, a::AbstractMatrix, b::AbstractMatrix) =
    sqrt!(_pairwise!(r, MeanSqDeviation(), a, b))
function _pairwise!(r::AbstractMatrix, dist::NormRMSDeviation, a::AbstractMatrix, b::AbstractMatrix)
    _pairwise!(r, RMSDeviation(), a, b)
    @views for (i, j) in zip(axes(r, 1), axes(a, 2))
        amin, amax = extrema(a[:,j])
        r[i,:] ./= amax - amin
    end
    return r
end

function _pairwise!(r::AbstractMatrix, dist::MeanSqDeviation, a::AbstractMatrix)
    _pairwise!(r, SqEuclidean(), a)
    # TODO: Replace by rdiv!(r, size(a, 1)) once julia compat ≥v1.2
    s = size(a, 1)
    for I in eachindex(r)
        r[I] /= s
    end
    return r
end
_pairwise!(r::AbstractMatrix, dist::RMSDeviation, a::AbstractMatrix) =
    sqrt!(_pairwise!(r, MeanSqDeviation(), a))
function _pairwise!(r::AbstractMatrix, dist::NormRMSDeviation, a::AbstractMatrix)
    _pairwise!(r, RMSDeviation(), a)
    @views for (i, j) in zip(axes(r, 1), axes(a, 2))
        amin, amax = extrema(a[:,j])
        r[i,:] ./= amax - amin
    end
    return r
end

# CosineDist

function _pairwise!(r::AbstractMatrix, ::CosineDist, a::AbstractMatrix, b::AbstractMatrix)
    require_one_based_indexing(r, a, b)
    m, na, nb = get_pairwise_dims(r, a, b)
    inplace = promote_type(eltype(r), typeof(oneunit(eltype(a))'oneunit(eltype(b)))) === eltype(r)
    R = inplace ? mul!(r, a', b) : a'b
    ra = norm_percol(a)
    rb = norm_percol(b)
    for j = 1:nb
        for i = 1:na
            r[i, j] = max(1 - R[i, j] / (ra[i] * rb[j]), 0)
        end
    end
    r
end
function _pairwise!(r::AbstractMatrix, ::CosineDist, a::AbstractMatrix)
    require_one_based_indexing(r, a)
    m, n = get_pairwise_dims(r, a)
    inplace = promote_type(eltype(r), typeof(oneunit(eltype(a))'oneunit(eltype(a)))) === eltype(r)
    R = inplace ? mul!(r, a', a) : a'a
    ra = norm_percol(a)
    for j = 1:n
        for i = 1:(j - 1)
            r[i, j] = r[j, i]
        end
        r[j, j] = zero(eltype(r))
        for i = j + 1:n
            r[i, j] = max(1 - R[i, j] / (ra[i] * ra[j]), 0)
        end
    end
    r
end

# CorrDist
# This part of codes is accelerated because:
# 1. It calls the accelerated `_pairwise` specilization for CosineDist
# 2. pre-calculated `_centralize_colwise` avoids four times of redundant computations
#    of `_centralize` -- ~4x speed up
_centralize_colwise(x::AbstractMatrix) = x .- mean(x, dims=1)
_pairwise!(r::AbstractMatrix, ::CorrDist, a::AbstractMatrix, b::AbstractMatrix) =
    _pairwise!(r, CosineDist(), _centralize_colwise(a), _centralize_colwise(b))
_pairwise!(r::AbstractMatrix, ::CorrDist, a::AbstractMatrix) =
    _pairwise!(r, CosineDist(), _centralize_colwise(a))


text:
"""
    sample(A::AbstractArray{<:Real,N}, rn::Real) -> CartesianIndex{N}

Select a cartesian index of given weight array `A` and random number `rn`,
where `A` should be an `AbstractArray` with element type of `Real`.

# Example

```jldoctest
julia> A = reshape(1:4, 2, 2)
2×2 reshape(::UnitRange{Int64}, 2, 2) with eltype Int64:
 1  3
 2  4

julia> EvolutionaryModelingTools.sample(A, 1)
CartesianIndex(1, 1)

julia> EvolutionaryModelingTools.sample(A, 2)
CartesianIndex(2, 1)

julia> EvolutionaryModelingTools.sample(A, 4)
CartesianIndex(1, 2)

julia> EvolutionaryModelingTools.sample(A, 7)
CartesianIndex(2, 2)
```
"""
sample
# sample for number
sample(::Real, ::Real) = CartesianIndex()
# sample for scalar
sample(::AbstractArray{<:Real,0}, ::Real) = CartesianIndex()
# sample for vector
sample(V::AbstractVector{<:Real}, rn::Real) = CartesianIndex(_sample(V, rn))
# sample for multi-dimensional array
function sample(A::AbstractArray{<:Real}, rn::Real)
    ind = _sample(vec(A), rn)
    return CartesianIndex(ind2sub(A, ind))
end

# convert linear index `ind` of  `A` to cartesian index `sub`
ind2sub(A::AbstractArray, ind::Integer) = Base._ind2sub(axes(A), ind)

# Select a index in `1:length(V)` with weight `V` and random number `rn`
# This method is a modification of `sample([rng], wv::AbstractWeights)` of
# `StatsBase.jl` [MIT License](https://github.com/JuliaStats/StatsBase.jl)
# Copyright (c) 2012-2016: Dahua Lin, Simon Byrne, Andreas Noack, Douglas Bates,
# John Myles White, Simon Kornblith, and other contributors.
function _sample(V, rn::Real) # sample with a given random number
    n = length(V)
    i = 1
    cw = V[1]
    while cw < rn && i < n
        i += 1
        @inbounds cw += V[i]
    end
    return i
end


text:
function cholesky(mat::Matrix)
    # if the matrix is no a square matrix, return one with zeros
    if size(mat, 1) != size(mat, 2)
        return zeros(size(mat, 1), size(mat, 1))
    end
    # since verifing that the matrix is positive definite is so
    # time-consuming, the step is skipped
    # it's assumed that the matrix has $n$ rows and columns
    n = size(mat, 1)
    # $L$ is the matrix to be returned
    L = zeros(n, n)
    # the formula is proved by induction, basically
    for itr1 in 1 : n, itr2 in itr1 : n
        if itr2 == itr1
            L[itr2, itr1] = mat[itr2, itr1]
            for itr3 in 1 : (itr1 - 1)
                L[itr2, itr1] -= L[itr2, itr3] ^ 2
            end
            L[itr2, itr1] = sqrt(L[itr2, itr1])
        else
            L[itr2, itr1] = mat[itr2, itr1]
            for itr3 in 1 : (itr1 - 1)
                L[itr2, itr1] -= L[itr2, itr3] * L[itr1, itr3]
            end
            L[itr2, itr1] /= L[itr1, itr1]
        end
    end
    return L
end

max_stars_repo_name: Shimushushushu/Etude-Rhapsody

max_issues_repo_name: Shimushushushu/Etude-Rhapsody

max_forks_repo_name: Shimushushushu/Etude-Rhapsody


text:
try
    EcologicalNetworks._check_species_validity(::Type{MangalNode}) = nothing
    EcologicalNetworks._check_species_validity(::Type{MangalReferenceTaxon}) = nothing
catch e
    @info "Compatibility with EcologicalNetworks 0.3 will be removed soon"
    EcologicalNetworks.check_species_validity(::Type{MangalNode}) = nothing
    EcologicalNetworks.check_species_validity(::Type{MangalReferenceTaxon}) = nothing
end

function get_all_interactions(n::MangalNetwork, query::Pair...)
    page_size = 250
    network_interactions = MangalInteraction[]
    interactions_to_get = count(MangalInteraction, n, query...)
    pages_to_do = convert(Int64, ceil(interactions_to_get/page_size))
    for page in 1:pages_to_do
        paging_query = ["page" => "$(page-1)", "count" => page_size, query...]
        append!(network_interactions, Mangal.interactions(n, paging_query...))
    end
    return network_interactions
end

import Base.convert

"""
    convert(::Type{EcologicalNetworks.UnipartiteNetwork}, n::MangalNetwork, query::Pair...)

TODO
"""
function convert(::Type{EcologicalNetworks.UnipartiteNetwork}, n::MangalNetwork, query::Pair...)
    all_interactions = Mangal.get_all_interactions(n, query...)
    return convert(EcologicalNetworks.UnipartiteNetwork, all_interactions)
end

"""
    convert(::Type{EcologicalNetworks.UnipartiteNetwork}, n::MangalNetwork, query::Pair...)

TODO
"""
function convert(::Type{EcologicalNetworks.UnipartiteQuantitativeNetwork}, n::MangalNetwork, query::Pair...)
    all_interactions = Mangal.get_all_interactions(n, query...)
    return convert(EcologicalNetworks.UnipartiteQuantitativeNetwork, all_interactions)
end

"""
    convert(::Type{UnipartiteNetwork}, interac::Vector{MangalInteraction})

TODO
"""
function convert(::Type{EcologicalNetworks.UnipartiteNetwork}, interac::Vector{MangalInteraction})

    all_object_nodes = MangalNode[]

    for i in interac
        append!(all_object_nodes, [i.from, i.to])
    end

    object_nodes = unique(all_object_nodes)
    S = length(object_nodes)
    A = zeros(Bool, (S,S))
    N = EcologicalNetworks.UnipartiteNetwork(A, object_nodes)
    for i in interac
        if i.strength != 0
            N[i.from, i.to] = true
        end
    end

    return N
end

"""
    convert(::Type{UnipartiteNetwork}, interac::Vector{MangalInteraction})

TODO
"""
function convert(::Type{EcologicalNetworks.UnipartiteQuantitativeNetwork}, interac::Vector{MangalInteraction})

    all_object_nodes = MangalNode[]

    for i in interac
        append!(all_object_nodes, [i.from, i.to])
    end

    object_nodes = unique(all_object_nodes)
    S = length(object_nodes)
    A = zeros(Float64, (S,S))
    N = EcologicalNetworks.UnipartiteQuantitativeNetwork(A, object_nodes)
    for i in interac
        if i.strength != 0
            N[i.from, i.to] = convert(Float64, i.strength)
        end
    end

    return N
end

"""
missing are dropped
"""
function taxonize(N::T) where {T <: EcologicalNetworks.UnipartiteNetwork}
    @assert eltype(species(N)) == MangalNode
    unique_ref_taxa = unique([s.taxon for s in EcologicalNetworks.species(N)])
    @warn "This function really should inform of dropped nodes"
    filter!(!ismissing, unique_ref_taxa)
    u_refs = convert(Vector{MangalReferenceTaxon}, unique_ref_taxa)
    S = length(u_refs)
    A = zeros(Bool, (S,S))
    K = EcologicalNetworks.UnipartiteNetwork(A, u_refs)
    for s1 in EcologicalNetworks.species(N)
        if !ismissing(s1.taxon)
            for s2 in EcologicalNetworks.species(N)
                if !ismissing(s2.taxon)
                    K[s1.taxon,s2.taxon] = N[s1,s2]
                end
            end
        end
    end
    return EcologicalNetworks.simplify(K)
end


text:
set_global_client()

@testset "SUBSCRIBE" begin
    channels = ["first", "second", "third"]
    publisher = Client()
    subscriber = Client()
    messages = []
    @test subscriber.is_subscribed == false

    @async subscribe(channels...; client=subscriber) do msg
        push!(messages, msg)
    end

    wait_until_subscribed(subscriber)
    @test subscriber.is_subscribed == true
    @test subscriber.subscriptions == Set{String}(channels)

    @test publish("first", "hello"; client=publisher) == 1
    @test publish("second", "world"; client=publisher) == 1
    @test publish("something", "else"; client=publisher) == 0

    @test length(messages) == 2
    @test messages[1] == ["message", "first", "hello"]
    @test messages[2] == ["message", "second", "world"]

    @test_throws RedisError set("already", "subscribed"; client=subscriber)
    @test_throws RedisError subscribe("alreadsubscribed"; client=subscriber) do msg end

    unsubscribe("first"; client=subscriber)

    wait_until_channel_unsubscribed(subscriber, "first")
    @test subscriber.is_subscribed == true
    @test subscriber.subscriptions == Set{String}(["second", "third"])

    @test publish("first", "not subscribed anymore"; client=publisher) == 0

    @test length(messages) == 2

    unsubscribe(; client=subscriber) # unsubscribe from everything

    wait_until_unsubscribed(subscriber)
    @test subscriber.is_subscribed == false
    @test isempty(subscriber.subscriptions)

    stop_fn(msg) = msg[end] == "close subscription"

    @async subscribe(channels...; stop_fn=stop_fn, client=subscriber) do msg end

    wait_until_subscribed(subscriber)
    @test subscriber.is_subscribed == true
    @test subscriber.subscriptions == Set{String}(channels)

    @test publish("first", "close subscription"; client=publisher) == 1

    wait_until_unsubscribed(subscriber)
    @test subscriber.is_subscribed == false
    @test isempty(subscriber.subscriptions)
end

@testset "PSUBSCRIBE" begin
    patterns = ["first*", "second*", "third*"]
    publisher = Client()
    subscriber = Client()
    messages = []
    @test subscriber.is_subscribed == false

    @async psubscribe(patterns...; client=subscriber) do msg
        push!(messages, msg)
    end

    wait_until_subscribed(subscriber)
    @test subscriber.is_subscribed == true
    @test subscriber.psubscriptions == Set{String}(patterns)

    @test publish("first_pattern", "hello"; client=publisher) == 1
    @test publish("second_pattern", "world"; client=publisher) == 1
    @test publish("something", "else"; client=publisher) == 0

    @test length(messages) == 2
    @test messages[1] == ["pmessage", "first*", "first_pattern", "hello"]
    @test messages[2] == ["pmessage", "second*", "second_pattern", "world"]

    @test_throws RedisError set("already", "subscribed"; client=subscriber)
    @test_throws RedisError psubscribe("alreadsubscribed"; client=subscriber) do msg end

    punsubscribe("first*"; client=subscriber)

    wait_until_pattern_unsubscribed(subscriber, "first*")
    @test subscriber.is_subscribed == true
    @test subscriber.psubscriptions == Set{String}(["second*", "third*"])

    @test publish("first_pattern", "not subscribed anymore"; client=publisher) == 0

    @test length(messages) == 2

    punsubscribe(; client=subscriber) # unsubscribe from everything

    wait_until_unsubscribed(subscriber)
    @test subscriber.is_subscribed == false
    @test isempty(subscriber.psubscriptions)

    stop_fn(msg) = msg[end] == "close subscription"

    @async psubscribe(patterns...; stop_fn=stop_fn, client=subscriber) do msg end

    wait_until_subscribed(subscriber)
    @test subscriber.is_subscribed == true
    @test subscriber.psubscriptions == Set{String}(patterns)

    @test publish("first_pattern", "close subscription"; client=publisher) == 1

    wait_until_unsubscribed(subscriber)
    @test subscriber.is_subscribed == false
    @test isempty(subscriber.subscriptions)
end

flushall()


text:
using Distributions
using Klara
using PGUManifoldMC

DATADIR = "../../data"
SUBDATADIR = "smmala"

nchains = 10
nmcmc = 110000
nburnin = 10000

covariates, = dataset("swiss", "measurements");
ndata, npars = size(covariates);

covariates = (covariates.-mean(covariates, 1))./repmat(std(covariates, 1), ndata, 1);

outcome, = dataset("swiss", "status");
outcome = vec(outcome);

function ploglikelihood(p::Vector{Float64}, v::Vector)
  Xp = v[2]*p
  dot(Xp, v[3])-sum(log(1+exp(Xp)))
end

plogprior(p::Vector{Float64}, v::Vector) = -0.5*(dot(p, p)/v[1]+npars*log(2*pi*v[1]))

pgradlogtarget(p::Vector{Float64}, v::Vector) = v[2]'*(v[3]-1./(1+exp(-v[2]*p)))-p/v[1]

function ptensorlogtarget(p::Vector{Float64}, v::Vector)
  r = 1./(1+exp(-v[2]*p))
  broadcast(*, r.*(1-r), v[2])'*v[2]+(eye(npars)/v[1])
end

p = BasicContMuvParameter(
  :p,
  loglikelihood=ploglikelihood,
  logprior=plogprior,
  gradlogtarget=pgradlogtarget,
  tensorlogtarget=ptensorlogtarget,
  nkeys=4
)

model = likelihood_model([Hyperparameter(:λ), Data(:X), Data(:y), p], isindexed=false)

sampler = SMMALA(0.02)

mcrange = BasicMCRange(nsteps=nmcmc, burnin=nburnin)

outopts = Dict{Symbol, Any}(:monitor=>[:value], :diagnostics=>[:accept])

times = Array(Float64, nchains)
stepsizes = Array(Float64, nchains)
i = 1

while i <= nchains
  v0 = Dict(:λ=>100., :X=>covariates, :y=>outcome, :p=>rand(Normal(0, 3), npars))

  job = BasicMCJob(
    model,
    sampler,
    mcrange,
    v0,
    tuner=AcceptanceRateMCTuner(0.7, score=x -> logistic_rate_score(x, 3.), verbose=false),
    outopts=outopts
  )

  tic()
  run(job)
  runtime = toc()

  chain = output(job)
  ratio = acceptance(chain)

  if 0.65 < ratio < 0.75
    writedlm(joinpath(DATADIR, SUBDATADIR, "chain"*lpad(string(i), 2, 0)*".csv"), chain.value, ',')
    writedlm(joinpath(DATADIR, SUBDATADIR, "diagnostics"*lpad(string(i), 2, 0)*".csv"), vec(chain.diagnosticvalues), ',')

    times[i] = runtime
    stepsizes[i] = job.sstate.tune.step

    println("Iteration ", i, " of ", nchains, " completed with acceptance ratio ", ratio)
    i += 1
  end
end

writedlm(joinpath(DATADIR, SUBDATADIR, "times.csv"), times, ',')
writedlm(joinpath(DATADIR, SUBDATADIR, "stepsizes.csv"), stepsizes, ',')


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    Condition c;

    void foo() {
	c = new Condition();
    }

    public static void main(String[] args) {}
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module RunnerTests
using Test

include("../../../src/julia/shared/load.jl")
include("../../../src/julia/runner/load.jl")
import ADPerfTest
using GMMData
import TestLoader
import Benchmark

mockgmmcalcname = "Mock" * "GMM"

dir = @__DIR__
if !(dir ∈ LOAD_PATH)
    push!(LOAD_PATH, dir)
end

mutable struct Counter
    count::Int
    Counter() = new(0)
end

@testset "Runner Tests" begin
    @testset "Module Loading" begin
        @test isa(TestLoader.get_gmm_test(mockgmmcalcname), ADPerfTest.Test{GMMInput,GMMOutput})
        @test_throws ArgumentError TestLoader.get_ba_test(mockgmmcalcname)
        test = TestLoader.get_gmm_test(mockgmmcalcname)
        output = empty_gmm_output()
        test.output!(output, test.context)
        @test output.objective == 1.0
        test = TestLoader.get_gmm_test(mockgmmcalcname * "2")
        output = empty_gmm_output()
        test.output!(output, test.context)
        @test output.objective == 2.0
    end

    @testset "Time Limit" begin
        # Run_count guarantees total time greater than the time_limit
        run_count = 100
        time_limit = 0.1
        # Execution_time should be more than minimum_measurable_time
        # because we can expect find_repeats_for_minimum_measurable_time
        # to call calculate_objective function only once in that case.
        minimum_measurable_time = 0.0
        execution_time = 0.01
        ctx = Counter()
        objective(ctx, ::Int) = begin ctx.count += 1; sleep(execution_time) end
        Benchmark.measure_shortest_time!(ctx, minimum_measurable_time, run_count, time_limit, objective)
        # Number of runs should be less then run_count variable because total_time will be reached.
        @test ctx.count ≤ ceil(time_limit / execution_time)
    end

    @testset "Number of Runs Limit" begin
        # Run_count guarantees total time lesser than the time_limit
        run_count = 10
        time_limit = 10.0
        # Execution_time should be more than minimum_measurable_time
        # because we can expect find_repeats_for_minimum_measurable_time
        # to call calculate_objective function only once in that case.
        minimum_measurable_time = 0.0
        execution_time = 0.01
        ctx = Counter()
        objective(ctx, ::Int) = begin ctx.count += 1; sleep(execution_time) end
        Benchmark.measure_shortest_time!(ctx, minimum_measurable_time, run_count, time_limit, objective)
        # Number of runs should be equal to run_count variable because total_time won't be reached.
        @test ctx.count == run_count
    end

    @testset "Time Measurement" begin
        # Run_count guarantees total time lesser than the time_limit
        run_count = 10
        time_limit = 100000.0
        # Execution_time should be more than minimum_measurable_time
        # because we can expect find_repeats_for_minimum_measurable_time
        # to call calculate_objective function only once in that case.
        minimum_measurable_time = 0.0
        execution_time = 0.01
        ctx = Counter()
        objective(ctx, ::Int) = begin ctx.count += 1; sleep(execution_time) end
        shortest_time = Benchmark.measure_shortest_time!(ctx, minimum_measurable_time, run_count, time_limit, objective)
        # Number of runs should be equal to run_count variable because total_time won't be reached.
        @test ctx.count == run_count
        @test shortest_time ≥ execution_time
    end

    @testset "Search for Repeats" begin
        assumed_repeats = 16
        execution_time = 0.01
        minimum_measurable_time = execution_time * assumed_repeats
        ctx = Counter()
        objective(ctx, repeats) = begin sleep(repeats * execution_time) end
        repeats, _, _ = Benchmark.find_repeats_for_minimum_measurable_time!(ctx, minimum_measurable_time, objective)
        @test repeats ≠ Benchmark.measurable_time_not_achieved
        @test repeats ≤ assumed_repeats
    end

    @testset "Repeats Not Found" begin
        minimum_measurable_time = 1000.0
        ctx = Counter()
        objective(ctx, repeats) = nothing
        repeats, _, _ = Benchmark.find_repeats_for_minimum_measurable_time!(ctx, minimum_measurable_time, objective)
        @test repeats == Benchmark.measurable_time_not_achieved
    end
end
end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
function TurnAroundRange(arg0::GroundStation, arg1::GroundStation, arg2::AbsoluteDate, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::ObservableSatellite)
    return TurnAroundRange((GroundStation, GroundStation, AbsoluteDate, jdouble, jdouble, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function get_master_station(obj::TurnAroundRange)
    return jcall(obj, "getMasterStation", GroundStation, ())
end

function get_slave_station(obj::TurnAroundRange)
    return jcall(obj, "getSlaveStation", GroundStation, ())
end


text:
# Parallel Map, copy/pasted from Julia with minor adaption
function pmaplm(f, lsts...; err_retry=true, err_stop=false, ncores = 1)
    @show ncores
    len = length(lsts)

    @compat results = Dict{Int,Any}()

    retryqueue = Any[]
    task_in_err = false
    is_task_in_error() = task_in_err
    set_task_in_error() = (task_in_err = true)

    nextidx = 0
    getnextidx() = (nextidx += 1)

    states = [start(lsts[idx]) for idx in 1:len]
    function getnext_tasklet()
        if is_task_in_error() && err_stop
            return nothing
        elseif !any(idx->done(lsts[idx],states[idx]), 1:len)
            nxts = [next(lsts[idx],states[idx]) for idx in 1:len]
            for idx in 1:len; states[idx] = nxts[idx][2]; end
            nxtvals = [x[1] for x in nxts]
            return (getnextidx(), nxtvals)
        elseif !isempty(retryqueue)
            return shift!(retryqueue)
        else
            return nothing
        end
    end

    @sync begin
        for wpid in workers()[1:ncores]
            @async begin
                tasklet = getnext_tasklet()
                while (tasklet != nothing)
                    (idx, fvals) = tasklet
                    try
                        result = remotecall_fetch(wpid, f, fvals...)
                        if isa(result, Exception)
                            ((wpid == myid()) ? rethrow(result) : throw(result))
                        else
                            results[idx] = result
                        end
                    catch ex
                        if err_retry
                            push!(retryqueue, (idx,fvals, ex))
                        else
                            results[idx] = ex
                        end
                        set_task_in_error()
                        break # remove this worker from accepting any more tasks
                    end

                    tasklet = getnext_tasklet()
                end
            end
        end
    end

    for failure in retryqueue
        results[failure[1]] = failure[3]
    end
    [results[x] for x in 1:nextidx]
end


text:
module GaussLaw1D

export Problem, solve, Constant

using OrdinaryDiffEq: Vern7
using ArgCheck: @argcheck
using BoundaryValueDiffEq: BoundaryValueDiffEq
const BV = BoundaryValueDiffEq
using Unitful: Unitful
const UF = Unitful

# export units for convenience
using Unitful: V, kV
export V, kV
using Unitful: pC, nC, μC, mC, C
export pC, nC, μC, mC, C
using Unitful: F
using Unitful: nm, μm, mm, cm, m
export nm, μm, mm, cm, m
using Unitful: ns, μs, ms, s
export ns, μs, ms, s
using Unitful: nGy, μGy, mGy, Gy
export nGy, μGy, mGy, Gy
using UnitfulRecipes # enable unitful recipes for convenience

const VACUUM_PERMITIVITY = 8.8541878128e-12*F * m^-1
const U_LENGTH = mm
const U_POT = V
const U_E = U_POT / U_LENGTH

const Length = typeof(1.0U_LENGTH)
const Potential = typeof(1.0U_POT)
const EFieldStrength = typeof(1.0U_E)

struct Problem{Fn}
    charge_density::Fn
    limits::NTuple{2, Length}
    shape::Symbol
    potential::NTuple{2, Potential}
    electric_permitivity::typeof(VACUUM_PERMITIVITY)
    radii::Vector{Length}
end

struct Constant{X}
    value::X
end
(c::Constant)(r) = c.value

function Problem(charge_density;
                 potential,
                 limits,
                 shape,
        electric_permitivity=VACUUM_PERMITIVITY,
        radii=collect(range(limits[1], stop=limits[2], length=1000)),
    )
    @argcheck issorted(limits)
    @argcheck limits[1] > 0.0mm
    @argcheck shape in (:parallel, :cylindrical, :spherical)
    r = limits[1]
    cd = charge_density(r)
    actual_cd_dim = UF.dimension(cd)
    expected_cd_dim = UF.dimension(1.0C/(mm^3))
    if  actual_cd_dim != expected_cd_dim
        msg = """Charge density $cd has the wrong unit. Expected unit is C*mm^-3.
        """
        throw(ArgumentError(msg))
    end
    Problem(charge_density, limits, shape, potential, electric_permitivity, radii)
end

function to_compute_unit(pot::Potential)
    UF.ustrip(U_POT, pot)
end

@noinline function bc!(residual, u, w, r)
    # @show residual
    # @show u
    # @show w
    # @show r
    # error()
    wanted_phi_left, wanted_phi_right = w.potential
    actual_phi_left = u[1][1]
    actual_phi_right = u[end][1]
    residual[1] = actual_phi_left  - UF.ustrip(U_POT, wanted_phi_left)
    residual[2] = actual_phi_right - UF.ustrip(U_POT, wanted_phi_right)
    residual
end

@noinline function step!(du, u, w, r0)
    # @show du
    # @show u
    # @show w
    # @show r0
    # error()
    r = r0*U_LENGTH
    phi0, E0 = u
    E = E0 * U_E
    rho = w.charge_density(r)
    dphi = -E0
    dEdr = compute_dEdr(E, rho, r, w)
    du[1] = dphi
    du[2] = UF.ustrip(U_E/U_LENGTH, dEdr)
    du
end

function compute_dEdr(E, rho, r, w)
    shape = w.shape
    ε = w.electric_permitivity
    ret = rho/ε
    if shape == :parallel
        ret -= 0E/r
    elseif shape == :cylindrical
        ret -= 1E/r
    elseif shape == :spherical
        ret -= 2E/r
    else
        error("This error should be unreachable. Unknown shape $shape")
    end
    ret
end

# count = 0
function initial_guess(w, r_strip)
    # @show r_strip
    r = U_LENGTH * r_strip
    phi_left, phi_right = w.potential
    r_left, r_right = w.limits
    E = (phi_right - phi_left) / (r_right - r_left)
    w_left = (r_right - r) / (r_right - r_left)
    w_right = (r - r_left) / (r_right - r_left)
    phi = w_left * phi_left + w_right * phi_right
    # global count += 1
    # @show count
    # if mod(count, 10000) == 0
    #     error()
    # end
    [UF.ustrip(U_POT, phi), UF.ustrip(U_E,E)]
end

function build_problem(w::Problem)
    r1, r2 = w.limits
    rspan = (UF.ustrip(U_LENGTH, r1)::Float64, UF.ustrip(U_LENGTH,r2)::Float64)
    ret = BV.BVProblem(step!, bc!, initial_guess, rspan, w)
    # ret = BV.TwoPointBVProblem(step!, bc!, initial_guess(w, w.limits[1]/U_LENGTH), rspan, w)
end

function solve(w::Problem,
        solver=BV.Shooting(Vern7()),
        ; kw...)
    prob = build_problem(w)
    sol = BV.solve(prob, solver; kw...)

    itp = sol.(UF.ustrip.(U_LENGTH, w.radii))
    potential = first.(itp) .* U_POT
    E = last.(itp) .* U_E
    return (E=E,
        potential=potential,
        rho=w.charge_density.(w.radii),
        radii=w.radii,
        raw=sol)
end

end # module


text:
module ConstructionBase

export getproperties
export setproperties
export constructorof

# Use markdown files as docstring:
for (name, path) in [
    :ConstructionBase => joinpath(dirname(@__DIR__), "README.md"),
    :constructorof => joinpath(@__DIR__, "constructorof.md"),
    :getproperties => joinpath(@__DIR__, "getproperties.md"),
    :setproperties => joinpath(@__DIR__, "setproperties.md"),
]
    # Don't fail when somehow importing docstrings doesn't work (we
    # don't loose any functionalities for that).  We use explicit
    # `@docs` in `docs/src/index.md` to make sure importing docstrings
    # succeeds in CI.
    try
        include_dependency(path)
        str = read(path, String)
        @eval @doc $str $name
    catch err
        @error "Failed to import docstring for $name" exception=(err, catch_backtrace())
    end
end

@generated function constructorof(::Type{T}) where T
    getfield(parentmodule(T), nameof(T))
end

constructorof(::Type{<:Tuple}) = tuple
constructorof(::Type{<:NamedTuple{names}}) where names =
    NamedTupleConstructor{names}()

struct NamedTupleConstructor{names} end

@generated function (::NamedTupleConstructor{names})(args...) where names
    quote
        Base.@_inline_meta
        $(NamedTuple{names, Tuple{args...}})(args)
    end
end

getproperties(o::NamedTuple) = o
getproperties(o::Tuple) = o
@generated function getproperties(obj)
    fnames = fieldnames(obj)
    fvals = map(fnames) do fname
        Expr(:call, :getproperty, :obj, QuoteNode(fname))
    end
    fvals = Expr(:tuple, fvals...)
    :(NamedTuple{$fnames}($fvals))
end

################################################################################
##### setproperties
################################################################################
function setproperties(obj; kw...)
    setproperties(obj, (;kw...))
end

setproperties(obj             , patch::Tuple      ) = setproperties_object(obj     , patch )
setproperties(obj             , patch::NamedTuple ) = setproperties_object(obj     , patch )
setproperties(obj::NamedTuple , patch::Tuple      ) = setproperties_namedtuple(obj , patch )
setproperties(obj::NamedTuple , patch::NamedTuple ) = setproperties_namedtuple(obj , patch )
setproperties(obj::Tuple      , patch::Tuple      ) = setproperties_tuple(obj      , patch )
setproperties(obj::Tuple      , patch::NamedTuple ) = setproperties_tuple(obj      , patch )

setproperties_namedtuple(obj, patch::Tuple{}) = obj
@noinline function setproperties_namedtuple(obj, patch::Tuple)
    msg = """
    setproperties(obj::NamedTuple, patch::Tuple) only allowed for empty Tuple. Got:
    obj = $obj
    patch = $patch
    """
    throw(ArgumentError(msg))
end
function setproperties_namedtuple(obj, patch)
    res = merge(obj, patch)
    validate_setproperties_result(res, obj, obj, patch)
    res
end
function validate_setproperties_result(
    nt_new::NamedTuple{fields}, nt_old::NamedTuple{fields}, obj, patch) where {fields}
    nothing
end
@noinline function validate_setproperties_result(nt_new, nt_old, obj, patch)
    O = typeof(obj)
    P = typeof(patch)
    msg = """
    Failed to assign properties $(fieldnames(P)) to object with fields $(fieldnames(O)).
    You may want to overload
    ConstructionBase.setproperties(obj::$O, patch::NamedTuple)
    ConstructionBase.getproperties(obj::$O)
    """
    throw(ArgumentError(msg))
end
function setproperties_namedtuple(obj::NamedTuple{fields}, patch::NamedTuple{fields}) where {fields}
    patch
end

setproperties_tuple(obj::Tuple, patch::NamedTuple{()}) = obj
@noinline function setproperties_tuple(obj::Tuple, patch::NamedTuple)
    msg = """
    setproperties(obj::Tuple, patch::NamedTuple) only allowed for empty NamedTuple. Got:
    obj  ::Tuple      = $obj
    patch::NamedTuple = $patch
    """
    throw(ArgumentError(msg))
end
function setproperties_tuple(obj::NTuple{N,Any}, patch::NTuple{N,Any}) where {N}
    patch
end
append(x,y) = (x..., y...)
@noinline function throw_setproperties_tuple_error(obj, patch)
    msg = """
    Cannot call `setproperties(obj::Tuple, patch::Tuple)` with `length(obj) < length(patch)`. Got:
    obj = $obj
    patch = $patch
    """
    throw(ArgumentError(msg))
end
function setproperties_tuple(obj::NTuple{N,Any}, patch::NTuple{K,Any}) where {N,K}
    if K > N
        throw_setproperties_tuple_error(obj, patch)
    end
    append(patch, after(obj, Val{K}()))
end
function after(xs::Tuple, ::Val{N}) where {N}
    after(Base.tail(xs), Val{N-1}())
end
function after(x::Tuple, ::Val{0})
    x
end

setproperties_object(obj, patch::Tuple{}) = obj
@noinline function setproperties_object(obj, patch::Tuple)
    msg = """
    setproperties(obj, patch::Tuple) only allowed for empty Tuple. Got:
    obj = $obj
    patch = $patch
    """
end
setproperties_object(obj, patch::NamedTuple{()}) = obj
function setproperties_object(obj, patch)
    nt = getproperties(obj)
    nt_new = merge(nt, patch)
    validate_setproperties_result(nt_new, nt, obj, patch)
    constructorof(typeof(obj))(Tuple(nt_new)...)
end

include("nonstandard.jl")
include("functions.jl")

end # module


text:
__precompile__(true)
module CrystalBalls

using StaticArrays

# abstract types, d is number of spacetime dimensions
abstract type AbstractLattice{d} end
# the type parameters give the spacetime dimensionality
abstract type AbstractField{d} end
abstract type AbstractGaugeField{d} <: AbstractField{d} end
export AbstractLattice, AbstractField, AbstractGaugeField

import Base.rand
import Base: size, getindex, sub2ind, dec

include("unitary.jl")
include("lattice.jl")
include("gaugefields.jl")
include("metropolis.jl")

end # module


text:
using SimradEK60
using SimradEK60TestData
filename = EK60_SAMPLE
ps = SimradEK60.load(filename)
ps38 = [p for p in ps if p.frequency == 38000]
Sv38 = Sv(ps38) # Volume backscatter
al38 = alongshipangle(ps38) # Split beam angle
at38 = athwartshipangle(ps38)
_R = R(ps38) # Range / depth
_t = filetime(ps38) # timestamps