text:
using Statistics
using Plots
using BenchmarkTools
using Profile
using FFTW
using Statistics
using Optim
using Measures
using Images, FileIO

#using Profile
using LinearAlgebra
using SparseArrays

# put the cwd on the path so Julia can find the module
push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils

@time fhash = fink_filter_hash(1, 8, nx=64, pc=1, wd=1)

function realspace_filter(Nx, f_i, f_v)

    zarr = zeros(ComplexF64, Nx, Nx)
    for i = 1:length(f_i)
        zarr[f_i[i]] = f_v[i] # filter*image in Fourier domain
    end
    filt = ifft(zarr)  # real space, complex
    return filt
end

Nx = 256
bar = rand(2500,100)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
#M = S2_iso_matrix(fhash, sparse=false)
#Ms = sparse(M)
Ms = fhash["S2_iso_mat"]
M = Array(Ms)
@benchmark M*bar
@benchmark Ms*bar

# M*bar test
# Nx         full   sparse
# 128       0.8 ms     4 μs
# 128x100   45  ms   500 μs
# 256       1.9 ms     7 μs
# 256x100   90  ms   800 μs

function myrealifft(fbar)
    (Nx, Ny)  = size(fbar)
    fbars = circshift(fbar[end:-1:1,end:-1:1],(1,1))
    myarg = fbar[1:Ny÷2+1,:]+conj(fbars[1:Ny÷2+1,:])
    irfft(myarg,Ny).* 0.5
end

function wst_S1S2_derivfast(image::Array{Float64,2}, filter_hash::Dict)
    #Works now.
    #Possible Bug: You assumed zero freq in wrong place and N-k is the wrong transformation.
    # Use 2 threads for FFT
    FFTW.set_num_threads(2)

    # array sizes
    (Nx, Ny)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate output array
    dS1dp  = zeros(Float64, Nx, Nx, Nf)
    dS2dp  = zeros(Float64, Nx, Nx, Nf, Nf)

    # allocate image arrays for internal use
    im_rdc_0_1 = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
    im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex

    ## 1st Order
    im_fd_0 = fft(image)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    f_ind_rev = [[CartesianIndex(mod1(Nx+2 - ci[1],Nx), mod1(Nx+2 - ci[2],Nx)) for ci in f_Nf] for f_Nf in f_ind]

    # tmp arrays for all tmp[f_i] = arr[f_i] .* f_v opns
    zarr1     = zeros(ComplexF64, Nx, Nx)
    fz_fψ1    = zeros(ComplexF64, Nx, Nx)
    fterm_a   = zeros(ComplexF64, Nx, Nx)
    fterm_ct1 = zeros(ComplexF64, Nx, Nx)
    fterm_ct2 = zeros(ComplexF64, Nx, Nx)
    #rterm_bt2 = zeros(ComplexF64, Nx, Nx)
    # make a FFTW "plan" for an array of the given size and type
    P = plan_ifft(im_fd_0)  # P is an operator, P*im is ifft(im)

    # Loop over filters
    for f1 = 1:Nf
        f_i1 = f_ind[f1]  # CartesianIndex list for filter1
        f_v1 = f_val[f1]  # Values for f_i1
        f_i1rev = f_ind_rev[f1]

        zarr1[f_i1] = f_v1 .* im_fd_0[f_i1] #F[z]
        zarr1_rd = P*zarr1 #for frzi
        fz_fψ1[f_i1] = f_v1 .* zarr1[f_i1]
        fz_fψ1_rd = P*fz_fψ1
        dS1dp[:,:,f1] = 2 .* real.(fz_fψ1_rd) #real.(conv(I_λ, ψ_λ))
        #CHECK; that this equals derivS1fast code

        #dS2 loop prep
        #ψ_λ1  = realspace_filter(Nx, f_i1, f_v1) #Slow
        #fcψ_λ1 = fft(conj.(ψ_λ1)) #Slow
        I_λ1 = sqrt.(abs2.(zarr1_rd))
        fI_λ1 = fft(I_λ1)
        rterm_bt1 = zarr1_rd./I_λ1 #Z_λ1/I_λ1_bar.
        rterm_bt2 = conj.(zarr1_rd)./I_λ1
        for f2 = 1:Nf
            f_i2 = f_ind[f2]  # CartesianIndex list for filter2
            f_v2 = f_val[f2]  # Values for f_i2

            fterm_a[f_i2] = fI_λ1[f_i2] .* (f_v2).^2 #fterm_a = F[I_λ1].F[ψ_λ]^2
            rterm_a = P*fterm_a                #Finv[fterm_a]

            fterm_bt1 = fft(rterm_a .* rterm_bt1) #F[Finv[fterm_a].*Z_λ1/I_λ1_bar] for T1
            fterm_bt2 = fft(rterm_a .* rterm_bt2) #F[Finv[fterm_a].*Z_λ1_bar/I_λ1_bar] for T2
            fterm_ct1[f_i1] = fterm_bt1[f_i1] .* f_v1    #fterm_b*F[ψ_λ]
            #fterm_ct2 = fterm_bt2 .* fcψ_λ1             #Slow
            #println(size(fterm_ct2[f_i1rev]),size(fterm_bt2[f_i1rev]),size(f_v1),size(conj.(f_v1)))
            fterm_ct2[f_i1rev] = fterm_bt2[f_i1rev] .* f_v1 #f_v1 is real
            #fterm_ct2slow = fterm_bt2 .* fcψ_λ1
            dS2dp[:, :, f1, f2] = real.(P*(fterm_ct1 + fterm_ct2))
            #println("Term2",fterm_ct2)

            #Reset f2 specific tmp arrays to 0
            fterm_a[f_i2] .= 0
            fterm_ct1[f_i1] .=0
            fterm_ct2[f_i1rev] .=0

        end
        # reset all reused variables to 0
        zarr1[f_i1] .= 0
        fz_fψ1[f_i1] .= 0

    end
    return dS1dp, dS2dp
end

function derivtestS1S2(Nx)
    eps = 1e-4
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
    im = rand(Float64, Nx, Nx).*0.1
    im[6,6]=1.0
    im1 = copy(im)
    im0 = copy(im)
    im1[2,3] += eps/2
    im0[2,3] -= eps/2
    dS1dp, dS2dp = wst_S1S2_derivfast(im, fhash)
    println(size(dS2dp))
    der0=DHC_compute(im0,fhash,doS2=true)
    der1=DHC_compute(im1,fhash,doS2=true)
    dS = (der1-der0) ./ eps

    Nf = length(fhash["filt_index"])
    i0 = 3+Nf
    blarg = dS2dp[2,3,:,:]
    diff = dS[i0:end]-reshape(blarg,Nf*Nf)
    println(dS[i0:end])
    println("and")
    println(blarg)
    println("stdev: ",std(diff))
    println()
    println(diff)

    return
end

derivtestS1S2(8)

println(1)

function derivtest2(Nx)
    eps = 1e-4
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
    im = rand(Float64, Nx, Nx).*0.1
    im[6,6]=1.0
    im1 = copy(im)
    im0 = copy(im)
    im1[2,3] += eps/2
    im0[2,3] -= eps/2
    dS20dp = wst_S20_deriv(im, fhash)

    der0=DHC_compute(im0,fhash,doS2=false,doS20=true)
    der1=DHC_compute(im1,fhash,doS2=false,doS20=true)
    dS = (der1-der0) ./ eps

    Nf = length(fhash["filt_index"])
    i0 = 3+Nf
    blarg = dS20dp[2,3,:,:]
    diff = dS[i0:end]-reshape(blarg,Nf*Nf)
    println(dS[i0:end])
    println("and")
    println(blarg)
    println("stdev: ",std(diff))
    println()
    println(diff)

    return
end

derivtest2(8)

Nx=128
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
im = zeros(Float64, Nx, Nx)
im[6,6]=1.0
@benchmark blarg = wst_S1_deriv(im, fhash)

# S1 deriv time, Jan 30 version (old) compared to Feb 14 version from NM
# Nx    old      Feb 14
#   32     17 ms   0.5 ms
#   64     34 ms   3.5 ms
#  128    115 ms    20 ms
#  256    520 ms    92 ms
#  512   2500 ms   720 ms   540 ms with 2 threads
# 1024   9500 ms  3300 ms  2500 ms with 2

Nx = 128
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
im = zeros(Float64, Nx, Nx)
im[6,6]=1.0
@benchmark blarg,blafb2 = wst_S1S2_derivfast(im, fhash)

# S20 vs S2 deriv time, Feb 15
# Nx     S20
#   8      1 ms    27 ms
#  16      7      104
#  32     50      300
#  64    400 ms   1.1 sec
# 128    3.3 sec  5.0
# 256   17.2 sec
# 512   ---

Nx = 32
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
im = zeros(Float64, Nx, Nx)
im[6,6]=1.0
@benchmark blarg = wst_S20_deriv(im, fhash)

Profile.clear()
@profile blarg = wst_S20_deriv(im, fhash)
Juno.profiler()

# S20 deriv time, Jan 30
# Nx     Jan 30  Feb 14
#   8     28 ms   1 ms
#  16    112      7
#  32    320     50
#  64   1000    400
# 128   5 sec     3.3 sec
# 256   ---      17.2 sec
# 512   ---

print(1)

# wst_S1_deriv agrees with brute force at 1e-10 level.
function derivtest(Nx)
    eps = 1e-4
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
    im = rand(Float64, Nx, Nx).*0.1
    im[6,6]=1.0
    im1 = copy(im)
    im0 = copy(im)
    im1[2,3] += eps/2
    im0[2,3] -= eps/2
    blarg = wst_S1_deriv(im, fhash)

    der0=DHC_compute(im0,fhash,doS2=false)
    der1=DHC_compute(im1,fhash,doS2=false)
    dS = (der1-der0) ./ eps

    diff = dS[3:end]-blarg[2,3,:]
    println(dS[3:end])
    println("and")
    println(blarg[2,3,:])
    println("stdev: ",std(diff))

    #plot(dS[3:end])
    #plot!(blarg[2,3,:])
    #plot(diff)
    return
end

derivtest(128)

Nx=32
im = rand(Nx,Nx)
@time fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
@benchmark Sarr = DHC_compute(im, fhash, doS2=false)  # 0.36 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=true)   # 6.4 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=false, doS12=true)  # 0.7 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=false, doS20=true)  # 2.1 ms

Nx=64
im = rand(Nx,Nx)
@time fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
@benchmark Sarr = DHC_compute(im, fhash, doS2=false)  # 1.3 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=true)   # 19 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=false, doS12=true)  # 2.8 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=false, doS20=true)  # 6.6 ms

Nx=256
im = rand(Nx,Nx)
@time fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
@benchmark Sarr = DHC_compute(im, fhash, doS2=false)  # 25 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=true)   # 375 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=false, doS12=true)  # 55 ms
@benchmark Sarr = DHC_compute(im, fhash, doS2=false, doS20=true)  # 130 ms

Nx = 64
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
im = rand(Nx,Nx)
S_targ = DHC_compute(im, fhash, doS2=false)

im = rand(Nx,Nx)

function wst_synth(im_init, fixmask)
    # fixmask -  0=float in fit   1=fixed

    function wst_synth_chisq(vec_in)

        thisim = copy(im_init)
        thisim[indfloat] = vec_in

        Sarr  = DHC_compute(thisim, fhash, doS2=false)
        diff  = (Sarr .- S_targ)[3:end]

        # should have some kind of weight here
        chisq = diff'*diff
        return chisq

    end

    function wst_synth_dchisq(storage, vec_in)

        thisim = copy(im_init)
        thisim[indfloat] = vec_in
        dSdp   = wst_S1_deriv(thisim, fhash)
        S1arr  = DHC_compute(thisim, fhash, doS2=false)
        diff  = (S1arr - S_targ)[3:end]

        # dSdp matrix * S1-S_targ is dchisq
        dchisq_im = (reshape(dSdp, Nx*Nx, Nf) * diff).*2
        dchisq = reshape(dchisq_im, Nx, Nx)[indfloat]

        storage .= dchisq
    end

    (Nx, Ny)  = size(im_init)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(fhash["filt_index"])

    # index list of pixels to float in fit
    indfloat = findall(fixmask .== 0)

    # initial values for floating pixels
    vec_init = im_init[indfloat]
    println(length(vec_init), " pixels floating")

    eps = zeros(size(vec_init))
    eps[1] = 1e-4
    chisq1 = wst_synth_chisq(vec_init+eps./2)
    chisq0 = wst_synth_chisq(vec_init-eps./2)
    brute  = (chisq1-chisq0)/1e-4

    clever = zeros(size(vec_init))
    _bar = wst_synth_dchisq(clever, vec_init)
    println("Brute:  ",brute)
    println("Clever: ",clever[1])

    # call optimizer
    res = optimize(wst_synth_chisq, wst_synth_dchisq, vec_init, BFGS())

    # copy results into pixels of output image
    im_synth = copy(im_init)
    im_synth[indfloat] = Optim.minimizer(res)
    println(res)

    return im_synth
end

Nx    = 128
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
im    = rand(Nx,Nx)
fixmask = im .> 0.5
S_targ = DHC_compute(im, fhash, doS2=false)

init = copy(im)
init[findall(fixmask .==0)] .= 0

foo = wst_synth(init, fixmask)

#Optim.minimizer(optimize(FOM, im, LBFGS()))

# using dchisq function
# size   t(BFGS) t(LBFGS) [sec]
# 16x16    1
# 32x32    4        9
# 64x64    20      52
# 128x128         153 (fitting 50% of pix)

# without dchisq function
# size   t(BFGS) t(LBFGS) [sec]
# 16x16   27
# 32x32  355      110
# 64x64          1645
# 128x128        est 28,000  (8 hrs)

function readdust()

    RGBA_img = load(pwd()*"/scratch_DF/t115_clean_small.png")
    img = reinterpret(UInt8,rawview(channelview(RGBA_img)))

    return img[1,:,:]
end

function wst_synthS20(im_init, fixmask, S_targ, S20sig; iso=false)
    # fixmask -  0=float in fit   1=fixed

    function wst_synth_chisq(vec_in)

        thisim = copy(im_init)
        thisim[indfloat] = vec_in

        M20 = fhash["S2_iso_mat"]

        S20 = DHC_compute(thisim, fhash, doS2=false, doS20=true, norm=false, iso=iso)
        i0 = 3+(iso ? N1iso : Nf)
        diff  = ((S20[i0:end] - S_targ)./S20sig)

        # should have some kind of weight here
        chisq = diff'*diff
        println(chisq)
        return chisq

    end

    function wst_synth_dchisq(storage, vec_in)

        thisim = copy(im_init)
        thisim[indfloat] = vec_in

        dS20dp = reshape(wst_S20_deriv(thisim, fhash), Nx*Nx, Nf*Nf)
        if iso
            M20 = fhash["S2_iso_mat"]
            dS20dp = dS20dp * M20'
        end
        i0 = 3+(iso ? N1iso : Nf)

        S20arr = (DHC_compute(thisim, fhash, doS2=false, doS20=true, norm=false, iso=iso))[i0:end]

        # put both factors of S20sig in this array to weight
        diff   = (S20arr - S_targ)./(S20sig.^2)
        #println("size of diff", size(diff))
        #println("size of dS20dp", size(reshape(dS20dp, Nx*Nx, Nf*Nf)))
        # dSdp matrix * S1-S_targ is dchisq
        #dchisq_im = (reshape(dS20dp, Nx*Nx, Nf*Nf) * diff).*2
        dchisq_im = (dS20dp * diff).*2
        dchisq = reshape(dchisq_im, Nx, Nx)[indfloat]

        storage .= dchisq
    end

    (Nx, Ny)  = size(im_init)
    if Nx != Ny error("Input image must be square") end
    (N1iso, Nf)    = size(fhash["S1_iso_mat"])

    # index list of pixels to float in fit
    indfloat = findall(fixmask .== 0)

    # initial values for floating pixels
    vec_init = im_init[indfloat]
    println(length(vec_init), " pixels floating")

    eps = zeros(size(vec_init))
    eps[1] = 1e-4
    chisq1 = wst_synth_chisq(vec_init+eps./2)
    chisq0 = wst_synth_chisq(vec_init-eps./2)
    brute  = (chisq1-chisq0)/1e-4

    clever = zeros(size(vec_init))
    _bar = wst_synth_dchisq(clever, vec_init)
    println("Brute:  ",brute)
    println("Clever: ",clever[1])

    # call optimizer
    #res = optimize(wst_synth_chisq, wst_synth_dchisq, copy(vec_init), BFGS())
    res = optimize(wst_synth_chisq, wst_synth_dchisq, copy(vec_init),
        ConjugateGradient(), Optim.Options(iterations=1000))

    # copy results into pixels of output image
    im_synth = copy(im_init)
    im_synth[indfloat] = Optim.minimizer(res)
    println(res)

    return im_synth
end

function S20_weights(im, fhash, Nsam=10; iso=iso)

    (Nx, Ny)  = size(im)
    if Nx != Ny error("Input image must be square") end
    (N1iso, Nf)    = size(fhash["S1_iso_mat"])

    # fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)

    S20   = DHC_compute(im, fhash, doS2=false, doS20=true, norm=false, iso=iso)
    Ns    = length(S20)
    S20arr = zeros(Float64, Ns, Nsam)
    for j=1:Nsam
        noise = rand(Nx,Nx)
        S20arr[:,j] = DHC_compute(im+noise, fhash, doS2=false, doS20=true, norm=false, iso=iso)
    end

    wt = zeros(Float64, Ns)
    for i=1:Ns
        wt[i] = std(S20arr[i,:])
    end

    return wt
end

# read dust map
dust = Float64.(readdust())
dust = dust[1:256,1:256]

Nx     = 16
doiso  = false
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(N1iso, Nf)    = size(fhash["S1_iso_mat"])
i0 = 3+(doiso ? N1iso : Nf)
img    = imresize(dust,(Nx,Nx))
#fixmask = rand(Nx,Nx) .< 0.1
fixmask = falses((Nx, Nx))

S_targ = DHC_compute(img, fhash, doS2=false, doS20=true, norm=false, iso=doiso)
S_targ = S_targ[i0:end]

init = copy(img)
floatind = findall(fixmask .==0)
init[floatind] .+= rand(length(floatind)).*50 .-25

S20sig = S20_weights(img, fhash, 100, iso=doiso)
S20sig = S20sig[i0:end]
foo = wst_synthS20(init, fixmask, S_targ, S20sig, iso=doiso)
S_foo = DHC_compute(foo, fhash, doS2=false, doS20=true, norm=false, iso=doiso)

plot_synth_QA(im, init, foo, fhash)

function plot_synth_QA(ImTrue, ImInit, ImSynth, fhash; fname="test2.png")

    # -------- define plot1 to append plot to a list
    function plot1(ps, image; clim=nothing, bin=1.0, fsz=16, label=nothing)
        # ps    - list of plots to append to
        # image - image to heatmap
        # clim  - tuple of color limits (min, max)
        # bin   - bin down by this factor, centered on nx/2+1
        # fsz   - font size
        marg   = 1mm
        nx, ny = size(image)
        nxb    = nx/round(Integer, 2*bin)

        # -------- center on nx/2+1
        i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
        i1 = min(nx,round(Integer, (nx/2)+nxb+1))
        lims = [i0,i1]
        subim = image[i0:i1,i0:i1]
        push!(ps, heatmap(image, aspect_ratio=:equal, clim=clim,
            xlims=lims, ylims=lims, size=(400,400),
            legend=false, xtickfontsize=fsz, ytickfontsize=fsz,#tick_direction=:out,
            rightmargin=marg, leftmargin=marg, topmargin=marg, bottommargin=marg))
        if label != nothing
            annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,:white,32))
        end
        return
    end

    # -------- initialize array of plots
    ps   = []
    clim  = (0,200)
    clim2 = (0,200).-100

    # -------- 6 panel QA plot
    plot1(ps, ImTrue, clim=clim, label="True")
    plot1(ps, ImSynth, clim=clim, label="Synth")
    plot1(ps, ImInit, clim=clim, label="Init")
    plot1(ps, ImInit-ImTrue, clim=clim2, label="Init-True")
    plot1(ps, ImInit-ImSynth, clim=clim2, label="Init-Synth")
    plot1(ps, ImSynth-ImTrue, clim=clim2, label="Synth-True")

    myplot = plot(ps..., layout=(3,2), size=(1400,2000))
    savefig(myplot, fname)
end

plot_synth_QA(im, init, foo, fhash)

# using dchisq function with S20
# size   t(BFGS) t(LBFGS) [sec]
# 8x8      30
# 16x16    90
# 32x32   189                           74
# 64x64  1637                          531      642 iso
# 128x128                              2 hrs

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
"""
initroot!(toCall::Stack{Function}, ::ILDSearch, model::CPModel, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection)

generic function to instantiate the research based on a specific Strategy <: SearchStrategy. The max discrepancy correspond to the number of branchable variables
at the beginning of the search. Calls to expandIlds! with a decreasing discrepancy is stacked in the toCall Stack.
"""
function initroot!(toCall::Stack{Function}, search::ILDSearch , model::CPModel, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection)

    # Note that toCall stack is a LIFO data structure, expandIlds with a discrepancy threshold of 0 will be the first one to execute (then with 1, 2, 3, etc.)
    for k in search.d:-1:1
        push!(toCall, (model) -> (restart_search!(model); expandIlds!(toCall,k, model, variableHeuristic, valueSelection)))
    end
    return expandIlds!(toCall, 0, model, variableHeuristic, valueSelection,nothing)
end

"""
        expandIlds!(toCall::Stack{Function}, discrepancy::Int64, previousdepth::Int64, direction::Union{Nothing, Symbol} , model::CPModel, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection, newConstraints=nothing)

This function fills the toCall Stack (LIFO) and perform a recursive Limited Discrepancy Search. Some procedures will contain a call to `expandIlds!` itself. Each `expandIlds!`
call is wrapped around a `saveState!` and a `restoreState!` to be able to backtrack thanks to the trailer.

This implementation is based on this paper : Limited Discrepancy Search - 1995 - William  D.  Harvey  and  Matthew  L.  Ginsberg. This method is not efficiant compared to the Korf approach but doesn't need any given max depth of the search
tree, which is unknown for CP search tree.
We should maybe look at this : https://www.researchgate.net/publication/220639800_Limited_discrepancy_search_revisited
"""
function expandIlds!(toCall::Stack{Function}, discrepancy::Int64, model::CPModel, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection, newConstraints=nothing; prunedDomains::Union{CPModification,Nothing}=nothing)
    # Dealing with limits
    model.statistics.numberOfNodes += 1
    model.statistics.numberOfNodesBeforeRestart += 1

    if !belowTimeLimit(model)
        return :TimeLimitStop
    end
    if !belowNodeLimit(model)
        return :NodeLimitStop
    end
    if !belowSolutionLimit(model)
        return :SolutionLimitStop
    end
    # Fix-point algorithm
    feasible, pruned = fixPoint!(model, newConstraints, prunedDomains)
    model.statistics.lastPruning = sum(map(x-> length(x[2]),collect(pruned)))

    if !feasible
        model.statistics.numberOfInfeasibleSolutions += 1
        model.statistics.numberOfInfeasibleSolutionsBeforeRestart += 1

        return :Infeasible
    end
    if solutionFound(model)
        #TODO understand this
        if (discrepancy == 0)
            triggerFoundSolution!(model)
            return :FoundSolution
        end
        return :alreadyFoundSolution
    end
    # Variable selection
    x = variableHeuristic(model)
    # Value selection
    v = valueSelection(DecisionPhase, model, x)

    if (discrepancy>0)
    push!(toCall, (model) -> (restoreState!(model.trailer); :BackTracking))
    push!(toCall, (model) -> (
        prunedDomains = CPModification();
        addToPrunedDomains!(prunedDomains, x, remove!(x.domain, v));
        expandIlds!(toCall, discrepancy-1, model, variableHeuristic, valueSelection, getOnDomainChange(x), prunedDomains=prunedDomains)
    ))
    push!(toCall, (model) -> (saveState!(model.trailer); :SavingState))
    end

    push!(toCall, (model) -> (restoreState!(model.trailer); :BackTracking))
    push!(toCall, (model) -> (
        prunedDomains = CPModification();
        addToPrunedDomains!(prunedDomains, x, assign!(x, v));
        expandIlds!(toCall, discrepancy, model, variableHeuristic, valueSelection, getOnDomainChange(x), prunedDomains=prunedDomains)
    ))
    push!(toCall, (model) -> (saveState!(model.trailer); :SavingState))
    return :Feasible
end

ext: jl


