text:
function fake_wordvectors(type_element::Type{T}=Float32) where T<:AbstractFloat
    vocab = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"]
    vectors = diagm(0 => T[1.0, 2, 4, 5, 12, 1, 2, 1, 0, 10])
    vocab_hash = Dict(v => i for (i,v) in enumerate(vocab))
    return WordVectors(vocab, vectors, vocab_hash)
end

@testset "write2disk: WordVectors" begin
    filepath, io = mktemp()
    for T in [Float32, Float64]
        for kind in [:text, :binary]
            # create data
            wv_orig = fake_wordvectors(T)
            # dump data
            write2disk(filepath, wv_orig, kind=kind)
            # reload data
            wv_loaded = wordvectors(filepath, T, kind=kind, normalize=false)
            # tests
            @test eltype(wv_orig.vectors) == eltype(wv_loaded.vectors)
            @test wv_orig.vectors == wv_loaded.vectors
            @test wv_orig.vocab == wv_loaded.vocab
            @test wv_orig.vocab_hash == wv_loaded.vocab_hash
        end
    end
end

@testset "write2disk: CompressedWordVectors" begin
    filepath, io = mktemp()
    for T in [Float32, Float64]
        for kind in [:text, :binary]
            # create data
            wv_orig = fake_wordvectors(T)
            # compress data
            wv_compressed = compress(wv_orig, sampling_ratio=1.0,
                                     k=2, m=5, method=:pq)
            # write compressed data & reload
            write2disk(filepath, wv_compressed, kind=kind)
            wv_loaded = compressedwordvectors(filepath, T, kind=kind)
            # start tests
            @test wv_compressed.vocab == wv_loaded.vocab
            @test wv_compressed.vocab_hash == wv_loaded.vocab_hash
            @test typeof(wv_compressed) == typeof(wv_loaded)
            @test all(wv_compressed.vectors.data .==
                      wv_loaded.vectors.data)
            @test wv_compressed.vectors.quantizer.rot ==
                wv_loaded.vectors.quantizer.rot
        end
    end
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
using ChainRulesCore: length
import DPFEHM
import GaussianRandomFields
import DifferentiableBackwardEuler
import Optim
import Random
import Zygote
import ChainRulesCore

using JLD2
using Statistics: mean, std

global losses_train = Float64[]
global losses_test = Float64[]
global rmses_train = Float64[]
global rmses_test = Float64[]
global train_time = Float64[]

@everywhere begin
    using Flux
    using ChainRulesCore
    using Random
    using GaussianRandomFields
    using DPFEHM

    n = 51
    ns = [n, n]
    steadyhead = 0e0
    sidelength = 200
    thickness  = 1.0
    mins = [-sidelength, -sidelength] #meters
    maxs = [sidelength, sidelength] #meters
    num_eigenvectors = 200
    sigma = 1.0
    lambda = 50
    cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma))
    x_pts = range(mins[1], maxs[1]; length=ns[1])
    y_pts = range(mins[2], maxs[2]; length=ns[2])
    grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, y_pts)

    pressure_target  = 1.0

    coords, neighbors, areasoverlengths, volumes = DPFEHM.regulargrid2d(mins, maxs, ns, thickness)
    dirichletnodes = Int[]
    dirichleths = zeros(size(coords, 2))
    for i = 1:size(coords, 2)
        if abs(coords[1, i]) == sidelength || abs(coords[2, i]) == sidelength
            push!(dirichletnodes, i)
            dirichleths[i] = steadyhead
        end
    end

    function getQs(Qs::Vector, rs::Vector)
        return sum(map(i->getQs(Qs[i], rs[i]), 1:length(Qs)))
    end

    function getQs(Q::Number, r::Number)#this splits the Q across all the nodes that are as close as possible to r (should be at least 4)
        dists = map(i->(sqrt(sum(coords[:, i] .^ 2)) - r) ^ 2, 1:size(coords, 2))
        mindist = minimum(dists)
        goodnodes = (dists .≈ mindist)
        Qs = Q * goodnodes / sum(goodnodes)
        return Qs
    end

    function solve_numerical(Qs, T, rs)
        logKs2Ks_neighbors(Ks) = exp.(0.5 * (Ks[map(p->p[1], neighbors)] .+ Ks[map(p->p[2], neighbors)]))
        Qs = getQs(Qs, rs)
        @assert length(T) == length(Qs)

        Ks_neighbors = logKs2Ks_neighbors(T)
        h_gw = DPFEHM.groundwater_steadystate(Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs)
        goodnode = div(size(coords, 2), 2) + 1
        @assert coords[:, goodnode] == [0, 0]#make sure we are looking at the right "good node"
        return h_gw[goodnode, end] - steadyhead
    end
end

model = Chain(
    Conv((3, 3), 1=>8, pad=(1,1), relu),
    x -> maxpool(x, (2,2)),
    Conv((3, 3), 8=>16, pad=(1,1), relu),
    x -> maxpool(x, (2,2)),
    Conv((3, 3), 16=>8, pad=(1,1), relu),
    x -> maxpool(x, (2,2)),
    flatten,
    Dense(288, 1),
) |> f64

# Make neural network parameters trackable by Flux
θ = params(model)

function loss(x)
    Ts = reshape(hcat(map(y->y[1], x)...), size(x[1][1], 1), size(x[1][1], 2), 1, length(x))
    targets = map(y->y[2], x)
    #@show size(Ts)
    Q1 = model(Ts)
    Qs = map(Q->[Q, Qinj], Q1)
    loss = sum(pmap(i->solve_numerical(Qs[i], Ts[:, :, 1, i], rs) - targets[i], 1:size(Ts, 4)).^2)
    return loss
end

function cb2()
    # Terminal output
    loss_train = sum(map(x->loss(x),data_train_batch))
    loss_test = sum(map(x->loss(x),data_test))
    println(string("callback: train loss: ", loss_train, " test loss: ", loss_test))
end

cb = Flux.throttle(cb2, 10)
opt = ADAM()

# Set ranges for P, T
# Fix pressure target for now. It should be removed, but when I removed it with MultiTheis, the PIML framework didn't work as well. Moral of the story, NN's are squirrelly things that don't always make sense.
PRES = (0.0,0.0) # Pressure target range [m]
TRANS = (-2.0,0.0) # log10 transmissivity [m^2/s]
# Assume that reservoir thickness is 100m (S = S_s*b)
STOR = (-4.0,-1.0) # log10 storativity

# Location to monitor pressure at
mon_well_crds = [[-75.0,-75.0]]
# Location of injection and extraction wells
well_crds = [[-50.0,-50.0],[50.0,50.0]]
# Injection rate
Qinj = 0.031688 # [m^3/s] (1 MMT water/yr)
# Training epochs
epochs = 1:2
# Calculate distance between extraction and injection wells
rs = [sqrt((sum((cs-mon_well_crds[1]).^2))) for cs in well_crds]

# Input/output normalization function
lerp(x, lo, hi) = x*(hi-lo)+lo

# Random T, S, and P generator
sample() = [10^lerp(rand(), TRANS...), 10^lerp(rand(), STOR...), lerp(rand(), PRES...)]

# batch 1:1 for one batch size
batch_size = 64
data_train_batch = [[(GaussianRandomFields.sample(grf), pressure_target) for i = 1:batch_size] for v in 1:4]
data_test = [[(GaussianRandomFields.sample(grf), pressure_target)] for i = 1:4]

println("The training has started..")
loss_train = sum(map(x->loss(x), data_train_batch))
rmse_train = sqrt(loss_train/(batch_size*length(data_train_batch)))
loss_test = sum(map(x->loss(x), data_test))
rmse_test = sqrt(loss_test/length(data_test))

println(string("epoch: 0 train rmse: ", rmse_train, " test rmse: ", rmse_test))
push!(losses_train, loss_train)
push!(losses_test, loss_test)
push!(rmses_train, rmse_train)
push!(rmses_test, rmse_test)
for epoch in epochs
    tt = @elapsed Flux.train!(loss, θ, data_train_batch, opt)#, cb = cb)
    push!(train_time, tt)
    loss_train = sum(map(x->loss(x), data_train_batch))
    rmse_train = sqrt(loss_train/(batch_size*length(data_train_batch)))
    loss_test = sum(map(x->loss(x), data_test))
    rmse_test = sqrt(loss_test/length(data_test))
    # Terminal output
    println(string("epoch: ", epoch, " train rmse: ", rmse_train, " test rmse: ", rmse_test))
    # Save convergence metrics
    push!(losses_train, loss_train)
    push!(losses_test, loss_test)
    push!(rmses_train, rmse_train)
    push!(rmses_test, rmse_test)
end
println("The training has finished!")

@save string("loss_data_2.jld2") epochs losses_train  rmses_train  losses_test  rmses_test  train_time
println("The data has beens saved!")

ext: jl


text:
# PS8 Dayuan Xie
# using Pkg
# Pkg.add("MultivariateStats")

using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using DataFramesMeta
using Distributions
using CSV
using MultivariateStats

function allwrap()
# Question 1
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS8-factor/nlsy.csv"
df = CSV.read(HTTP.get(url).body)

ols_lm = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr), df)
println(ols_lm)

# Question 2
matrix_asvab = hcat(df.asvabAR, df.asvabCS,df.asvabMK,df.asvabNO,df.asvabPC,df.asvabWK)
cor_asvab = cor(matrix_asvab)
println(cor_asvab)

# Question 3
ols_asvab_lm = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr +
                               asvabAR + asvabCS + asvabMK + asvabNO + asvabPC + asvabWK), df)
println(ols_asvab_lm)

#  It will be problematic to directly include  the six asvab variables in the regression
#  since these variables are highly correlated with each other.

# Question 4
asvabMat = (hcat(df.asvabAR, df.asvabCS,df.asvabMK,df.asvabNO,df.asvabPC,df.asvabWK))'
M = fit(PCA, asvabMat; maxoutdim=1)
asvabPCA = MultivariateStats.transform(M, asvabMat)
asvabPCA= vec(asvabPCA')
insert!(df,7, asvabPCA, :asvabPCA)
ols_PCA = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabPCA), df)
println(ols_PCA)

# Question 5
asvabMat = (hcat(df.asvabAR, df.asvabCS,df.asvabMK,df.asvabNO,df.asvabPC,df.asvabWK))'
M = fit(FactorAnalysis, asvabMat; maxoutdim=1)
asvabFactorAnalysis = MultivariateStats.transform(M, asvabMat)
asvabFactorAnalysis= vec(asvabFactorAnalysis')
insert!(df,8, asvabFactorAnalysis, :asvabFactorAnalysis)
ols_FactorAnalysis = lm(@formula(logwage ~ black + hispanic + female + schoolt + gradHS + grad4yr + asvabFactorAnalysis), df)
println(ols_FactorAnalysis)

# Question 6
function lgwt(N0::Integer,a::Real=-1,b::Real=1)
    N  = N0-1
    N1 = N+1
    N2 = N+2
    xu = range(-1,stop=1,length=N1)
    y = cos.((2*(0:N) .+ 1)*pi/(2*N .+ 2))  .+  ( 0.27/N1 ) .* sin.( pi .* xu .* N/N2 )
    L  = zeros(N1,N2)
    Lp = zeros(N1,N2)
    y0 = 2
    vareps = 2e-52
    i = 0
    tracker=0
    it_max=10

    while (norm(y.-y0,Inf)>vareps && tracker<=it_max)
        d=norm(y.-y0,Inf)

        L[:,1]  .= 1
        Lp[:,1] .= 0

        L[:,2] .= y

        for k=2:N1
            L[:,k+1] = ( (2*k-1)*y .* L[:,k] .- (k-1)*L[:,k-1] )/k
        end

        Lp = (N2)*( L[:,N1] .- y .* L[:,N2] )./(1 .- y.^2)
        y0 = y
        y  = y0 - L[:,N2]./Lp
        if norm(y.-y0,Inf)==d
            tracker+=1
        end
        i+=1

    end

    x = (a.*(1 .- y) .+ b .* (1 .+ y))./2
    w=(b-a)./((1 .- y.^2).*Lp.^2)*(N2/N1)^2

    return x,w
end

D= Normal(0,1)
ξ= rand(Normal(0,1),size(df,1),1)
M= hcat(df.asvabAR, df.asvabCS,df.asvabMK,df.asvabNO,df.asvabPC,df.asvabWK)
Xm= hcat(ones(size(df,1)), df.black, df.hispanic, df.female)
X= hcat(ones(size(df,1)), df.black, df.hispanic, df.female, df.schoolt, df.gradHS, df.grad4yr)
Y= df.logwage
L= zeros(size(df,1),1)
function likelihoodfunction(sigma,alpha,beta,gamma,delta)
    sigmaj=sigma[1:end-1]
    sigmaw=sigma[end]
    for i= 1: size(df,1)
        for j= 1:6
            part1vec=zeros(6,1)
            part1vec[j,1]=(1/sigma[1,j])*pdf.(D, (M[i,j]-Xm[i,:]*alpha-gamma*ξ[i,1])/1/sigma[1,j]))
            part1=part1vec[1,1]*part1vec[2,1]*part1vec[3,1]*part1vec[4,1]*part1vec[5,1]*part1vec[6,1]
        end
        part2=(1/sigmaw)*pdf.(D,(Y[i,1]-X[i,:]*beta-delta*ξ[i,1])/sigmaw)
        L[i,1]= part1*part2
    end
    return L
end