text:
using Test
using LinearAlgebra
using Random
using Statistics
using LinearAlgebra
using CalibrateEDMF.ModelTypes
using CalibrateEDMF.ReferenceModels
using CalibrateEDMF.ReferenceStats
using CalibrateEDMF.TurbulenceConvectionUtils
using CalibrateEDMF.HelperFuncs
using CalibrateEDMF.DistributionUtils

@testset "ReferenceStatistics" begin
    # Choose same SCM to speed computation
    data_dir = mktempdir()
    scm_dirs = repeat([joinpath(data_dir, "Output.Bomex.000000")], 2)
    # Reduce resolution and t_max to speed computation as well
    t_max = 4 * 3600.0
    dt_max = 30.0
    dt_min = 20.0
    io_frequency = 120.0
    namelist_args = [
        ("time_stepping", "t_max", t_max),
        ("time_stepping", "dt_max", dt_max),
        ("time_stepping", "dt_min", dt_min),
        ("stats_io", "frequency", io_frequency),
        ("grid", "dz", 150.0),
        ("grid", "nz", 20),
    ]

    kwargs_ref_model = Dict(
        :y_names => [["u_mean"], ["v_mean"]],
        :y_dir => scm_dirs,
        :scm_dir => scm_dirs,
        :case_name => repeat(["Bomex"], 2),
        :t_start => repeat([t_max - 2.0 * 3600], 2),
        :t_end => repeat([t_max], 2),
        :namelist_args => repeat([namelist_args], 2),
    )
    ref_models = construct_reference_models(kwargs_ref_model)
    run_reference_SCM.(ref_models, overwrite = false, run_single_timestep = false)

    # Test penalty behavior of ReferenceStatistics.get_profile
    filename = get_stats_path(scm_dirs[1])

    # if simulation lines up with ti, tf => no penalty
    y = ReferenceStats.get_profile(filename, ["u_mean", "v_mean"]; ti = 0.0, tf = t_max)
    @test maximum(y) < 1e4
    # if simulation end within dt margin => no penalty
    y = ReferenceStats.get_profile(filename, ["u_mean", "v_mean"]; ti = 0.0, tf = t_max + io_frequency)
    @test maximum(y) < 1e4
    # if simulation end outside of dt margin => penalty
    y = ReferenceStats.get_profile(filename, ["u_mean", "v_mean"]; ti = 0.0, tf = t_max + io_frequency + dt_max)
    @test all(y .> 1e4)
    # if simulation end after ti  => no penalty
    y = ReferenceStats.get_profile(filename, ["u_mean", "v_mean"]; ti = t_max - 0.01, tf = t_max)
    @test maximum(y) < 1e4
    # if simulation end before ti  => penalty
    y = ReferenceStats.get_profile(filename, ["u_mean", "v_mean"]; ti = t_max + 0.01, tf = t_max)
    @test all(y .> 1e4)

    # Test only tikhonov vs PCA and tikhonov
    pca_list = [false, true]
    norm_list = [false, true]
    mode_list = ["absolute", "relative"]
    dim_list = [false, true]
    for (perform_PCA, norm, tikhonov_mode, dim_scaling) in zip(pca_list, norm_list, mode_list, dim_list)
        ref_stats = ReferenceStatistics(
            ref_models;
            perform_PCA = perform_PCA,
            normalize = norm,
            variance_loss = 0.1,
            tikhonov_noise = 0.1,
            tikhonov_mode = tikhonov_mode,
            dim_scaling = dim_scaling,
            y_type = SCM(),
            Σ_type = SCM(),
        )

        @test pca_length(ref_stats) == size(ref_stats.Γ, 1)
        @test full_length(ref_stats) == size(ref_stats.Γ_full, 1)
        @test (pca_length(ref_stats) < full_length(ref_stats)) == perform_PCA
        for (ci, m) in enumerate(ref_models)
            y_case, _, _ = get_obs(m, SCM(), SCM(), norm, z_scm = get_z_obs(m))
            @test full_length(ref_stats, ci) == length(y_case)
            y_case_pca = ref_stats.pca_vec[ci]' * y_case
            @test pca_length(ref_stats, ci) == length(y_case_pca)
            case_full_inds = full_inds(ref_stats, ci)
            case_pca_inds = pca_inds(ref_stats, ci)
            @test ref_stats.y_full[case_full_inds] == y_case
            @test ref_stats.y[case_pca_inds] == y_case_pca
        end
        @test_throws ArgumentError pca_inds(ref_stats, 0)
        @test_throws ArgumentError pca_inds(ref_stats, 3)
        if perform_PCA == false
            # Tikhonov regularization results in variance inflation
            @test tr(ref_stats.Γ) > tr(ref_stats.Γ_full)
            # Since the same configuration is added twice,
            # the full covariance is singular,
            @test det(ref_stats.Γ_full) ≈ 0
            # But not the regularized covariance.
            @test det(ref_stats.Γ) > 0
            for ci in 1:2
                @test pca_length(ref_stats, ci) == full_length(ref_stats, ci)
                @test pca_inds(ref_stats, ci) == full_inds(ref_stats, ci)
            end
        end
    end

    # Verify that incorrect definitions throw error
    @test_throws AssertionError ReferenceStatistics(
        ref_models;
        perform_PCA = false,
        normalize = false,
        tikhonov_mode = "relative",
        y_type = SCM(),
        Σ_type = SCM(),
    )

    l2_reg = Dict("foo" => [0.1], "bar" => [0.3, 0.0, 0.4])
    reg_indices = flat_dict_keys_where(l2_reg, above_eps)
    @test length(reg_indices) == 3
    # Dict is not ordered
    @test reg_indices == [1, 3, 4]

    @testset "PCA" begin
        dofs = [5, 10, 25, 100, 250]
        ts = [30, 300, 200, 50, 40]
        fs = [0.2, 0.1, 0.05, 0.2, 0.1]
        for (dof, t, f) in zip(dofs, ts, fs)
            time_evol = rand(dof, t)
            μ = vcat(mean(time_evol, dims = 2)...)
            Γ = cov(time_evol, dims = 2)
            λ, eigvecs = eigen(Γ)
            λ_pca, P_pca = pca(Γ, f)
            len_pca = length(λ_pca)
            @test Diagonal(λ_pca) ≈ P_pca' * Γ * P_pca
            @test sum(λ_pca) < tr(Γ)
            @test sum(λ_pca) > (1 - f) * tr(Γ)
            @test sum(λ_pca) - λ[dof - len_pca + 1] < (1 - f) * tr(Γ)

            μ_pca, Γ_pca, P_pca = obs_PCA(μ, Γ, f)
            @test tr(Γ_pca) < tr(Γ)
            @test tr(Γ_pca) ≈ sum(λ_pca)
            @test tr(Γ_pca) > (1 - f) * tr(Γ)
            @test tr(Γ_pca) - λ[dof - len_pca + 1] < (1 - f) * tr(Γ)
        end
    end
end


text:
using TimeseriesPrediction, DelayEmbeddings
using Documenter, PyPlot, Literate, DocumenterMarkdown

PyPlot.ioff()
cd(@__DIR__)

# Expand Spatio-temporal examples using Literate:
function replace_includes(str)

    included = ["1Dfield_temporalprediction.jl",
    "2Dfield_crossprediction.jl", "2Dfield_temporalprediction.jl"]

    path = dirname(dirname(pathof(TimeseriesPrediction)))*"/examples/"

    for ex in included
        content = read(path*ex, String)
        str = replace(str, "include(\"$(ex)\")" => content)
    end
    return str
end
# Literate it:
Literate.markdown("src/stexamples.jl", "src/";
                  name = "stexamples", preprocess = replace_includes)
# Literate.notebook("src/tsprediction/stexamples.jl", "src/tsprediction/";
#                   name = "stexamples", preprocess = replace_includes)

makedocs(modules=[TimeseriesPrediction],
doctest=false, root = @__DIR__, format = Markdown())

close("all")

if !Sys.iswindows()
    deploydocs(
        deps   = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4",
        "python-markdown-math", "pygments", "pymdown-extensions"),
        repo   = "github.com/JuliaDynamics/TimeseriesPrediction.jl.git",
        target = "site",
        make = () -> run(`mkdocs build`)
    )
end


text:
"""
Abstract model type. All models are subtypes of this one, see *e.g.* [`AbstractAModel`](@ref)
"""
abstract type AbstractModel end

"""
Abstract structure that defines component models, which are structures that lists the processes that can be simulated for a given component.
For example [`LeafModels`](@ref) is a concrete implementation of an `AbstractComponentModel` for leaves.
"""
abstract type AbstractComponentModel <: AbstractModel end

"""
Assimilation (photosynthesis) abstract model. All photosynthesis models must be a subtype of
this.
"""
abstract type AbstractAModel <: AbstractModel end

"""
Stomatal conductance abstract model. All stomatal conductance models must be a subtype of
this.

An AbstractGsModel subtype struct must implement at least a g0 field.
"""
abstract type AbstractGsModel <: AbstractModel end

"""
Light interception abstract struct. All light interception models must be a subtype of this.
"""
abstract type AbstractInterceptionModel <: AbstractModel end

"""
Energy balance abstract struct. All energy balance models must be a subtype of this.
"""
abstract type AbstractEnergyModel <: AbstractModel end


text:
using CSV, DataFrames, Statistics
using AxisKeys, MonteCarloMeasurements, DimensionalData
using StatsBase, NamedTupleTools
using StanSample

df = CSV.read(joinpath(@__DIR__, "..", "..", "..", "data", "WaffleDivorce.csv"), DataFrame);
df.D = zscore(df.Divorce)
df.M = zscore(df.Marriage)
df.A = zscore(df.MedianAgeMarriage)
data = (N=size(df, 1), D=df.D, A=df.A, M= df.M)

stan5_1 = "
data {
    int < lower = 1 > N; // Sample size
    vector[N] D; // Outcome
    vector[N] A; // Predictor
}
parameters {
    real a; // Intercept
    real bA; // Slope (regression coefficients)
    real < lower = 0 > sigma;    // Error SD
}
transformed parameters {
    vector[N] mu;               // mu is a vector
    for (i in 1:N)
        mu[i] = a + bA * A[i];
}
model {
    a ~ normal(0, 0.2);         //Priors
    bA ~ normal(0, 0.5);
    sigma ~ exponential(1);
    D ~ normal(mu , sigma);     // Likelihood
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

m5_1s = SampleModel("m5.1s", stan5_1)
rc5_1s = stan_sample(m5_1s; data)

if success(rc5_1s)

    nt5_1s = read_samples(m5_1s, :particles)
    NamedTupleTools.select(nt5_1s, (:a, :bA, :sigma)) |> display

    ka_chns = read_samples(m5_1s, :keyedarray)
    ka_chns |> display

    da_chns = read_samples(m5_1s, :dimarray)
    da_chns |> display
end


text:
"""
    function temporalsubset(C::ClimGrid, startdate::Date, enddate::Date)

Returns the temporal subset of ClimGrid C. The temporal subset is defined by a start and end date.

"""
function temporalsubset(C::ClimGrid, datebeg::Tuple, dateend::Tuple)

    T = typeof(get_timevec(C)[1])
    timeV = get_timevec(C)
    idxtimebeg, idxtimeend = ClimateTools.timeindex(timeV, datebeg, dateend, T)

    startdate = ClimateTools.buildtimetype(datebeg, T)
    enddate = ClimateTools.buildtimetype(dateend, T)

    # some checkups
    @argcheck startdate <= enddate
    # @argcheck startdate >= C[1][Axis{:time}][1]
    # @argcheck enddate <= C[1][Axis{:time}][end]

    dataOut = C[1][Axis{:time}(idxtimebeg:idxtimeend)]

    # The following control ensure that a 1-timestep temporal subset returns a 3D Array with time information on the timestep. i.e. startdate == enddate
    if ndims(dataOut) == 2
        timeV = startdate
        latsymbol = Symbol(C.dimension_dict["lat"])
        lonsymbol = Symbol(C.dimension_dict["lon"])
        data2 = fill(NaN, (size(dataOut,1), size(dataOut, 2), 1))
        data2[:,:,1] = dataOut
        dataOut = AxisArray(data2, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]), Axis{:time}(C[1][Axis{:time}][:]))

    end

    return ClimGrid(dataOut, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency=C.frequency, experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)

end

"""
    resample(C::ClimGrid, startmonth::Int64, endmonth::Int64)

Return a resampled subset of ClimGrid C based on months.
"""
function resample(C::ClimGrid, startmonth::Int64, endmonth::Int64)

    @argcheck startmonth >= minimum(Dates.month.(C[1][Axis{:time}][:]))
    @argcheck startmonth <= maximum(Dates.month.(C[1][Axis{:time}][:]))

    if startmonth <= endmonth
        # Each matrix [:,:,i] represent data for a day
        datain = C.data.data
        # Date vector
        datevecin = C[1][Axis{:time}][:]
        # Where are the data between startmonth and endmonth
        index = (Dates.month.(datevecin) .<= endmonth) .&  (Dates.month.(datevecin) .>= startmonth)
        # Keep only data between startmonth and endmonth
        dataout = datain[:,:,index]
        datevecout = datevecin[index]
        # Create the ClimGrid output
        lonsymbol = Symbol(C.dimension_dict["lon"])
        latsymbol = Symbol(C.dimension_dict["lat"])
        axisout = AxisArray(dataout, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]),    Axis{:time}(datevecout))

    elseif endmonth < startmonth
        # Each matrix [:,:,i] represent data for a day
        datain = C.data.data
        # Date vector
        datevecin = C[1][Axis{:time}][:]
        # Where are the data between startmonth and endmonth
        index = (Dates.month.(datevecin) .<= endmonth) .|  (Dates.month.(datevecin) .>= startmonth)
        # Keep only data between startmonth and endmonth
        dataout = datain[:,:,index]
        datevecout = datevecin[index]
        # Create the ClimGrid output
        lonsymbol = Symbol(C.dimension_dict["lon"])
        latsymbol = Symbol(C.dimension_dict["lat"])
        axisout = AxisArray(dataout, Axis{lonsymbol}(C[1][Axis{lonsymbol}][:]), Axis{latsymbol}(C[1][Axis{latsymbol}][:]), Axis{:time}(datevecout))
    end

    return ClimGrid(axisout, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping,dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency=C.frequency, experiment=C.experiment, run=C.run, project=C.project,institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable,typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    resample(C::ClimGrid, season::String)

Return a resampled subset of ClimGrid C for a given season. Season options are: "DJF" (December-February), "MAM" (March-May), "JJA" (June-August), "SON" (September-November)
"""
function resample(C::ClimGrid, season::String)
    if lowercase(season) == "djf"
      D = resample(C, 12, 2)
    elseif lowercase(season) == "mam"
      D = resample(C, 3, 5)
    elseif lowercase(season) == "jja"
      D = resample(C, 6, 8)
    elseif lowercase(season) == "son"
      D = resample(C, 9, 11)
    else
      error("Wrong season name. Options are DJF (December-February; Winter), MAM (March-May; Spring), JJA (June-August; Summer) and SON (September-November; Fall)")
    end

    return D
end

"""
    periodmean(C::ClimGrid; startdate::Tuple, enddate::Tuple)

Mean of array data over a given period.
"""
function periodmean(C::ClimGrid; start_date::Tuple=(Inf, ), end_date::Tuple=(Inf,), level=1)

    if start_date != (Inf,) || end_date != (Inf,)
        C = temporalsubset(C, start_date, end_date)
    end

    datain = C.data.data

    # Mean and squeeze
    if ndims(datain) == 2
        dataout = datain
    elseif ndims(datain) == 3
        if size(datain, 3) == 1 # already an average on single value
            dataout = dropdims(datain, dims=3)
        else
            dataout = dropdims(Images.meanfinite(datain, 3), dims=3)
        end
    elseif ndims(datain) == 4
        if size(datain, 4) == 1
            dataout = dropdims(datain[:, :, level, :], dims = 3)
        else
            dataout = dropdims(Images.meanfinite(datain, 3), dims=3)
        end
    end

    # Build output AxisArray
    FD = buildarray_climato(C, dataout)

    # Return climGrid type containing the indice
    return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency=C.frequency, experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable="periodmean", typeofvar=C.typeofvar, typeofcal="climatology", varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    timeindex(timeVec, start_date, end_date, T)

Return the index of time vector specified by start_date and end_date. T is the DateTime type (see NCDatasets.jl documentation).
"""
function timeindex(timeV, datebeg::Tuple, dateend::Tuple, T)

    # Start Date
    if !isinf(datebeg[1])
        # Build DateTime type
        start_date = ClimateTools.buildtimetype(datebeg, T)
        # @argcheck start_date >= timeV[1]
        idxtimebeg = findfirst(timeV .>= start_date)[1]
    else
        idxtimebeg = 1
    end
    # End date
    if !isinf(dateend[1])
        # Build DateTime type
        end_date = ClimateTools.buildtimetype(dateend, T)
        # @argcheck end_date <= timeV[end]
        idxtimeend = findlast(timeV .<= end_date)[1]
    else
        idxtimeend = length(timeV)
    end

    if !isinf(datebeg[1]) && !isinf(dateend[1])
        @argcheck start_date <= end_date
    end
    return idxtimebeg, idxtimeend
end

# function timeindex(C::ClimGrid, datebeg::Tuple, dateend::Tuple)
#
#     T = typeof(get_timevec(C)[1])
#
#     # Start Date
#     if !isinf(datebeg[1])
#
#         # Build DateTime type
#         start_date = ClimateTools.buildtimetype(datebeg, T)
#
#         # Check
#         @argcheck start_date >= timeV[1]
#
#         idxtimebeg = findfirst(timeV .>= start_date)[1]
#     else
#         idxtimebeg = 1
#     end
#
#     # End date
#     if !isinf(dateend[1])
#
#         # Build DateTime type
#         end_date = ClimateTools.buildtimetype(dateend, T)
#
#         @argcheck end_date <= timeV[end]
#
#         idxtimeend = findlast(timeV .<= end_date)[1]
#     else
#         idxtimeend = length(timeV)
#     end
#
#     if !isinf(datebeg[1]) && !isinf(dateend[1])
#         @argcheck start_date <= end_date
#     end
#
#     return idxtimebeg, idxtimeend
# end

"""
    buildtimetype(datetuple, f)

Returns the adequate DateTime for temporal subsetting using DateType *f*
"""
function buildtimetype(date_tuple, f)

    if length(date_tuple) == 1
        dateout = f(date_tuple[1], 01, 01)
    elseif length(date_tuple) == 2
        dateout = f(date_tuple[1], date_tuple[2], 01)
    elseif length(date_tuple) == 3
        dateout = f(date_tuple[1], date_tuple[2], date_tuple[3])
    elseif length(date_tuple) == 4
        dateout = f(date_tuple[1], date_tuple[2], date_tuple[3], date_tuple[4], 00, 00)
    elseif length(date_tuple) == 5
        dateout = f(date_tuple[1], date_tuple[2], date_tuple[3], date_tuple[4], date_tuple[5], 00)
    elseif length(date_tuple) == 6
        dateout = f(date_tuple[1], date_tuple[2], date_tuple[3], date_tuple[4], date_tuple[5], date_tuple[6])
    end

    return dateout
end

"""
    daymean(C::ClimGrid)

Returns the daily average given a sub-daily ClimGrid.
"""
function daymean(C::ClimGrid)

    datain = C[1].data
    timevec   = get_timevec(C)
    T = typeof(timevec[1])

    nbdays = unique(yearmonthday.(timevec))
    nbdays_len = length(nbdays)
    dataout = Array{typeof(datain[1])}(undef, (size(C[1], 1), size(C[1], 2), nbdays_len))
    newtime = Array{T}(undef, nbdays_len)

    # loop over year-month-days
    Threads.@threads for iday in 1:nbdays_len

        daytmp = nbdays[iday]
        datefind = T(daytmp[1], daytmp[2], daytmp[3])
        idx = findall(x -> Dates.year(x) == Dates.year(datefind) && Dates.month(x) == Dates.month(datefind) && Dates.day(x) == Dates.day(datefind), timevec)

        dataout[:, :, iday] = Statistics.mean(datain[:, :, idx], dims=3)
        newtime[iday] = datefind

    end

    # Build output AxisArray
    FD = buildarray_resample(C, dataout, newtime)

    return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="day", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)

end

"""
    daysum(C::ClimGrid)

Returns the daily sum given a sub-daily ClimGrid C.
"""
function daysum(C::ClimGrid)

    datain = C[1].data
    timevec   = get_timevec(C)
    T = typeof(timevec[1])

    nbdays = unique(yearmonthday.(timevec))
    nbdays_len = length(nbdays)
    dataout = Array{typeof(datain[1])}(undef, (size(C[1], 1), size(C[1], 2), nbdays_len))
    newtime = Array{T}(undef, nbdays_len)

    # loop over year-month-days
    Threads.@threads for iday in 1:nbdays_len

        daytmp = nbdays[iday]
        datefind = T(daytmp[1], daytmp[2], daytmp[3])
        idx = findall(x -> Dates.year(x) == Dates.year(datefind) && Dates.month(x) == Dates.month(datefind) && Dates.day(x) == Dates.day(datefind), timevec)

        dataout[:, :, iday] = Statistics.sum(datain[:, :, idx], dims=3)
        newtime[iday] = datefind

    end

    # Build output AxisArray
    FD = buildarray_resample(C, dataout, newtime)

    return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="day", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)

end

"""
    monthmean(C::ClimGrid)

Returns monthly means of ClimGrid C.
"""
function monthmean(C::ClimGrid)

    datain = C[1].data
    timevec = get_timevec(C)
    T = typeof(timevec[1])

    nbmonth = unique(yearmonth.(timevec))
    nbmonth_len = length(nbmonth)
    dataout = Array{typeof(datain[1])}(undef, (size(C[1], 1), size(C[1], 2), nbmonth_len))
    newtime = Array{T}(undef, nbmonth_len)

    # loop over year-month-days
    Threads.@threads for imonth in 1:nbmonth_len

        monthtmp = nbmonth[imonth]
        datefind = T(monthtmp[1], monthtmp[2])
        idx = findall(x -> Dates.year(x) == Dates.year(datefind) && Dates.month(x) == Dates.month(datefind), timevec)

        dataout[:, :, imonth] = Statistics.mean(datain[:, :, idx], dims=3)
        newtime[imonth] = datefind

    end

    # Build output AxisArray
    FD = buildarray_resample(C, dataout, newtime)

    return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="day", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)

end

"""
    monthsum(C::ClimGrid)

Returns monthly sums of ClimGrid C.
"""
function monthsum(C::ClimGrid)

    datain = C[1].data
    timevec = get_timevec(C)
    T = typeof(timevec[1])

    nbmonth = unique(yearmonth.(timevec))
    nbmonth_len = length(nbmonth)
    dataout = Array{typeof(datain[1])}(undef, (size(C[1], 1), size(C[1], 2), nbmonth_len))
    newtime = Array{T}(undef, nbmonth_len)

    # loop over year-month-days
    Threads.@threads for imonth in 1:nbmonth_len

        monthtmp = nbmonth[imonth]
        datefind = T(monthtmp[1], monthtmp[2])
        idx = findall(x -> Dates.year(x) == Dates.year(datefind) && Dates.month(x) == Dates.month(datefind), timevec)

        dataout[:, :, imonth] = Statistics.sum(datain[:, :, idx], dims=3)
        newtime[imonth] = datefind

    end

    # Build output AxisArray
    FD = buildarray_resample(C, dataout, newtime)

    return ClimGrid(FD, longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency="day", experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)

end

"""
    timeresolution(timevec::Array{N,1} where N)

Return the time resolution of the vector timevec.

"""
function timeresolution(timevec::Array{N,1} where N)

    # timevec = (NetCDF.ncread(str, "time"))
    if length(timevec) > 1
        timediff = diff(timevec)[2]
        if timediff == 1. || timediff == 1
            return "24h"
        elseif round(timediff, digits=5) == round(12/24, digits=5)
            return "12h"
        elseif round(timediff, digits=5) == round(6/24, digits=5)
            return "6h"
        elseif round(timediff, digits=5) == round(3/24, digits=5)
            return "3h"
        elseif round(timediff, digits=5) == round(1/24, digits=5)
            return "1h"
        elseif round(timediff, digits=5) == 365.0 || round(timediff, digits=5) == 366.0 || round(timediff, digits=5) == 360.0
            return "Yearly"
        elseif timediff == 31.0 || timediff == 30.0 || timediff == 28.0 || timediff == 29.0
            return "Monthly"
        end
    else
        return "N/A"
    end
end

"""
    function pr_timefactor(rez::String)

Return the time factor that should be applied to precipitation to get accumulation for resolution "rez"

"""
function pr_timefactor(rez::String)

    if rez == "24h"
        return 86400.0
    elseif rez == "12h"
        return 43200.0
    elseif rez == "6h"
        return 21600.0
    elseif rez == "3h"
        return 10800.0
    elseif rez == "1h"
        return 3600.0
    elseif rez == "N/A"
        return 1.0
    end

end

"""
    function daymean_factor(rez::String)

Return the time factor that should be applied to precipitation to get accumulation for resolution "rez"

"""
function daymean_factor(rez::String)

    if rez == "24h" || rez == "day" || rez == "daily"
        return 1
    elseif rez == "12h"
        return 2
    elseif rez == "6h"
        return 4
    elseif rez == "3h"
        return 8
    elseif rez == "1h"
        return 24
    elseif rez == "N/A"
        return 1
    end

end

"""
    dropfeb29(C::ClimGrid)

Removes february 29th. Needed for bias correction.
"""
function dropfeb29(C::ClimGrid)
    date_vec = get_timevec(C) # [1][Axis{:time}][:]
    f = typeof(date_vec[1])
    feb29th = (Dates.month.(date_vec) .== Dates.month(f(2000, 2, 2))) .& (Dates.day.(date_vec) .== Dates.day(29))
    dataout = C[1].data[:, :, .!feb29th] # drop for data
    date_vec = date_vec[.!feb29th] # drop for calendar

    datevec_noleap = Array{DateTimeNoLeap}(undef, length(date_vec))

    for id = 1:size(date_vec, 1)
        datevec_noleap[id] = reinterpret(DateTimeNoLeap, date_vec[id])
    end

    dataout = buildarray_resample(C, dataout, datevec_noleap)

    return ClimGrid(dataout; longrid=C.longrid, latgrid=C.latgrid, msk=C.msk, grid_mapping=C.grid_mapping, dimension_dict=C.dimension_dict, timeattrib=C.timeattrib, model=C.model, frequency=C.frequency, experiment=C.experiment, run=C.run, project=C.project, institute=C.institute, filename=C.filename, dataunits=C.dataunits, latunits=C.latunits, lonunits=C.lonunits, variable=C.variable, typeofvar=C.typeofvar, typeofcal=C.typeofcal, varattribs=C.varattribs, globalattribs=C.globalattribs)
end

"""
    leapyears(datevec)

Returns the leap years contained into datevec.
"""
function leapyears(datevec)

    years = unique(Dates.year.(datevec))
    lyrs = years[Dates.isleapyear.(years)]

    return lyrs

end

function find_julianday_idx(julnb, ijulian, window)
    if ijulian <= window
        idx = @. (julnb >= 1) & (julnb <= (ijulian + window)) | (julnb >= (365 - (window - ijulian))) & (julnb <= 365)
    elseif ijulian > 365 - window
        idx = @. (julnb >= 1) & (julnb <= ((window-(365-ijulian)))) | (julnb >= (ijulian-window)) & (julnb <= 365)
    else
        idx = @. (julnb <= (ijulian + window)) & (julnb >= (ijulian - window))
    end
    return idx
end

"""
    yearmonthdayhour(dt::AbstractCFDateTime) -> (Int64, Int64, Int64, Int64)
Simultaneously return the year, month, day and hour parts of `dt`.
Author: Alexander-Barth (Github)
"""
yearmonthdayhour(dt::DT) where DT <: Dates.TimeType = (Dates.year(dt),Dates.month(dt), Dates.day(dt), Dates.hour(dt))


text:
for trial = 1:10
    lnum = length(labels)
    R_stats = zeros(lnum,3)
    seed_sets = spzeros(n,lnum)
    for lab = 1:length(labels)
        label = labels[lab]
        T = findnz(LabelMatrix[:,label])[1]
        nT = length(T)
        println("$lab \t $nT \t $label")

        # Generate a new seed set
        seednum = round(Int64,nT/20)
        grownum = round(Int64,min(nT*2,n/2-seednum))
        p = randperm(nT)
        Rstart = T[p[1:seednum]]
        OneHop = get_immediate_neighbors(H,Ht,Rstart)
        Rmore = BestNeighbors(H,d,Rstart,OneHop,grownum)
        R = union(Rmore,Rstart)
        Rs = findall(x->in(x,Rstart),R)     # Force seed nodes to be in output set
        prr, rer, f1r = PRF(T,R)
        R_stats[lab,:] = [prr, rer, f1r]
        seed_sets[Rstart,lab] .= 1
    end
# matwrite("Seed_sets_stack45_$trial.mat",Dict("R_stats"=>R_stats,"seed_sets"=>seed_sets))
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
# This file was generated, do not modify it. # hide
y = G / Diagonal(diag(G))

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
if !@isdefined(Erdos)
    include("../src/Erdos.jl")
    using Erdos
end

using BenchmarkTools
using Dates
import JLD: load, save

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 40.
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 0.05
BenchmarkTools.DEFAULT_PARAMETERS.samples = 50000

VERS = true ? "v0.7" : "v0.6"
bench_dir = Base.source_dir()
res_dir = joinpath(bench_dir, "results", VERS)
par_dir = joinpath(bench_dir, "parameters")
res_dir6 = joinpath(bench_dir, "results", "v0.6")
res_dir7 = joinpath(bench_dir, "results", "v0.7")

### ADD BENCHMARKS  ###############
suite = BenchmarkGroup()
GLIST = [Graph{Int64}, Network]
DGLIST = [DiGraph{Int64}, DiNetwork]
GROUPS = [
            "core",
            "generators",
            "flow",
            "centrality",
            "dismantling",
            "connectivity",
            "persistence",
            "shortestpaths",
            "traversals",
            "matching"
         ]

for group in GROUPS
    include("$group/$group.jl")
end

####  LOADING / SAVING PARAMS
function savepars(suite, groups=GROUPS)
    for group in groups
        d = joinpath(par_dir, group)
        !isdir(d) && mkdir(d)
        path = joinpath(d, "$(Date(now())).jld")
        println("saving $path")
        # BenchmarkTools.save(path, group, params(suite[group]))
        save(path, group, params(suite[group]))
    end
end

function loadpars!(suite, groups=GROUPS)
    for group in groups
        d = joinpath(par_dir,group)
        !isdir(d) && mkdir(d)
        files = readdir(d)
        if length(files) > 0
            dates = map(x -> Date(split(x, ['.'])[1]), files)
            f = joinpath(d, "$(maximum(dates)).jld")
            println("loading $f")
            # loadparams!(suite[group], BenchmarkTools.load(f, group), :evals, :samples)
            loadparams!(suite[group], load(f, group), :evals, :samples)
        end
    end
end

function saveres(res)
    for group in GROUPS
        d = joinpath(res_dir, group)
        !isdir(d) && mkdir(d)
        path = joinpath(d, "$(Date(now())).jld")
        # BenchmarkTools.save(path, group, res[group])
        save(path, group, res[group])
        println("Benchmarks' result saved in $path")
    end
end

loadres(prev::Int=0) = loadres(res_dir, prev)
function loadres(rdir::String, prev::Int=0)
    res = BenchmarkGroup()
    for group in GROUPS
        d = joinpath(rdir, group)
        !isdir(d) && mkdir(d)
        files = readdir(d)
        if length(files) > 0
            dates = sort(map(x -> Date(split(x, ['.'])[1]), files))
            f = joinpath(d, "$(dates[end-prev]).jld")
            # res[group] = BenchmarkTools.load(f, group)
            res[group] = load(f, group)
        end
    end
    return res
end

"""
example: myjudge("core", "edges")
"""
function myjudge(names::String...)
    s = res
    sold = resold
    for n in names
        s = s[n]
        sold = sold[n]
    end
    return judge(median(s),median(sold))
end

function runbench()
    res = run(suite, verbose=true)
    return res
end

comparebench(res,prev::Int=0) = comparebench(loadres(prev), res)

function comparebench67()
    comparebench(loadres(res_dir6),loadres(res_dir7))
end

function comparebench(resold, res)
    has_regressions = false
    has_improves = true
    for group in GROUPS
        println("GROUP $group")
        m = median(res[group])
        !haskey(resold, group) && continue
        mold = median(resold[group])
        judgement = judge(m, mold)
        println(m)
        regr = regressions(judgement)
        if length(regr) > 0
            has_regressions = false
            print_with_color(:red, "REGRESSIONS FOUND:\n")
            @show regr
            print_with_color(:red, "******************\n")
        end
        improvs = improvements(judgement)
        if length(improvs) > 0
            has_improves = true
            print_with_color(:green, "IMPROVEMENTS FOUND:\n")
            @show improvs
            print_with_color(:green, "******************\n")
        end
    end
end

loadpars!(suite)

println("Retune the benchmarks and save the parameters with `tune!(suite); savepars(suite)`")
println("Run the benchmarks with `res = runbench()`")
println("Compare to previous results with `comparebench(res)`")
println("Save the results with `saveres(res)`")


text:
################################################################################
#                              ACCESSER METHODS
################################################################################
"""
    collection_domains(domain::AbstractInfiniteDomain)

Return the array of domains associated with a `CollectionDomain`. Error if the input
domain is not a `CollectionDomain`.
"""
function collection_domains(domain::AbstractInfiniteDomain)
    type = typeof(domain)
    error("`collection_domains` not defined for infinite domain of type $type.")
end
function collection_domains(domain::CollectionDomain{S}
                         )::Vector{S} where {S <: InfiniteScalarDomain}
    return domain.domains
end

################################################################################
#                               BASIC EXTENSIONS
################################################################################
Base.length(s::CollectionDomain)::Int = length(collection_domains(s))
Base.length(s::MultiDistributionDomain)::Int = length(s.distribution)
Base.length(s::InfiniteScalarDomain)::Int = 1

################################################################################
#                           SUPPORT VALIDITY METHODS
################################################################################
"""
    supports_in_domain(supports::Union{Real, Vector{<:Real}, Array{<:Real, 2}},
                    domain::AbstractInfiniteDomain)::Bool

Used to check if `supports` are in the domain of `domain`. Returns `true` if
`supports` are in domain of `domain` and returns `false` otherwise.
This is primarily an internal method for performing checks but can be extended
for user-defined domain types. Extending this is optional, but recommended where
possible. Note by fallback, this returns `true` for unrecognized domain types such
that an error won't be thrown.
"""
function supports_in_domain end

# IntervalDomain
function supports_in_domain(supports::Union{Real, Vector{<:Real}},
                         domain::IntervalDomain)::Bool
    min_support = minimum(supports)
    max_support = maximum(supports)
    if min_support < JuMP.lower_bound(domain) || max_support > JuMP.upper_bound(domain)
        return false
    end
    return true
end

# UnivariateDistribution domain
function supports_in_domain(supports::Union{Real, Vector{<:Real}},
                         domain::UniDistributionDomain)::Bool
    return all(Distributions.insupport(domain.distribution, supports))
end

# MultivariateDistribution domain
function supports_in_domain(
    supports::Array{<:Real},
    domain::MultiDistributionDomain{<:Distributions.MultivariateDistribution{S}}
    )::Bool where {S <: Distributions.ValueSupport}
    if length(domain.distribution) != size(supports, 1)
        error("Support dimensions does not match distribution dimensions.")
    end
    return all(Distributions.insupport(domain.distribution, supports))
end

# CollectionDomain
function supports_in_domain(supports::Array{<:Real, 2}, domain::CollectionDomain)::Bool
    domains = collection_domains(domain)
    if length(domains) != size(supports, 1)
        error("Support dimensions does not match CollectionDomain dimensions.")
    end
    for i in eachindex(domains)
        if !supports_in_domain(supports[i, :], domains[i])
            return false
        end
    end
    return true
end

# Fallback
function supports_in_domain(supports, domain::AbstractInfiniteDomain)::Bool
    return true
end

################################################################################
#                              LOWER BOUND METHODS
################################################################################
"""
    JuMP.has_lower_bound(domain::AbstractInfiniteDomain)::Bool

Return `Bool` indicating if `domain` has a lower bound that can be determined. This
should be extended for user-defined infinite domains. It defaults to `false`
for unrecognized domain types.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> domain = InfiniteDomain(0, 1);

julia> has_lower_bound(domain)
true
```
"""
function JuMP.has_lower_bound(domain::AbstractInfiniteDomain)::Bool # fallback
    return false
end

# IntervalDomain
JuMP.has_lower_bound(domain::IntervalDomain)::Bool = true

# DistributionDomain (Univariate)
JuMP.has_lower_bound(domain::UniDistributionDomain)::Bool = true

# CollectionDomain
function JuMP.has_lower_bound(domain::CollectionDomain)::Bool
    for s in collection_domains(domain)
        if !JuMP.has_lower_bound(s)
            return false
        end
    end
    return true
end

"""
    JuMP.lower_bound(domain::AbstractInfiniteDomain)::Union{Real, Vector{<:Real}}

Return the lower bound of `domain` if one exists. This should be extended for
user-defined infinite domains if appropriate. Errors if `JuMP.has_lower_bound`
returns `false`. Extensions are enabled by `JuMP.has_lower_bound(domain)` and
`JuMP.lower_bound(domain)`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> domain = InfiniteDomain(0, 1);

julia> lower_bound(domain)
0.0
```
"""
function JuMP.lower_bound(domain::AbstractInfiniteDomain) # fallback
    type = typeof(domain)
    error("`JuMP.lower_bound` not defined for infinite domain of type $type.")
end

# IntervalDomain
JuMP.lower_bound(domain::IntervalDomain)::Float64 = domain.lower_bound

# DistributionDomain (Univariate)
function JuMP.lower_bound(domain::UniDistributionDomain)::Real
    return minimum(domain.distribution)
end

# CollectionDomain
function JuMP.lower_bound(domain::CollectionDomain)::Vector{<:Real}
    return [JuMP.lower_bound(i) for i in collection_domains(domain)]
end

"""
    JuMP.set_lower_bound(domain::AbstractInfiniteDomain,
                         lower::Union{Real, Vector{<:Real}})::AbstractInfiniteDomain

Set and return the lower bound of `domain` if such an operation makes sense. Errors
if the type of `domain` does not support this operation or has not been extended.
User-defined domain types should extend this if appropriate.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> domain = InfiniteDomain(0, 1);

julia> set_lower_bound(domain, 0.5)
[0.5, 1]
```
"""
function JuMP.set_lower_bound(domain::AbstractInfiniteDomain,
                              lower::Union{Real, Vector{<:Real}}) # fallback
    type = typeof(domain)
    error("`JuMP.set_lower_bound` not defined for infinite domain of type $type.")
end

# IntervalDomain
function JuMP.set_lower_bound(domain::IntervalDomain, lower::Real)::IntervalDomain
    return IntervalDomain(lower, domain.upper_bound)
end

# DistributionDomain
function JuMP.set_lower_bound(domain::Union{UniDistributionDomain,
                                         MultiDistributionDomain},
                              lower::Union{Real, Vector{<:Real}})
    error("Cannot set the lower bound of a distribution, try using " *
          "`Distributions.Truncated` instead.")
end

# CollectionDomain
function JuMP.set_lower_bound(domain::CollectionDomain, lower::Vector{<:Real})::CollectionDomain
    domains = collection_domains(domain)
    new_domains = [JuMP.set_lower_bound(domains[i], lower[i]) for i in eachindex(domains)]
    return CollectionDomain(new_domains)
end

################################################################################
#                              UPPER BOUND METHODS
################################################################################
"""
    JuMP.has_upper_bound(domain::AbstractInfiniteDomain)::Bool

Return `Bool` indicating if `domain` has a upper bound that can be determined. This
should be extended for user-defined infinite domains. It defaults to `false`
for unrecognized domain types.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> domain = InfiniteDomain(0, 1);

julia> has_upper_bound(domain)
true
```
"""
function JuMP.has_upper_bound(domain::AbstractInfiniteDomain)::Bool # fallback
    return false
end

# IntervalDomain
JuMP.has_upper_bound(domain::IntervalDomain)::Bool = true

# DistributionDomain (Univariate)
JuMP.has_upper_bound(domain::UniDistributionDomain)::Bool = true

# CollectionDomain
function JuMP.has_upper_bound(domain::CollectionDomain)::Bool
    for i in collection_domains(domain)
        if !JuMP.has_upper_bound(i)
            return false
        end
    end
    return true
end

"""
    JuMP.upper_bound(domain::AbstractInfiniteDomain)::Union{Real, Vector{<:Real}}

Return the upper bound of `domain` if one exists. This should be extended for
user-defined infinite domains if appropriate. Errors if `JuMP.has_upper_bound`
returns `false`. Extensions are enabled by `JuMP.has_upper_bound(domain)` and
`JuMP.upper_bound(domain)`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> domain = InfiniteDomain(0, 1);

julia> upper_bound(domain)
1.0
```
"""
function JuMP.upper_bound(domain::AbstractInfiniteDomain) # fallback
    type = typeof(domain)
    error("`JuMP.upper_bound` not defined for infinite domain of type $type.")
end

# IntervalDomain
JuMP.upper_bound(domain::IntervalDomain)::Float64 = domain.upper_bound

# DistributionDomain (Univariate)
function JuMP.upper_bound(domain::UniDistributionDomain)::Real
    return maximum(domain.distribution)
end

# CollectionDomain
function JuMP.upper_bound(domain::CollectionDomain)::Vector{<:Real}
    return [JuMP.upper_bound(i) for i in collection_domains(domain)]
end

"""
    JuMP.set_upper_bound(domain::AbstractInfiniteDomain,
                         upper::Real)::AbstractInfiniteDomain

Set and return the upper bound of `domain` if such an aoperation makes sense. Errors
if the type of `domain` does not support this operation or has not been extended.
User-defined domain types should extend this if appropriate.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> domain = InfiniteDomain(0, 1);

julia> set_upper_bound(domain, 0.5)
[0, 0.5]
```
"""
function JuMP.set_upper_bound(domain::AbstractInfiniteDomain, upper::Real) # fallback
    type = typeof(domain)
    error("`JuMP.set_lower_bound` not defined for infinite domain of type $type.")
end

# IntervalDomain
function JuMP.set_upper_bound(domain::IntervalDomain, upper::Real)::IntervalDomain
    return IntervalDomain(domain.lower_bound, upper)
end

# DistributionDomain
function JuMP.set_upper_bound(domain::Union{UniDistributionDomain,
                                         MultiDistributionDomain}, lower::Real)
    error("Cannot set the upper bound of a distribution, try using " *
          "`Distributions.Truncated` instead.")
end

# CollectionDomain
function JuMP.set_upper_bound(domain::CollectionDomain, lower::Vector{<:Real})::CollectionDomain
    domains = collection_domains(domain)
    new_domains = [JuMP.set_upper_bound(domains[i], lower[i]) for i in eachindex(domains)]
    return CollectionDomain(new_domains)
end

################################################################################
#                        SUPPORT AND LABEL GENERATION
################################################################################
"""
    AbstractSupportLabel

An abstract type for support label types. These are used to distinguish different
kinds of supports that are added to infinite parameters.
"""
abstract type AbstractSupportLabel end

"""
    All <: AbstractSupportLabel

This support label is unique in that it isn't associated with a particular set of
supports, but rather is used used to indicate that all supports should be used.
"""
struct All <: AbstractSupportLabel end

# Filler label for NoGenerativeSupports
struct _NoLabel <: AbstractSupportLabel end

"""
    PublicLabel <: AbstractSupportLabel

An abstract type used to denote that labels that should be given to the user by
default.
"""
abstract type PublicLabel <: AbstractSupportLabel end

"""
    UserDefined <: PublicLabel

A support label for supports that are supplied by the user directly to an infinite
parameter.
"""
struct UserDefined <: PublicLabel end

"""
    UniformGrid <: PublicLabel

A support label for supports that are generated uniformly accross a given interval.
"""
struct UniformGrid <: PublicLabel end

"""
    SampleLabel <: PublicLabel

An abstract type for labels of supports that are generated via some sampling technique.
"""
abstract type SampleLabel <: PublicLabel end

"""
    MCSample <: SampleLabel

A support label for supports that are generated via Monte Carlo Sampling.
"""
struct MCSample <: SampleLabel end

"""
    WeightedSample <: SampleLabel

A support label for supports that are generated by sampling from a statistical
distribution.
"""
struct WeightedSample <: SampleLabel end

"""
    Mixture <: PublicLabel

A support label for multi-dimensional supports that are generated from a variety
of methods.
"""
struct Mixture <: PublicLabel end

"""
    UniqueMeasure{S::Symbol} <: PublicLabel

A support label for supports that are provided from the `DiscreteMeasureData`
associated with a measure where a unique label is generated to distinguish those
supports. This is done by invoking [`generate_unique_label`](@ref).
"""
struct UniqueMeasure{S} <: PublicLabel end

"""
    MeasureBound <: PublicLabel

A support label for supports that are generated using the upper and lower bounds
for `FunctionalDiscreteMeasureData`.
"""
struct MeasureBound <: PublicLabel end

"""
    InternalLabel <: AbstractSupportLabel

An abstract type for support labels that are associated with supports that should
not be reported to the user by default.
"""
abstract type InternalLabel <: AbstractSupportLabel end

"""
    generate_unique_label()::Type{UniqueMeasure}

Generate and return a unique support label for measures.
"""
function generate_unique_label()::DataType
    return UniqueMeasure{gensym()}
end

# Define default values of sig_digits and num_supports keywords
const DefaultSigDigits = 12
const DefaultNumSupports = 10

# a user interface of generate_support_values
"""
    generate_supports(domain::AbstractInfiniteDomain
                      [method::Type{<:AbstractSupportLabel}];
                      [num_supports::Int = DefaultNumSupports,
                      sig_digits::Int = DefaultSigDigits]
                      )::Tuple{Array{<:Real}, DataType}

Generate `num_supports` support values with `sig_digits` significant digits in
accordance with `domain` and return them along with the correct generation label(s).
`IntervalDomain`s generate supports uniformly with label `UniformGrid` and
distribution domains generate them randomly accordingly to the
underlying distribution. Moreover, `method` indicates the generation method that
should be used. These `methods` correspond to parameter support labels. Current
labels that can be used as generation methods include (but may not be defined
for certain domain types):
- [`MCSample`](@ref): Uniformly distributed Monte Carlo samples.
- [`WeightedSample`](@ref): Monte Carlo samples that are weighted by an underlying PDF.
- [`UniformGrid`](@ref): Samples that are generated uniformly over the domain.

Extensions that employ user-defined infinite domain types and/or methods
should extend [`generate_support_values`](@ref) to enable this. Errors if the
`domain` type and /or methods are unrecognized. This is intended as an internal
method to be used by methods such as [`generate_and_add_supports!`](@ref).
"""
function generate_supports(domain::AbstractInfiniteDomain;
                           num_supports::Int = DefaultNumSupports,
                           sig_digits::Int = DefaultSigDigits
                           )::Tuple
    return generate_support_values(domain, num_supports = num_supports,
                                   sig_digits = sig_digits)
end

# 2 arguments
function generate_supports(domain::AbstractInfiniteDomain,
                           method::Type{<:AbstractSupportLabel};
                           num_supports::Int = DefaultNumSupports,
                           sig_digits::Int = DefaultSigDigits
                           )::Tuple
    return generate_support_values(domain, method,
                                   num_supports = num_supports,
                                   sig_digits = sig_digits)
end

"""
    generate_support_values(domain::AbstractInfiniteDomain,
                            [method::Type{MyMethod} = MyMethod];
                            [num_supports::Int = DefaultNumSupports,
                            sig_digits::Int = DefaultSigDigits]
                            )::Tuple{Array{<:Real}, Symbol}

A multiple dispatch method for [`generate_supports`](@ref). This will return
a tuple where the first element are the supports and the second is their
label. This can be extended for user-defined infinite domains and/or generation
methods. When defining a new domain type the default method dispatch should
make `method` an optional argument (making it the default). Otherwise, other
method dispatches for a given domain must ensure that `method` is positional
argument without a default value (contrary to the definition above). Note that the
`method` must be a subtype of either [`PublicLabel`](@ref) or [`InternalLabel`](@ref).
"""
function generate_support_values(domain::AbstractInfiniteDomain,
                                 args...; kwargs...)
    if isempty(args)
        error("`generate_support_values` has not been extended for infinite domains " *
              "of type `$(typeof(domain))`. This automatic support generation is not " *
              "implemented.")
    else
        error("`generate_support_values` has not been extended for infinite domains " *
              "of type `$(typeof(domain))` with the generation method `$(args[1])`. " *
              "This automatic support generation is not implemented.")
    end
end

# IntervalDomain and UniformGrid
function generate_support_values(domain::IntervalDomain,
                                 method::Type{UniformGrid} = UniformGrid;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits,
                                 )::Tuple{Vector{<:Real}, DataType}
    lb = JuMP.lower_bound(domain)
    ub = JuMP.upper_bound(domain)
    new_supports = round.(range(lb, stop = ub, length = num_supports),
                          sigdigits = sig_digits)
    return new_supports, method
end

# IntervalDomain and MCSample
function generate_support_values(domain::IntervalDomain,
                                 method::Type{MCSample};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits,
                                 )::Tuple{Vector{<:Real}, DataType}
    lb = JuMP.lower_bound(domain)
    ub = JuMP.upper_bound(domain)
    dist = Distributions.Uniform(lb, ub)
    new_supports = round.(Distributions.rand(dist, num_supports),
                          sigdigits = sig_digits)
    return new_supports, method
end

# UniDistributionDomain and MultiDistributionDomain (with multivariate only)
function generate_support_values(
    domain::Union{UniDistributionDomain, MultiDistributionDomain{<:Distributions.MultivariateDistribution}},
    method::Type{WeightedSample} = WeightedSample;
    num_supports::Int = DefaultNumSupports,
    sig_digits::Int = DefaultSigDigits
    )::Tuple{Array{<:Real}, DataType}
    dist = domain.distribution
    new_supports = round.(Distributions.rand(dist, num_supports),
                          sigdigits = sig_digits)
    return new_supports, method
end

# UniDistributionDomain and MCSample
function generate_support_values(
    domain::UniDistributionDomain,
    method::Type{MCSample};
    num_supports::Int = DefaultNumSupports,
    sig_digits::Int = DefaultSigDigits
    )::Tuple{Vector{Float64}, DataType}
    return generate_support_values(domain, WeightedSample; num_supports = num_supports,
                                   sig_digits = sig_digits)[1], method # TODO use an unwieghted sample...
end

# MultiDistributionDomain (matrix-variate distribution)
function generate_support_values(
    domain::MultiDistributionDomain{<:Distributions.MatrixDistribution},
    method::Type{WeightedSample} = WeightedSample;
    num_supports::Int = DefaultNumSupports,
    sig_digits::Int = DefaultSigDigits
    )::Tuple{Array{Float64, 2}, DataType}
    dist = domain.distribution
    raw_supports = Distributions.rand(dist, num_supports)
    new_supports = Array{Float64}(undef, length(dist), num_supports)
    for i in 1:size(new_supports, 2)
        new_supports[:, i] = round.(reduce(vcat, raw_supports[i]),
                                    sigdigits = sig_digits)
    end
    return new_supports, method
end

# Generate the supports for a collection domain
function _generate_collection_supports(domain::CollectionDomain, num_supports::Int,
                                       sig_digits::Int)::Array{Float64, 2}
    domains = collection_domains(domain)
    # build the support array transpose to fill in column order (leverage locality)
    trans_supports = Array{Float64, 2}(undef, num_supports, length(domains))
    for i in eachindex(domains)
        @inbounds trans_supports[:, i] = generate_support_values(domains[i],
                                                   num_supports = num_supports,
                                                   sig_digits = sig_digits)[1]
    end
    return permutedims(trans_supports)
end

function _generate_collection_supports(domain::CollectionDomain,
                                       method::Type{<:AbstractSupportLabel},
                                       num_supports::Int,
                                       sig_digits::Int)::Array{Float64, 2}
    domains = collection_domains(domain)
    # build the support array transpose to fill in column order (leverage locality)
    trans_supports = Array{Float64, 2}(undef, num_supports, length(domains))
    for i in eachindex(domains)
        @inbounds trans_supports[:, i] = generate_support_values(domains[i],
                                                   method,
                                                   num_supports = num_supports,
                                                   sig_digits = sig_digits)[1]
    end
    return permutedims(trans_supports)
end

# CollectionDomain (IntervalDomains)
function generate_support_values(domain::CollectionDomain{IntervalDomain},
                                 method::Type{UniformGrid} = UniformGrid;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(domain, num_supports, sig_digits)
    return new_supports, method
end

function generate_support_values(domain::CollectionDomain{IntervalDomain},
                                 method::Type{MCSample};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(domain, method, num_supports, sig_digits)
    return new_supports, method
end

# CollectionDomain (UniDistributionDomains)
function generate_support_values(domain::CollectionDomain{<:UniDistributionDomain},
                                 method::Type{WeightedSample} = WeightedSample;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(domain, num_supports, sig_digits)
    return new_supports, method
end

# CollectionDomain (InfiniteScalarDomains)
function generate_support_values(domain::CollectionDomain,
                                 method::Type{Mixture} = Mixture;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(domain, num_supports, sig_digits)
    return new_supports, method
end

# CollectionDomain (InfiniteScalarDomains) using purely MC sampling
# this is useful for measure support generation
function generate_support_values(domain::CollectionDomain,
                                 method::Type{MCSample};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(domain, method, num_supports, sig_digits)
    return new_supports, method
end

# For label All: dispatch to default methods
function generate_support_values(domain::AbstractInfiniteDomain, ::Type{All};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits)
    return generate_support_values(domain, num_supports = num_supports,
                                   sig_digits = sig_digits)
end


text:
struct BoolVarFactor <: AbstractFactor{UInt8}
    name::String
    boolvar::AbstractBoolVariate
end

getlevels(factor::BoolVarFactor) = ["True", "False"]

Base.length(var::BoolVarFactor) = length(var.boolvar)

function slice(factor::BoolVarFactor, fromobs::Integer, toobs::Integer, slicelength::Integer)
    f = (x::Bool -> x ? UInt8(1) : UInt8(2))
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(factor.boolvar, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, UInt8)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
function SimpleTimeStampedTableParser(arg0::jint, arg1::SimpleTimeStampedTableParser_RowConverter)
    return SimpleTimeStampedTableParser((jint, SimpleTimeStampedTableParser_RowConverter), arg0, arg1)
end

function parse(obj::SimpleTimeStampedTableParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", List, (InputStream, JString), arg0, arg1)
end


text:
################################################################################
#
#  Simplify
#
################################################################################

function simplify(K::NfRel{nf_elem}; cached::Bool = true, prec::Int = 100)
  Kabs, mK, mk = absolute_field(K, cached = false)
  OK = maximal_order(K)
  new_basis = Vector{nf_elem}(undef, degree(Kabs))
  B = pseudo_basis(OK)
  ideals = Dict{NfOrdIdl, Vector{nf_elem}}()
  for i = 1:length(B)
    I = B[i][2].num
    if !haskey(ideals, I)
      bas = lll_basis(I)
      ideals[I] = nf_elem[mK\(K(x)) for x in bas]
    end
  end
  ind = 1
  for i = 1:degree(OK)
    I = B[i][2]
    bI = ideals[I.num]
    el = mK\(B[i][1])
    for j = 1:length(bI)
      new_basis[ind] = divexact(el*bI[j], I.den)
      ind += 1
    end
  end
  O = NfOrd(new_basis)
  O.ismaximal = 1
  O.disc = absolute_discriminant(OK)
  if prec == 100
    OLLL = lll(O)
  else
    OLLL = lll(O, prec = prec)
  end
  el = _simplify(OLLL, mK)
  pel = mK(el)
  f = minpoly(pel)
  @assert degree(f) == degree(K)
  Ks = number_field(f, cached = cached, check = false)[1]
  mKs = hom(Ks, K, pel)
  return Ks, mKs
end

#Finds a small elements given by small combinations of the basis of O
#that generates the extension given by mK
function _simplify(O::NfOrd, mK::NfToNfRel)
  L = nf(O)
  #First, we choose the candidates
  B = basis(O, L)
  a = gen(L)
  Bnew = nf_elem[]
  nrep = min(3, degree(L))
  for i = 1:length(B)
    push!(Bnew, B[i])
    for j = 1:nrep
      push!(Bnew, B[i]+B[j])
      push!(Bnew, B[i]-B[j])
    end
  end

  #Now, we test for primitiveness.
  Lrel = codomain(mK)
  OLrel = maximal_order(Lrel)
  K = base_field(Lrel)
  OK = maximal_order(K)

  n = degree(Lrel)

  P, d = _find_prime(Lrel)
  p = minimum(P, copy = false)
  abs_deg = degree(P)*d
  FP, mFP = ResidueField(OK, P)
  mFP1 = extend_easy(mFP, K)
  #First, we search for elements that are primitive using block systems
  F = FlintFiniteField(p, abs_deg, "w", cached = false)[1]
  emb = find_embedding(FP, F)
  rt = roots(map_coeffs(emb, map_coeffs(mFP1, Lrel.pol)))

  indices = Int[]
  for i = 1:length(Bnew)
    if isone(denominator(Bnew[i]))
      continue
    end
    if _is_primitive_via_block(mK(Bnew[i]), rt, mFP1, emb)
      push!(indices, i)
    end
  end
  #Now, we select the one of smallest T2 norm
  I = t2(a)
  for i = 1:length(indices)
    t2n = t2(Bnew[indices[i]])
    if t2n < I
      a = Bnew[indices[i]]
      I = t2n
    end
  end
  return a
end

function find_embedding(F::FqFiniteField, K::FqFiniteField)
  f = defining_polynomial(F)
  rt = roots(f, K)
  img = rt[1]
  return (x -> sum([coeff(x, i)*img^i for i = 0:degree(F)-1]))
end

function _is_primitive_via_block(a::NfRelElem{nf_elem}, rt::Vector{fq}, mF, emb)
  n = degree(parent(a))
  pol = data(a)
  polF = map_coeffs(emb, map_coeffs(mF, pol))
  nconjs = 1
  conjs = Set{fq}([evaluate(polF, rt[1])])
  for i = 2:length(rt)
    ev = evaluate(polF, rt[i])
    if ev in conjs
      return false
    end
    push!(conjs, ev)
    nconjs += 1
    if nconjs > div(n, 2)
      return true
    end
  end
  error("Something went wrong")
end

function _block(a::NfRelElem{nf_elem}, rt::Vector{fq}, mF, emb)
  pol = data(a)
  polF = map_coeffs(emb, map_coeffs(mF, pol))
  evs = fq[evaluate(polF, x) for x in rt]
  b = Vector{Int}[]
  a = BitSet()
  i = 0
  n = length(rt)
  while i < length(evs)
    i += 1
    if i in a
      continue
    end
    z = evs[i]
    push!(b, findall(x-> evs[x] == z, 1:n))
    for j in b[end]
      push!(a, j)
    end
  end
  return b
end

function _find_prime(L::NfRel{nf_elem})
  p = 2^10
  K = base_field(L)
  OK = maximal_order(K)
  OL = maximal_order(L)

  n_attempts = min(degree(L), 10)
  candidates = Vector{Tuple{NfOrdIdl, Int}}(undef, n_attempts)
  i = 1
  f = L.pol
  threshold = degree(f)^2
  while i < n_attempts+1
    p = next_prime(p)
    if isindex_divisor(OK, p)
      continue
    end
    lp = prime_decomposition(OK, p)
    P = lp[1][1]
    if isindex_divisor(OL, P)
      continue
    end
    F, mF = ResidueField(OK, P)
    mF1 = extend_easy(mF, K)
    fF = map_coeffs(mF1, f)
    if degree(fF) != degree(f) || !issquarefree(fF)
      continue
    end
    FS = factor_shape(fF)
    d = lcm(Int[x for (x, v) in FS])
    if d < threshold
      candidates[i] = (P, d)
      i += 1
    end
  end
  res =  candidates[1]
  for j = 2:n_attempts
    if candidates[j][2] < res[2]
      res = candidates[j]
    end
  end
  return res[1], res[2]
end

function _sieve_primitive_elements(B::Vector{T}) where T <: NumFieldElem
  K = parent(B[1])
  n = absolute_degree(K)
  B1 = typeof(B)()
  k = div(n, 2)
  for x in B
    c = conjugates_arb(x, 16)
    isprimitive = true
    for i = 2:k+1
      for j = 1:i-1
        if overlaps(c[i], c[j])
          isprimitive = false
          break
        end
      end
      if !isprimitive
        break
      end
    end
    if isprimitive
      push!(B1, x)
    end
  end
  return B1
end

function _find_prime(L::NfRelNS{nf_elem})
  p = 2^10
  K = base_field(L)
  OK = maximal_order(K)
  OL = maximal_order(L)
  dL = numerator(discriminant(L, FlintQQ))

  n_attempts = min(degree(L), 10)
  candidates = Vector{Tuple{NfOrdIdl, Int}}(undef, n_attempts)
  i = 1
  pols = L.pol
  threshold = absolute_degree(L)^2
  polsR = Vector{fq_poly}(undef, length(pols))
  while i < n_attempts+1
    p = next_prime(p)
    if isindex_divisor(OK, p) || divisible(dL, p)
      continue
    end
    lp = prime_decomposition(OK, p)
    P = lp[1][1]
    @assert !isindex_divisor(OL, P)
    F, mF = ResidueField(OK, P)
    mF1 = extend_easy(mF, K)
    is_proj = true
    for j = 1:length(pols)
      fF = isunivariate(map_coeffs(mF1, pols[j]))[2]
      if degree(fF) != total_degree(pols[j]) || !issquarefree(fF)
        is_proj = false
        break
      end
      polsR[j] = fF
    end
    if !is_proj
      continue
    end
    d = 1
    for j = 1:length(polsR)
	    FS = factor_shape(polsR[j])
      d1 = lcm(Int[x for (x, v) in FS])
      d = lcm(d, d1)
    end
    if d < threshold
      candidates[i] = (P, d)
      i += 1
    end
  end
  res = candidates[1]
  for j = 2:n_attempts
    if candidates[j][2]*degree(candidates[j][1]) < res[2]*degree(res[1])
      res = candidates[j]
    end
  end
  return res[1], res[2]
end

function _sieve_primitive_elements(B::Vector{NfRelNSElem{nf_elem}})
  Lrel = parent(B[1])
  #First, we choose the candidates
  Bnew = NfRelNSElem{nf_elem}[]
  nrep = min(3, absolute_degree(Lrel))
  for i = 1:length(B)
    push!(Bnew, B[i])
    for j = 1:nrep
      if i != j
        push!(Bnew, B[i]+B[j])
        push!(Bnew, B[i]-B[j])
      end
    end
  end
  #Now, we test for primitiveness.
  K = base_field(Lrel)
  OK = maximal_order(K)
  Zx = ZZ["x"][1]

  n = absolute_degree(Lrel)

  P, d = _find_prime(Lrel)
  p = minimum(P, copy = false)
  abs_deg = d*degree(P)
  #First, we search for elements that are primitive using block systems
  Fp = GF(p, cached = false)
  Fpx = PolynomialRing(Fp, cached = false)[1]
  F = FlintFiniteField(p, abs_deg, "w", cached = false)[1]
  Fx = PolynomialRing(F, cached = false)[1]
  rt_base_field = roots(Zx(K.pol), F)
  rt = Dict{fq, Vector{Vector{fq}}}()
  Rxy = PolynomialRing(F, ngens(Lrel), cached = false)[1]
  tmp = Fpx()
  for r in rt_base_field
    vr = Vector{Vector{fq}}()
    for f in Lrel.pol
      g = isunivariate(f)[2]
      coeff_gF = fq[]
      for i = 0:degree(g)
        nf_elem_to_gfp_fmpz_poly!(tmp, coeff(g, i))
        push!(coeff_gF, evaluate(tmp, r))
      end
      gF = Fx(coeff_gF)
      push!(vr, roots(gF))
    end
    rt[r] = vr
  end
  rt1 = Dict{fq, Vector{Vector{fq}}}()
  ind = 1
  nconjs_needed = div(n, 2)+1
  for (r, v) in rt
    rtv = Vector{Vector{fq}}()
    it = cartesian_product_iterator([1:length(v[i]) for i in 1:length(v)])
    for i in it
      push!(rtv, [v[j][i[j]] for j = 1:length(v)])
      ind += 1
      if ind > nconjs_needed
        break
      end
    end
    rt1[r] = rtv
    if ind > nconjs_needed
      break
    end
  end
  indices = Int[]
  for i = 1:length(Bnew)
    if length(vars(Bnew[i].data)) < ngens(Lrel)
      continue
    end
    if _is_primitive_via_block(Bnew[i], rt1, Rxy, tmp)
      push!(indices, i)
    end
  end
  return Bnew[indices]
end

function _is_primitive_via_block(a::NfRelNSElem{nf_elem}, rt::Dict{fq, Vector{Vector{fq}}}, Rxy, tmp)
  n = degree(parent(a))
  pol = data(a)
  conjs = Set{fq}()
  for (r, vr) in rt
    ctx = MPolyBuildCtx(Rxy)
    for (c, v) in zip(coeffs(pol), exponent_vectors(pol))
      nf_elem_to_gfp_fmpz_poly!(tmp, c)
      push_term!(ctx, evaluate(tmp, r), v)
    end
    g = finish(ctx)
    for i = 1:length(vr)
      ev = evaluate(g, vr[i])
      if ev in conjs
        return false
      end
      push!(conjs, ev)
    end
  end
  return true
end

function simplified_absolute_field(L::NfRelNS; cached = false)
  OL = maximal_order(L)
  B = lll_basis(OL)
  B1 = _sieve_primitive_elements(B)
  a = B1[1]
  I = t2(a)
  for i = 2:min(50, length(B1))
    J = t2(B1[i])
    if J < I
      a = B1[i]
      I = J
    end
  end
  f = absolute_minpoly(a)
  @assert degree(f) == absolute_degree(L)
  K = number_field(f, check = false, cached = cached)[1]
  mp = hom(K, L, a)
  return K, mp
end

function simplified_absolute_field(L::NfRel{nf_elem}; cached::Bool = false)
  OL = maximal_order(L)
  B = lll_basis(OL)
  B1 = _sieve_primitive_elements(B)
  a = B1[1]
  I = t2(a)
  for i = 2:min(50, length(B1))
    J = t2(B1[i])
    if J < I
      a = B1[i]
      I = J
    end
  end
  f = absolute_minpoly(a)
  @assert degree(f) == absolute_degree(L)
  K = number_field(f, check = false, cached = cached)[1]
  mp = hom(K, L, a)
  imp = inv(mp)
  return K, mp, hom(base_field(L), K,  imp(L(gen(base_field(L)))))
end


text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Jan 24th, 2020

"Get the optimized value of decision vector."
function value_vec(vec_x::Union{Array{VariableRef,1}, VariableRef})
    if isa(vec_x, VariableRef)
        vec_value = value(vec_x)
    else
        vec_value = [value(vec_x[i]) for i = 1:length(vec_x)]
    end

    return vec_value
end

"Get the optimized value of decision matrix."
function value_mat(mat_x::Array{VariableRef,2})
  m, n = size(mat_x)
  mat_value = [value(mat_x[i, j]) for i=1:m, j=1:n]

  return mat_value
end

"Get the optimized value of a 3-D decision matrix."
function value_mat3(mat3_x::Array{VariableRef,3})
  m, n, o = size(mat3_x)
  mat3_value =
  mat3_value = [value(mat3_x[i, j, k]) for i=1:m, j=1:n, k=1:o]

  return mat3_value
end

"Get the optimized value of a 4-D decision matrix."
function value_mat4(mat4_x::Array{VariableRef,4})
  m, n, o, p = size(mat4_x)
  mat4_value = [value(mat4_x[i, j, k, l]) for i=1:m, j=1:n, k=1:o, l=1:p]

  return mat4_value
end

"Get the optimized value of decision matrix."
function value_mat_extra1(mat_x)
  m, n = size(mat_x)
  mat_value = zeros(m, n+1)
  for i=1:m, j=2:(n+1)
    mat_value[i, j] = value(mat_x[i, j])
  end

  return mat_value
end

"Get the optimized value of a 3-D decision matrix."
function value_mat3_extra1(mat3_x)
  m, n, o = size(mat3_x)
  mat3_value = zeros(m, n+1, o)
  for i=1:m, j=2:(n+1), k=1:o
    mat3_value[i, j, k] = value(mat3_x[i, j, k])
  end

  return mat3_value
end

"Get the optimized value of a 4-D decision matrix."
function value_mat4_extra1(mat4_x)
  m, n, o, p = size(mat4_x)
  mat4_value = zeros(m, n+1, o, p)
  for i=1:m, j=2:(n+1), k=1:o, l=1:p
    mat4_value[i, j, k, l] = value(mat4_x[i, j, k, l])
  end

  return mat4_value
end


text:
## Low-level _Message type for sending/receiving small ZMQ messages directly, without the complications
## and overhead of the Message object for sharing buffers between Julia and libzmq.

# Low-level message type, matching the declaration of
# zmq_msg_t in the header: char _[64];
primitive type _Message 64 * 8 end

const _MessageOrRef = Union{_Message,Base.RefValue{_Message}}

function msg_init()
    zmsg = Ref{_Message}()
    rc = ccall((:zmq_msg_init, libzmq), Cint, (Ref{_Message},), zmsg)
    rc != 0 && throw(StateError(jl_zmq_error_str()))
    return zmsg
end

function msg_init(nbytes::Int)
    zmsg = Ref{_Message}()
    rc = ccall((:zmq_msg_init_size, libzmq), Cint, (Ref{_Message}, Csize_t), zmsg, nbytes % Csize_t)
    rc != 0 && throw(StateError(jl_zmq_error_str()))
    return zmsg
end

# note: no finalizer for _Message, so we need to call close manually!
function Base.close(zmsg::_MessageOrRef)
    rc = ccall((:zmq_msg_close, libzmq), Cint, (Ref{_Message},), zmsg)
    rc != 0 && throw(StateError(jl_zmq_error_str()))
    return nothing
end

Base.length(zmsg::_MessageOrRef) = ccall((:zmq_msg_size, libzmq), Csize_t, (Ref{_Message},), zmsg) % Int
Base.unsafe_convert(::Type{Ptr{UInt8}}, zmsg::_MessageOrRef) =
    ccall((:zmq_msg_data, libzmq), Ptr{UInt8}, (Ref{_Message},), zmsg)

# isbits data, vectors thereof, and strings can be converted to/from _Message

function _MessageRef(x::T) where {T}
    isbitstype(T) || throw(MethodError(_MessageRef, (x,)))
    n = sizeof(x)
    zmsg = msg_init(n)
    @preserve zmsg unsafe_store!(Ptr{T}(Base.unsafe_convert(Ptr{UInt8}, zmsg)), x)
    return zmsg
end

function _MessageRef(x::Vector{T}) where {T}
    isbitstype(T) || throw(MethodError(_MessageRef, (x,)))
    n = sizeof(x)
    zmsg = msg_init(n)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{UInt8}, Ptr{T}, Csize_t), zmsg, x, n)
    return zmsg
end

function _MessageRef(x::String)
    n = sizeof(x)
    zmsg = msg_init(n)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), zmsg, x, n)
    return zmsg
end

function unsafe_copy(::Type{Vector{T}}, zmsg::_MessageOrRef) where {T}
    isbitstype(T) || throw(MethodError(unsafe_copy, (T, zmsg,)))
    n = length(zmsg)
    len, remainder = divrem(n, sizeof(T))
    iszero(remainder) || error("message length $n not a multiple of sizeof($T)")
    a = Array{T}(undef, len)
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{T}, Ptr{UInt8}, Csize_t), a, zmsg, n)
    return a
end

function unsafe_copy(::Type{T}, zmsg::_MessageOrRef) where {T}
    isbitstype(T) || throw(MethodError(unsafe_copy, (T, zmsg,)))
    n = length(zmsg)
    n == sizeof(T) || error("message length $n ≠ sizeof($T)")
    return @preserve zmsg unsafe_load(Ptr{T}(Base.unsafe_convert(Ptr{UInt8}, zmsg)))
end

function unsafe_copy(::Type{String}, zmsg::_MessageOrRef)
    n = length(zmsg)
    return @preserve zmsg unsafe_string(Base.unsafe_convert(Ptr{UInt8}, zmsg), n)
end

unsafe_copy(::Type{IOBuffer}, zmsg::_MessageOrRef) = IOBuffer(unsafe_copy(Vector{UInt8}, zmsg))


text:
function solveNSLAE(D)
    basic_numbers = []
    rows, columns = size(D)

    if D[rows, columns] != 0
        for i in (columns-1):-1:1
            if D[rows, i] != 0
                break
            end

            if i == 1
                return nothing
            end
        end
    end

    for i in 1:rows
        for j in (columns - 1):-1:1
            if j == 1
                push!(basic_numbers, 1)
                break
            end
            if D[i, j] == 0
                push!(basic_numbers, j-1)
                break
            end
        end
    end

    D[rows, :] = D[rows, :]/D[rows,basic_numbers[1]]
    rows -= 1
    for k in 2:lastindex(basic_numbers)
        for i in 1:k
            D[rows, i:columns] = D[rows, i:columns] - D[rows, i] * D[rows + i, i:columns]
        end
        D[rows, :] = D[rows, :] / D[rows, basic_numbers[k]]
        rows -= 1
    end
    x = D[:, columns]
    return x
end