nodes, weights = lgwt(7,-4,4)
likelihood=zeros(size(df,1),1)
for i= 1: size(df,1)
likelihood[i,1]=log(sum(weights.* L[i,1].*pdf.(D,nodes)))
end

sigma_hat = optimize(sigma -> -sum(likelihood), rand(1,7), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(sigma_hat.minimizer)

end

allwrap()


text:
using JLD2
using DataFrames
using BenchmarkTools

file_names = ["benchmarks/runs/conv_N1000.jld2", "benchmarks/runs/deconv_N1000.jld2"]

run = jldopen(file_names[1])
results_conv = run["results"]

run = jldopen(file_names[2])
results_deconv = run["results"]

results = DataFrame()
insertcols!(results, 1, :type => results_conv[!, :type], :conv => results_conv[!, :trial], :deconv => results_deconv[!, :trial])

# io = IOContext(stdout, :logbins => true)
# for row in eachrow(results)
#     show(io, MIME("text/plain"), row[2])
# end


text:
function update_scope_conf!(scope_conf, scope::TDS2002B, refresh_cnt, base)
     # Get Channel Volt per Div
     #refresh_cnt==base*1  && (scope_conf.CH1_Volt_div = string(get_vertical_scale(scope, "CH1")))
     #refresh_cnt==base*2 && (scope_conf.CH2_Volt_div = string(get_vertical_scale(scope, "CH2")))
     ##refresh_cnt==base*3 && (scope_conf.CH1_Offset = string(get_ch_position(scope, "CH1")))
     ##refresh_cnt==base*4 && (scope_conf.CH2_Offset = string(get_ch_position(scope, "CH2")))
     #Get Time per Div
     #refresh_cnt==base*3 && (scope_conf.Time_div = string(get_horizontal_scale(scope, "Time")))
     #Get Triger Seting
     #refresh_cnt==base*4 && (scope_conf.Trigger_source = get_trig_data(scope, "Source"))
     #refresh_cnt==base*5 && (scope_conf.Trigger_level = get_trig_data(scope, "Level"))
     #refresh_cnt==base*6 && (scope_conf.Trigger_mode = get_trig_data(scope, "Mode"))
     #Get Meas - not working yet
     #(refresh_cnt==base*3) && (scope_conf.Measurement_Value1 = string(get_meas_data(scope, "Meas_Nr1")))
     #(refresh_cnt==base*4) && (scope_conf.Measurement_Value2 = string(get_meas_data(scope, "Meas_Nr2")))
     #(refresh_cnt==base*5) && (scope_conf.Measurement_Value3 = string(get_meas_data(scope, "Meas_Nr3")))
	 #(refresh_cnt==base*6) && (scope_conf.Measurement_Value4 = string(get_meas_data(scope, "Meas_Nr4")))
     #(refresh_cnt==base*7) && (scope_conf.Measurement_Value5 = string(get_meas_data(scope, "Meas_Nr5")))

	 # conf acq chan1 and acquire
	 (refresh_cnt==base*1) && (conf_acq_ch(scope, "CH1"))
	 (refresh_cnt==base*5) && ((scope_conf.t, scope_conf.y1) = Trigger_Aquistion(scope, "CH1"))

	 # conf acq chan2 and acquire
	 (refresh_cnt==base*10) && (conf_acq_ch(scope, "CH2"))
	 (refresh_cnt==base*15) && ((scope_conf.t, scope_conf.y2) = Trigger_Aquistion(scope, "CH2"))

end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
# Cross-correlation and convolution
function FFTConvolution( image::Array{<:Real,N}, kernel::Array{<:Real,N}; typ=Float32 ) where {N}
	return FFTConvolution!( image, kernel, zeros( typ, size(image).+size(kernel).-1) );
end

function FFTConvolution!( image::Array{<:Real,N}, kernel::Array{<:Real,N}, conv::Array{T,N} ) where {T,N}

	@assert all( size(conv) .== size(image) .+ size(kernel) .- 1 );

	paddedkernel = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedkernel[ 1:size(kernel,1), 1:size(kernel,2), 1:size(kernel,3) ] .= convert.( Complex{T}, kernel )

	paddedimage  = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedimage[ 1:size(image,1), 1:size(image,2), 1:size(image,3) ] .= convert.( Complex{T}, image )

	FFTW.fft!( paddedimage  );
	FFTW.fft!( paddedkernel );

	@inbounds @simd for idx in 1:length(paddedimage)
		paddedimage[idx] = paddedimage[idx] * paddedkernel[idx]
	end

	FFTW.ifft!( paddedimage );

	@inbounds @simd for idx in 1:length(paddedimage)
		conv[idx] = real( paddedimage[idx] )
	end

	return conv
end

function FFTConvolution_crop( image::Array{<:Real,N}, kernel::Array{<:Real,N}; typ=Float32 ) where {N}
	return FFTConvolution_crop!( image, kernel, zeros( typ, size(image) ) );
end

function FFTConvolution_crop!( image::Array{<:Real,N}, kernel::Array{<:Real,N}, conv::Array{T,N} ) where {T,N}

	@assert all( size(conv) .== size(image) );

	paddedkernel = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedkernel[ 1:size(kernel,1), 1:size(kernel,2), 1:size(kernel,3) ] .= convert.( Complex{T}, kernel )

	paddedimage  = zeros( Complex{T}, size(image) .+ size(kernel) .- 1 );
	paddedimage[ 1:size(image,1), 1:size(image,2), 1:size(image,3) ] .= convert.( Complex{T}, image )

	FFTW.fft!( paddedimage  );
	FFTW.fft!( paddedkernel );

	@inbounds @simd for idx in 1:length(paddedimage)
		paddedimage[idx] = paddedimage[idx] * paddedkernel[idx]
	end

	FFTW.ifft!( paddedimage );

	offs = div.( size( kernel ), 2 );
	zoff = ( N == 2 ) ? 0 : offs[3];

	for zet in 1+zoff:size(image,3)+zoff, col in 1+offs[2]:size(image,2)+offs[2]
		@simd for row in 1+offs[1]:size(image,1)+offs[1]
			conv[row-offs[1],col-offs[2],zet-zoff] = real( paddedimage[row,col,zet] )
	end end

	return conv
end

function FFTConvolution_b( image::Array{T1,2}, kernel::Array{T2,2}; typ=Float32 ) where {T1<:Real,T2<:Real}

	paddedkernel = zeros( typ, size(image) .+ size(kernel) .- 1 );
	paddedkernel[ 1:size(kernel,1), 1:size(kernel,2) ] .= convert.( typ, kernel )

	paddedimage  = zeros( typ, size(image) .+ size(kernel) .- 1 );
	paddedimage[ 1:size(image,1), 1:size(image,2) ] .= convert.( typ, image )

	Ff = FFTW.fft( paddedimage  );
	Fg = FFTW.fft( paddedkernel );

	Fcorr =  Ff .* Fg;

	return real.( FFTW.bfft( Fcorr ) )
end

function paddedConvolution( image::Array{T1,2}, kernel::Array{T2,2}; typ=Float32 ) where {T1<:Real,T2<:Real}

	N1, N2 = size( image  );
	M1, M2 = size( kernel );

	paddedImage = similar(image, N1+2*M1-2, N2+2*M2-2);
	paddedImage[ M1:(M1+N1-1), M2:(M2+N2-1) ] .= image;

	corrmatrix = zeros( typ, N1 + M1 - 1, N2 + M2 - 1 )

	for col in 1:( N2 + M2 - 1)
		for row in 1:( N1 + M1 - 1 )
			corrmatrix[ row, col ] = sum( paddedImage[ row:row+M1-1, col:col+M2-1 ] .* kernel )
		end
	end
	return corrmatrix
end

max_stars_repo_name: BrucePBon/ImageAnalysis

max_issues_repo_name: BrucePBon/ImageAnalysis

max_forks_repo_name: BrucePBon/ImageAnalysis


text:
"""
    $(SIGNATURES)

Multiply various full and partial dimension proposal densities.

DevNotes
- FIXME consolidate partial and full product AMP API, relates to #1010
- TODO better consolidate with full dimension product
- TODO -- reuse memory rather than rand here
"""
function prodmultiplefullpartials(dens::Vector{BallTreeDensity},
                                  partials::Dict{Int, Vector{BallTreeDensity}},
                                  Ndims::Int,
                                  N::Int,
                                  manis::Tuple )
  #
  # calculate products over all dimensions, legacy proposals held in `dens` vector
  pGM = AMP.manifoldProduct(dens, manis, Niter=1) |> getPoints

  # whats up with this?
  for (dimnum,pp) in partials
    push!(pp, AMP.manikde!(pGM[dimnum:dimnum,:], (manis[dimnum],) ))
  end

  # do each partial dimension individually
  for (dimnum,pp) in partials
    pGM[dimnum,:] = AMP.manifoldProduct(pp, (manis[dimnum],), Niter=1) |> getPoints
  end

  return pGM
end

"""
    $SIGNATURES

Take product of `dens` accompanied by optional `partials` proposal belief densities.

Notes
-----
- `d` dimensional product approximation
- `partials` are treated as one dimensional
- Incorporate ApproxManifoldProducts to process variables in individual batches.
"""
function productbelief( dfg::AbstractDFG,
                        vertlabel::Symbol,
                        dens::Vector{<:BallTreeDensity},
                        partials::Dict{Int, <:AbstractVector{<:BallTreeDensity}},
                        N::Int;
                        dbg::Bool=false,
                        logger=ConsoleLogger()  )
  #
  vert = DFG.getVariable(dfg, vertlabel)
  manis = getVariableType(vert) |> getManifolds
  pGM = Array{Float64,2}(undef, 0,0)
  lennonp, lenpart = length(dens), length(partials)
  denspts = getPoints(getBelief(vert))
  Ndims = size(denspts,1)
  with_logger(logger) do
    @info "[$(lennonp)x$(lenpart)p,d$(Ndims),N$(N)],"
  end

  dens_ = 0 < lennonp ? dens : [deepcopy(getBelief(vert));]

  if 0 < lennonp # || (lennonp == 0 && 0 < lenpart)
    # multiple non-partials
    pGM = prodmultiplefullpartials(dens_, partials, Ndims, N, manis)
  elseif lennonp == 0 && 0 < lenpart
    # only partials, must get other existing values for vertlabel from dfg
    pGM = deepcopy(denspts)
      # do each partial dimension individually
    for (dimnum,pp) in partials
      pGM[dimnum,:] = AMP.manifoldProduct(pp, (manis[dimnum],), Niter=1) |> getPoints
    end
  else
    with_logger(logger) do
      @warn "Unknown density product on variable=$(vert.label), lennonp=$(lennonp), lenpart=$(lenpart)"
    end
    pGM = Array{Float64,2}(undef, 0,1)
  end

  return pGM
end

"""
    $SIGNATURES

Calculate the proposals and products on `destvert` using `factors` in factor graph `dfg`.

Notes
- Returns tuple of product and whether full dimensional (=true) or partial (=false).
"""
function predictbelief( dfg::AbstractDFG,
                        destvert::DFGVariable,
                        factors::Vector{<:DFGFactor};
                        solveKey::Symbol=:defaul,
                        needFreshMeasurements::Bool=true,
                        N::Int=0,
                        dbg::Bool=false,
                        logger=ConsoleLogger(),
                        dens = Array{BallTreeDensity,1}(),
                        partials = Dict{Int, Vector{BallTreeDensity}}()  )
  #

  # determine number of particles to draw from the marginal
  nn = N != 0 ? N : size(getVal(destvert),2)

  # get proposal beliefs
  destvertlabel = destvert.label
  inferdim = proposalbeliefs!(dfg, destvertlabel, factors, dens, partials, N=nn, dbg=dbg)

  # take the product
  pGM = productbelief(dfg, destvertlabel, dens, partials, nn, dbg=dbg, logger=logger )

  return pGM, sum(inferdim)
end

function predictbelief( dfg::AbstractDFG,
                        destvertsym::Symbol,
                        factorsyms::AbstractVector{Symbol};
                        solveKey::Symbol=:default,
                        needFreshMeasurements::Bool=true,
                        N::Int=0,
                        dbg::Bool=false,
                        logger=ConsoleLogger(),
                        dens = Array{BallTreeDensity,1}(),
                        partials = Dict{Int, Vector{BallTreeDensity}}()  )
  #
  factors = getFactor.(dfg, factorsyms)
  vert = getVariable(dfg, destvertsym)

  # determine the number of particles to draw from the marginal
  nn = N != 0 ? N : size(getVal(vert, solveKey=solveKey),2)

  # do the belief prediction
  predictbelief(dfg, vert, factors, solveKey=solveKey, needFreshMeasurements=needFreshMeasurements, N=nn, dbg=dbg, logger=logger, dens=dens, partials=partials)
end

function predictbelief( dfg::AbstractDFG,
                        destvertsym::Symbol,
                        ::Colon;
                        solveKey::Symbol=:default,
                        needFreshMeasurements::Bool=true,
                        N::Int=0,
                        dbg::Bool=false,
                        logger=ConsoleLogger(),
                        dens = Array{BallTreeDensity,1}(),
                        partials = Dict{Int, Vector{BallTreeDensity}}() )
  #
  predictbelief(dfg, destvertsym, getNeighbors(dfg, destvertsym), solveKey=solveKey, needFreshMeasurements=needFreshMeasurements, N=N, dbg=dbg, logger=logger, dens=dens, partials=partials )
end

"""
    $(SIGNATURES)

Using factor graph object `dfg`, project belief through connected factors
(convolution with likelihood) to variable `sym` followed by a approximate functional product.

Return: product belief, full proposals, partial dimension proposals, labels
"""
function localProduct(dfg::AbstractDFG,
                      sym::Symbol;
                      solveKey::Symbol=:default,
                      N::Int=100,
                      dbg::Bool=false,
                      logger=ConsoleLogger() )
  #
  # vector of all neighbors as Symbols
  lb = getNeighbors(dfg, sym)

  # # get proposal beliefs
  dens = Array{BallTreeDensity,1}()
  partials = Dict{Int, Vector{BallTreeDensity}}()
  pGM, sinfd = predictbelief(dfg, sym, lb, solveKey=solveKey, logger=logger, dens=dens, partials=partials)

  # make manifold belief from product
  vari = getVariable(dfg, sym)
  pp = AMP.manikde!(pGM, getVariableType(vari) |> getManifolds )

  return pp, dens, partials, lb, sinfd
end
localProduct(dfg::AbstractDFG, lbl::AbstractString; solveKey::Symbol=:default, N::Int=100, dbg::Bool=false) = localProduct(dfg, Symbol(lbl), solveKey=solveKey, N=N, dbg=dbg)

"""
    $SIGNATURES

Basic wrapper to take local product and then set the value of `sym` in `dfg`.

Notes
- returns `::Tuple{BallTreeDensity, Float64, Vector{Symbol}}`

DevNotes:
- Unknown issue first occurred here near IIF v0.8.4 tag, recorded case at 2020-01-17T15:26:17.673
"""
function localProductAndUpdate!(dfg::AbstractDFG,
                                sym::Symbol,
                                setkde::Bool=true,
                                logger=ConsoleLogger() )
  #
  # calculate new points for sym using existing structure around sym in dfg
  newPts, dens, parts, lbl, infdim = localProduct(dfg, sym, N=getSolverParams(dfg).N, logger=logger)
  # maybe update dfg sym with newly calculated points
  setkde && 0 < size(getPoints(newPts),2) ? setValKDE!(dfg, sym, newPts, false, infdim) : nothing

  return newPts, infdim, lbl
end

#


text:
function kvadratnajna(a,b,c)
    if (b*b - 4*a*c) == 0
        return -b / 2 * a , -b / 2 * a
    end
    return (-b - sqrt(b*b - 4*a*c)) / 2 * a, (-b + sqrt(b*b - 4*a*c)) / 2 * a

end


text:
include("solve_problem.jl")

# Path of the file containing the dataset
file = "dataset/multi-depot/stations/ceuze-1.txt"
# Number of vehicules available
nb_vehicles = 3
# Wether using BnP-opt or BnP-approx in the case of a multi-depot dataset
approx = true
# Choose to display or not Branch-and-Price steps
verbose = true

solve_problem(file, nb_vehicles, approx, verbose)

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
export DoubleSymLayer,getDoubleSymLayer

"""
 Implementation of symmetric double layer model

 Y(theta,Y0) = K(th1)'(activation( K(th1)\*Y0 + trafo.Bin\*th2))) + trafo.Bout\*th3
"""
mutable struct DoubleSymLayer{T, TK <: AbstractConvKernel{T}, TN <: Union{batchNormNN{T}, normLayer{T}}} <: AbstractMeganetElement{T}
    activation  :: Function   # activation function
    activation! :: Function # in place activation function
    K           :: TK   # Kernel model, e.g., convMod
    nLayer      :: TN   # normalization layer
    Bin         :: Array{T}   # Bias inside the nonlinearity
    Bout        :: Array{T}   # bias outside the nonlinearity
end

function getDoubleSymLayer(TYPE::Type,K,nLayer::AbstractMeganetElement{T};
                           Bin=zeros(nFeatOut(K),0),Bout=zeros(nFeatIn(K),0),
                           activation=tanhActivation,activation_inplace=tanhActivation!) where {T <: Number}
    BinT = convert.(T, Bin)
    BoutT = convert.(T, Bout)
    return DoubleSymLayer(activation,activation_inplace,K,nLayer,BinT,BoutT);

end

function splitWeights(this::DoubleSymLayer{T},theta::Array{T}) where {T<:Number}

    th1 = theta[1:nTheta(this.K)::Int]
    cnt = length(th1)
    th2 = theta[cnt+(1:size(this.Bin,2))]
    cnt = cnt + length(th2)
    th3 = theta[cnt+(1:size(this.Bout,2))]
    cnt = cnt + length(th3)

    th4 = theta[cnt+1:end];

    return th1, th2, th3, th4
end

function apply(this::DoubleSymLayer{T},theta::Array{T},Yin::Array{T,2},tmp,doDerivative=true) where {T<:Number}
    if isempty(tmp)
        tmp = Array{Any}(2)
        tmp[1] = Array{Any}(0)
        tmp[2] = Array{Any}(0)
    end
    #QZ = []
    nex = div(length(Yin),nFeatIn(this))::Int
    Y   = reshape(Yin,:,nex)

    theta1,theta2,theta3,theta4 = splitWeights(this,theta)
    Kop    = getOp(this.K,theta1)
    KY     = Kop*Y # TODO: Look into making convolution in place

    KY,dummy,tmp[1] = apply(this.nLayer,theta4,KY,tmp[1],doDerivative)
    Yt     = KY
    if !isempty(theta2)
     Yt .+= this.Bin*theta2
    end

    if doDerivative
        if isempty(tmp[2])
            tmp[2] = copy(Yt)
        else
            tmp2 = tmp[2]
            tmp2 .= Yt
        end
    end

    Z::Array{T,2},      = this.activation!(Yt,[],false)
    Z      = -(Kop'*Z)
    if !isempty(theta3)
        Z  .+= this.Bout*theta3
    end
    return Z, Z, tmp
end

function nTheta(this::DoubleSymLayer)
    return nTheta(this.K) + size(this.Bin,2)+ size(this.Bout,2) + nTheta(this.nLayer)
end

function nFeatIn(this::DoubleSymLayer)
    return nFeatIn(this.K)
end

function nFeatOut(this::DoubleSymLayer)
    return nFeatIn(this.K)
end

function nDataOut(this::DoubleSymLayer)
    return nFeatIn(this)
end

function initTheta(this::DoubleSymLayer{T})  where {T<:Number}
    theta = [vec(initTheta(this.K));
    T(0.01)*ones(T,size(this.Bin,2),1);
    T(0.01)*ones(T,size(this.Bout,2),1);
    initTheta(this.nLayer)];
    return theta
end

function Jthetamv(this::DoubleSymLayer{T},dtheta::Array{T},theta::Array{T},Y::Array{T},tmp)  where {T<:Number}

    A,dA = this.activation(tmp[2],true)
    th1, th2,th3,th4    = splitWeights(this,theta)
    dth1,dth2,dth3,dth4 = splitWeights(this,dtheta)

	Kop    = getOp(this.K,th1)
    dKop   = getOp(this.K,dth1)
    dY     = dKop*Y

    dY = Jmv(this.nLayer,dth4,dY,th4,Kop*Y,copy(tmp[1]))[2]
    dY = dY .+ this.Bin*dth2

    dY = -(Kop'*(dA.*dY) + dKop'*A) .+ this.Bout*dth3
    return dY, dY
end

function JYmv(this::DoubleSymLayer{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp)  where {T<:Number}
    #TODO: Look into why this activation cannot be done in place (tests fail)
    dA = this.activation(tmp[2],true)[2]

    nex = div(length(dY),nFeatIn(this))
    dY  = reshape(dY,:,nex)
    Y   = reshape(Y,:,nex)
    th1, th2,th3,th4  = splitWeights(this,theta)

    Kop = getOp(this.K,th1)
    dY = Kop*dY
    dY = JYmv(this.nLayer,dY,th4,Kop*Y,copy(tmp[1]))[2]
    dZ = -(Kop'*(dA.*dY))
    return dZ, dZ
end

function Jmv(this::DoubleSymLayer{T},dtheta::Array{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp)  where {T<:Number}
    A,dA = this.activation(copy(tmp[2]),true)
    nex = div(length(Y),nFeatIn(this))

    th1, th2,th3,th4    = splitWeights(this,theta)
    dth1,dth2,dth3,dth4 = splitWeights(this,dtheta)

    Kop    = getOp(this.K,th1)
    dKop   = getOp(this.K,dth1)
    if length(dY)>1
        dY  = reshape(dY,:,nex)
        KdY = Kop*dY
    else
        KdY = 0
    end
    dY = dKop*Y+KdY
    dY = Jmv(this.nLayer,dth4,dY,th4,Kop*Y,tmp[1])[2]

    dY = reshape(dY,:,nex)
    if !isempty(dth2)
        dY .+= this.Bin*dth2
    end

    dY = -(Kop'*(dA.*dY) + dKop'*A)
    if !isempty(dth3)
        dth3 .+= this.Bout*dth3
    end

    return dY, dY
end

function JthetaTmv(this::DoubleSymLayer{T},Z::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp)  where {T<:Number}

    nex       = div(length(Y),nFeatIn(this))
    Z         = reshape(Z,:,nex)
    th1,th2,th3,th4  = splitWeights(this,theta)
    Kop       = getOp(this.K,th1)
    A,dA      = this.activation(tmp[2],true)

    dth3      = vec(sum(this.Bout'*Z,2))
    dAZ       = dA.*(Kop*Z)
    dth2      = vec(sum(this.Bin'*dAZ,2))

    dth4,dAZ  = JTmv(this.nLayer,dAZ,zeros(T,0),th4,Kop*Y,tmp[1])

    dth1      = JthetaTmv(this.K,A,zeros(T,0),Z)
    dth1     += JthetaTmv(this.K,dAZ,zeros(T,0),Y)
    dtheta    = [-vec(dth1); -vec(dth2); vec(dth3); -vec(dth4)]
    return dtheta
end

function JYTmv(this::DoubleSymLayer{T},Zin::Array{T},dummy::Array{T},theta::Array{T},Y::Array{T},tmp)  where {T<:Number}

    nex       = div(length(Y),nFeatIn(this))
    Z         = reshape(Zin,:,nex)
    th1,th2,th3,th4  = splitWeights(this,theta)
    Kop       = getOp(this.K,th1)
    A,dA      = this.activation(tmp[2],true)

    dAZ       = dA.*(Kop*Z)
    dAZ       = JYTmv(this.nLayer,dAZ,(T)[],th4,Kop*Y,tmp[1])
    dAZ_out       = reshape(dAZ,:,nex)
    dY  = -(Kop'*dAZ_out)
    return dY
end

function JTmv(this::DoubleSymLayer{T}, Zin::Array{T}, dummy::Array{T},
                                            theta::Array{T}, Yin::Array{T}, tmp) where {T<:Number}

    nex  = div(length(Yin),nFeatIn(this))
    Z         = reshape(Zin, :, nex)
    Yt        = reshape(tmp[2]::Array{T,2},:,nex)
    Y         = reshape(Yin,:,nex)
    th1, th2, th3, th4  = splitWeights(this,theta)
    #Kop       = getOp(this.K,th1)
    A::Array{T,2}, dA::Array{T,2}    = this.activation(Yt,true)

    dth3      = vec(sum(this.Bout'*Z,2))

    KopZ = Amv(this.K, th1, Z)
    dAZ1       = dA.*KopZ

    dth2      = vec(sum(this.Bin'*dAZ1,2))
    KopY      = Amv(this.K, th1, Y)
    dth4, dAZ2  = JTmv(this.nLayer,dAZ1,zeros(T,0),th4,KopY,tmp[1])
    dth1      = JthetaTmv(this.K,dAZ2,zeros(T,0),Y)
    dth1 += JthetaTmv(this.K,A,(T)[],Z)
    dtheta    = [-vec(dth1); -vec(dth2); vec(dth3);-vec(dth4)]

    dAZ_out = reshape(dAZ2,:,nex)
    KopTdAZ = ATmv(this.K, th1, dAZ_out)
    dY  = -KopTdAZ
    return dtheta, dY
end


text:
"""
    AsTable(cols)

A type used for selection operations to signal that the columns selected by the
wrapped selector should be passed as a `NamedTuple` to the function.
"""
struct AsTable
    cols
end

Base.broadcastable(x::AsTable) = Ref(x)

function make_unique!(names::Vector{Symbol}, src::AbstractVector{Symbol};
                      makeunique::Bool=false)
    if length(names) != length(src)
        throw(DimensionMismatch("Length of src doesn't match length of names."))
    end
    seen = Set{Symbol}()
    dups = Int[]
    for i in 1:length(names)
        name = src[i]
        if in(name, seen)
            push!(dups, i)
        else
            names[i] = src[i]
            push!(seen, name)
        end
    end

    if length(dups) > 0
        if !makeunique
            dupstr = join(string.(':', unique(src[dups])), ", ", " and ")
            msg = "Duplicate variable names: $dupstr. Pass makeunique=true " *
                  "to make them unique using a suffix automatically."
            throw(ArgumentError(msg))
        end
    end

    for i in dups
        nm = src[i]
        k = 1
        while true
            newnm = Symbol("$(nm)_$k")
            if !in(newnm, seen)
                names[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end

    return names
end

function make_unique(names::AbstractVector{Symbol}; makeunique::Bool=false)
    make_unique!(similar(names), names, makeunique=makeunique)
end

"""
    gennames(n::Integer)

Generate standardized names for columns of a DataFrame.
The first name will be `:x1`, the second `:x2`, etc.
"""
function gennames(n::Integer)
    res = Vector{Symbol}(undef, n)
    for i in 1:n
        res[i] = Symbol(@sprintf "x%d" i)
    end
    return res
end

function funname(f)
    n = nameof(f)
    String(n)[1] == '#' ? :function : n
end

if isdefined(Base, :ComposedFunction) # Julia >= 1.6.0-DEV.85
    using Base: ComposedFunction
else
    using Compat: ComposedFunction
end

funname(c::ComposedFunction) = Symbol(funname(c.outer), :_, funname(c.inner))

# Compute chunks of indices, each with at least `basesize` entries
# This method ensures balanced sizes by avoiding a small last chunk
function split_indices(len::Integer, basesize::Integer)
    len′ = Int64(len) # Avoid overflow on 32-bit machines
    @assert len′ > 0
    @assert basesize > 0
    np = Int64(max(1, len ÷ basesize))
    return split_to_chunks(len′, np)
end

function split_to_chunks(len::Integer, np::Integer)
    len′ = Int64(len) # Avoid overflow on 32-bit machines
    np′ = Int64(np)
    @assert len′ > 0
    @assert 0 < np′ <= len′
    return (Int(1 + ((i - 1) * len′) ÷ np):Int((i * len′) ÷ np) for i in 1:np)
end

if VERSION >= v"1.4"
    function _spawn_for_chunks_helper(iter, lbody, basesize)
        lidx = iter.args[1]
        range = iter.args[2]
        quote
            let x = $(esc(range)), basesize = $(esc(basesize))
                @assert firstindex(x) == 1

                nt = Threads.nthreads()
                len = length(x)
                if nt > 1 && len > basesize
                    tasks = [Threads.@spawn begin
                                 for i in p
                                     local $(esc(lidx)) = @inbounds x[i]
                                     $(esc(lbody))
                                 end
                             end
                             for p in split_indices(len, basesize)]
                    foreach(wait, tasks)
                else
                    for i in eachindex(x)
                        local $(esc(lidx)) = @inbounds x[i]
                        $(esc(lbody))
                    end
                end
            end
            nothing
        end
    end
else
    function _spawn_for_chunks_helper(iter, lbody, basesize)
        lidx = iter.args[1]
        range = iter.args[2]
        quote
            let x = $(esc(range))
                for i in eachindex(x)
                    local $(esc(lidx)) = @inbounds x[i]
                    $(esc(lbody))
                end
            end
            nothing
        end
    end
end

"""
    @spawn_for_chunks basesize for i in range ... end

Parallelize a `for` loop by spawning separate tasks
iterating each over a chunk of at least `basesize` elements
in `range`.

A number of task higher than `Threads.nthreads()` may be spawned,
since that can allow for a more efficient load balancing in case
some threads are busy (nested parallelism).
"""
macro spawn_for_chunks(basesize, ex)
    if !(isa(ex, Expr) && ex.head === :for)
        throw(ArgumentError("@spawn_for_chunks requires a `for` loop expression"))
    end
    if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @spawn_for_chunks"))
    end
    return _spawn_for_chunks_helper(ex.args[1], ex.args[2], basesize)
end

function _nt_like_hash(v, h::UInt)
    length(v) == 0 && return hash(NamedTuple(), h)

    h = hash((), h)
    for i in length(v):-1:1
        h = hash(v[i], h)
    end

    return xor(objectid(Tuple(propertynames(v))), h)
end

_findall(B) = findall(B)

@inline _blsr(x) = x & (x-1)

# findall returning a range when possible (all true indices are contiguous), and optimized for B::BitVector
# the main idea is taken from Base.findall(B::BitArray)
function _findall(B::BitVector)::Union{UnitRange{Int}, Vector{Int}}
    nnzB = count(B)
    nnzB == 0 && return 1:0
    nnzB == length(B) && return 1:length(B)
    local I
    Bc = B.chunks
    Bi = 1 # block index
    i1 = 1 # index of current block beginng in B
    i = 1  # index of the _next_ one in I
    c = Bc[1] # current block

    start = -1 # the begining of ones block
    stop = -1  # the end of ones block

    @inbounds while true # I not materialized
        if i > nnzB # all ones in B found
            Ir = start:start + i - 2
            @assert length(Ir) == nnzB
            return Ir
        end

        if c == 0
            if start != -1 && stop == -1
                stop = i1 - 1
            end
            while c == 0 # no need to return here as we returned above
                i1 += 64
                Bi += 1
                c = Bc[Bi]
            end
        end
        if c == ~UInt64(0)
            if stop != -1
                I = Vector{Int}(undef, nnzB)
                for j in 1:i-1
                    I[j] = start + j - 1
                end
                break
            end
            if start == -1
                start = i1
            end
            while c == ~UInt64(0)
                if Bi == length(Bc)
                    Ir = start:length(B)
                    @assert length(Ir) == nnzB
                    return Ir
                end

                i += 64
                i1 += 64
                Bi += 1
                c = Bc[Bi]
            end
        end
        if c != 0 # mixed ones and zeros in block
            tz = trailing_zeros(c)
            lz = leading_zeros(c)
            co = c >> tz == (one(UInt64) << (64 - lz - tz)) - 1 # block of countinous ones in c
            if stop != -1  # already found block of ones and zeros, just not materialized
                I = Vector{Int}(undef, nnzB)
                for j in 1:i-1
                    I[j] = start + j - 1
                end
                break
            elseif !co # not countinous ones
                I = Vector{Int}(undef, nnzB)
                if start != -1
                    for j in 1:i-1
                        I[j] = start + j - 1
                    end
                end
                break
            else # countinous block of ones
                if start != -1
                    if tz > 0 # like __1111__ or 111111__
                        I = Vector{Int}(undef, nnzB)
                        for j in 1:i-1
                            I[j] = start + j - 1
                        end
                        break
                    else # lz > 0, like __111111
                        stop = i1 + (64 - lz) - 1
                        i += 64 - lz

                        # return if last block
                        if Bi == length(Bc)
                            Ir = start:stop
                            @assert length(Ir) == nnzB
                            return Ir
                        end

                        i1 += 64
                        Bi += 1
                        c = Bc[Bi]
                    end
                else # start == -1
                    start = i1 + tz

                    if lz > 0 # like __111111 or like __1111__
                        stop = i1 + (64 - lz) - 1
                        i += stop - start + 1
                    else # like 111111__
                        i += 64 - tz
                    end

                    # return if last block
                    if Bi == length(Bc)
                        Ir = start:start + i - 2
                        @assert length(Ir) == nnzB
                        return Ir
                    end

                    i1 += 64
                    Bi += 1
                    c = Bc[Bi]
                end
            end
        end
    end
    @inbounds while true # I materialized, process like in Base.findall
        if i > nnzB # all ones in B found
            @assert nnzB == i - 1
            return I
        end

        while c == 0 # no need to return here as we returned above
            i1 += 64
            Bi += 1
            c = Bc[Bi]
        end

        while c == ~UInt64(0)
            for j in 0:64-1
                I[i + j] = i1 + j
            end
            i += 64
            if Bi == length(Bc)
                @assert nnzB == i - 1
                return I
            end
            i1 += 64
            Bi += 1
            c = Bc[Bi]
        end

        while c != 0
            tz = trailing_zeros(c)
            c = _blsr(c) # zeros last nonzero bit
            I[i] = i1 + tz
            i += 1
        end
    end
    @assert false "should not be reached"
end


text:
module Parser #JSON

using Compat

#Define ordered dictionary from DataStructures if present
_HAVE_DATASTRUCTURES = try
    using DataStructures
    import DataStructures.OrderedDict
    true
catch
    false
end

const TYPES = Any # Union(Dict, Array, AbstractString, Number, Bool, Nothing) # Types it may encounter
const KEY_TYPES = Union(AbstractString) # Types it may encounter as object keys

export parse

type ParserState{T<:AbstractString}
    str::T
    s::Int
    e::Int
end
ParserState(str::AbstractString,s::Int,e::Int) = ParserState(str, s, e)

charat{T<:AbstractString}(ps::ParserState{T}) = ps.str[ps.s]
incr(ps::ParserState) = (ps.s += 1)
hasmore(ps::ParserState) = (ps.s < ps.e)

# UTILITIES

# Eat up spaces starting at s.
function chomp_space{T<:AbstractString}(ps::ParserState{T})
    c = charat(ps)
    while isspace(c) && hasmore(ps)
        incr(ps)
        c = charat(ps)
    end
end

# Run past the separator in a name : value pair
function skip_separator{T<:AbstractString}(ps::ParserState{T})
    while (charat(ps) != ':') && hasmore(ps)
        incr(ps)
    end
    (charat(ps) != ':') && error("Separator not found ", ps)
    incr(ps)
    nothing
end

# Used for line counts
function _count_before{T<:AbstractString}(haystack::T, needle::Char, _end::Int)
    count = 0
    i = 1
    while i < _end
        haystack[i]==needle && (count += 1)
        i += 1
    end
    count
end

# Prints an error message with an indicator to the source
function _error(message::AbstractString, ps::ParserState)
    lines = _count_before(ps.str, '\n', ps.s)
    # Replace all special multi-line/multi-space characters with a space.
    strnl = replace(ps.str, r"[\b\f\n\r\t\s]", " ")
    li = (ps.s > 20) ? ps.s - 9 : 1 # Left index
    ri = min(ps.e, ps.s + 20)       # Right index
    error(message *
      "\nLine: " * string(lines) *
      "\nAround: ..." * strnl[li:ri] * "..." *
      "\n           " * (" " ^ (ps.s - li)) * "^\n"
    )
end

# PARSING

function parse_array{T<:AbstractString}(ps::ParserState{T}, ordered::Bool)
    incr(ps) # Skip over the '['
    _array = TYPES[]
    chomp_space(ps)
    charat(ps)==']' && (incr(ps); return _array) # Check for empty array
    while true # Extract values from array
        v = parse_value(ps, ordered) # Extract value
        push!(_array, v)
        # Eat up trailing whitespace
        chomp_space(ps)
        c = charat(ps)
        if c == ','
            incr(ps)
            continue
        elseif c == ']'
            incr(ps)
            break
        else
            _error("Unexpected char: " * string(c), ps)
        end
    end
    return _array
end

function parse_object{T<:AbstractString}(ps::ParserState{T}, ordered::Bool)
    if ordered
        parse_object(ps, ordered, OrderedDict{KEY_TYPES,TYPES}())
    else
        parse_object(ps, ordered, Dict{KEY_TYPES,TYPES}())
    end
end

function parse_object{T<:AbstractString}(ps::ParserState{T}, ordered::Bool, obj)
    incr(ps) # Skip over opening '{'
    chomp_space(ps)
    charat(ps)=='}' && (incr(ps); return obj) # Check for empty object
    while true
        chomp_space(ps)
        _key = parse_string(ps)           # Key
        skip_separator(ps)
        _value = parse_value(ps, ordered) # Value
        obj[_key] = _value                             # Building object
        chomp_space(ps)
        c = charat(ps) # Find the next pair or end of object
        if c == ','
            incr(ps)
            continue
        elseif c == '}'
            incr(ps)
            break
        else
            _error("Unexpected char: " * string(c), ps)
        end
    end
    return obj
end

if VERSION <= v"0.3-"
    utf16_is_surrogate(c::Uint16) = (c & 0xf800) == 0xd800
    utf16_get_supplementary(lead::Uint16, trail::Uint16) = char((lead-0xd7f7)<<10 + trail)
else
    const utf16_is_surrogate = Base.utf16_is_surrogate
    const utf16_get_supplementary = Base.utf16_get_supplementary
end

# TODO: Try to find ways to improve the performance of this (currently one
#       of the slowest parsing methods).
function parse_string{T<:AbstractString}(ps::ParserState{T})
    str = ps.str
    s = ps.s
    e = ps.e

    str[s]=='"' || _error("Missing opening string char", ps)
    s = nextind(str, s) # Skip over opening '"'
    b = IOBuffer()
    found_end = false
    while s <= e
        c = str[s]
        if c == '\\'
            s = nextind(str, s)
            c = str[s]
            if c == 'u' # Unicode escape
                u = unescape_string(str[s - 1:s + 4]) # Get the string
                c = u[1]
                if utf16_is_surrogate(@compat(UInt16(c)))
                    if str[s+5] != '\\' || str[s+6] != 'u'
                        _error("Unmatched UTF16 surrogate", ps)
                    end
                    u2 = unescape_string(str[s + 5:s + 10])
                    c = utf16_get_supplementary(@compat(UInt16(c)),@compat(UInt16(u2[1])))
                    # Skip the additional 6 characters
                    for _ = 1:6
                        s = nextind(str, s)
                    end
                end
                write(b, c)
                # Skip over those next four characters
                for _ = 1:4
                    s = nextind(str, s)
                end
            elseif c == '"'  write(b, '"' )
            elseif c == '\\' write(b, '\\')
            elseif c == '/'  write(b, '/' )
            elseif c == 'b'  write(b, '\b')
            elseif c == 'f'  write(b, '\f')
            elseif c == 'n'  write(b, '\n')
            elseif c == 'r'  write(b, '\r')
            elseif c == 't'  write(b, '\t')
            else _error("Unrecognized escaped character: " * string(c), ps)
            end
        elseif c == '"'
            found_end = true
            s = nextind(str, s)
            break
        else
            write(b, c)
        end
        s = nextind(str, s)
    end
    ps.s = s
    found_end || _error("Unterminated string", ps)
    takebuf_string(b)
end

function parse_simple{T<:AbstractString}(ps::ParserState{T})
    c = charat(ps)
    if c == 't' && ps.str[ps.s + 3] == 'e'     # Looks like "true"
        ps.s += 4
        ret = true
    elseif c == 'f' && ps.str[ps.s + 4] == 'e' # Looks like "false"
        ps.s += 5
        ret = false
    elseif c == 'n' && ps.str[ps.s + 3] == 'l' # Looks like "null"
        ps.s += 4
        ret = nothing
    else
        _error("Unknown simple: " * string(c), ps)
    end
    ret
end

function parse_value{T<:AbstractString}(ps::ParserState{T}, ordered::Bool)
    chomp_space(ps)
    (ps.s > ps.e) && return nothing # Nothing left

    ch = charat(ps)
    if ch == '"'
        ret = parse_string(ps)
    elseif ch == '{'
        ret = parse_object(ps, ordered)
    elseif (ch >= '0' && ch <= '9') || ch=='-' || ch=='+'
        ret = parse_number(ps)
    elseif ch == '['
        ret = parse_array(ps, ordered)
    elseif ch == 'f' || ch == 't' || ch == 'n'
        ret = parse_simple(ps)
    else
        _error("Unknown value", ps)
    end
    return ret
end

if VERSION < v"0.4.0-dev+3874"
    function tryparse(T::Type{Float64},str::AbstractString)
        temp = Array(Float64,1)
        if float64_isvalid(str,temp)
            return Nullable(temp[1])
        else
            return Nullable{Float64}()
        end
    end
end

function parse_number{T<:AbstractString}(ps::ParserState{T})
    str = ps.str
    p = ps.s
    e = ps.e
    is_float = false

    c = str[p]
    if c=='-' || c=='+' # Look for sign
        p += 1
        (p <= e) ? (c = str[p]) : _error("Unrecognized number", ps)  # Something must follow a sign
    end

    if c == '0' # If number begins with 0, it must be int(0) or a floating point
        p += 1
        if p <= e
            if str[p] == '.'
                is_float = true
                p += 1
            end
        end
    elseif '0' < c <= '9' # Match more digits
        while '0' <= c <= '9'
            p += 1
            (p <= e) ? (c = str[p]) : break
        end
        if (p <= e) && (c == '.')
            is_float = true
            p += 1
        end
    else
        _error("Unrecognized number", ps)
    end

    if p <= e
        c = str[p]

        if is_float # Match digits after decimal
            while '0' <= c <= '9'
                p += 1
                (p <= e) ? (c = str[p]) : break
            end
        end

        if c == 'E' || c == 'e' || c == 'f' || c == 'F'
            is_float = true
            p += 1
            (p > e) && _error("Unrecognized number", ps)
            c = str[p]
            if c == '-' || c == '+' # Exponent sign
                p += 1
                (p > e) && _error("Unrecognized number", ps)
                c = str[p]
            end
            while '0' <= c <= '9' # Exponent digits
                p += 1
                (p <= e) ? (c = str[p]) : break
            end
        end
    end

    vs = SubString(ps.str, ps.s, p-1)
    ps.s = p
    if is_float
        tmp64 = tryparse(Float64,vs)
        isnull(tmp64) ? error("Invalid floating point number", ps) : return get(tmp64)
    else
        return Base.parse(Int, vs)
    end
end

function parse(str::AbstractString; ordered::Bool=false)
    pos::Int = 1
    len::Int = endof(str)
    len < 1 && return
    ordered && !_HAVE_DATASTRUCTURES && error("DataStructures package required for ordered parsing: try `Pkg.add(\"DataStructures\")`")
    parse_value(ParserState(str, pos, len), ordered)
end

end #module Parser


text:
module Ensembles

export Ensemble, workheatlist, erasure_rate, loadparas, save_ensemble, load_ensemble

import YAML
import Base.+
import JLD2

using ..Force
using ..UpdateMethods
using ..Configurations
using ..Trajectories: lastx

struct Ensemble{T <: AbstractForce}
    nsample :: Integer
    force :: T
    heun :: Heun
    initlist :: Vector{<:AbstractFloat}
    configurations :: Vector{Configuration}
end

function +(e1::Ensemble{T}, e2::Ensemble{T}) where T
    @assert e1.force == e2.force
    @assert e1.heun == e2.heun
    Ensemble(
        e1.nsample + e2.nsample,
        e1.force,
        e1.heun,
        [e1.initlist; e2.initlist],
        [e1.configurations; e2.configurations]
    )
end

function load_ens(files::Vector{String})::Ensemble
    sum(load_configuration(f) for f in files)
end

function workheatlist(ens::Ensemble)
    worklist = [conf.workheat.work for conf in ens.configurations]
    heatlist = [conf.workheat.heat for conf in ens.configurations]
    return worklist, heatlist
end

function erasure_rate end

function erasure_rate(ens::Ensemble)::Real
    endlist = [lastx(conf.traj) for conf in ens.configurations]
    erased_num = (x->erased(ens.force, x)).(endlist) |> sum
    erased_num / ens.nsample
end

function loadparas(filename="para.yml")
    paras = nothing
    open(filename) do file
        paras = YAML.load(file)
    end
    paras
end

function save_ensemble(ens::Ensemble, filename)
    @JLD2.save filename ens
end

function load_ensemble(filename)::Ensemble
    @JLD2.load filename ens
    return ens
end

end


text:
module CodecZstd

export
    ZstdCompressor,
    ZstdCompressorStream,
    ZstdDecompressor,
    ZstdDecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs

using Libdl
const libzpath = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(libzpath)
    error("CodecZlib.jl is not installed properly, run Pkg.build(\"CodecZlib\") and restart Julia.")
end
include(libzpath)
check_deps()

include("libzstd.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
@testset "Reproducing Naveau Results" begin
    @testset "Model (i)" begin

        data = Extremes.dataset("rain")

        # Fit of the first model by maximum likelihood
        fm = EGPpowerfit(data.Rainfall[data.Rainfall .> 0])

        # Parameter estimates
        θ̂ = [params(fm)[1]; params(fm)[2]; params(fm)[3]]

        # Parameter estimates with Naveau (R)
        θ = [4.56372682168789; 0.223166954694046; 1.19212081334056]

        @test θ̂ ≈ θ rtol = 0.01
    end

    @testset "Model (iii)" begin

        data = Extremes.dataset("rain")

        # Fit of the first model by maximum likelihood
        fm = EGPbetafit(data.Rainfall[data.Rainfall .> 0], initialvalues = [log(5), 0.15, log(10)])

        # Parameter estimates
        θ̂ = [params(fm)[1]; params(fm)[2]; params(fm)[3]]

        # Parameter estimates with Naveau (R)
        θ = [5.30852235690306; 0.167910382987591; 29.9938206633414]

        @test θ̂ ≈ θ rtol = 0.01
    end

    @testset "Model (iv)" begin

        data = Extremes.dataset("rain")

        # Fit of the first model by maximum likelihood
        fm = EGPbetapowerfit(data.Rainfall[data.Rainfall .> 0], initialvalues = [log(5), 0.15, log(1.2), log(2)])

        # Parameter estimates
        θ̂ = [params(fm)[1]; params(fm)[2]; params(fm)[3]; params(fm)[4]]

        # Parameter estimates with Naveau (R)
        θ = [3.15761254; 0.21744742; 0.00861047; 1.20754036]

        @test θ̂ ≈ θ rtol = 0.01
    end
end


text:
module NamedTupleUtilities

"""
    select(a::NamedTuple, v::Val{n})
Select a field `n` from `a` if it is in `a`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.select((a=1,b=2,c=3),Val(:a))
(a = 1,)
```
"""
@generated function select(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if i == bn)...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    remove(a::NamedTuple, v::Val{n})
Remove a field `n` from the `a` if it is in `a`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.remove((a=1,b=2,c=3),Val(:c))
(a = 1, b = 2)
```
"""
@generated function remove(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if i != bn)...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    range(a::NamedTuple, b::Val{n}, c::Val{n})
Return a NamedTuple which retains the fields from `b` to `c` in `a`.
If `b` is not in `a`, then it will return the empty NamedTuple.
If `c` is not in `a`, then it will return everything starting with `b`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.range((a=1,b=2,c=3),Val(:a),Val(:b))
(a = 1, b = 2)
```
"""
@generated function range(a::NamedTuple{an}, ::Val{bn}, ::Val{cn}) where {an, bn, cn}
    rangeStarted = false
    names = Symbol[]
    for n in an
        if n == bn
            rangeStarted = true
        end
        if rangeStarted
            push!(names, n)
        end
        if n == cn
            rangeStarted = false
            break
        end
    end
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :( NamedTuple{$(names...,),$types}(($(vals...),)) )
end

"""
    rename(a::NamedTuple, b::Val{n}, c::Val{n})
Return a NamedTuple derived from `a` in which the the field from `b` is renamed to `c`.
If `b` is not in `a`, then it will return the original NamedTuple.
If `c` is in `a`, then `ERROR: duplicate field name in NamedTuple: "c" is not unique` will occur.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:a),Val(:d))
(d = 1, b = 2, c = 3)
julia> QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:m),Val(:d))
(a = 1, b = 2, c = 3)
julia> QueryOperators.NamedTupleUtilities.rename((a = 1, b = 2, c = 3),Val(:a),Val(:c))
ERROR: duplicate field name in NamedTuple: "c" is not unique
```
"""
@generated function rename(a::NamedTuple{an}, ::Val{bn}, ::Val{cn}) where {an, bn, cn}
    names = Symbol[]
    typesArray = DataType[]
    vals = Expr[]
    for n in an
        if n == bn
            push!(names, cn)
        else
            push!(names, n)
        end
        push!(typesArray, fieldtype(a, n))
        push!(vals, :(getfield(a, $(QuoteNode(n)))))
    end
    types = Tuple{typesArray...}
    return :(NamedTuple{$(names...,),$types}(($(vals...),)))
end

"""
    startswith(a::NamedTuple, b::Val{n})
Return a NamedTuple which retains the fields with names starting with `b` in `a`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.startswith((abc=1,bcd=2,cde=3),Val(:a))
(abc = 1,)
```
"""
@generated function startswith(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if Base.startswith(String(i), String(bn)))...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    not_startswith(a::NamedTuple, b::Val{n})
Return a NamedTuple which retains the fields with names that do not start with `b` in `a`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.not_startswith((abc=1,bcd=2,cde=3),Val(:a))
(bcd = 2, cde = 3)
```
"""
@generated function not_startswith(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if !Base.startswith(String(i), String(bn)))...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    endswith(a::NamedTuple, b::Val{n})
Return a NamedTuple which retains the fields with names that do not end with `b` in `a`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.not_endswith((abc=1,bcd=2,cde=3),Val(:d))
(abc = 1, cde = 3)
```
"""
@generated function endswith(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if Base.endswith(String(i), String(bn)))...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    not_endswith(a::NamedTuple, b::Val{n})
Return a NamedTuple which retains the fields with names ending with `b` in `a`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.endswith((abc=1,bcd=2,cde=3),Val(:d))
(bcd = 2,)
```
"""
@generated function not_endswith(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if !Base.endswith(String(i), String(bn)))...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    occursin(a::NamedTuple, b::Val{n})
Return a NamedTuple which retains the fields with names containing `b` as a substring.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.occursin((abc=1,bcd=2,cde=3),Val(:d))
(bcd = 2, cde = 3)
```
"""
@generated function occursin(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if Base.occursin(String(bn), String(i)))...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    not_occursin(a::NamedTuple, b::Val{n})
Return a NamedTuple which retains the fields without names containing `b` as a substring.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.not_occursin((abc=1,bcd=2,cde=3),Val(:d))
(abc = 1,)
```
"""
@generated function not_occursin(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if !Base.occursin(String(bn), String(i)))...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    oftype(a::NamedTuple, b::DataType)
Returns a NamedTuple which retains the fields whose elements have type `b`.
```jldoctest
julia> QueryOperators.NamedTupleUtilities.oftype((a = [4,5,6], b = [3.,2.,1.], c = ["He","llo","World!"]), Val(Int64))
(a = [4, 5, 6],)
julia> QueryOperators.NamedTupleUtilities.oftype((a = [4,5,6], b = [3.,2.,1.], c = ["He","llo","World!"]), Val(Number))
(a = [4, 5, 6], b = [3., 2., 1.])
julia> QueryOperators.NamedTupleUtilities.oftype((a = [4,5,6], b = [3.,2.,1.], c = ["He","llo","World!"]), Val(Float32))
NamedTuple()
```
"""
@generated function oftype(a::NamedTuple{an}, ::Val{b}) where {an, b}
    names = ((i for i in an if fieldtype(a, i) <: b)...,)
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

end


text:
using Test
using Printf: @printf
using QuantumLattices.Essentials.QuantumOperators
using QuantumLattices.Interfaces: id, value, rank, add!, sub!, mul!, div!
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Prerequisites.Combinatorics: Combinations
using QuantumLattices.Prerequisites.Traits: contentnames, parameternames, parametertype, isparameterbound

import QuantumLattices.Interfaces: ⊗, ⋅, permute

struct AID{O<:Real, S<:Real} <: SingularID
    orbital::O
    nambu::S
end

@testset "ID" begin
    @test ID{SingularID, 1}|>rank == 1
    @test ID|>rank == Any
    @test promote_type(ID{AID{Int, Int}, 1}, ID{AID{Int, Int}}) == ID{AID{Int, Int}}
    @test promote_type(ID{AID{Int, Int}}, ID{AID{Int, Int}, 1}) == ID{AID{Int, Int}}
    @test promote_type(ID{AID{Int, Int}}, ID{AID{Int, Float}}) == ID{AID{Int, <:Real}}
    @test promote_type(Tuple{}, ID{AID{Int, Int}}) == ID{AID{Int, Int}}
    @test promote_type(Tuple{}, ID{AID{Int, Int}, 2}) == ID{AID{Int, Int}}

    cid = AID(2, 1) * AID(1, Inf)
    @test cid == ID(AID(2, 1), AID(1, Inf))
    @test cid|>string == "ID(AID(2, 1), AID(1, Inf))"
    @test cid|>eltype == AID{Int, <:Real}
    @test cid|>rank == 2
    @test cid|>typeof|>rank == 2

    sid = AID(1, 1)
    @test sid * cid == ID(AID(1, 1), AID(2, 1), AID(1, Inf))
    @test cid * sid == ID(AID(2, 1), AID(1, Inf), AID(1, 1))
    @test cid * cid == ID(AID(2, 1), AID(1, Inf), AID(2, 1), AID(1, Inf))

    @test deepcopy(ID(AID(1, 1))) == ID(AID(1, 1))
    @test isless(SingularID, ID(AID(1, 2)), ID(AID(1, Inf)))
    @test isless(SingularID, ID(AID(2, 1)), ID(AID(1, 2), AID(1, Inf)))

    cid = AID(2, 1) * AID(3, 4.0)
    @test cid == ID(AID, (2, 3), (1, 4.0))
    @test cid|>propertynames == (:orbitals, :nambus)
    @test cid.orbitals == (2, 3)
    @test cid.nambus == (1, 4.0)
end

struct AOperator{V<:Number, I<:ID{SingularID}} <: OperatorProd{V, I}
    value::V
    id::I
end
⊗(m₁::AOperator, m₂::AOperator) = m₁ * m₂
⋅(m₁::AOperator, m₂::AOperator) = m₁ * m₂
Base.show(io::IO, opt::AOperator) = @printf io "AOperator(value=%s, id=%s)" opt.value opt.id
Base.repr(opt::AOperator) = "AOperator($(opt.value), $(opt.id))"

function permute(id₁::AID, id₂::AID)
    @assert id₁ ≠ id₂ "permute error: permuted ids should not be equal to each other."
    if (id₁.nambu == 3-id₂.nambu) && (id₁.orbital == id₂.orbital)
        if id₁.nambu == 2
            return (AOperator(1, ID()), AOperator(1, ID(id₂, id₁)))
        else
            return (AOperator(-1, ID()), AOperator(1, ID(id₂, id₁)))
        end
    else
        return (AOperator(1, ID(id₂, id₁)),)
    end
end

@testset "QuantumOperator" begin
    @test parameternames(QuantumOperator) == (:id, :value)
    @test isparameterbound(QuantumOperator, :value, Any) == true
    @test isparameterbound(QuantumOperator, :id, ID) == true
    @test isparameterbound(QuantumOperator, :id, ID{AID{Int, Int}, 2}) == false

    @test contentnames(OperatorProd) == (:value, :id)
    @test parameternames(OperatorProd) == (:value, :id)
    @test isparameterbound(OperatorProd, :value, Any) == false
    @test valtype(OperatorProd) == parametertype(OperatorProd, :value) == parametertype(OperatorProd, 1) == Any
    @test valtype(OperatorProd{Int}) == parametertype(OperatorProd{Int}, :value) == parametertype(OperatorProd{Int}, 1) == Int
    @test idtype(OperatorProd{<:Number}) == parametertype(OperatorProd{<:Number}, :id) == parametertype(OperatorProd{<:Number}, 2) == ID{SingularID}
    @test idtype(OperatorProd{<:Number, ID{AID}}) == parametertype(OperatorProd{<:Number, ID{AID}}, :id) == parametertype(OperatorProd{<:Number, ID{AID}},2) == ID{AID}
    @test promote_type(AOperator{Int}, AOperator) == AOperator
    @test promote_type(AOperator, AOperator{Int}) == AOperator
    @test promote_type(AOperator{Int}, AOperator{Float}) == AOperator{Float}
    @test promote_type(AOperator{Int, ID{AID{Int, Int}, 2}}, AOperator{Float, ID{AID{Int, Int}, 2}}) == AOperator{Float, ID{AID{Int, Int}, 2}}
    @test promote_type(AOperator{Int, ID{AID}}, AOperator{Float, ID{AID}}) == AOperator{Float, <:ID{AID}}
    @test promote_type(AOperator{Complex{Int}, ID{AID}}, Float, Val(:*)) == promote_type(Float, AOperator{Complex{Int}, ID{AID}}, Val(:*)) == AOperator{Complex{Float}, ID{AID}}

    opt = AOperator(1.0, ID(AID(1, 1)))
    @test value(opt) == 1.0
    @test id(opt) == ID(AID(1, 1))
    @test opt|>deepcopy == opt
    @test isequal(opt|>deepcopy, opt)
    @test isapprox(opt, replace(opt, value=opt.value+10^-6); atol=10^-5)
    @test opt|>valtype == opt|>typeof|>valtype == parametertype(opt|>typeof, :value) == Float
    @test opt|>idtype == opt|>typeof|>idtype == parametertype(opt|>typeof, :id) == ID{AID{Int, Int}, 1}
    @test opt|>rank == opt|>typeof|>rank == 1
    @test opt[1] == AOperator(1, ID(AID(1, 1)))
    @test length(opt) == 1
    @test firstindex(opt) == 1
    @test lastindex(opt) ==1
    @test opt|>typeof|>one == opt|>one == AOperator(1.0, ID())
    @test +opt == opt
    @test -opt == AOperator(-1.0, ID(AID(1, 1)))
    @test opt*2 == 2*opt == opt*AOperator(2, ID()) == AOperator(2, ID())*opt == AOperator(2.0, ID(AID(1, 1)))
    @test opt/2 == opt/AOperator(2, ID()) == AOperator(0.5, ID(AID(1, 1)))
    @test opt^2 == opt*opt

    @test convert(AOperator{Float, Tuple{}}, AOperator(2, ID())) == AOperator(2.0, ID())
    @test convert(AOperator{Float, Tuple{}}, 3) == AOperator(3.0, ID())
    @test AOperator(1, ID())+1 == 1+AOperator(1, ID()) == AOperator(1, ID())+AOperator(1, ID()) == AOperator(2, ID())
    @test opt+1 == 1+opt == opt+AOperator(1, ID()) == AOperator(1, ID())+opt
    @test AOperator(1, ID())-2 == 1-AOperator(2, ID()) == AOperator(1, ID())-AOperator(2, ID()) == AOperator(-1, ID())
    @test AOperator(2, ID())*AOperator(3, ID()) == AOperator(6, ID())

    opt₁ = AOperator(1.0, ID(AID(1, 1)))
    opt₂ = AOperator(2.0, ID(AID(1, 2)))
    opts = OperatorSum(opt₁, opt₂)
    optp₁ = promote_type(typeof(opts), OperatorSum{NTuple{2, AID{Float, Float}}, AOperator{Complex{Int}, NTuple{2, AID{Float, Float}}}})
    optp₂ = OperatorSum{Tuple{AID, Vararg{AID{Float, Float}}}, AOperator{Complex{Float}, <:Tuple{AID, Vararg{AID{Float, Float}}}}}
    @test optp₁ == optp₂
    @test string(opts) == "OperatorSum with 2 AOperator:\n  AOperator(value=2.0, id=ID(AID(1, 2)))\n  AOperator(value=1.0, id=ID(AID(1, 1)))\n"
    @test repr(opts) == "OperatorSum with 2 AOperator:\n  AOperator(2.0, ID(AID(1, 2)))\n  AOperator(1.0, ID(AID(1, 1)))"
    @test opts|>deepcopy == opts
    @test isapprox(opts, deepcopy(opts)) && isapprox(opts, opts+10^-6; atol=10^-5) && isapprox(opts+10^-6, opts; atol=10^-5)
    @test !isapprox(opts, opts+10^-6) && !isapprox(opts+10^-6, opts)
    @test add!(deepcopy(opts)) == opts
    @test sub!(deepcopy(opts)) == opts
    @test mul!(deepcopy(opts), 2.0) == mul!(deepcopy(opts), AOperator(2, ID())) == opts*2 == opts*AOperator(2, ID())
    @test div!(deepcopy(opts), 2.0) == div!(deepcopy(opts), AOperator(2, ID())) == opts/2 == opts/AOperator(2, ID())
    @test +opts == opts
    @test opt₁+opt₂ == opts
    @test opt₁+opts == opts+opt₁ == OperatorSum(opt₁*2, opt₂)
    @test opts+opts == OperatorSum(opt₁*2, opt₂*2)
    @test -opts == OperatorSum(-opt₁, -opt₂)
    @test 1.0-opts == AOperator(1.0, ID())-opts == -opts+1.0
    @test opts-1.0 == opts-AOperator(1.0, ID())
    @test opt-1.0 == opt-AOperator(1.0, ID()) == opt+AOperator(-1.0, ID())
    @test 1.0-opt == AOperator(1.0, ID())-opt == 1.0+(-opt)
    @test opts-opt₁ == OperatorSum(opt₂)
    @test opt₁-opts == OperatorSum(-opt₂)
    @test opts-opts == zero(opts) == zero(typeof(opts))
    @test 2*opts == opts*2 == OperatorSum(opt₁*2, opt₂*2) == AOperator(2, ID())*opts
    @test opts/2 == OperatorSum(opt₁/2, opt₂/2)
    @test opts^2 == opts*opts

    temp = OperatorSum{ID{AID{Int, Int}}, AOperator{Float}}(opts)
    @test add!(deepcopy(temp), 1) == add!(deepcopy(temp), AOperator(1, ID())) == opts+1 == 1+opts == opts+AOperator(1, ID()) == AOperator(1, ID())+opts
    @test sub!(deepcopy(temp), 1) == sub!(deepcopy(temp), AOperator(1, ID())) == opts-1 == opts-AOperator(1, ID())

    opt₃ = AOperator(3, ID(AID(1, 3)))
    @test opt₁*opt₂ == opt₁⊗opt₂ == opt₁⋅opt₂ == AOperator(2.0, ID(AID(1, 1), AID(1, 2)))
    @test opts*opt₃ == opts⊗opt₃ == opts⋅opt₃ == OperatorSum(AOperator(3.0, ID(AID(1, 1), AID(1, 3))), AOperator(6.0, ID(AID(1, 2), AID(1, 3))))
    @test opt₃*opts == opt₃⊗opts == opt₃⋅opts == OperatorSum(AOperator(3.0, ID(AID(1, 3), AID(1, 1))), AOperator(6.0, ID(AID(1, 3), AID(1, 2))))
    @test opts*OperatorSum(opt₃) == opts⊗OperatorSum(opt₃) == opts⋅OperatorSum(opt₃) == OperatorSum(AOperator(3.0, ID(AID(1, 1), AID(1, 3))), AOperator(6.0, ID(AID(1, 2), AID(1, 3))))

    table = Dict(AID(1, 1)=>1, AID(1, 2)=>2)
    opt = AOperator(2.0, ID(AID(1, 1), AID(1, 2)))
    @test sequence(opt, table) == (1, 2)
    @test split(opt) == (2.0, AOperator(1.0, ID(AID(1, 1))), AOperator(1.0, ID(AID(1, 2))))

    opt₁ = AOperator(1, ID(AID(1, 1)))
    opt₂ = AOperator(2, ID(AID(1, 2)))
    opts = OperatorSum(opt₁, opt₂)
    opt = AOperator(3, ID())
    @test opt₁//3 == opt₁//opt == AOperator(1//3, ID(AID(1, 1)))
    @test opt₂//3 == opt₂//opt == AOperator(2//3, ID(AID(1, 2)))
    @test opts//3 == opts//opt == OperatorSum(AOperator(1//3, ID(AID(1, 1))), AOperator(2//3, ID(AID(1, 2))))
end

@testset "replace" begin
    id₁, id₂ = AID(1, 1), AID(1, 2)
    opt = AOperator(1.5, ID(id₁, id₂))
    id₃, id₄ = AID(2, 1), AID(2, 2)
    opt₁ = AOperator(2.0, ID(id₃))+AOperator(3.0, ID(id₄))
    opt₂ = AOperator(2.0, ID(id₃))-AOperator(3.0, ID(id₄))
    @test replace(opt, id₁=>opt₁, id₂=>opt₂) == 1.5*opt₁*opt₂
    @test replace(OperatorSum(opt), id₂=>opt₂, id₁=>opt₁) == 1.5*opt₁*opt₂
end

@testset "permute" begin
    id₁, id₂ = AID(1, 1), AID(1, 2)
    opt₁ = AOperator(1.5, ID(id₁, id₂))
    opt₂ = AOperator(1.5, ID(id₂, id₁))
    opt₀ = AOperator(1.5, ID())
    @test permute(opt₁, Dict(id₁=>1, id₂=>2)) == opt₂-opt₀
    @test permute(opt₁+opt₂, Dict(id₁=>1, id₂=>2)) == 2*opt₂-opt₀
    @test permute(opt₂, Dict(id₁=>2, id₂=>1)) == opt₁+opt₀
    @test permute(opt₁+opt₂, Dict(id₁=>2, id₂=>1)) == 2*opt₁+opt₀
end

@testset "Transformation" begin
    m = AOperator(1, ID(AID(1, 1)))
    s = OperatorSum(m)

    i = Identity()
    @test i==deepcopy(i) && isequal(i, deepcopy(i))
    @test valtype(i, m) == valtype(typeof(i), typeof(m)) == typeof(m)
    @test valtype(i, s) == valtype(typeof(i), typeof(s)) == typeof(s)
    @test i(m)==m && i(s)==s

    n = Numericalization{Float64}()
    @test valtype(n) == valtype(typeof(n)) == Float64
    @test valtype(n, m) == valtype(typeof(n), typeof(m)) == AOperator{Float64, ID{AID{Int, Int}, 1}}
    @test valtype(n, s) == valtype(typeof(n), typeof(s)) == OperatorSum{ID{AID{Int, Int}, 1}, AOperator{Float64, ID{AID{Int, Int}, 1}}}
    @test n(m)==replace(m, value=1.0) && n(s)==OperatorSum(replace(m, value=1.0))
end


text:
## GENERATE SYNTHETIC DATA:

function noisy(label; λ=0.0)
    if label
        return 1 - λ*rand()
    else
        return λ*rand()
    end
end

labels = rand(Bool, 200);

scores(λ) = map(labels) do label
    noisy(label, λ=λ)
end

## COMPARE AUC

using ROC

roc_good = roc(scores(0.6), labels, true);
roc_bad = roc(scores(1.0), labels, true);
area_good = AUC(roc_good)
area_bad =  AUC(roc_bad)

## PLOT THE ROC's

using Plots
pyplot()

plot(roc_good, label="good");
plot!(roc_bad, label="bad")

savefig("rocs.png")


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export iseps
"""
  `iseps(x)`

Checks if `x` is the epsilon label. Defined only for `String` and `Int`, where `<eps>` and `0` are the epsilon symbols. Extend this method if you want to create WFST with a user defined type.
"""
iseps(x::T) where {T<:Integer} = x == zero(T)
iseps(x::String) = x == "<eps>"
iseps(x::Char) = x == 'ϵ'

export get_eps
"""
  `get_eps(x::Type)`

Returns the epsilon label for a given `Type`. Defined only for `String` and `Int`, where `<eps>` and `0` are the epsilon symbols. Extend this method if you want to create WFST with a user defined type.
"""
get_eps(x::Type{T}) where {T<:Integer} = T(0)
get_eps(x::Type{String}) = "<eps>"
get_eps(x::Type{Char}) = 'ϵ'


text:
module LCMGL

depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(depsjl) ? include(depsjl) : error("LCMGL not properly ",
    "installed. Please run\nPkg.build(\"LCMGL\")")

import Base: unsafe_convert
using Libdl
export LCM, LCMGLClient,
	switch_buffer,
	begin_mode,
	end_mode,
	vertex,
	color,
	normal,
	scale_axes,
	point_size,
	line_width,
	translate,
	rotate,
	push_matrix,
	pop_matrix,
	load_identity,
	sphere,
	draw_axes

POINTS         = 0x0000
LINES          = 0x0001
LINE_LOOP      = 0x0002
LINE_STRIP     = 0x0003
TRIANGLES      = 0x0004
TRIANGLE_STRIP = 0x0005
TRIANGLE_FAN   = 0x0006
QUADS          = 0x0007
QUAD_STRIP     = 0x0008
POLYGON        = 0x0009

mutable struct LCM
    pointer::Ptr{Cvoid}

    LCM() = begin
        lc = new(ccall((:lcm_create, liblcm), Ptr{Cvoid}, (Ptr{UInt8},), ""))
		finalizer(close, lc)
        lc
    end
end
unsafe_convert(::Type{Ptr{Cvoid}}, lc::LCM) = lc.pointer

function close(lcm::LCM)
	if lcm.pointer != C_NULL
	    ccall((:lcm_destroy, liblcm), Cvoid, (Ptr{Cvoid},), lcm)
		lcm.pointer = C_NULL
	end
end

function LCM(func::Function)
    lc = LCM()
    try
        func(lc)
    finally
        close(lc)
    end
end

struct Clcmgl
    lcm::Ptr{Cvoid}
    name::Ptr{Cchar}
    channel_name::Ptr{Cchar}
    scene::Int32
    sequence::Int32
    data::Ptr{UInt8}
    datalen::Cint
    data_alloc::Cint
    texture_count::UInt32
end

mutable struct LCMGLClient
    lcm::LCM
	name::AbstractString
    pointer::Ptr{Clcmgl}

    LCMGLClient(lcm::LCM, name::AbstractString) = begin
        gl = new(lcm, name,
                 ccall((:bot_lcmgl_init, libbot2_lcmgl_client),
				       Ptr{Clcmgl}, (Ptr{Cvoid}, Ptr{UInt8}), lcm, name))
        finalizer(close, gl)
        gl
    end
end
unsafe_convert(::Type{Ptr{Clcmgl}}, gl::LCMGLClient) = gl.pointer

function datalen(lcmgl::LCMGLClient)
    (lcmgl.pointer == C_NULL) && error("LCMGLClient is unininitialized")
    cgl = unsafe_load(lcmgl.pointer)
    cgl.datalen
end

function close(lcmgl::LCMGLClient)
	if lcmgl.pointer != C_NULL
		ccall((:bot_lcmgl_destroy, libbot2_lcmgl_client),
		      Cvoid, (Ptr{Clcmgl},), lcmgl)
		lcmgl.pointer = C_NULL
	end
end

LCMGLClient(name::AbstractString) = LCMGLClient(LCM(), name)

function LCMGLClient(func::Function, name::AbstractString, automatically_switch_buffer::Bool=true)
    LCM() do lc
        LCMGLClient(func, lc, name, automatically_switch_buffer)
    end
end

function LCMGLClient(func::Function, lcm::LCM, name::AbstractString, automatically_switch_buffer::Bool=true)
    gl = LCMGLClient(lcm, name)
    try
        func(gl)
        if automatically_switch_buffer && datalen(gl) > 0
            switch_buffer(gl)
        end
    finally
		close(gl)
    end
end

switch_buffer(gl::LCMGLClient) = ccall((:bot_lcmgl_switch_buffer, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl},), gl)

# begin and end are reserved keywords in Julia, so I've renamed
# them to begin_mode and end_mode
begin_mode(gl::LCMGLClient, mode::Integer) = ccall((:bot_lcmgl_begin, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cuint), gl, mode)
end_mode(gl::LCMGLClient) = ccall((:bot_lcmgl_end, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl},), gl)

vertex(gl::LCMGLClient, x, y) = ccall((:bot_lcmgl_vertex2d, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cdouble, Cdouble), gl, x, y)
vertex(gl::LCMGLClient, x, y, z) = ccall((:bot_lcmgl_vertex3d, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cdouble, Cdouble, Cdouble), gl, x, y, z)

color(gl::LCMGLClient, red, green, blue) = ccall((:bot_lcmgl_color3f, libbot2_lcmgl_client),
    Cvoid, (Ptr{Clcmgl}, Cfloat, Cfloat, Cfloat), gl.pointer, red, green, blue)
color(gl::LCMGLClient, red, green, blue, alpha) = ccall((:bot_lcmgl_color4f, libbot2_lcmgl_client),
    Cvoid, (Ptr{Clcmgl}, Cfloat, Cfloat, Cfloat, Cfloat), gl.pointer, red, green, blue, alpha)
normal(gl::LCMGLClient, x, y, z) = ccall((:bot_lcmgl_normal3f, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cfloat, Cfloat, Cfloat), gl, x, y, z)
scale_axes(gl::LCMGLClient, x, y, z) = ccall((:bot_lcmgl_scalef, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cfloat, Cfloat, Cfloat), gl, x, y, z)

point_size(gl::LCMGLClient, size) = ccall((:bot_lcmgl_point_size, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cfloat), gl, size)
line_width(gl::LCMGLClient, width) = ccall((:bot_lcmgl_line_width, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cfloat), gl, width)

translate(gl::LCMGLClient, v0, v1, v2) = ccall((:bot_lcmgl_translated, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cdouble, Cdouble, Cdouble), gl, v0, v1, v2)
rotate(gl::LCMGLClient, angle, x, y, z) = ccall((:bot_lcmgl_rotated, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl}, Cdouble, Cdouble, Cdouble, Cdouble), gl, angle, x, y, z)
push_matrix(gl::LCMGLClient) = ccall((:bot_lcmgl_push_matrix, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl},), gl)
pop_matrix(gl::LCMGLClient) = ccall((:bot_lcmgl_pop_matrix, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl},), gl)

sphere(gl::LCMGLClient, origin, radius, slices, stacks) = ccall((:bot_lcmgl_sphere, libbot2_lcmgl_client),
    Cvoid, (Ptr{Clcmgl}, Ptr{Cdouble}, Cdouble, Cint, Cint), gl, origin, radius, slices, stacks)

draw_axes(gl::LCMGLClient) = ccall((:bot_lcmgl_draw_axes, libbot2_lcmgl_client), Cvoid, (Ptr{Clcmgl},), gl)

function __init__()
    @static if Sys.islinux()
        Libdl.dlopen(liblcm, Libdl.RTLD_GLOBAL)
    end
end

end

import LCMGL


text:
using DashBootstrapComponents, DashHtmlComponents

list_group = html_div([
    dbc_listgroup([
        dbc_listgroupitem("Internal link", href = "/l/components/list_group"),
        dbc_listgroupitem("External link", href = "https://google.com"),
        dbc_listgroupitem("Disabled link", href = "https://google.com", disabled = true),
        dbc_listgroupitem("Button", id = "button-item", n_clicks = 0, action = true),
    ]),
    html_p(id = "counter"),
]);

callback!(app, Output("counter", "children"), Input("button-item", "n_clicks")) do n
    return "Button clicked $n times"
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
 BenchmarkDefaults =
     BenchmarkTools.Parameters(
       seconds = 4.0,
       samples = 7_500,
       evals = 1,
       overhead = BenchmarkTools.estimate_overhead(),
       gctrial = true,
       gcsample = false,
       time_tolerance = 1.25e-11,
       memory_tolerance = 1/4096, # was 1/1024,
     );

BenchmarkTools.DEFAULT_PARAMETERS.overhead = BenchmarkTools.estimate_overhead()
BenchmarkTools.DEFAULT_PARAMETERS.evals = 1
BenchmarkTools.DEFAULT_PARAMETERS.samples = 7_500
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 1.35e-11
BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = 1/4096

# introduce `@refd` benchmarking assist (`@refd @btime ..`)

walk(x, inner, outer) = outer(x)
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))
MacroTools.postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

function _refd(expr::Expr)
    if expr.head == :$
        :($(Expr(:$, :(Ref($(expr.args...)))))[])
    else
        expr
    end
end
_refd(x) = x

macro refd(expr)
    out = postwalk(_refd, expr) |> esc
end


text:
using PyDSTool, Test

@time begin
@time @testset "Linear Solve" begin include("linear.jl") end
@time @testset "Calcium Direct" begin include("calcium.jl") end
@time @testset "Calcium Parameterized" begin include("pf_calcium.jl") end
end


text:
# abs, iszero, isone, sign, and sqrt specialize Base's functions
# specially, float masks to_double
export compare,
       integral_division,
       inverse,
       is_negative,
       is_positive,
       is_square,
       square,
       to_interval,
       unit_part

const NT = RT
export NT

if !iscxxtype(NT)
    for F ∈ (:abs, :sqrt, :float, :iszero, :isone, :sign)
        @eval Base.$F(x::CxxBaseRef{NT}) = $F(x[])
    end
end

"""
    compare(x::NT, y::NT)

The function [`compare()`](@ref) compares the first argument with respect to the
second, i.e. it returns [`LARGER`](@ref) if ``x`` is larger than ``y``.
"""
compare(x::Real, y::Real) = compare(convert(NT, x), convert(NT, y))

"""
    integral_division(x::NT, y::NT)

The function [`integral_division()`](@ref) (a.k.a. exact division or division
without remainder) maps ring elements ``(x, y)`` to ring element ``z`` such that
``x = yz`` if such a ``z`` exists (i.e. if ``x`` is divisible by ``y``).
"""
integral_division(x::Real, y::Real) = integral_division(convert(NT, x),
                                                        convert(NT, y))

"""
    inverse(x::NT)

The function [`inverse()`](@ref) returns the inverse element with respect to
multiplication.

!!! info "Precondition"

    ``x ≠ 0``.
"""
inverse(x::Real) = inverse(convert(NT, x))

"""
    is_negative(x::NT)

The function [`is_negative()`](@ref) determines if a value is negative or not.
"""
is_negative(x::Real) = is_negative(convert(NT, x))

"""
    is_positive(x::NT)

The function [`is_positive()`](@ref) determines if a value is positive or not.
"""
is_positive(x::Real) = is_positive(convert(FT, x))

"""
    is_square(x::NT, y::Ref{NT} = Ref(0))

A ring element ``x`` is said to be a square iff there exists a ring element
``y`` such that ``x = y * y``.
"""
is_square(x::Real, y::Ref{Real} = Ref(0)) = is_square(convert(NT, x), Ref(convert(NT, y[])))

"""
    to_interval(x::NT)

The function [`to_interval()`](@ref) computes for a given number ``x`` a double
interval containing ``x``.
"""
to_interval(x::Real) = to_interval(convert(NT, x))

"""
    unit_part(x::NT)

The function [`unit_part()`](@ref) computes the unit part of a given ring
element.
"""
unit_part(x::Real) = unit_part(convert(NT, x))


text:
using Dispersal, Test

init = [1.0 4.0 7.0;
        2.0 5.0 8.0;
        3.0 6.0 9.0]

output = ArrayOutput(init; tspan=1:2)

rule = Ruleset(AlleeExtinction(minfounders = 8.0))

sim!(output, rule)
@test output[1] == [1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

@test output[2] == [0.0 0.0 0.0;
                    0.0 0.0 8.0;
                    0.0 0.0 9.0]


text:
# http://web.mit.edu/puzzle/www/2012/puzzles/watson_2_0/r_e_s_p_e_c_t/
@testset "R.E.S.P.E.C.T" begin
    f = best_feature(String.(shuffle.(collect.(wordlist("""
        ABMNOT
        AENORTY
        BCEKLORSTU
        BFLU
        CDEILNOTU
        CIK
        GIOPS
        ABCEKNO
        ABLMNOR
        ACEILMPR
        ACEILTUV
        HINOPSY
        ACDEIKLNST
        ACDEMY
        ACEMZ
        ACHIMNST
        DEIO
        GLNOTU
        ABCEGILORTY
        ABEINORTX
        ACRTUY
        AGHTUY
        CDEIMNOPU
        ACDILMOPT
        ACHILP
        BCEKUY
        MNOTU
        """)))), 5)
    @test f.description == "has a 1-letter transaddition"
end


text:
function getLabel(args...; Tree::Array{SubString{String},1} = [],
	Start::Int = Int[], stop_char::Array{String,1} = [])

  label = String[]
	u = 0
	while countnz(Tree[u+Start] .== stop_char) == 0
		u = u + 1
		label = join(Tree[Start:(u - 1 + Start)])
	end
  labelStart = (label,u + Start)
	labelStart = Dict(zip(["label","end"], labelStart))
	return labelStart
end


text:
# function to build the fast application of the dual Lippmann-Schwinger equation
# we suppose u(x) = \int_D G(x,y) \sigma(y) dy
# and we solve ( \Delta + omega^2(1+q)) u = - omega^2q u^{inc}

type FastMDual
    # type to encapsulate the fast application of M = I + omega^2G*spadiagm(nu)
    GFFT :: Array{Complex128,2}
    nu :: Array{Float64,1}
    # number of points in the extended domain
    ne :: Int64
    me :: Int64
    # number of points in the original domain
    n  :: Int64
    m  :: Int64
    # frequency
    omega :: Float64
    quadRule :: String
    function FastMDual(GFFT,nu,ne,me,n,m,k; quadRule::String = "trapezoidal")
      return new(GFFT,nu,ne,me,n, m, k, quadRule)
    end
end

# we need to overload these functions for gmres
import Base.*
import Base.A_mul_B!
import Base.size

function size(M::FastMDual, dim::Int64)
  # function to returns the size of the underliying matrix (M.m*M.n)^2
  return M.m*M.n
end

function *(M::FastMDual, b::Array{Complex128,1})
    # function to overload the applyication of
    # M using a Toeplitz reduction via a FFT
    # dummy function to call fastconvolution
    return fastconvolution(M,b)
end

function A_mul_B!(Y,
                  M::FastMDual,
                  V)
    # in place matrix matrix multiplication
    @assert(size(Y) == size(V))
    # print(size(V))
    for ii = 1:size(V,2)
        Y[:,ii] = M*V[:,ii]
    end
end

@inline function fastconvolution(M::FastMDual, b::Array{Complex128,1})
    # function to overload the applyication of
    # M using a Toeplitz reduction via a FFT

    # computing omega^2 nu G*(b)
    B = M.omega^2*(M.nu.*FFTconvolution(M,b))

    # returning -b + omega^2(nu G*b)
    return (-b + B[:])
end

@inline function FFTconvolution(M::FastMDual, b::Array{Complex128,1})
    # function to overload the applyication of
    # convolution of b times G

    if M.quadRule == "trapezoidal"

      #obtaining the middle index
      indMiddle = round(Integer, M.n)

      # Allocate the space for the extended B
      BExt = zeros(Complex128,M.ne, M.me);
      # Apply spadiagm(nu) and ented by zeros
      BExt[1:M.n,1:M.m]= reshape(M.nu.*b,M.n,M.m) ;

      # Fourier Transform
      BFft = fft(BExt)
      # Component-wise multiplication
      BFft = M.GFFT.*BFft
      # Inverse Fourier Transform
      BExt = ifft(BFft)

      # multiplication by omega^2
      B = (BExt[M.n:M.n+M.n-1, M.m:M.m+M.m-1]);

    elseif M.quadRule == "Greengard_Vico"
      # for this we use the Greengard Vico method in the
      # frequency domain

      # Allocate the space for the extended B
      BExt = zeros(Complex128,M.ne, M.me);
      # Apply spadiagm(nu) and ented by zeros
      BExt[1:M.n,1:M.m]= reshape(M.nu.*b,M.n,M.m) ;

      # Fourier Transform
      BFft = fftshift(fft(BExt))
     # Component-wise multiplication
      BFft = M.GFFT.*BFft
      # Inverse Fourier Transform
      BExt = ifft(ifftshift(BFft))

      # multiplication by omega^2
      B = BExt[1:M.n, 1:M.m];
    end
    return  B[:]
end

function buildFastConvolutionDual(x::Array{Float64,1},y::Array{Float64,1},
                              h::Float64,k,nu::Function; quadRule::String = "trapezoidal")

  if quadRule == "trapezoidal"

    (ppw,D) = referenceValsTrapRule();
    D0      = D[round(Int,k*h)];
    (n,m) = length(x), length(y)
    Ge    = buildGConv(x,y,h,n,m,D0,k);
    GFFT  = fft(Ge);
    X = repmat(x, 1, m)[:]
    Y = repmat(y', n,1)[:]

    return FastMDual(GFFT,nu(X,Y),2*n-1,2*m-1,n, m, k);

  elseif quadRule == "Greengard_Vico"

      Lp = 4*(abs(x[end] - x[1]) + h)
      L  =   (abs(x[end] - x[1]) + h)*1.5
      (n,m) = length(x), length(y)
      X = repmat(x, 1, m)[:]
      Y = repmat(y', n,1)[:]

      # this is depending if n is odd or not
      if mod(n,2) == 0
        kx = (-(2*n):1:(2*n-1));
        ky = (-(2*m):1:(2*m-1));

        KX = (2*pi/Lp)*repmat(kx, 1, 4*m);
        KY = (2*pi/Lp)*repmat(ky', 4*n,1);

        S = sqrt.(KX.^2 + KY.^2);

        GFFT = Gtruncated2D(L, k, S)
        return FastMDual(GFFT, nu(X,Y), 4*n, 4*m,
                     n, m, k , quadRule="Greengard_Vico");
      else
        # kx = (-2*(n-1):1:2*(n-1) )/4;
        # ky = (-2*(m-1):1:2*(m-1) )/4;

        # KX = (2*pi/Lp)*repmat(kx, 1, 4*m-3);
        # KY = (2*pi/Lp)*repmat(ky', 4*n-3,1);

        # S = sqrt.(KX.^2 + KY.^2);

        # GFFT = Gtruncated2D(L, k, S)

        # return FastM(GFFT,nu(X,Y),4*n-3,4*m-3,
        #              n,m, k,quadRule = "Greengard_Vico");

        kx = (-2*n:1:2*n-1);
        ky = (-2*m:1:2*m-1);

        KX = (2*pi/Lp)*repmat( kx, 1,4*m);
        KY = (2*pi/Lp)*repmat(ky',4*n,  1);

        S = sqrt.(KX.^2 + KY.^2);

        GFFT = Gtruncated2D(L, k, S)

        return FastMDual(GFFT,nu(X,Y),4*n,4*m,
                     n,m, k,quadRule = "Greengard_Vico");

    end
  end
end

##########################################################################
########  Function fo build the preconditioner for the  ##################
########        Dual Lippmann-Schwinger equation    ######################
##########################################################################

function buildGSparseAandG(k::Float64,X::Array{Float64,1},Y::Array{Float64,1},
                          D0::Complex128, n::Int64 ,m::Int64, nu::Function ; method::String = "normal")

    nuArray = nu(X,Y)
    Ind = reshape(collect(1:n*m),n,m);

    if method=="normal"
      (Indices, Values) = entriesSparseA(k,X,Y,D0, n ,m);
    elseif method == "randomized"
      (Indices, Values) = entriesSparseARand(k,X,Y,D0, n ,m);
    end

    Entries = entriesSparseG(k,X,Y,D0, n ,m);

    (rowA, colA, valA) = createIndices(Ind[2:end-1,2:end-1][:],
                                       Indices[1][:], Values[1][:]);

    (rowG, colG, valG) = createIndices(Ind[2:end-1,2:end-1][:],
                                       Indices[1][:],Values[1][:], Entries[1],nuArray);

    (Row, Col, Val) = createIndices(Ind[1,2:end-1][:],
                                    Indices[2][:], Values[2][:]);

    rowA = vcat(rowA,Row);
    colA = vcat(colA,Col);
    valA = vcat(valA,Val);

    (Row, Col, Val) = createIndices(Ind[end,2:end-1][:],
                                    Indices[3][:], Values[3][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    (Row, Col, Val) = createIndices(Ind[2:end-1,1][:],
                                    Indices[4][:], Values[4][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    (Row, Col, Val) = createIndices(Ind[2:end-1,end][:],
                                    Indices[5][:], Values[5][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    (Row, Col, Val) = createIndices(Ind[1,1],
                                    Indices[6][:], Values[6][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    (Row, Col, Val) = createIndices(Ind[end,1],
                                    Indices[7][:], Values[7][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    (Row, Col, Val) = createIndices(Ind[1,end],
                                    Indices[8][:], Values[8][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    (Row, Col, Val) = createIndices(Ind[end,end],
                                    Indices[9][:], Values[9][:]);

    rowA = vcat(rowA,Row)
    colA = vcat(colA,Col)
    valA = vcat(valA,Val)

    A = sparse(rowA,colA,valA);

    (Row, Col, Val) = createIndices(Ind[1,2:end-1][:],
                                    Indices[2][:],Values[2][:], Entries[2],nuArray);

    rowG = vcat(rowG,Row);
    colG = vcat(colG,Col);
    valG = vcat(valG,Val);

    (Row, Col, Val) = createIndices(Ind[end,2:end-1][:],
                                    Indices[3][:], Values[3][:], Entries[3],nuArray);
    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    (Row, Col, Val) = createIndices(Ind[2:end-1,1][:],
                                    Indices[4][:], Values[4][:], Entries[4],nuArray);

    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    (Row, Col, Val) = createIndices(Ind[2:end-1,end][:],
                                    Indices[5][:], Values[5][:], Entries[5],nuArray);
    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    (Row, Col, Val) = createIndices(Ind[1,1],
                                    Indices[6][:], Values[6][:], Entries[6],nuArray);
    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    (Row, Col, Val) = createIndices(Ind[end,1],
                                    Indices[7][:], Values[8][:], Entries[8],nuArray);

    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    (Row, Col, Val) = createIndices(Ind[1,end],
                                    Indices[8][:], Values[8][:], Entries[7],nuArray);

    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    (Row, Col, Val) = createIndices(Ind[end,end],
                                    Indices[9][:], Values[9][:], Entries[9],nuArray);

    rowG = vcat(rowG,Row)
    colG = vcat(colG,Col)
    valG = vcat(valG,Val)

    ADGs = sparse(rowG,colG,valG);

    return (A, ADGs)

end

function createIndices(row::Array{Int64,1}, col::Array{Int64,1},
                       valA::Array{Complex128,1},
                       valG::Array{Complex128,2}, nu::Array{Float64,1})
  # function to create the indices for a sparse matrix
  @assert length(col) == length(valA)
  nn = length(col);
  mm = length(row);

  ll = length(valA)

  Row = kron(row, ones(Int64, nn));
  Col = kron(ones(Int64,mm), col) + Row;
  Val = nu[Col].*(kron(ones(Int64,mm), valA))

  # multiplying by the blocks of G.
  for ii = 1:mm
    Val[(ii-1)*ll+1:ii*ll] = (Val[(ii-1)*ll+1:ii*ll].')*valG
  end

  return (Row,Col,Val)
end

function createIndices(row::Int64, col::Array{Int64,1},
                       valA::Array{Complex128,1},
                       valG::Array{Complex128,2}, nu::Array{Float64,1})

  @assert length(col) == length(valA)
  nn = length(col);
  mm = 1;
  ll = length(valA)

  Row = kron(row, ones(Int64, nn));
  Col = kron(ones(Int64,mm), col) + Row;
  Val = nu[Col].*(kron(ones(Int64,mm), valA))

  print(size(Val))

  # multiplying by the blocks of G.
  for ii = 1:mm
    Val[(ii-1)*ll+1:ii*ll] = (Val[(ii-1)*ll+1:ii*ll].')*valG
  end

  return (Row,Col,Val)
end

# function sparsifiedMatrix()

#     As = buildSparseA(k,X,Y,D0, n ,m);

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer